Summary
该协议从系外行星的光曲线中提取信息，并构造其表面图。它使用地球的光曲线，作为代理系行星，来演示这种方法。
Abstract
从单点观测中对系外行星要素进行空间解析对于评估系外行星的潜在可居住性至关重要。该协议的最终目标是确定这些行星世界是否包含地质特征和/或气候系统。我们提出了一种从多波长单点光曲线中提取信息并检索表面贴图的方法。它使用奇异值分解 （SVD） 来分离导致光曲线变化的源，并推断部分多云气候系统的存在。通过分析从SVD获得的时间序列，可以推断出主要组件（PC）的物理归属，而无需任何光谱属性的假设。结合查看几何体，如果发现其中一台 PC 包含曲面信息，则重建曲面贴图是可行的。退化源于像素几何体的卷积，频谱信息决定了重建表面图的质量，这就要求引入正则化。为了演示该协议，分析了作为代理系行星的地球多波长光曲线。提出了结果与地面真实性的比较，以表明协议的性能和局限性。这项工作为系外行星应用的未来泛化提供了基准。
Introduction
确定宜居世界是天体生物学1的最终目标之一。自第一次探测2以来，4000多颗系外行星已被确认为3号行星，其中有许多地球类比（例如，TRAPPIST-1e）4。这些行星的轨道和行星性质与地球相似，因此具有潜在的可居住性。在这种情况下，从有限的观察中评估它们的可居住性至关重要。基于对地球生命知识的了解，地质和气候系统对可居住性至关重要，因此，可居住性可以作为生物特征。原则上，即使一个行星在空间上不能比单个点更好地解析，也可以从远处观测到这些系统的特征。在这种情况下，在评估系外行星的可居住性时，从单点光曲线确定地质特征和气候系统至关重要。这些系外行星的表面测绘变得紧迫。
尽管查看几何体和光谱要素之间有卷积，但系外行星表面的信息仍包含在其时间解析的单点光曲线中，该曲线可从远处获取，并获得足够的观测结果。然而，由于云层的影响，可能适合地球的系外行星的二维（2D）表面测绘具有挑战性。检索二D地图的方法已经开发和测试使用模拟光曲线和已知的光谱5，6，7，8，但这些方法尚未应用于实际观测。此外，在分析现在和近期的系外行星观测时，当行星表面成分没有受到良好约束时，特征光谱的假设可能会引起争议。
本文演示了一种地球外行星的表面测绘技术。我们使用 SVD 来评估和分离来自多波长光曲线中包含的不同来源的信息，而无需对任何特定光谱进行假设。结合查看几何，我们使用及时解析但空间复杂的表面信息来呈现曲面图的重建。为了演示这种方法，分析了深空气候观测站/地球多色成像相机（DSCOVR/EPIC;www.nesdis.noaa.gov/DSCOVR/spacecraft.html）获得的两年多波长单点地球观测结果。我们使用地球作为代理系系外行星来评估这种方法，因为目前对系外行星的观测是不够的。我们以纸张附加代码为例。它是在 python 3.7 下开发的，具有 anaconda 和 healpy 封装，但协议的数学也可以在其他编程环境（例如 IDL 或 MATLAB）中完成。
Subscription Required. Please recommend JoVE to your librarian.
Protocol
1. 编程设置
- 为附加代码设置编程环境。需要一台具有 Linux 操作系统的计算机，因为 Windows 上没有修复包。代码的计算费用并不昂贵，因此普通个人计算机可以处理协议。
- 按照说明 （https://docs.anaconda.com/anaconda/install/linux/） 将 Anaconda 与 Python 3.7 安装到系统中，然后在终端中使用以下命令来设置编程环境：
$ conda 创建 - 命名 myenv python = 3.7
$ conda 激活 myenv
$ 康达安装阿纳康达
$ conda 安装愈合
注：这些步骤可能需要数十分钟，具体取决于硬件和 Internet 速度。前两个命令行中的环境名称"myenv"可以更改为任何其他字符串。
2. 从观测中获取多波长光曲线和查看几何图形
- 在查看几何体中，包括子恒星的经度和纬度以及每个相应时间帧的子观察点。
若要使用以下附加代码，请确保这两个文件的格式与LightCurve、.csv几何.csv。 - 运行 PlotTimeSeries.py 可视化数据并检查其质量。将创建 两.png 光 曲线.png 几何学（补充图1-2）。如果应用于不同的观测值，可能需要调整此参数和以下绘图代码中的参数。
$ python PlotTimeSeries.py曲线
$ python PlotTimeSeries.py几何
3. 从光曲线中提取表面信息
- 将系外行星的多波长反照率光曲线居中，并在每个波长上通过相应的标准偏差使其规范化。这将导致每个通道的同等重要性。
其中R't，k和R t，k分别是t-th时间步长和k-th波长的缩放和观测反照率;μk和 σk是 k-th 波长反照率时间序列的均值和标准偏差。- 运行Normalize.py使光曲线（Rt，k）正常化。输出保存在规范化LightCurve.csv。
$ 巨蛇Normalize.py
- 运行Normalize.py使光曲线（Rt，k）正常化。输出保存在规范化LightCurve.csv。
- 运行PlotTimeSeries.py以可视化规范化的光曲线。将创建一个.png光曲线（补充图3）。
$ python PlotTimeSeries.py归一化光曲线 - 在缩放的图照率光曲线上应用 SVD 以查找占统治地位的 PC 及其相应的时间序列。
在左侧，T 和 K 是时间步长和观测波长的总数; R' 是缩放反照率观测的矩阵，其 （t，k） - th 元素为 R't，k。在右侧，V 列是 PC， 正畸矢量定义空间 SVD 项目到; • 是 一个对角线矩阵，其 （k，k）-th 元素是沿 V 的 k-th 列定义的 k-th 轴缩放光曲线的标准 偏差。U 列 是 V 中每台 PC 的相应时间 序列。- 运行SingularValueDecomposition.py分解 R'。生成的U、 +， VT分别保存在输出文件中 U.csv、单.csv和V_T.csv中。
$ 巨蛇SingularValueDecomposition.py
- 运行SingularValueDecomposition.py分解 R'。生成的U、 +， VT分别保存在输出文件中 U.csv、单.csv和V_T.csv中。
- 使用PlotTimeSeries.py PlotSVD.py可视化SVD 结果。将创建三.png U.png、V_T.png和数字（补充图 4-6）。
$ 巨PlotTimeSeries.py U
$ 巨蛇PlotSVD.py - 分析 PC 的贡献和相应的时间序列，以确定包含表面信息的 PC。
- 比较对角线 + 的 奇异值。一颗类似地球的部分多云的系外行星预计将有两个可比较的占优势奇异值。
注： + 可能 包含较少或超过两个占主导地位的奇异值，如下文讨论。 - 比较两种主流 PC 的时间序列模式。包含表面信息的 PC 往往具有比另一个更常规的形状。由于纵向不对称和表面在连续两天内发生小变化的再现，相应的时间序列往往具有大致恒定的每日变化。
- 使用 Lomb-Scargle 周期图9、10计算两台主流 PC 的周期，以确认 PC 的选择。包含表面信息的 PC 往往具有与功率密度频谱中旋转周期对应的较高峰值。
- 运行Periodogram.py，获取每台 PC 的时间序列的功率光谱。功率光谱保存在周期图中.csv。
$ 巨蛇Periodogram.py - 运行 PlotPeriodogram.py 以可视化这些周期图并确认 PC 的选择。将创建 图.png 图（补充图 7）。当前绘图代码以虚线添加表示年、半年、日间和半日周期的虚线，供参考，当应用于其他观测值时，可能需要更改这些虚线。
$ 巨蛇PlotPeriodogram.py - 选择包含曲面信息及其时间序列 u j的 PC vj。
其中 V[：，j] 和 U[：，j] 分别是 V 和 U 的 j-th列;j 是步骤 3.3 中推断的包含表面信息的 PC 索引。
- 比较对角线 + 的 奇异值。一颗类似地球的部分多云的系外行星预计将有两个可比较的占优势奇异值。
4. 构建行星表面图
- 使用分层等值面积等值度像素化 （HEALPix）11 方法对检索映射进行像素化。它将行星的球面划分为具有相同面积和均匀分布的像素。将 p-th 像素的未知值表示为 xp。
- 运行 HEALPixRandom.py 可视化像素化方法。将创建 一个.png 的"HEALPixRandom"（补充图8）。可以更改第 17 行的参数 N 侧，以换不同的分辨率。此步骤可能需要几秒钟到几分钟的时间，具体取决于分辨率。
$ 巨蛇HEALPixRandom.py
- 运行 HEALPixRandom.py 可视化像素化方法。将创建 一个.png 的"HEALPixRandom"（补充图8）。可以更改第 17 行的参数 N 侧，以换不同的分辨率。此步骤可能需要几秒钟到几分钟的时间，具体取决于分辨率。
- 使用观察几何计算在 t-th 时间步长 wt，p 的观测值中 p-th像素的重量。
其中α、t、p、β、p是太阳和 航天器天顶角，在 t-th 时间步长时为 p-th 像素;ct 是 t-th 观测值的规范化术语，因此每个时间步数的总重量总和是统一的。
注：假定几何学在此步骤中已知，也可以从其他分析派生，如下文讨论。- 运行ComputeWeight.py计算 wt，p.更改第 23行的 N 侧的值，以换取检索到的地图的其他分辨率。输出由于其大小而保存为 W.npz。
$ 巨蛇ComputeWeight.py
- 运行ComputeWeight.py计算 wt，p.更改第 23行的 N 侧的值，以换取检索到的地图的其他分辨率。输出由于其大小而保存为 W.npz。
- 使用 PlotWeight.py 可视化这些权重。将在一个文件夹"权重"中创建多个数字，每个时间步数一 个。将它们合并结果为 补充视频 1，它显示了每个像素的重量如何随时间变化。由于可视化效果大量，此步骤可能需要数小时才能完成。
$ 巨蛇PlotWeight.py - 组合几何和观测结果，以达到线性回归问题。
其中 P 是检索像素的总数; W 是权重矩阵，w t、p 为 （t，p）-th 元素; x 由 xp 作为 p-th 元素组成，这是在此问题中要解决的数量。
用L-2规范的正则化解决线性回归问题。
其中 I 是标识矩阵，+ 是正则化参数。
注：当T=10 4 和 P=3*103 时，10 -3是 Α 的一个良好值。应通过比较正方形误差 e 中的两个术语的值来调整它们，如下所示。
- 运行LinearRegression.py解决此线性回归问题。x 的结果保存在文件像素值.csv。更改行 16 处的 α 值，以换行化的不同强度。
$ 巨蛇LinearRegression.py
- 运行LinearRegression.py解决此线性回归问题。x 的结果保存在文件像素值.csv。更改行 16 处的 α 值，以换行化的不同强度。
- 根据 HEALPix 的映射规则将 x 转换为 2D 曲面贴图。
- 运行PlotMap.py使用不同的正则参数构造检索到的映射。将用Map_设置Map_-2.png、Map_-3.png和Map_-4.png三个数字（补充图9）。 HEALPix 文档 11 中描述了像素索引与其在地图上的位置之间的关系。此步骤需要数十秒。
$ 巨蛇PolotMap.py
- 运行PlotMap.py使用不同的正则参数构造检索到的映射。将用Map_设置Map_-2.png、Map_-3.png和Map_-4.png三个数字（补充图9）。 HEALPix 文档 11 中描述了像素索引与其在地图上的位置之间的关系。此步骤需要数十秒。
5. 估计检索到的地图的不确定性
- 将步骤 4.3 的线性回归问题重写为 x 的"真实值"，将观测噪声ε 。
- 假设 ε 遵循高斯分布 N （0，σ2I[T*T]）并估计其协方差。T-P 是检索到的地图固定时 uj 从观测的自由度。
- 在步骤 4.4 和 5.1 中合并方程。它导致 x 的高斯矢 量。
- 计算 x 的期望值和协方差 矩阵。
- 获取 x 中每个元素的不确定性，作为Cov = x + 对角线上相应元素的平方根。
其中 ep 是 xp 的不确定性; Diag=Cov=x=p 是 Cov = x = 对角线上的 p-th元素。 - 运行 Covariance.py计算 x 的协方差 矩阵。由于其大小，结果 保存在协方差.npz中。此步骤需要几十秒到几分钟，具体取决于 W 的大小。
$ 巨 蛇Covariance.py
- 假设 ε 遵循高斯分布 N （0，σ2I[T*T]）并估计其协方差。T-P 是检索到的地图固定时 uj 从观测的自由度。
- 根据 HEALPix 的映射规则将 ep 转换为检索到的 2D 地图。
- 运行PlotCovariance.py以可视化Cov[x]，然后将不确定性 ep映射到检索到的地图。将创建两个数字协.png和.png（补充图10-11）。
$ 巨蛇PlotCovariance.py
- 运行PlotCovariance.py以可视化Cov[x]，然后将不确定性 ep映射到检索到的地图。将创建两个数字协.png和.png（补充图10-11）。
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
我们使用地球的多波长单点光曲线来演示协议，并将结果与地面真相进行比较，以评估表面测绘的质量。此处使用的观测结果由DSCOVR/EPIC获得，这是一颗位于地球和太阳之间的第一个拉格朗日点（L1）附近的卫星，拍摄地球太阳表面十个波长的图像。两年（2016年和2017年）的观测值用于此演示，与江等人（2018年）12 和范等人（2019年）13的观测结果相同，其中提供了有关观测的更多细节。图 9：27 UTC 的样本观测值，2017 年 2 月 8 日，如图 1 所示。地球的图像被整合到单个点，以模拟外星人，遥远的观察者获得的光曲线观测，谁不能比一个像素更好的空间解析地球。因此，生成具有 ±10，000 个时间步长的多波长单点系外行星光曲线，这是该协议的输入数据。
按照步骤 3，我们发现多波长光曲线中的两台主要 PC，第二台 PC （PC2） 包含表面信息。PC2 的时间序列派生为步骤 3.5，显示更规律的形态，每天变化大致恒定，其功率谱显示比第一台 PC 更强的日周期（PC1， 图 2）。因此，此代理系外行星的表面图按照步骤4（图3a）进行构建，该图由每个像素的PC2值组成。与地球的地面真相（图3b）相比，重建后的地图恢复了所有主要大陆，尽管在南半球存在一些分歧，云层部分地阻止观测地表信息。根据步骤 5 （图3c） 获得的每个像素值的不确定性在检索地图中为 10% 左右，表明曲面映射质量良好，结果为正数。
图1：地球阳光半球的反射图像。
观测结果由DSCOVR/EPIC在10个波长和9：27 UTC，2017年2月8日进行。 请单击此处查看此图的较大版本。
图2：两种主流PC的时间序列和功率光谱。
（a） PC1 的时间序列。每日最大值和最小值用黑线表示。（b） PC1 时间序列的功率频谱。年度、半年、日间和半日周期表示为黑色虚线。（c） 和 （d） 与 （ a ）和（ b ）相同， 但与 PC2 相对应.这个数字取自范等人（2019）13。请单击此处查看此图的较大版本。
图3：地球表面的重建。
（a） 地球表面图，作为代理系行星，由多波长光曲线重建。地图中的颜色是每个像素的 PC2 值。中值的轮廓表示为黑线。（b） 地球表面图的地面真相。（c） （ a ） 中所示的重建地图的不确定性。这个数字修改自范等人 （2019）13。 请单击此处查看此图的较大版本。
图4：找到最佳正则化参数的结果。
当重建的 + 2（实线）达到其最小值时，正则化参数 = 的最佳值为10 -3.153（虚线）。请单击此处查看此图的较大版本。
图5：观测噪声的灵敏度测试。
（a） 从具有不同信噪比（S/Ns）的观测结果得出的视场中PC2与陆地分数（实线）之间的相关系数。与无噪声的光曲线的原始相关性显示为虚线。（b） 每台 PC 在具有不同观测值的 S/N 得出的土地分数上的重要性。重要性是使用梯度提升回归树 （GBRT） 模型计算的，如 Fan 等人 （2019）13 中所述。 请单击此处查看此图的较大版本。
补充图1：地球十波长反射的时间序列。请点击这里下载此图。
补充图2：（a）子观察点纬度的时间序列。（b） 与 （a） 相同，但经度相同。（c） 和 （d） 与 （a） 和 （b） 相同，但对应于子恒星点。请点击这里下载此图。
补充图3：地球十波长的规范化反射时间序列。请点击这里下载此图。
补充图4：十台 PC 的时间序列，U.请点击这里下载此图。
补充图 5：对应于每个 PC 的奇异值，对角线元素的 α。请点击这里下载此图。
补充图6：10台PC的规范化反射光谱，V列。请点击这里下载此图。
补充图7：10台PC时间序列的功率密度光谱。请点击这里下载此图。
补充图 8：检索地图的像素化，填充随机像素值。请点击这里下载此图。
补充图9：使用（a）10-2、（b）10-3和 （c）10-4的不同正则化参数重建的地球地图。请点击这里下载此图。
补充图10：x的协方差矩阵。请点击这里下载此图。
补充图11：x协方差矩阵对角线元素的平方根，映射到检索到的曲面贴图上。请点击这里下载此图。
视频 S1：2016 和 2017 年每个时间帧的观测值像素权重。
补充文件。请点击这里下载这些文件。
Subscription Required. Please recommend JoVE to your librarian.
Discussion
该协议的一个关键要求是从光曲线中提取表面信息的可行性，这取决于云的覆盖范围。在步骤 3.5.1 中，在系外行星之间，PC 的相对值可能不同。就地球而言，前两台PC主导光曲线变化，对应于与表面无关的云层和表面（Fan等人，2019年）13。它们具有可比较的奇异值，因此可以按照步骤 3.5.2 和 3.5.3 分隔曲面信息。对于未来对系外行星的观测，在极端情况下，无论是完全多云的还是无云的系外行星，SVD中只有一个占主导地位的PC出现在步骤3.3中。在这种情况下，光谱分析对于解释此 PC 的含义是必要的，因为云和表面的组成是不同的。如果主导 PC 对应于表面，则仍可以遵循步骤 4 和步骤 5;如果主要 PC 与表面相对应，则仍可以遵循步骤 4 和步骤 5。如果它对应于云，可以得出一个结论，即表面信息被云阻挡，因此不能使用给定波长的光曲线提取。在这种情况下，曲面映射是不可行的。第三个甚至第四个可比的主导PC也可能存在，它可以对应于另一层云或大规模水文过程，并且不会使方法的以下步骤失效，只要提取表面信息。
几何和光谱卷积导致的退化是限制检索地图质量的主要因素，如科万海峡（2013年）14和藤井等人（2017年）15。由于主流 PC 的时间序列只覆盖了 PC 平面的一小部分，因此在空间和光谱变化之间始终存在权衡。换句话说，检索到的地图（图2a）即使使用无限数量的时间步长和完美的观测值，只要使用相同的波长使用光曲线，也不可能得到太大的改善。我们引入正规化，以部分缓解退化。步骤 4.4 中正则化项 α 的最佳值是使用由地面真实点合成的观测结果确定的，其中观测到的 uj被视场 （FOV） 中的加权和缩放土地分数所取代。为了生成合成观测，我们使用步骤4.3中的方程，用每个像素的地真理陆地分数y代替x。y使用观测的PC2 u 2和平均 FOV 陆地分数 13 之间的强线性相关性，y与 x 缩放到相同的范围。由于退化，y无法从步骤4.4中的线性回归中完美恢复，因此我们通过查找最小值=2，平方残差按每个像素的方差缩放来确定 =的最佳值。后者按每个像素的绝对值进行估计。这与川原和藤井（2011）中的L曲线标准类似。在T=9739和P=3072的本论文中，α的最佳值为10-3.153（图4）。
观测噪声是影响制图质量的另一个因素，在实际操作中可能会破坏光曲线的 SVD 分析。我们通过将不同级别的观测噪声引入原始光曲线来测试协议的鲁棒性。它们被假定包含所有噪声源（例如天空背景、暗电流和读出噪声），并遵循高斯分布。在没有噪声的原始光曲线中，PC2 显示与 FOV 地位分数13的强 （r2±0.91） 线性相关性，因此其时间序列用于曲面映射。随着噪声水平的增加，PC2 与表面之间的相关性变弱（图 5）。当信噪比 （S/N） 小于 10 （图 5a） 时，相关系数 r2将低于0.5，尽管 PC2 的重要性仍然占主导地位（图 5b）。我们建议在将来泛化中自信地应用协议，最低为 30 的 S/N。值得注意的是，S/N是系外行星信号与观测噪声的比率，母星信号被移除。
在步骤 4.2 中假定查看几何体是已知的，因为从系外行星观测结果精确派生查看几何体超出了这项工作的范围。除了可以从光曲线观测和功率密度光谱（图2b和2d）中得出的轨道元素外，表面测绘只需要两个数量，夏/冬至和倾斜。夏/冬至通常与表面对应PC的时间序列的极值重合，只要北半球和南半球之间存在明显的不对称。从系外行星对光曲线17、18的振幅和频率的影响可以推断出外行星的弯曲性。所有这些推导都需要观测采样频率至少高于行星自转的频率，而行星自转目前很少满足系外行星的观测采样频率。
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
作者没有什么可透露的。
Acknowledgments
这项工作部分得到加州理工学院喷气推进实验室的支持，该实验室与美国宇航局签订合同。YLY 感谢华盛顿大学虚拟行星实验室的支持。
Materials
| Name | Company | Catalog Number | Comments |
| Python 3.7 with anaconda and healpy packages | Other programming environments (e.g., IDL or MATLAB) also work. |
References
- Schwieterman, E. W., et al. Exoplanet Biosignatures: A Review of Remotely Detectable Signs of Life. Astrobiology. 18 (6), 663-708 (2018).
- Campbell, B., Walker, G. A. H., Yang, S. A Search for Substellar Companions to Solar-type Stars. The Astrophysical Journal. 331, 902 (1988).
- NASA. NASA Exoplanet Archive (2019) Confirmed Planets Table. , (2019).
- Gillon, M., et al. Seven temperate terrestrial planets around the nearby ultracool dwarf star TRAPPIST-1. Nature. 542 (7642), 456-460 (2017).
- Kawahara, H., Fujii, Y. Global Mapping of Earth-like Exoplanets from Scattered Light Curves. The Astrophysical Journal. 720 (2), 1333 (2010).
- Fujii, Y., Kawahara, H. Mapping Earth Analogs from Photometric Variability: Spin-Orbit Tomography for Planets in Inclined Orbits. The Astrophysical Journal. 755 (2), 101 (2012).
- Cowan, N. B., Fujii, Y. Mapping Exoplanets. Handbook of Exoplanets. , Springer, Cham. (2018).
- Farr, B., Farr, W. M., Cowan, N. B., Haggard, H. M., Robinson, T. exocartographer: A Bayesian Framework for Mapping Exoplanets in Reflected Light. The Astronomical Journal. 156 (4), 146 (2018).
- Lomb, N. R. Least-Squares Frequency Analysis of Unequally Spaced Data. Astrophysics and Space Science. 39 (2), 447 (1976).
- Scargle, J. D. Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data. The Astrophysical Journal. 263, 835 (1982).
- Górski, K. M., et al. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. The Astrophysical Journal. 622 (2), 759 (2005).
- Jiang, J. H., et al. Using Deep Space Climate Observatory Measurements to Study the Earth as an Exoplanet. The Astronomical Journal. 156 (1), 26 (2018).
- Fan, S., et al. Earth as an Exoplanet: A Two-dimensional Alien Map. The Astrophysical Journal Letters. 882 (1), 1 (2019).
- Cowan, N. B., Strait, T. E. Determining Reflectance Spectra of Surfaces and Clouds on Exoplanets. The Astrophysical Journal Letters. 765 (1), 17 (2013).
- Fujii, Y., Lustig-Yaeger, J., Cowan, N. B. Rotational Spectral Unmixing of Exoplanets: Degeneracies between Surface Colors and Geography. The Astronomical Journal. 154 (5), 189 (2017).
- Kawahara, H., Fujii, Y. Mapping Clouds and Terrain of Earth-like Planets from Photomertic Variability: Demonstration with Planets in Face-on Orbits. The Astrophysical Journal Letters. 739 (2), 62 (2011).
- Kawahara, H. Frequency Modulation of Directly Imaged Exoplanets: Geometric Effect as a Probe of Planetary Obliquity. The Astrophysical Journal. 822 (2), 112 (2016).
- Schwartz, J. C., Sekowski, C., Haggard, H. M., Pall ́e, E., Cowan, N. B. Inferring planetary obliquity using rotational and orbital photometry. Monthly Notices of the Royal Astronomical Society. 457 (1), 926-938 (2016).
