0 Research please specify Bahador Saket is with Georgia Tech. E-mail:
<EMAIL_ADDRESS>Dominik Moritz, Halden Lin, and Jeffrey Heer are with
University of Washington. E-mail: domoritz<EMAIL_ADDRESS>Çağatay
Demiralp is with Columbia University. E-mail<EMAIL_ADDRESS>Victor
Dibia is with IBM Research. E-mail<EMAIL_ADDRESS>Saket and Moritz et
al.: Beyond Heuristics: Learning Visualization Design

# Beyond Heuristics: Learning Visualization Design

Bahador Saket Bahador Saket and Dominik Moritz contributed equally to this
work. Dominik Moritz∗ Halden Lin  Victor Dibia  Çağatay Demiralp  Jeffrey Heer

###### Abstract

In this paper, we describe a research agenda for deriving design principles
directly from data. We argue that it is time to go beyond manually curated and
applied visualization design guidelines. We propose learning models of
visualization design from data collected using graphical perception studies
and build tools powered by the learned models. To achieve this vision, we need
to 1) develop scalable methods for collecting training data, 2) collect
different forms of training data, 3) advance interpretability of machine
learning models, and 4) develop adaptive models that evolve as more data
becomes available.

###### keywords:

Automated visualization design, machine learning, design guidelines,
visualization recommendation, feature engineering, visualization
recommendation systems

K.6.1Management of Computing and Information SystemsProject and People
ManagementLife Cycle; K.7.mThe Computing ProfessionMiscellaneousEthics

Introduction

The demand for data visualization has significantly grown in recent years with
the increasing availability and complexity of digitized data across everyday
domains. By visually encoding data properties, visualizations aim to enhance
data understanding by leveraging human visual perception, which has evolved
for fast pattern detection and recognition. Understanding the effectiveness of
a given visualization in achieving this goal is a fundamental pursuit in data
visualization research and has important implications in practice.
Visualization researchers regularly conduct empirical studies to investigate
how people decode and interpret visualizations (e.g., [6, 12, 27, 23, 26]).
Such empirical studies are important for understanding and improving the
effectiveness of visualizations. Indeed, design guidelines and heuristics that
we use today in data visualization are an accumulation of what we have learned
from such empirical studies over decades.

It is not, however, always possible to distill guidelines through analysis of
data collected by empirical studies due to various confounding factors. Even
when this is possible, derived guidelines might be inadequate for conveying
the subtleties present in user data or might be marred by those providing the
guidelines. Moreover, design guidelines provided by empirical studies often
make their way to visualization tools slowly for two main reasons. First, our
design knowledge is continually evolving as visualization researchers
regularly publish results from empirical studies and provide new sets of
guidelines for designing effective visualizations. Second, designers of
visualization tools need to spend a significant amount of time and effort to
manually incorporate these guidelines.

In recent years, there has been an increasing trend to publish data collected
from empirical studies with increased awareness of the importance of
replicability and reproducibility. Researchers have made large datasets of
experimental data of visualizations’ effectiveness publicly available (e.g.,
[10, 26, 12, 23]). We advocate learning models from this data and building
tools that automatically apply the learned model. We believe machine learning
models provide a practical opportunity to implicitly engineer the insights
provided by empirical user performance data into visualization systems. The
main advantage is that new guidelines can be applied in practice faster in an
unbiased and reproducible manner.

In this paper, we discuss how we might automatically derive visualization
design principles from experimental data in a scalable way by using machine
learning. We further categorize and discuss the existing approaches for
learning visualization designs to illustrate feasibility and how future
systems fit into these categories. We argue that visualization design
principles used today are also derived from data. However, these principles
had to be abstracted by a visualization researcher, taught by a teacher, and
applied manually by the designer. We propose that the next step for the
research community is to curate a knowledge base of design principles that can
be applied automatically. Next, we should aim to learn from data both design
principles and how to trade off among them. As more data becomes available, we
may one day be able to learn visualization design end-to-end. In the
following, we describe each of these steps, and future research directions to
achieve this vision.

##  1 Visualization Recommendation Systems

In this section, we discuss existing and previous work on automated
visualization design and recommendation engines that incorporate guidelines
derived from graphical perception studies. We categorize these systems into
knowledge-based, data-driven, and hybrid visualization design tools. These
categories are used in recommender systems [2] and represent where the
knowledge that the system uses to recommend visualizations comes from. Table 1
gives an overview of existing machine learning and knowledge-based systems.

###  1.1 Knowledge-Based Systems

A large body of existing automated visualization design tools focuses on
suggesting visualizations based on user-defined constraints (such as which
data attributes to visualize) and design constraints. Mackinlay’s APT [16]
leverages a compositional algebra to enumerate the space of visualizations. It
then uses _expressiveness_ and _effectiveness_ criteria based on the work by
Bertin [5] and Cleveland & McGill [6] to prune and rank visualizations.
_Expressiveness_ refers to the ability of a visualization to convey all and
only the facts in the data. _Effectiveness_ refers to the ability of a
visualization when the information it conveys is more readily perceived than
with other visualizations. The SAGE project [18] extends APT by taking into
account additional data properties such as cardinality, uniqueness, and
functional dependencies. Tableau’s Show Me [17] introduces a set of heuristic
rules to recommend visualizations. Following this line of work, CompassQL [28]
also uses similar heuristic rules to develop expressiveness and effectiveness
criteria to evaluate visualization options. However, CompassQL extends the
earlier work by using a set of hand-tuned scores to specify criteria such as
space efficiency and legibility based on visualization and data properties.

Many of the automated design tools discussed above prune and rank candidate
visualizations based on a set of manually curated design guidelines derived
from previous empirical studies. Designers of these tools often spent a
considerable amount of time and effort to incorporate the findings of the
previous and current empirical studies.

System | Recommender Type | Modeling Approach | Learning | Model | Input Features | Design Space
---|---|---|---|---|---|---
VizDeck | Data Driven | Basic Features | Yes | Linear | 5 Data properties per field | 9+ Types
Kopol | Data Driven | Basic Features | No | Tree | Data properties, Task | 5 Types
Kim et al. | Data Driven | Basic Features | No | Linear | Data properties, Task | 12 Scatterplots
Data2Vis | Data Driven | Raw Features | Yes | Deep | Raw Dataset | Vega-Lite
APT | Knowledge Based | Hand Tuned | No | Rules | Data Types, Importance ordering | APT
CompassQL | Knowledge Based | Hand Tuned | No | Linear | Partial Specification | Vega-Lite
Draco | Hybrid | Learning Trade-Offs | Yes | Linear | Partial Specification, Task | Vega-Lite
Table 1: Comparison of different approaches to model automated visualization
design by the type of recommender(section 1), the modeling approach based on
the kind of features that are used (section 2), whether the approach uses
machine learning with genralization, the type of model, the input to the
model, and the design space.

###  1.2 Data-Driven Systems

The visualization literature is no stranger to data-driven models elicited
through human-subject experiments. For example, low-level perceptual models
such as Weber-Fechner Law [9], Stevens’ Power Law [25], and perceptual kernels
[7] are all based on fitting parametric and non-parametric models to empirical
user data, informing low-level visual encoding design. While data-driven
models are prevalent, using data-driven models for automated visualization
design is a nascent area and only a handful of papers exist to date.

With advances in machine learning, more researchers in the visualization
community started taking initial steps towards developing machine learning
models to recommend visualizations. A machine learning model tries to best
predict what visualization is most appropriate based on the given inputs
(e.g., tasks, data attributes, etc.). Developers of visualization tools need
to hand-craft informative, discriminating, and independent features for
learning such models. To the best of our knowledge, VizDeck [11] was the first
attempt at learning a recommendation model for high-level visual design.
VizDeck is a web-based visualization recommender tool designed for exploratory
data analysis of unorganized data. Using users’ up and down votes on a gallery
of visualization, VizDeck learns to recommend charts that the user is most
likely going to vote up. It does so by learning a linear scoring function over
statistical properties of the dataset. Visualizations are picked from a corpus
of possible candidate visualizations.

Saket et al. [23] recently conducted a crowdsourced study to evaluate the
effectiveness of five visualization types (Table, Line Chart, Bar Chart,
Scatterplot, and Pie Chart) across 10 different visual analysis tasks and from
two different datasets. Based on their findings, they developed
Kopol111https://kopoljs.github.io/, a mini JavaScript prototype visualization
recommender that uses a decision tree to suggest visualizations based on the
given tasks and data attribute types. Kim et al. [12] also recently developed
a model for 12 scatterplot encodings, the task type, and the cardinality and
entropy of some data fields. During their crowd-sourced experiment, Kim et al.
had users perform tasks for a combination of features and used results to
create a ranking.

Luo et al. [15] conducted a study in which 100 participants annotated 33,412
visualizations as good/bad, and provided 285,236 pairwise comparisons between
visualizations. They used 42 public datasets to create the visualizations for
their experiment. Luo et al. then developed DeepEye [15], a visualization
recommender system that uses a binary classifier to decide if a visualization
is good or bad, and a supervised learning-to-rank model to rank the
visualizations based on users’ preferences data collected in their experiment.

A more recent trend in machine learning is to learn models on all available
features instead of hand-crafting good features—a labor intensive and often
biased process. For example, Data2Vis [8] is a neural translation model that
generates visualization specifications in the Vega-Lite grammar from
descriptions of a dataset. The model was trained on thousands of example pairs
of data and visualizations recommended by CompassQL.

###  1.3 Hybrid Systems

In knowledge-based systems, the system designer provides the knowledge about
visualization design. In data-driven systems the system learns a
recommendation model from data. Hybrid recommender systems are both knowledge-
based and data-driven. The system designer has full control over the
recommendation process but can augment the knowledge base with machine
learning. Recently, machine learning experts argued that learning with
structure must be a top priority for AI to achieve human-like abilities [3].

Draco [19] is a formal model that represents (1) visualizations as sets of
logical facts and (2) design principles as a collection of hard and soft
constraints over these facts. Using constraints, we can take theoretical
design knowledge and express it in a concrete, extensible, and testable form.
Compared to VizDeck and Kopol, Draco’s visualization model covers a larger set
of visualizations supported by Vega-Lite [24]. Similar to CompassQL [28],
Draco’s default model uses rules that are informed by empirical studies but
formalized by experts. To avoid ad-hoc choices and support an evolving
knowledge base, Draco can learn how to trade off among competing design rules
using a simple linear model. Draco’s learning algorithm uses learning to rank,
a machine learning method that enables it to learn a preference model from
ordered pairs of complete visualizations. Using this method, Draco can learn
without the need to normalize scores between perceptual studies with different
methods and conditions.

##  2 Features for Modeling Visualization Design

In this section, we discuss the different approaches to automated
visualization design with respect to the features that they use. In
particular, we discuss building models from basic features (visualization
type, data properties, task), feature engineering, learning design rules, and
learning without feature engineering. Table 1 provides an overview over the
systems discussed here.

In the discussion below, we assume that a machine learning model for
visualization design takes as input a set of features that can include some
specification of the visualization, data, and task and outputs a corresponding
score. This score can represent how preferred a visualization is (based on
effectiveness, how easy it to read the visualization, etc.). The magnitude of
the score may not be meaningful, but it provides a rank ordering of the
feature vectors (i.e., possible designs). An automated visualization design
system can use any model of this form by enumerating possible designs and
recommending the one with the highest score.

###  2.1 Modeling Using Basic Features

The simplest and fastest way to model a score is to use the results of
experimental studies of effectiveness. These models use data such as the type
of visualization, data properties, and task type as features. For example,
Kopol by Saket et al. [23] and Kim et al. [12] are examples of this approach.
They use data such as tasks, visualization types, and data attribute types as
features. VizDeck [11] is another example of a system that uses basic
features. Instead of learning from effectiveness studies, VizDeck [11] uses
users’ up and down votes to change the scoring of each visualization
alternative. In addition to the visualization type, VizDeck’s model uses
statistics of the data as input.

While it is often simple and fast to use the results of experimental studies
to learn models (similar to Kopol [23] and Kim et al. [12]), the design space
of these models is small since they only support a small set of visualization
types and a basic set of features. In practice, visualization designers might
want to consider a larger design space. Going forward, the visualization
community needs to develop models that cover a much broader design space. To
effectively discriminate data in a larger design space, models need to use
more expressive features. In the next section, we discuss a method for
learning generalizable models with expressive features for visualization
design.

###  2.2 Using Visualization Design Rules as Features

Another method for automating visualization design is to use the violations of
design rules as features for learning models. A design rule is a predicate
over properties of a visualization, the user’s task, and the data. For
example, “Do not use the same field for the x and the y channel of a
scatterplot” describes the properties scatterplot and the same field on x and
y should not occur together. Assuming a set of design rules, the feature
vector representation of a visualization design is the number of times a
design rule is being violated by the design. For example, a design that does
not violate any design rules can be represented as a feature vector with only
zeros. We can use these feature vectors to learn a model. Since each feature
corresponds to a design rule, the weight of each feature in the learned model
is a measure of the relevance of its corresponding rule. Many of the design
rules that we use today are not always prescriptive. Today’s visualization
designers need to decide what rules they prefer for the specific visualization
they are working on. Machine learning models use statistical models to handle
uncertainty and noise in training data and are thus well equipped to handle
design rules that are not prescriptive. A model that uses design rules as
features learns the trade-offs among competing rules from data and is more
deterministic than human designers when applying them.

Recently, Moritz et al. proposed Draco [19], a method to derive feature
vectors for a machine learning model from Vega-Lite specifications [24]. This
allows for more sophisticated features to be composed, through feature-
engineering, from the underlying grammar. The main idea in Draco’s learning
approach is to use the violations of _design rules as features_ as outlined
above.

The main advantage of this approach over pure data-driven learning is that the
model can incorporate existing knowledge into the algorithm and thus learn a
model that generalizes from little data. This approach is also a
generalization of learning from basic features (subsection 2.1), as every
feature about the visualization type or the data can be written as a predicate
over the structural representation of the visualization (e.g., in Vega-Lite)
or the data schema. For example, to support tasks, we may just add rules for
each task. However, one of the main challenges with this approach is the
limited availability of design rules encoded in the system. Researchers should
work on systematically enumerating design rules and encode them in a machine
readable form so that they can be used as features in machine learning
systems. In the next section, we discuss how machine learning may support them
in this endeavour.

###  2.3 Learning Features from Data

Systems that use design rules for feature engineering and then learn a linear
scoring function over a vector of violations are limited by these rules.
Design rules relate properties of a visualization, and learning these
relations is known in machine learning as structure learning. For example,
inductive logic programming methods can infer logic programs from specific
instances [21]. However, many algorithms assume no or very little noise.
Design rules, however, are not prescriptive and have exceptions. For example,
the design rule that quantitative scales should start at zero does not make
sense when the data has no meaningful zero such as temperature data. However,
that does not mean that the zero rule should not be used. Markov logic
networks [22] are statistical models that support Bayesian reasoning and thus
can handle this uncertainty well. Their structure can be learned [13]. It
remains an open question whether these learning methods produce reasonable
design rules from the experimental data that is available today. Moreover,
more work is needed to investigate how to design an experiment to collect data
that results in reliable rules. One approach may be to learn rules from data
but have experts confirm them.

The key to learning structured rules is the availability of large high quality
datasets of examples. However, as more data becomes available, we may be able
to skip the feature engineering step and learn an end-to-end model, as we
discuss in the next section.

###  2.4 Learning without Feature Engineering

A recent trend in machine learning is to learn models on all available
features instead of hand-crafting good features. In particular, deep learning
models shift the burden of feature engineering to the system by automatically
learning abstract features and complex relationships that would otherwise be
difficult to capture by humans. This also means models can be more flexible.
Machine learning on the full data has led to impressive results in computer
vision and machine translation among other areas. However, deep learning
models in particular are extremely data hungry often requiring millions of
training examples. The data also needs to cover a large fraction of the design
space. For example, to learn a visualization design system that synthesizes
visualizations not just from templates but a grammar, the model needs to first
learn the grammar. For instance, Data2Vis [8] uses $215\,000$ training
examples to learn a subset of the Vega-Lite grammar and recommendations from
CompassQL. For this reason, this approach is not practical yet and more work
is needed to collect enough high quality training data. A main concern with
the quality of the training data is also whether the data is representative of
the true distribution of visualizations in the wild. Machine learning models
are probabilistic models that rank examples higher if they have seen many
similar examples in their training corpus. Consequently, if the training data
is biased, the model may only recommend a single visualization type.
Understanding how the bias in training data affects neural models is an open
area of research.

##  3 Next Steps in Learning Visualization Models

We view the existing body of work as the first step towards moving beyond
heuristics and learning visualization design from data. Multiple avenues for
future work lie in designing more interpretable machine learning models,
developing scalable methods for collecting different forms of training data,
and designing adaptive and evolving models.

###  3.1 Tooling to Help Designers Understand Models

Machine learning models can remove the manual effort of writing and tweaking
rules by building rules on the fly. Moreover, as we are getting more data, it
is easy to retrain the machine learning models rather quickly and frequently,
thus improving the accuracy of the models. Despite the advantages of learning
models, potential downsides of incorporating machine learning models include
an extra layer of complexity and diminished transparency [1], both of which
can make it difficult for both developers and end users to understand how the
system works. Unlike cases where the designers apply visualization design
guidelines manually, incorporating machine learning models might result in
losing the ability to look at the designed guidelines and their underlying
rationale.

Going forward, designers incorporating the empirical data to learn models
should provide methods to better investigate the underlying rationale and
convey an understanding of how the trained models will behave in the future.
Such visualization systems should communicate a variety of technical
information, such as the most representative features used in training the
model, the level of correlation among different features, and others. This
would help designers better understand the underlying logic behind rules
extracted by the model. Moreover, we envision visualization systems combining
machine learning models with state-of-the-art human-computer interface
techniques capable of translating models into understandable and useful
explanation dialogues for end users. Such systems should explain to the end
users such things as: why does the learned model recommend a specific set of
visualizations? What factors does the learned model use to prune and rank
visualizations? How do user interactions with the system affect the
recommendations?

###  3.2 Collecting Training Data

In order to learn the design of effective visualizations from user data,
access to high quality data with sufficient variety and context is critical.

Perceptual studies measuring effectiveness as defined in a controlled setting
can provide training data to learn design guidelines. However, current data
from graphical perception studies are stored in different destinations all
over the web (e.g., different repositories, personal webpages, etc.). These
inconsistencies make it harder for researchers and designers to track and
access available datasets. Better organizing data collected by graphical
perception studies can make deriving models easier. Going forward, one
solution might be to create a single data repository where the community can
share data collected from these empirical studies, thus improving the
accessibility of the collected data.

Graphical perception experiments often have specific research questions they
set out to answer. These studies are typically conducted under different
conditions, varying sample sizes, varying datasets, and for a diverse set of
tasks. As such, they may provide useful but inherently partial and incomplete
data to be used as training data for a generalizable models. Going forward, we
need large-scale, focused data collection efforts to provide the desired
consistency, size, context, and variation required to train generalizable
models with large learning capacities. Active learning methods can guide the
data collection.

An alternative to collecting new data is to use existing corpora of
visualizations on the web, such as Tableau Public [20], Beagle [4], and the
Chartmaker directory222http://chartmaker.visualisingdata.com/ or figures in
papers from Viziometrics [14]. The design principles learned from this data
would reflect the reality of the kinds of visualizations scientists, data
analysts, and journalists create. However, the data may be confounded by the
tools used in a particular community, as well as network effects.

###  3.3 Adapting and Evolving Models

All systems have only partial knowledge of context and user intent. For
example, the analyst’s goals often depend on seeing a rendered chart with real
data before they realize what needs to be adjusted. As a specific example, an
analyst may decide to use a log scale upon seeing that the spread of data is
very large. Thus, it is crucial that recommendation systems support iterative
refinement of users’ intent. Moreover, individual preferences may mean that
the same model is not optimal for everyone. A core concern in machine learning
for visualization design should be how models can be adapted to the group or
individual using it. As such, we need to develop models that take into account
user feedback during visual data exploration. Ideally, the accuracy of such
models should increase as users interact with the system.

Even a model that adapts to users will never provide perfect recommendation,
as the models are limited. Current models are restricted by the number of
features they use and the data they were trained on. Models need to evolve and
expand as more data becomes available and researchers find new design rules.
However, this requires designers to spend a tremendous amount of time and
effort to combine existing and new data since the data are collected in
different formats. One possible next step towards creating adapting and
evolving models is to find a common language/format and destination to collect
the results of the empirical studies. Ideally, we can develop systems that
incorporate the incoming data and update the model automatically.

##  4 Conclusion

In the past, visualization recommendation systems have incorporated
visualization design guidelines through a manual process of curation and
application. In this paper, we argue it is time to move beyond this laborious
process that is hard to scale. We imagine a future in which these systems are
machine learning models learned from experimental data. To achieve this
future, steps must be taken in engineering robust and adaptable models,
developing tools for interpretability of the created models, and consolidating
data. Once achieved, however, new guidelines and data collected from graphical
perception studies can be applied in practice faster and in an unbiased and
reproducible manner.

### Acknowledgements

We thank the reviewers and Kevin Hu for their feedback.

## References

  * [1] A. Abdul, J. Vermeulen, D. Wang, B. Y. Lim, and M. Kankanhalli.  Trends and trajectories for explainable, accountable and intelligible systems: An hci research agenda.  CHI ’18, 2018. doi: 10 . 1145/3173574 . 3174156
  * [2] C. C. Aggarwal et al.  Recommender systems.  Springer, 2016.
  * [3] P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, C. Gulcehre, F. Song, A. Ballard, J. Gilmer, G. Dahl, A. Vaswani, K. Allen, C. Nash, V. Langston, C. Dyer, N. Heess, D. Wierstra, P. Kohli, M. Botvinick, O. Vinyals, Y. Li, and R. Pascanu.  Relational inductive biases, deep learning, and graph networks, 2018\.  arXiv:1806.01261v1.
  * [4] L. Battle, P. Duan, Z. Miranda, D. Mukusheva, R. Chang, and M. Stonebraker.  Beagle: Automated extraction and interpretation of visualizations from the web.  In Proceedings of CHI. ACM, 2018.
  * [5] J. Bertin.  Semiology of graphics: diagrams, networks, maps.  1983\.
  * [6] W. S. Cleveland and R. McGill.  Graphical perception: Theory, experimentation, and application to the development of graphical methods.  Journal of the American Statistical Association, 1984.
  * [7] Ç. Demiralp, M. Bernstein, and J. Heer.  Perceptual kernels for visualization design.  IEEE Vis (Proc. InfoVis), 2014.
  * [8] V. Dibia and Ç. Demiralp.  Data2Vis: Automatic generation of data visualizations using sequence to sequence recurrent neural networks.  CoRR, abs/1804.03126, 2018.
  * [9] G. Fechner.  Elements of psychophysics.  1966\.
  * [10] J. Heer, N. Kong, and M. Agrawala.  Sizing the horizon: The effects of chart size and layering on the graphical perception of time series visualizations.  In ACM Human Factors in Computing Systems (CHI), 2009.
  * [11] A. Key, B. Howe, D. Perry, and C. R. Aragon.  Vizdeck: self-organizing dashboards for visual analytics.  In Proceedings of the ACM SIGMOD International Conference on Management of Data, SIGMOD 2012, Scottsdale, AZ, USA, May 20-24, 2012, 2012\. doi: 10 . 1145/2213836 . 2213931
  * [12] Y. Kim and J. Heer.  Assessing effects of task and data distribution on the effectiveness of visual encodings.  Proc. EuroVis, 2018.
  * [13] S. Kok and P. Domingos.  Learning the structure of markov logic networks.  In Proceedings of the 22Nd International Conference on Machine Learning, ICML ’05. ACM, New York, NY, USA, 2005.
  * [14] P. Lee, J. West, and B. Howe.  Viziometrix: A platform for analyzing the visual information in big scholarly data.  In BigScholar Workshop (co-located at WWW), 2016.
  * [15] Y. Luo, X. Qin, N. Tang, G. Li, and X. Wang.  Deepeye: Creating good data visualizations by keyword search.  In Proceedings of the 2018 International Conference on Management of Data, SIGMOD ’18, pp. 1733–1736. ACM, New York, NY, USA, 2018\.
  * [16] J. Mackinlay.  Automating the design of graphical presentations of relational information.  ACM Transactions on Graphics, 1986.
  * [17] J. D. Mackinlay, P. Hanrahan, and C. Stolte.  Show Me: Automatic presentation for visual analysis.  IEEE Vis (Proc. InfoVis), 13, 2007.
  * [18] V. O. Mittal, G. Carenini, J. D. Moore, and S. Roth.  Describing complex charts in natural language: A caption generation system.  Computational Linguistics, 1998.
  * [19] D. Moritz, C. Wang, G. Nelson, H. Lin, A. Smith, B. Howe, and J. Heer.  Formalizing visualization design knowledge as constraints: Actionable and extensible models in draco.  In IEEE Vis (Proc. InfoVis), 2018.
  * [20] K. Morton, M. Balazinska, D. Grossman, R. Kosara, and J. Mackinlay.  Public data and visualizations: How are many eyes and tableau public used for collaborative analytics?  SIGMOD Rec., 2014.
  * [21] J. R. Quinlan.  Learning logical definitions from relations.  Machine Learning, 1990. doi: 10 . 1007/BF00117105
  * [22] M. Richardson and P. Domingos.  Markov logic networks.  Machine Learning, 62, 2006. doi: 10 . 1007/s10994-006-5833-1
  * [23] B. Saket, A. Endert, and Ç. Demiralp.  Task-based effectiveness of basic visualizations.  IEEE Vis (Proc. InfoVis), 2018.
  * [24] A. Satyanarayan, D. Moritz, K. Wongsuphasawat, and J. Heer.  Vega-Lite: A grammar of interactive graphics.  IEEE Vis (Proc. InfoVis), 2017.
  * [25] S. S. Stevens.  On the psychophysical law.  Psychological review, 1957.
  * [26] D. A. Szafir.  Modeling color difference for visualization design.  IEEE Vis (Proc. InfoVis), Jan 2018. doi: 10 . 1109/TVCG . 2017 . 2744359
  * [27] J. Talbot, J. Gerth, and P. Hanrahan.  Arc length-based aspect ratio selection.  IEEE Trans. Visualization & Comp. Graphics, 2011.
  * [28] K. Wongsuphasawat, D. Moritz, A. Anand, J. Mackinlay, B. Howe, and J. Heer.  Voyager: Exploratory analysis via faceted browsing of visualization recommendations.  IEEE Vis (Proc. InfoVis), 2016.


# Photometric Survey of Neptune’s Trojan Asteroids I: The Color Distribution

Larissa Markwardt Department of Physics, University of Michigan
450 Church Street
Ann Arbor, MI 48109-1107, USA  Hsing Wen Lin ( 林省文) Department of Physics,
University of Michigan
450 Church Street
Ann Arbor, MI 48109-1107, USA  David Gerdes Department of Astronomy,
University of Michigan
1085 South University Avenue
Ann Arbor, MI 48109-1107, USA  Department of Physics, University of Michigan
450 Church Street
Ann Arbor, MI 48109-1107, USA  Fred C. Adams Department of Astronomy,
University of Michigan
1085 South University Avenue
Ann Arbor, MI 48109-1107, USA  Department of Physics, University of Michigan
450 Church Street
Ann Arbor, MI 48109-1107, USA Larissa Markwardt<EMAIL_ADDRESS>
###### Abstract

In 2018, Jewitt identified the “The Trojan Color Conundrum”, namely that
Neptune’s Trojan asteroids (NTs) had no ultra-red members, unlike the the
nearby Kuiper Belt. Since then, numerous ultra-red NTs have been discovered,
seemingly resolving this conundrum (Lin et al., 2019; Bolin et al., 2023).
However, it is still unclear whether or not the Kuiper Belt has a color
distribution consistent with the NT population, as would be expected if it
were the source population. In this work, we present a new photometric survey
of 15 out of 31 NTs. We utilized the Sloan
$g^{\prime}$$r^{\prime}$$i^{\prime}$$z^{\prime}$ filters on the IMACS f/4
instrument which is mounted on the 6.5m Baade telescope. In this survey, we
identify four NTs as being ultra-red using a Principal Component Analysis
(PCA). This result brings the ratio of red to ultra-red NTs to 7.75:1, more
consistent with the corresponding Trans-Neptunian Object (TNO) ratio of
4-11:1. We also identify three targets as being blue (nearly Solar) in color.
Such objects may be C-type surfaces, but we see more of these blue NTs than
has been observed in the Kuiper Belt (Seccull et al., 2018). Finally, we show
that there are hints of a color-absolute magnitude (H) correlation, with
larger H (smaller sized, lower albedo) tending to be more red, but more data
is needed to confirm this result. The origin of such a correlation remains an
open question which will be addressed by future observations of the surface
composition of these targets and their rotational properties.

Neptune trojans (1097) – Multi-color photometry (1077) – CCD photometry (208)

##  1 Introduction

Trojan asteroids are planetary companions that reside in the asymmetric 1:1
mean-motion resonance of planets; these asteroids librate at the planet-Sun L4
and L5 Lagrange points, meaning that they have the same orbit as the planet
but librate about a point 60∘ ahead of (L4) or behind (L5) the planet.
Numerical simulations show that orbits of Trojan asteroids can be quite
stable, on order the age of the Solar System (Ćuk et al., 2012; Gomes &
Nesvorný, 2016; Lykawka et al., 2011). Therefore, the stable members of these
populations are likely relatively undisturbed remnants of our primordial
planetary disk. The physical properties of these populations can thus give us
a window into the early Solar System.

However, Neptune’s Trojan asteroids are not thought to have formed in-situ.
Rather, this population likely grew through capture of planetesimals during
the epoch of planetary migration, during which the outer planets migrated from
the location of their formation to their present day locations (Fernandez &
Ip, 1984; Malhotra, 1993, 1995; Hahn & Malhotra, 1999). Assuming Neptune
migrated significantly in its early evolution, the Lagrange points must have
also migrated with it (Kortenkamp et al., 2004) Therefore, the NT population
can be used to constrain migratory models (Gomes & Nesvorný, 2016; Nesvorný et
al., 2013, 2018; Pike et al., 2017a). Such migration would have occurred in
the first several hundred Myr in the history of the Solar System, so while
these objects may not have formed in-situ, they still are remnants of the very
early Solar System.

Such models show that primordial Jupiter Trojan populations do not survive
this planetary migration, indicating they must have originated from elsewhere
in the Solar System. (Roig & Nesvorný, 2015). Similarly, since the dynamics of
planetary migration likely dispersed any primordial NTs as well, from where
did the current population of NTs originate? The most likely source is the
nearby Kuiper Belt. If that were the case, one would expect these two
populations to be similar in size and color (surface composition). Regarding
the color of the KBOs, the bimodality of red ($g-i$ < 1.2) vs. ultra-red
($g-i$ > 1.2) members has been well established (Sheppard, 2010; Schwarz et
al., 2011; Hainaut et al., 2012; Peixinho et al., 2012; Sheppard, 2012;
Lacerda et al., 2014; Peixinho et al., 2015; Pike et al., 2017b; Wong & Brown,
2017; Schwamb et al., 2019). Similarly, the centaur population, small bodies
which orbit between Jupiter and Neptune, are thought to be fed by
planetesimals escaping the NT region (Horner & Lykawka, 2010). These objects
are also red/ultra-red in color (Peixinho et al., 2012, 2015).

Through 2018, no ultra-red NTs had been found, making their color distribution
distinctly different than their expected origins or offshoots. Termed the
“Trojan Color Conundrum”, this tension is not easy to resolve (Jewitt, 2018).
One explanation is that some sort of resurfacing has happened to the NT
population specifically that affected neither the centaurs or KBOs. Jupiter’s
Trojan population is also devoid of ultra-red members which is thought to be
due to thermal resurfacing (Luu & Jewitt, 1996; Jewitt, 2002). However, the
temperatures at the distance of Neptune are too cold for such a scenario to be
valid (Jewitt, 2018). Another potential explanation is collisional
resurfacing, which could strip the ultra-red crust off of the surfaces of
these bodies revealing a bluer surface underneath. One source of such
collisions could be Plutinos, 3:2 resonators with Neptune, which have
significant orbital overlap with the NT population (Almeida et al., 2009).
Such collisions are expected to occur when Plutinos have high libration
amplitudes, high eccentricities, and low inclinations; therefore, we would
expect the color distribution of NTs to be inclination-dependent as well,
where high inclination NTs avoid these collisions and retain their ultra-red
surfaces. Finally, this discrepancy could be due to a primordial boundary
between red/ultra-red bodies that was subsequently mixed by Neptune’s
migration (DeMeo & Carry, 2014; Neveu & Vernazza, 2019). Based on the exact
nature of the epochs of radial mixing, mass removal, and planet migration, the
resulting NT population could be devoid of ultra-red members while the Centaur
population is not (Neveu & Vernazza, 2019), but specific simulations of these
two populations have not been conducted. This hypothesis has been supported by
the discovery of two Trans-Neptunian Object (TNO)-like (red) objects all the
way in the asteroid belt (Hasegawa et al., 2021).

In 2019, the first ultra-red NT, 2013 VX30, was discovered (Lin et al., 2019),
and additional ultra-red NTs have been discovered since then (Bolin et al.,
2023). On the surface, these discoveries seem to resolve the conundrum.
However, the color distribution of NTs still appears distinct from that of
other TNO populations (Bolin et al., 2023). Further observations of the NT
population are needed to determine whether or not these distributions are
truly distinct.

The structure of this paper is as follows: Section 2 describes the design of
our photometric survey. Section 3 outlines our data reduction process. Section
4 presents the results of our survey. Section 5 discuss the meaning of our
results. Section 6 is conclusions drawn from these results.

##  2 Survey Design

The goal of this paper is to measure the optical colors of currently known NTs
in order to better understand the physical characteristics of their surfaces.
The targets are listed in Table 1. All of our targets have been previously
observed but not by the same survey. All of our targets, except 2015 VU207,
were already known to be stable for $\sim$Gyr (Lin et al., 2021, 2022).
Following the methods of Lin et al. (2022), we find that 2015 VU207 is also
stable for Gyr in our simulations.

We used the IMACS f/4 instrument on the 6.5m Baade telescope at Las Campanas
Observatory on 4 unique nights to observe this population. IMACS was most
suitable for this task with its optical wavelength coverage ($\sim$400 - 900
nm) and large FOV to account for the positional uncertainty of the targets.
The Sloan $g^{\prime}$$r^{\prime}$$i^{\prime}$$z^{\prime}$ filters were used
for our photometric measurements. In order to account for any variation due to
a target’s rotational period, we observed each target with “bounding”
$r^{\prime}$-band observations (i.e. each observation in a different filter
was preceded and followed by an observation in $r^{\prime}$). . We chose
$r^{\prime}$ to be the bounding observations since this filter reaches the
highest SNR in the shortest amount of time. The fast readout mode with 2x2
binning was used.

Table 1: NT targets of this survey. Columns: (1) Object Designation, has
previous color measurements taken from 1: Sheppard (2012), 2: Parker et al.
(2013), 3: Jewitt (2018), 4: Lin et al. (2019), 5: Bolin et al. (2023); (2)
Lagrange Point; (3) Eccentricity; (4) Inclination (∘); (5) Absolute Magnitude;
(6) Dates observed; (7) Measured ave. SDSS r-band magnitude; (8) Measured SDSS
g-r (mag); (9) Measured SDSS r-i (mag); (10) Measured SDSS r-z (mag); (11)
Color classification determined based on the Principal Component Analysis (see
Sec. 4.2).

Name | L4/L5 | e | i | H | Date Observed | Ave. r | $g-r$ | $r-i$ | $r-z$ | Color Class.
---|---|---|---|---|---|---|---|---|---|---
2006 RJ1031,2 | L4 | 0.03 | 8.2 | 7.56 | 113021 | 21.97 | 0.59 $\pm$ 0.045 | 0.16 $\pm$ 0.035 | 0.17 $\pm$ 0.058 | red
|  |  |  |  | 120222 | 21.88 | — | — | 0.24 $\pm$ 0.055 | indeterminate
2007 VL3051,2 | L4 | 0.07 | 28.1 | 8.51 | 113021 | 22.60 | 0.60 $\pm$ 0.054 | 0.25 $\pm$ 0.038 | -0.15 $\pm$ 0.109 | red
|  |  |  |  | 120222 | 22.60 | — | — | 0.30 $\pm$ 0.047 | indeterminate
2010 TS1913 | L4 | 0.05 | 6.6 | 8.07 | 113021 | 22.39 | 0.61 $\pm$ 0.029 | 0.30 $\pm$ 0.029 | 0.64 $\pm$ 0.078 | red
2011 SO2773 | L4 | 0.01 | 9.6 | 7.76 | 113021 | 22.43 | 0.60 $\pm$ 0.067 | — | —- | indeterminate
|  |  |  |  | 120222 | 22.53 | — | 0.57 $\pm$ 0.050 | 0.82 $\pm$ 0.047 | ultra-red
2012 UD1855 | L4 | 0.04 | 28.3 | 7.59 | 113021 | 22.32 | 0.61 $\pm$ 0.033 | 0.37 $\pm$ 0.045 | 0.12 $\pm$ 0.081 | red
2012 UV1775 | L4 | 0.07 | 20.8 | 9.28 | 113021 | 23.76 | 0.71 $\pm$ 0.058 | 0.23 $\pm$ 0.051 | — | red
2013 RL1245 | L4 | 0.03 | 10.1 | 8.83 | 113021 | 23.37 | 0.38 $\pm$ 0.075 | 0.54 $\pm$ 0.086 | 0.67 $\pm$ 0.128 | red
2013 TZ1875 | L4 | 0.07 | 13.1 | 8.19 | 113021 | 23.27 | 0.90 $\pm$ 0.053 | 0.30 $\pm$ 0.057 | — | ultra-red
2013 VX304,5 | L4 | 0.09 | 31.2 | 8.31 | 113021 | 22.60 | 1.01 $\pm$ 0.043 | 0.44 $\pm$ 0.043 | 0.86 $\pm$ 0.049 | ultra-red
|  |  |  |  | 091122 | 22.96 | 0.70 $\pm$ 0.104 | 0.47 $\pm$ 0.048 | 0.73 $\pm$ 0.045 | ultra-red
2014 RO745 | L4 | 0.05 | 29.5 | 8.39 | 120222 | 23.34 | 0.65 $\pm$ 0.052 | 0.42 $\pm$ 0.064 | 1.42 $\pm$ 0.069 | ultra-red
2014 SC3745 | L4 | 0.10 | 33.7 | 8.18 | 113021 | 23.24 | 0.43 $\pm$ 0.066 | 0.12 $\pm$ 0.081 | — | blue
2014 YB92${5}$ | L4 | 0.10 | 30.8 | 8.62 | 091222 | 23.41 | 0.46 $\pm$ 0.187 | 0.07 $\pm$ 0.100 | 0.36 $\pm$ 0.090 | blue
2015 VU2075 | L4 | 0.03 | 38.9 | 7.28 | 080922 | 22.23 | 0.31 $\pm$ 0.034 | 0.24 $\pm$ 0.031 | 0.40 $\pm$ 0.024 | red
|  |  |  |  | 091122 | 22.10 | 0.47 $\pm$ 0.052 | 0.09 $\pm$ 0.068 | 0.35 $\pm$ 0.028 | blue
2015 VV1655 | L4 | 0.09 | 16.8 | 9.02 | 113021 | 23.32 | 0.87 $\pm$ 0.049 | 0.32 $\pm$ 0.055 | — | ultra-red
2015 VW1655 | L4 | 0.05 | 5.0 | 8.39 | 113021 | 22.89 | 0.45 $\pm$ 0.032 | 0.36 $\pm$ 0.048 | — | red
|  |  |  |  | 120222 | 22.93 | — | — | 0.61 $\pm$ 0.060 | indeterminate

##  3 Photometric Reduction

###  3.1 Calibration

To calibrate the photometry of our IMACS observations, we cross-matched the
in-frame background stars against PS1 sources (Magnier et al., 2013). We first
converted the PS1 griz photometry to the SDSS system using the transformation
equations in Tonry et al. (2012), and then selected the sources with $g-r$
between 0.25 and 2.0, $r-i$ between 0.0 and 0.8 as the reference sources. By
solving the equation below using the apparent magnitude of the reference
sources, we determined the photometric zeropoint of each frame:

$m_{sdss}=m_{ins}+2.5\log_{10}(\tau_{exp})+m_{0},$ (1)

where $m_{sdss}$ is the apparent magnitude of a specific band of the cross-
matched reference sources, $m_{ins}$ is the instrumental magnitude of that
specific band measured from the IMACS image, $\tau_{exp}$ is the exposure
time, and $m_{0}$ is the photometric zeropoint of that frame.

After we determined the zeropoints of each frame, we used every cross-matched
star in every frame to evaluate the linear color conversions between the IMACS
and SDSS photmetric system by solving the following equation:

$m_{M}=m_{sdss}+a~{}(g-r)_{sdss}+b,$ (2)

where $m_{M}$ and $m_{sdss}$ are the IMACS and SDSS magnitude, respectively,
and a, b are the coefficients of the linear conversion. The results are:

$\begin{split}g_{M}=g_{sdss}-0.078(g-r)_{sdss}+0.069\\\
r_{M}=r_{sdss}-0.024(g-r)_{sdss}+0.024\\\
r_{M}=r_{sdss}-0.038(r-i)_{sdss}+0.015\\\
i_{M}=i_{sdss}-0.188(r-i)_{sdss}+0.134\\\
z_{M}=z_{sdss}-0.026(g-r)_{sdss}+0.031\end{split}$ (3)

With the photometric zeropoints and the color conversion equations, we are
able to measure the griz colors of targets in SDSS photmetry system.

###  3.2 PSF Modeling

To accurately measure the flux and apparent magnitude of NTs, we select stars
around the target NT to model the local PSF. Several popular analytical
functions are considered for modeling the PSF, such as Moffat (Moffat, 1969)
and the sum of 2D Gaussians (Bendinelli et al., 1990). Both functions can
adequately model the “wing” of PSF. However, considering our PSF can be
asymmetric (not round, see Figure 1), we model the PSF by using the
superposition of n asymmetric 2D Gaussians. The flux of the PSF at any point
in the $(x^{\prime},y^{\prime})$ orthogonal coordinate system is:

$PSF(x^{\prime},y^{\prime})=b(x^{\prime},y^{\prime})+\sum_{i=1}^{n}\mathrm{A_{i}}~{}\Big{(}\texttt{exp}\big{[}-(\frac{x^{\prime
2}}{2\sigma^{2}_{x^{\prime}i}}+\frac{y^{\prime
2}}{2\sigma^{2}_{y^{\prime}i}})\big{]}\big{)},$ (4)

where $b(x^{\prime},y^{\prime})$ is the background flux at that point, n is a
small number, Ai is the amplitude of individual Gaussian,
$\sigma_{x^{\prime}i}$ and $\sigma_{y^{\prime}i}$ are the widths on
$x^{\prime}$ and $y^{\prime}$ axes of individual Gaussian, respectively. This
equation can be rotated to the image reference frame $(x,y)$ with a position
angle $\theta$ and translating the centroid to $(x_{0},y_{0})$ such that

$\begin{pmatrix}x^{\prime}\\\
y^{\prime}\end{pmatrix}=\begin{bmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{bmatrix}\begin{pmatrix}x-x_{0}\\\
y-y_{0}\end{pmatrix}.$

Therefore, the Gaussian functions share the same center, position angle, and
ellipticity but have unequal contribution and different width. To proper chose
‘n’, the number of Gaussians we should use, we calculate Bayesian information
criterion (BIC) for each n we use. The BIC is defined as:

$\mathrm{BIC}=-2\>\mathrm{ln}(\hat{\mathcal{L}})+k\>\mathrm{ln(m)},$ (5)

where $\hat{\mathcal{L}}$ is the maximum likelihood of the model, $k$ is the
number of parameters estimated by the model, and m is number of data points we
use to fit the model. The models with lower BIC values are generally
preferred, and it penalizes the model with larger $k$ automatically. Since the
multiple Gaussian PSF model can be linearized by taking logarithm and assuming
that the errors are normally-distributed, the $\hat{\mathcal{L}}$ is
equivalent to the least squares estimation. Thus, BIC can be written as a
function of error variance $\hat{\sigma_{e}^{2}}$:

$\mathrm{BIC}=m\>\mathrm{ln}(\hat{\sigma_{e}^{2}})+k\>\mathrm{ln(m)},$ (6)

In other words, the model with lower residual and fewer parameters is
preferred. We find that the model with n = 1, a single 2D Gaussian, always has
highest BIC. On the other hand, the models with n = 2 and n = 3 generally have
similar BICs, therefore we conclude that using any model with n $>$ 3 is
redundant.

Finally, we use the PSF model with n = 2 or 3, depending on which one has
lower BIC. Once all of the parameters are measured via modeling the stars, the
target NT can be modeled by refitting the center and amplitude of the PSF. The
flux is the sum of the final model. Figure 1 demonstrates that both the star
and the NT can be properly subtracted by the PSF model.

Figure 1: PSF modeling and subtraction. top-left: A star with the PSF model
contour. bottom-left: The image of NT. middle: The model of the star (top) and
the NT (bottom). right: the images after subtraction of the model.

###  3.3 Rotation Curve Correction

The observed magnitudes, and the resulting colors we are trying to measure,
are subject to rotational variations on the surface of these objects. To
approximately account for this, we use a model that exhibits a linear
variation in source brightness (its $r^{\prime}$-band magnitudes) and constant
$g^{\prime}-r^{\prime}$, $r^{\prime}-i^{\prime}$, $r^{\prime}-z^{\prime}$
colors (to convert each measurement to and $r^{\prime}$-band magnitude). This
model was then fit using a least-squares approach (see Fig. 2). The resulting
colors have been converted to SDSS magnitudes ($griz$; see Eq. 3), which are
reported in Table 1.

###  3.4 Reddening Line

Taken from Hainaut & Delsanti (2002), the reddening, or the spectral index,
can be expressed as the percent of reddening per 100 nm:

$S(\lambda_{1},\lambda_{2})=100*\frac{R(\lambda_{2})-R(\lambda_{1})}{(\lambda_{2}-\lambda_{1})/100}$
(7)

where $R(\lambda)$ is taken from Jewitt & Meech (1986):

$R(\lambda)=10^{-0.4(m(\lambda)-m_{\odot}(\lambda))}$ (8)

such that $m(\lambda)$ and $m_{\odot}(\lambda)$ are the magnitude of the
object and the Sun, respectively, at a particular wavelength, $\lambda$.
Setting the reddening line to pass through the color of the Sun (i.e. for
$S(\lambda_{1},\lambda_{2})$ = 0,
$m(\lambda_{1})-m(\lambda_{2})=m_{\odot}(\lambda_{1})-m_{\odot}(\lambda_{2})$),
we can derive the following equation, assuming
$m(\lambda_{1})=m_{\odot}(\lambda_{1})$:

$m(\lambda_{2})=-2.5log[1-10^{-4}S(\lambda_{1},\lambda_{2})(\lambda_{1}-\lambda_{2})]+m_{\odot}(\lambda_{2})$
(9)

Assuming $S(\lambda_{1},\lambda_{2})$ varies from -10% to 80%, we can plot the
reddening line for $g-r$ vs $r-i$ and $g-r$ vs $r-z$ in Fig. 3 and Fig. 4
respectively. Note that our targets generally fall along the reddening line,
as has been observed for small bodies in the outer Solar System previously
(Hainaut & Delsanti, 2002). Objects that fall above/below the reddening line
must exhibit emission/absorption lines at those particular wavelengths causing
them to deviate from a flat spectral index.

Figure 2: This figure shows our least-squares approach to fitting
$r^{\prime}$-band lightcurves and colors for an example NT target, 2013 VX30.
Each observation is shown as a colored point (blue downward triangle for
$g^{\prime}$, green square for $r^{\prime}$, yellow diamond for $i^{\prime}$,
and orange sideways triangle for $z^{\prime}$

. We then used a constant, but free parameter, color term to convert each
observation to an $r^{\prime}$-band observations; these point are shown as
black circles. The solid line is our least-squares fit to the
$r^{\prime}$-band (black) points. The dotted lines show our 1$\sigma$
deviation from this fit.

##  4 Results

###  4.1 Color-Color Results

In Fig. 3, we show the $g-r$ and $r-i$ colors measured for our NT targets.
Similar to the scattered TNOs, our targets exhibit wide range in this color
space; while most targets fall within the “red” zone (principal component
<1.75; see Sec. 4.2), there are three firm and one potential NTs in the
“ultra-red” zone (principal component >1.75). Of these objects, we identified
two new ultra-red NTs, 2013 TZ187 and 2015 VV165, which were also
independently found and reported in Bolin et al. (2023). The potential “ultra-
red” NT, 2011 SO277 has varying results from different observations (Jewitt
(2018), Lin et al. (2019), and this work); see more discussion of this object
in Sec. 4.4.

With extra ultra-red colored NTs, the red to ultra-red ratio for our sample is
3.75:1, or 7.75:1 for the entire known population. This ratio is much more
consistent with the dynamically excited KBO ratio of between 4-11 : 1 (Schwamb
et al., 2019). However, comparing these ratios is not sufficient to determine
if the NT and KBO population come from the same source distribution (see Sec.
4.2). We also show the kernal density estimations (KDEs) of g-r and r-i color
in Fig. 3. Unlike the results from previous works, which claimed that the NTs
and JTs have very similar color distributions, our new results show that the
KDEs of NTs are closer the the KDEs of scattered TNOs. Further analysis is
presented in Sec. 4.2.

Figure 3: Measured $g-r$ vs $r-i$ of the NT population. Blue points are colors
of scattered TNOs and orange triangles are JTs, both taken from the literature
(Hainaut et al., 2012). Light blue x’s are previously observed colors of NTs
which the ”Trojan Color Conundrum” was based on (Sheppard & Trujillo, 2006;
Sheppard, 2012; Parker et al., 2013; Jewitt, 2018), while the blue plus signs
are more recently observed NT colors which bring this conundrum into question
(Lin et al., 2019; Bolin et al., 2023). Targets observed in this paper are
shown as green squares. Solar color and the reddening line (see 3.4) are
depicted as an yellow star and orange dotted line respectively. Objects that
have multiple observations in this paper are connected by a dot-dashed line.
NTs that have been previously observed in the literature are connected by a
dashed line. The yellow line marks values where the PCA yields values equal to
our cutoff of 1.75 (see Fig. 6 and Sec. 4.2). Objects in the yellow region are
above this cutoff and considered ultra-red in this paper. The blue line marks
values where the PCA yields values equal to our cutoff of -1.25 (see Fig. 6
and Sec. 4.2). Objects in the blue region are blue this cutoff and considered
blue in this paper. The top and right inset plots show the kernel density
estimation (KDE) of the g-r and r-i distributions respectively of the included
sub-populations.

In Fig. 4, we show the $g-r$ and $r-z$ colors measured for our NT targets. All
of our targets are consistent with the scattered/hot TNO populations. This
result is expected as NTs are thought to have originated from
scattered/unstable TNOs. The physical cause of this $z$-band colorization of
the cold TNO population is not currently clear, but must be due to some
absorption feature around 900 nm based on the displacement from reddening
line. Spectroscopic information, such as will be taken with JWST (Markwardt et
al., 2021), will shed further light on chemical links between these
populations.

Figure 4: Measured $g-r$ vs $r-z$ of the NT population. Navy upward triangles,
green downward triangles, and blue circles are measurements taken from the
literature of TNOs (scattered, cold, and hot respectively) (Schwamb et al.,
2019). Teal plus signs are colors of NTs taken from the literature (Lin et
al., 2019). Targets observed in this paper are shown as orange squares. Solar
color and the reddening line (see 3.4) are depicted as an yellow star and
orange dotted line respectively. Objects that have observations taken in this
paper and from the literature are connected with a dashed line. Objects that
have multiple observations in this paper are connected by a dot-dashed line.
The green ellipse demarcates the region of color-color space occupied only by
cold TNOs. The top and right inset plots show the kernel density estimation
(KDE) of the g-r and r-z distributions respectively of the included sub-
populations.

###  4.2 Comparison to Previous Observations

All of the targets in our sample have previous observations (though not all
from the same survey). Therefore, we compare the difference between our
measurements and those from the literature to our computed errors, shown in
Fig. 5, to determine if there is any systematic offset in our observations. We
find that the observed differences in g-r are mostly within our observational
errors, meaning our observations are roughly consistent with previous
literature. However, previous observations are split between being slightly
systematically larger in r-i and systematically lower than our measurements.
Further investigation indicated that the larger group has smaller offset in
the order of 0.05, and the lower group has larger offset about -0.15. We also
find an instrument dependency on the groups; the smaller offset samples were
mostly measured with Gemini and Dark Energy Survey, which both have proper
photometry transformation equations to SDSS system. On the other hand, the
larger offset samples were mostly measured by using the R and I filters or
without proper photometry transformation equations. Therefore it is likely
that the different photometry systems mostly contribute such systematic
offsets. In every case, this did not change the result much on the following
Principal Component Analysis (PCA) analysis, since the g-r axis is the
dominant element on our Principal Component.

Figure 5: The differences in observed color between NTs in this paper and the
literature as compared to the average error on our observations. The
differences in g-r, r-i, and r-z observations are shown as blue, orange, and
green histograms respectively. The average g-r, r-i, and r-z errors are shown
as blue dotted, green dot-dashed, and orange dashed lines respectively

###  4.3 Comparison to Other Populations

The ultimate goal of this work is to determine how similar the NT colors are
to other populations in the Solar System. A simple statistical test to measure
the likelihood that two distributions are drawn from the same underlying
distribution is the Kolmogorov-Smirnov (KS) test (Darling, 1957). Although the
KS test can be generalized to more than a single dimension, the interpretation
becomes complicated. For simplicity, we reduce the dimension of our data and
use the traditional statistical test. Specifically, we performed a Principal
Component Analysis (PCA) of our data, using the scikit-learn python package
(Pedregosa et al., 2011). Fig. 6 demonstrates that the PCA is able to
successfully reduce the g-r vs. r-i color-color plot to a 1-D parameter that
still distinguishes between the red and ultra-red populations of TNOs and the
whole JT population (which is comprised of only red objects). The principal
component value (PC1) which separates these populations is 1.75 (shown as a
dotted line in Fig. 6). We use this definition to classify our NT targets as
red or ultra-red; the corresponding region in g-r vs r-i space is shown in
Fig. 3 as a yellow shaded region. We then applied this PCA model to other
populations in the Solar System, including JTs and previous observations of
NTs, the results of which are shown in Fig. 7. By eye, the JT population is
clearly unique in that it is nearly devoid of any ultra-red members (i.e.
targets with a PC1 >1.75). Also of note, about 25% of the NT targets presented
in this paper occupy a unique region of PC1 $\sim-1$. This region corresponds
to blue objects that are not frequently present in the outer Solar System
populations (see Sec. 4.4 for a more in-depth discussion of these objects).

Figure 6: The results of running a Principal Component Analysis (PCA) with the g-r and r-i colors of certain Solar System populations. The green histogram corresponds to the JTs (taken from Hainaut et al. (2012)). The blue and orange histograms correspond to the red and ultra-red subpopulations of the scattered TNOs, taken from Hainaut et al. (2012); the classification of red vs ultra-red was determined by using a clustering algorithm (DBSCAN; Pedregosa et al. (2011)) which separated the TNOs into two sub-populations. Figure 7: Cumulative distributions of the Principal Component (see Sec. 4.2) values of populations in the Solar System. The cut-off between red and ultra-red as defined by this PCA is shown as a black dashed line (see Fig. 6). The cut-off between red and blue objects is similarly shown as a dot-dashed line. The JT and scattered TNO results are shown as orange and navy histograms respectively. The NTs observations from previous literature are shown as a blue histogram. The NT observations from this work are shown as a green histogram.  KS Test P-value | NTs (This Work) | NTs (Pre-2019) | NTs (Post-2019) | TNOs | JTs
---|---|---|---|---|---
NTs (This Work) | 1 | 0.020 | 0.61 | 0.56 | 0.003
NTs (Pre-2019) | 0.020 | 1 | 0.15 | 0.03 | 0.27
NTs (Post-2019) | 0.61 | 0.15 | 1 | 0.14 | 0.05
TNOs | 0.56 | 0.03 | 0.14 | 1 | 0.0002
JTs | 0.003 | 0.27 | 0.05 | 0.0002 | 1
Table 2: The resulting p-values of the KS Test on each combination of sub-
populations considered in this work.

We then ran a KS test for each combination of these Solar System populations
to determine the likelihood that they came from the same underlying
distribution; the results of these tests are recorded in Table 2. We conclude
that the compared populations are from different distributions if they have a
p-value of $\leq$ 0.05, corresponding to a 95% confidence level to reject the
null hypothesis. Therefore, we find that the population observed in this work
is not consistent with being drawn from the same distribution as the JTs, but
is instead more consistent with the TNO population. This result is the
opposite of what was found pre-2019, where the NTs were more consistent with
the JT population. The results from post-2019 data also show that the NT
population is more consistent with the TNO population, but this work shores up
this result significantly. Further observations of members of the NT
population in particular could also increase the statistical significance of
this result. However, we feel confident in claiming that our results show NTs
and TNOs are consistent with coming from the same underlying distribution
based on their optical colors with the greatest confidence to date.

###  4.4 Color-Absolute Magnitude Relations

In Fig. 8, we plot the Principal Component for our targets as a function of
absolute magnitude (H). We look for any significant clustering or correlations
in these plots which would indicate that the color classification of NTs is
dependent on their size.

To search for clustering in our datasets, we run a Mean Shift Clustering
algorithm (Pedregosa et al., 2011), which does not need a number of clusters
as an input parameter (just a bandwith which can be initialized with the
estimate_bandwith function). To test the significance of clustering we
calculate the Cluster Index. The Cluster Index from the SigClust evaluation
tool is defined as (Ferland et al., 2013):

$CI=\frac{\sum_{k=1}^{N}\sum_{i\in
C_{k}}\parallel\boldsymbol{x_{i}}-\boldsymbol{\bar{x}}^{(k)}\parallel^{2}}{\sum_{i=1}^{n}\parallel\boldsymbol{x_{i}}-\boldsymbol{\bar{x}}\parallel^{2}}$
(10)

where $\boldsymbol{\bar{x}}^{(k)}$ represents the mean of the kth cluster for
k = 1, 2, … N for N clusters and $\boldsymbol{\bar{x}}$ represents the overall
mean. The CI provides a p-value for the significance of the cluster between
these two clusters. To test if our data was correlated, we used the Pearson
Correlation Coefficient (Kirch, 2008) which is defined as:

$r=\frac{\sum(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum(x_{i}-\bar{x})^{2}\sum(y_{i}-\bar{y})^{2}}}$
(11)

where $x_{i}$ and $y_{i}$ are the data points and $\bar{x}$ and $\bar{y}$ are
the respective means. We calculated each of these values for all of the plots
shown in Fig. 8. To determine whether or not these values could be obtained
from random noise we generated 1000 sets of points with the same number of
objects as our observation within the same region of Principal Component vs H
space and ran the same analysis on those sets. These results are shown in the
inset histograms in Fig. 8.

We found that the cluster is consistent with random noise and should not be
considered significant. This result also suggest that the color of NTs are
distributed continuously from blue to ultra-red rather than bimodally. The
positive correlation with size is intriguing and may point to primordial
differences in objects of different sizes in the outer Solar System. However,
H is not a direct correlation to size as the object’s albedo must be taken
into account. Such observations do not currently exist for the NT population
and will be necessary to establish a color-size correlation. Indeed,
photometric observations of the rest of the NT population are necessary to
confirm this slight correlation.

Figure 8: NT colors as a function of absolute magnitude. Grey points are taken
from the literature (Sheppard & Trujillo, 2006; Sheppard, 2012; Parker et al.,
2013; Jewitt, 2018; Lin et al., 2019; Schwamb et al., 2019). Colored squares
were measured in this paper. Duplicate observations of the same object are
connected by dashed lines. The inset plots contain histograms of the Cluster
Indices and Pearson Correlation Coefficients of a random distribution colors
and absolute magnitude (see Sec. 4.1). Each grey dashed line in the inset
plots shows the corresponding value calculated for the observed distribution.

###  4.5 Unique Targets

While most of our targets are consistent with previous color measurements, one
object, 2011 SO277 is classified here as ultra-red while its previous
observations place it firmly with in the red zone. Based on our other
observations, we consider our results to be roughly consistent with previous
literature (see Fig. 5), so this result is indeed unexpected. One explanation
as to why this object has such different colors in independent observations is
that its surface is not homogeneous. To test this hypothesis, a more in depth
study of the rotational properties of the surface of this object is necessary,
which will be upcoming in our next work on the lightcurves of NTs.

Three of our targets, 2014 SC375, 2014 YB92, and 2015 VU207, are much bluer,
nearly solar in color, as compared to the other NTs or KBOs. Bolin et al.
(2023) also reported 2014 YB92 and 2015 VU207 have blue, near solar color. In
fact, these objects are as blue as the blue B/C-type asteroids, such as 3200
Phaetheon (Tabeshian et al., 2019; Lisse & Steckloff, 2022). A similarly blue
TNO has been observed, which appears to be covered ferric oxides and
phyllosilicates (Seccull et al., 2018). This TNO has a highly eccentric and
inclined orbit, suggesting it may have a common origin with C-type asteroids
and has since been implanted into trans-Neptunian space. It is possible that
these NTs originated elsewhere in the Solar System, but their current orbits
are stable for > Gyrs (see Sec. 2), implying that they were captured just
after Neptune’s migration. However, based on these results the blue ratio for
NTs is currently much higher than that of the TNO population. This result may
suggest that inner Solar System material may be more efficiently transferred
to NT orbits which have a smaller perihelion than the Kuiper Belt. Future
spectral observations would be necessary to reveal any compositional
differences this target may have as compared to the rest of the NT population.

##  5 Why were the ultra-red NTs rare before 2019?

Prior to 2019, the ultra-red NTs were very rare; none of the 13 NT samples in
Jewitt (2018) are ultra-red NTs, which led to the claim of a “Trojan Color
Conundrum”. Here we propose two possibilities to explain this inconsistency:

  1. 1.

Small number statistics: Small number statistics could generate such a
surprising result. If we assume a 7.75:1 apparent red to ultra-red ratio of
NTs, the chance to randomly select 13 objects without picking up any ultra-red
one is about 18%, which is very likely. If we use a 3.75:1 apparent red to
ultra-red ratio, the chance is now 0.5%. While it is not impossible, we may
also consider alternative explanations.

  2. 2.

Selection effect: Since bigger objects are easier to detect and obtain color
measurements for, the 13 objects in Jewitt (2018) trend to be large; 10 of 13
have H $\leqslant 8$. Moreover, many NTs have been discovered by deeper (Lin
et al., 2021) or wider (Bernardinelli et al., 2022, 2020; Lin et al., 2019)
surveys since 2018, which included many high-inclination objects. Thus the
Jewitt (2018) sample appears to be biased toward bigger sized and lower
inclination objects. In fact, 8 of 13 NTs in the Jewitt (2018) sample have
orbital inclination $<10^{\circ}$; 9 of the 31 currently known NTs have
inclination $<10^{\circ}$, meaning that 8 of the 9 total low-inclination NTs
were included in the Jewitt (2018) sample. Such objects has very similar red
color (see Figure 8). Therefore, the possible color-orbit-size correlation in
NT population could be at least partially explain why the “Trojan Color
Conundrum” was observed, especially when there were some selection biases in
that sample.

##  6 Conclusions

In this paper, we measure the griz colors for 15 of the 24 known NTs. We used
the IMACS f/4 instrument on the 6.5m Baade telescope with Sloan g’r’i’z’
filters to conduct our photometric survey. We confirm that 2013 VX30 is ultra-
red in color, and identify three NTs as ultra-red. This result brings the red
to ultra-red ratio of NTs to 7.75:1, much more consistent with the
corresponding TNO ratio and resolving the “Trojan Color Conundrum”. Moreover,
the color distribution of NTs is now indistinguishable from the scattered
population of TNOs and different from the Jovian Trojans. We also find three
targets which have solar color, the origin of which is unclear; the most
likely explanation is that these objects originated from the inner solar
system. For the entire NT population, we find that color of NTs may correlated
to their absolute magnitude, and the objects with larger H trend to have
redder color. The explanation behind this correlation remains an open question
that is difficult to address with current data. More discoveries of NTs
(especially around L5) are clearly needed. The L5 point has historically been
difficult to study due to its overlap with the galactic plane, but the NT L5
region is moving away from this high stellar density region, making now the
perfect time to start studying this population. The true degree of asymmetry
between the L4 and L5 clouds will be an important to distinguishing different
formation scenarios for the NT population. Moreover, our ongoing work to
measure the rotational period and specific composition of these small bodies
directly will be vital to understanding the true origin of the NT population.

This paper includes data gathered with the 6.5 meter Magellan Telescopes
located at Las Campanas Observatory, Chile. This material is based upon work
supported by the National Aeronautics and Space Administration under grant No.
NNX17AF21G issued through the SSO Planetary Astronomy Program and by the
National Science Foundation under grant No. AST-2009096. This research was
supported in part through computational resources and services provided by
Advanced Research Computing at the University of Michigan, Ann Arbor.

## References

  * Almeida et al. (2009) Almeida, A. J. C., Peixinho, N., & Correia, A. C. M. 2009, A&A, 508, 1021, doi: 10.1051/0004-6361/200911943
  * Bendinelli et al. (1990) Bendinelli, O., Parmeggiani, G., Zavatti, F., & Djorgovski, S. 1990, AJ, 99, 774, doi: 10.1086/115373
  * Bernardinelli et al. (2020) Bernardinelli, P. H., Bernstein, G. M., Sako, M., et al. 2020, ApJS, 247, 32, doi: 10.3847/1538-4365/ab6bd8
  * Bernardinelli et al. (2022) —. 2022, ApJS, 258, 41, doi: 10.3847/1538-4365/ac3914
  * Bolin et al. (2023) Bolin, B. T., Fremling, C., Morbidelli, A., et al. 2023, MNRAS, 521, L29, doi: 10.1093/mnrasl/slad018
  * Ćuk et al. (2012) Ćuk, M., Hamilton, D. P., & Holman, M. J. 2012, MNRAS, 426, 3051, doi: 10.1111/j.1365-2966.2012.21964.x
  * Darling (1957) Darling, D. A. 1957, The Annals of Mathematical Statistics, 28, 823.  http://www.jstor.org/stable/2237048
  * DeMeo & Carry (2014) DeMeo, F. E., & Carry, B. 2014, Nature, 505, 629, doi: 10.1038/nature12908
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137.  https://arxiv.org/abs/1302.4485
  * Fernandez & Ip (1984) Fernandez, J. A., & Ip, W. H. 1984, Icarus, 58, 109, doi: 10.1016/0019-1035(84)90101-5
  * Gomes & Nesvorný (2016) Gomes, R., & Nesvorný, D. 2016, A&A, 592, A146, doi: 10.1051/0004-6361/201527757
  * Hahn & Malhotra (1999) Hahn, J. M., & Malhotra, R. 1999, AJ, 117, 3041, doi: 10.1086/300891
  * Hainaut et al. (2012) Hainaut, O. R., Boehnhardt, H., & Protopapa, S. 2012, A&A, 546, A115, doi: 10.1051/0004-6361/201219566
  * Hainaut & Delsanti (2002) Hainaut, O. R., & Delsanti, A. C. 2002, A&A, 389, 641, doi: 10.1051/0004-6361:20020431
  * Hasegawa et al. (2021) Hasegawa, S., Marsset, M., DeMeo, F. E., et al. 2021, ApJ, 916, L6, doi: 10.3847/2041-8213/ac0f05
  * Horner & Lykawka (2010) Horner, J., & Lykawka, P. S. 2010, MNRAS, 402, 13, doi: 10.1111/j.1365-2966.2009.15702.x
  * Jewitt (2018) Jewitt, D. 2018, AJ, 155, 56, doi: 10.3847/1538-3881/aaa1a4
  * Jewitt & Meech (1986) Jewitt, D., & Meech, K. J. 1986, ApJ, 310, 937, doi: 10.1086/164745
  * Jewitt (2002) Jewitt, D. C. 2002, AJ, 123, 1039, doi: 10.1086/338692
  * Kirch (2008) Kirch, W., ed. 2008, Pearson’s Correlation Coefficient (Dordrecht: Springer Netherlands), 1090–1091, doi: 10.1007/978-1-4020-5614-7_2569
  * Kortenkamp et al. (2004) Kortenkamp, S. J., Malhotra, R., & Michtchenko, T. 2004, Icarus, 167, 347, doi: 10.1016/j.icarus.2003.09.021
  * Lacerda et al. (2014) Lacerda, P., Fornasier, S., Lellouch, E., et al. 2014, ApJ, 793, L2, doi: 10.1088/2041-8205/793/1/L2
  * Lin et al. (2022) Lin, H.-W., Markwardt, L., Napier, K. J., Adams, F. C., & Gerdes, D. W. 2022, Research Notes of the American Astronomical Society, 6, 79, doi: 10.3847/2515-5172/ac6752
  * Lin et al. (2019) Lin, H. W., Gerdes, D. W., Hamilton, S. J., et al. 2019, Icarus, 321, 426, doi: 10.1016/j.icarus.2018.12.006
  * Lin et al. (2021) Lin, H. W., Chen, Y.-T., Volk, K., et al. 2021, Icarus, 361, 114391, doi: 10.1016/j.icarus.2021.114391
  * Lisse & Steckloff (2022) Lisse, C., & Steckloff, J. 2022, Icarus, 381, 114995, doi: https://doi.org/10.1016/j.icarus.2022.114995
  * Luu & Jewitt (1996) Luu, J., & Jewitt, D. 1996, AJ, 112, 2310, doi: 10.1086/118184
  * Lykawka et al. (2011) Lykawka, P. S., Horner, J., Jones, B. W., & Mukai, T. 2011, MNRAS, 412, 537, doi: 10.1111/j.1365-2966.2010.17936.x
  * Magnier et al. (2013) Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20, doi: 10.1088/0067-0049/205/2/20
  * Malhotra (1993) Malhotra, R. 1993, Nature, 365, 819, doi: 10.1038/365819a0
  * Malhotra (1995) —. 1995, AJ, 110, 420, doi: 10.1086/117532
  * Markwardt et al. (2021) Markwardt, L., Adams, F., Gerdes, D., et al. 2021, The First Near-IR Spectroscopic Survey of Neptune Trojans, JWST Proposal. Cycle 1, ID. #2550
  * Moffat (1969) Moffat, A. F. J. 1969, A&A, 3, 455
  * Nesvorný et al. (2018) Nesvorný, D., Vokrouhlický, D., Bottke, W. F., & Levison, H. F. 2018, Nature Astronomy, 2, 878, doi: 10.1038/s41550-018-0564-3
  * Nesvorný et al. (2013) Nesvorný, D., Vokrouhlický, D., & Morbidelli, A. 2013, ApJ, 768, 45, doi: 10.1088/0004-637X/768/1/45
  * Neveu & Vernazza (2019) Neveu, M., & Vernazza, P. 2019, ApJ, 875, 30, doi: 10.3847/1538-4357/ab0d87
  * Parker et al. (2013) Parker, A. H., Buie, M. W., Osip, D. J., et al. 2013, AJ, 145, 96, doi: 10.1088/0004-6256/145/4/96
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Peixinho et al. (2015) Peixinho, N., Delsanti, A., & Doressoundiram, A. 2015, A&A, 577, A35, doi: 10.1051/0004-6361/201425436
  * Peixinho et al. (2012) Peixinho, N., Delsanti, A., Guilbert-Lepoutre, A., Gafeira, R., & Lacerda, P. 2012, A&A, 546, A86, doi: 10.1051/0004-6361/201219057
  * Pike et al. (2017a) Pike, R. E., Lawler, S., Brasser, R., et al. 2017a, AJ, 153, 127, doi: 10.3847/1538-3881/aa5be9
  * Pike et al. (2017b) Pike, R. E., Fraser, W. C., Schwamb, M. E., et al. 2017b, AJ, 154, 101, doi: 10.3847/1538-3881/aa83b1
  * Roig & Nesvorný (2015) Roig, F., & Nesvorný, D. 2015, AJ, 150, 186, doi: 10.1088/0004-6256/150/6/186
  * Schwamb et al. (2019) Schwamb, M. E., Fraser, W. C., Bannister, M. T., et al. 2019, ApJS, 243, 12, doi: 10.3847/1538-4365/ab2194
  * Schwarz et al. (2011) Schwarz, G. J., Ness, J.-U., Osborne, J. P., et al. 2011, ApJS, 197, 31, doi: 10.1088/0067-0049/197/2/31
  * Seccull et al. (2018) Seccull, T., Fraser, W. C., Puzia, T. H., Brown, M. E., & Schönebeck, F. 2018, ApJ, 855, L26, doi: 10.3847/2041-8213/aab3dc
  * Sheppard (2010) Sheppard, S. S. 2010, AJ, 139, 1394, doi: 10.1088/0004-6256/139/4/1394
  * Sheppard (2012) —. 2012, AJ, 144, 169, doi: 10.1088/0004-6256/144/6/169
  * Sheppard & Trujillo (2006) Sheppard, S. S., & Trujillo, C. A. 2006, Science, 313, 511, doi: 10.1126/science.1127173
  * Tabeshian et al. (2019) Tabeshian, M., Wiegert, P., Ye, Q., et al. 2019, The Astronomical Journal, 158, 30, doi: 10.3847/1538-3881/ab245d
  * Tonry et al. (2012) Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99, doi: 10.1088/0004-637X/750/2/99
  * Wong & Brown (2017) Wong, I., & Brown, M. E. 2017, AJ, 153, 145, doi: 10.3847/1538-3881/aa60c3