Zahira Nakhli

Ben Fafa

Martine Hatira

Patrick Pithioux

Chabrand

Khemais

Fafa Ben Hatira

Martine Pithioux

Patrick Chabrand

Khemais Saanouni

using a quasi-brittle continuum damage finite element model

Keywords: Vertebrae fracture, damage law, elastic modulus, finite element analysis

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The motivation of this research is to provide valuable improvement for the surgical diagnosis and aims at reducing healthy and osteoporotic vertebrae fractures, as well. The originality of this paper is the quasi-brittle law simplicity describing the bone damage in such cases like cars' crash or severe fall.

In this study, numerical models are developed and implemented in order to identify vertebral fracture over expected experimentally compression load ranges. Thanks to this study, answers are given to questions about the damage affecting vertebral bodies, such as the ultimate load magnitude as well as the damage state of all vertebrae.

The osteoporosis that is affecting more and more people (man and woman) is a disease which characterized by the bones 'weakness that increases the risk of their breaking. It is an age-related disease which causes progressive skeletal disorder characterized by a low bone mass and microarchitectural alterations increasing bone fragility and susceptibility to fracture [START_REF]Consensus development conference: "diagnosis, prophylaxis, and treatment of osteoporosis[END_REF].This causes pain and a significant loss of mobility. Vertebral fractures are the most [START_REF] Kanis | Diagnosis of osteoporosis and assessment of fracture risk[END_REF] common type of osteoporotic fracture and they are associated to the decreasing quality of life [START_REF] Lips | Quality of life in patients with osteoporosis[END_REF] and to the substantial morbidity [START_REF] Kendler | Vertebral Fractures: Clinical Importance and Management[END_REF].

Clinically speaking, osteoporosis is diagnosed using the dual-energy X-Ray absorptiometry (DXA) [START_REF] Kanis | Diagnosis of osteoporosis and assessment of fracture risk[END_REF][START_REF] Marshall | Meta-analysis of how well measures of bone mineral density predict occurrence of osteoporotic fractures[END_REF][START_REF] Schuit | Fracture incidence and association with bone mineral density in elderly men and women: the Rotterdam Study[END_REF]. The system measures the bone mineral density (BMD). The fracture risk is assumed to be highly correlated with the value measured by the DXA. Osteoporosis is defined for bones with a BMD between 833 and 648 mg/cm 2 (equivalent to a T-score between -1.0 and -2.5), whereas if a BMD is above 833 mg/cm 2 (T-score above -1.0), they are considered as healthy. Experts give osteoporosis medications with a BMD less than 648 mg/cm 2 (T-score less than -2.5). It is also noteworthy that DXA is submitted to restrictions for the assessment of bone fragility and osteoporosis diagnosis [START_REF] Mccreadie | Biomechanics of Fracture: "Is Bone Mineral Density Sufficient to Assess Risk?[END_REF]. DXA is based on 2D and not 3D measurements and used only for porous bone information and not for compact bone. It is known that microarchitecture alteration is currently included in the definition of osteoporosis [START_REF]Consensus development conference: "diagnosis, prophylaxis, and treatment of osteoporosis[END_REF]. From a numerical point of view, several FE model describing damage and fractures for bone structure are proposed in the literature. Most of these works introduce linear and non-linear isotropic and anisotropic finite element (FE) models. Each work presents a new model with some answers about the best way to predict bone structures damage. Doblaré et al [START_REF] Doblare | Anisotropic bone remodeling model based on a continuum damage-repair theory[END_REF] present a study of the proximal femoral extremity remodeling based on an anisotropic bone model and state that the evolution of bone microstructure internal variables can be formulated following the Continuum Damage Mechanics (CDM) principles.

They show that a good similarity is obtained between numerical and experimental results. The same authors present [START_REF] Doblare | On Numerical Modeling of Growth, Differentiation and Damage in Structural Living Tissues[END_REF] another work for both hard and soft tissues. A damage-bone remodeling theory presented in [START_REF] Ramtani | Computer simulation of an adaptive damage-bone remodeling law applied to three unit-bone bars structure[END_REF] shows that under constant load, the three unit bone bars structure exhibit inhomogeneous strains which strongly depend upon the manner in which the microdamage is distributed

To emphasize the distinct damage behavior in tension and compression, Garcia & al presented in [START_REF] Garcia | A three-dimensional elastic plastic damage constitutive law for bone tissue[END_REF] 3D anisotropic constitutive equations. For the sake of simplicity, Hambli et al ( [START_REF] Hambli | A quasi-brittle continuum damage finite element model of the human proximal femur based on element deletion[END_REF], [START_REF] Hambli | Finite element prediction of proximal femur fracture pattern based on orthotropic behavior law coupled to quasi-brittle damage[END_REF]) presented in these two works respectively a simple finite element (FE) model coupled with quasi-brittle damage model for the proximal femoral fracture prediction based on orthotropic behavior and a numerical FE-based model with an element deletion to describe the cracks propagation. Many Quantitative computed tomography (QCT)-based finite element (FE) voxel models have been presented to predict ultimate force of human lumbar vertebrae under axial compression. As mentioned in [START_REF] Schileoa | Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and fracture location on human femurs tested in vitro[END_REF], no clear choice of the failure criterion to adopt for the bone tissue has been found to predict fracture risk of bones. In that work, a combined experimentalnumerical approach has been used to build up a subject-specific finite element models which are able to accurately predict failure patterns of bones. In [START_REF] Ph | A 3D Model for Trabecular bone based on fabric tensors[END_REF] a theoretical model formulated within the framework of continuum damage mechanics and based on fabric tensor is proposed. A good correlation between experimental and numerical results has been found. The authors of [START_REF] Mengoni | A generic anisotropic continuum damage model integration scheme adaptable to both ductile damage and biological damage-like situations[END_REF] proposed an original bone remodeling law coupled to trabecular bone plasticity for the simulation of orthodontic tooth movements by using a phenomenological approach of anisotropic Continuum Damage Mechanics. Over the last decades, developments in (3D) have provided possibilities for measuring a variety of structural indices to characterize bone microarchitecture. However, standard mechanical tests did not provide easily the elastic constants due to the small size of specimens of human bone. The dependency of microarchitecture on bone's mechanical properties has been presented in [START_REF] Ulrich | The ability of three-dimensional structural indices to reflect mechanical aspects of trabecular bone[END_REF][START_REF] Kowalczykp | Elastic properties of cancellous bone derived from finite element models of parameterized microstructure cells[END_REF]. The authors showed that the elastic constants, specifically the young's modules are clearly correlated with the tissue morphology bone coefficient BV/TV. The results presented in [START_REF] Morgan | Trabecular bone modulusdensity relationships depend on anatomic site[END_REF] indicate that there is no universal modulus-density relationship for on-axis loading. The authors suggest that the sitespecificity in apparent modulus-density relationships may be attributed to the differences in architecture. In [START_REF] Ph | A 3D Model for Trabecular bone based on fabric tensors[END_REF], authors showed that there is a strong correlation between cumulated permanent strain and elastic modulus which indicate that these variables affect the damage process. For isotropic cellular materials, a basic model consists of a power relationship of the Lamé constants with respect to volume fraction [START_REF] Gibson | The Mechanical behavior of calcellous bone[END_REF]. For anisotropic materials, an averaging method was developed by Cowin et al in [START_REF] Cowin | Thermodynamic restriction of the elastic constants of bone[END_REF] to identify the type of elastic symmetry and express the dependency of the elastic constants toward the volume fraction in [START_REF] Yang | The Anisotropic Hooke's Law for Cancellous Bone and Wood[END_REF].

In general, most of the articles found in the literature show how it is important to provide a precise diagnosis of fracture state (magnitude and localization) and how to give different answers for different types of bone (healthy or not), and how to establish behavior model to quantify the bone fracture. However, they showed controversies in different conclusions on the same topic.

As it has been reported in [START_REF] Hambli | A quasi-brittle continuum damage finite element model of the human proximal femur based on element deletion[END_REF][START_REF] Hambli | Finite element prediction of proximal femur fracture pattern based on orthotropic behavior law coupled to quasi-brittle damage[END_REF] a quasi-brittle damage model coupled to an orthotropic behavior law can predict the proximal femur fracture. In [START_REF] Mirzaei | On prediction of the strength levels and failure patterns of human vertebrae using quantitative computed tomography (QCT)-based finite element method[END_REF], the authors present a study to predict the failure of human vertebra based only the strength level and failure patterns.

Our model based on a simple isotropic damage law gives a unified and reliable answer to different quality of bones and precisely predict the ultimate fracture load as well as the damage localization for vertebra.

Studies found in the literature followed different approaches to estimate the Young modulus. For example, authors of [START_REF] Ph | A 3D Model for Trabecular bone based on fabric tensors[END_REF][START_REF] Clouthier | Finite element analysis predicts experimental failure patterns in vertebral bodies loaded viaintervertebral discs up to large deformation[END_REF] reported that the mechanical behavior of trabecular bone is mainly governed by its tissue modulus and morphology, i.e. bone volume fraction (BV/TV) especially in the case of small strain. In studies [START_REF] Linde | The effect of specimen geometry on the mechanical behavior of trabecular bone specimens[END_REF][START_REF] Carter | The compressive behavior of bone as a two-phase porous structure[END_REF][START_REF] Snyder | Estimation of mechanical properties of cortical bone by computed tomography[END_REF][START_REF] Kaneko | Mechanical properties, density and quantitative CT scan data of trabecular bone with and without metastases[END_REF] , authors adopted elasticity-apparent density relationships to build FE models.

In the present work, three simple methods proposed in different works in the literature computing the Young's modulus E (MPa) that have been implemented in the numerical model in order to compare them and to show their performance toward reproducing the initiation and the propagation of the damage. The first method is based on the" grey value ". The second method has been proposed by Yang et al in [START_REF] Yang | The Anisotropic Hooke's Law for Cancellous Bone and Wood[END_REF]. Finally, the third one is presented by Hernandez in [START_REF] Hernandez | The Influence of Bone Volume Fraction and Ash Fraction on Bone Strength and Modulus[END_REF]. These three methods are going to be presented in the next section.

 Method 1: Many researchers such as Jacobs et al [START_REF] Jacobs | Numerical Simulation of Bone Adaptation to Mechanical Loading[END_REF], Keyaketal [START_REF] Keyak | Automated threedimensional finite element modeling of bone: A new method[END_REF], Helgason et al [START_REF] Helgason | Mathematical relationships between bone density and mechanical properties : A literature review[END_REF], Madi et al. [START_REF] Madi | Finite element simulations of the deformation of fused-cast refractories based on X-ray computed tomography[END_REF], Tawara et al [START_REF] Tawara | Mechanical evaluation by patient-specific finite element analyses demonstrates therapeutic effects for osteoporotic vertebrae[END_REF], and Jovanovici et al [START_REF] Jovanovic | Finite element modeling of the vertebra with geometry and material properties retrieved from CT-Scan Data[END_REF], Gislanson et al [START_REF] Gislason | Finite element modeling of the femur bone of a subject suffering from motor neuron lesion subjected to electrical stimulation[END_REF] showed that the Young's modulus for each finite element depends on the apparent density for each bone voxel (E = f ()), which in its turn was calculated on the basis of CT images. Relationship between the apparent density and CT values (HU: Hounsfield Units or GV: gray values) is also calibrated. The chosen formula can be expressed as:

E = 3500 ρ 2.2 (1)
 Method2: In [START_REF] Yang | The Anisotropic Hooke's Law for Cancellous Bone and Wood[END_REF], Yang et al propose a dependent orthotropic Hooke's model for cancellous bone based on the identification of the elastic constants of trabecular bone in function of volume fraction (BV/TV).

(

) 2
Where E T is the slope obtained from the experimental curves (see Figure .3)  Method 3: In [START_REF] Hernandez | The Influence of Bone Volume Fraction and Ash Fraction on Bone Strength and Modulus[END_REF], Hernandez et al present the used parameters to predict bone material properties with power law functions of the form:

(

) 3
Where y is the strength or elastic modulus, x a parameter (i.e apparent density) and a and b are empirical constants derived from experimental data. They propose the following relationship where the x on the preceding power function is the bone volume / total volume ratio (BV/TV): [START_REF] Kendler | Vertebral Fractures: Clinical Importance and Management[END_REF] In what follows the material, the experimental method and the numerical models with the obtained results are going to be detailed. The first one introduces the procedure to obtain data and geometry of the six CT scanned lumbar vertebra, followed by the experimental curves obtained for the vertebrae tested under compression till fracture. As pointed out previously, Continuum Damage Mechanics 'CDM' framework is chosen since it is the most appropriate constitutive framework to reproduce the bone structure failure as it will be detailed in section 2.In the third part, the boundary and load conditions are presented. The simulations of different cases are detailed in section 4 for the sake of a comparison between theoretical and experimental results. ). These parameters were also used in the finite element computations.

Material and Method

In Figure 1, the flow chart summarizes the protocol established to create FE models from CT data using material properties is presented following the steps detailed previously.

Experimental Mechanical Compression test

One of the most used classifications of type of trauma is the one introduced by Magerl [START_REF] Magerl | A comprehensive classification of thoracic and lumbar injuries[END_REF] based on the morphological analysis of lesion. The called type A is attributed to the compression injury. It is subdivided into three types (Group A1 fracture on the top of the vertebra; A2 Group fracture separation; A3 Group, burst-fracture [START_REF] Ch | Comment nous traitons les fractures récentes du rachis thoraco-lombaire[END_REF]. To reproduce the trauma produced by cars' crash or severe fall, compression of lumbar vertebra investigations were performed on six lumbar vertebrae obtained from four human cadavers. An INSTRON 5566 device was used for the compression test. Vertebrae with embedded epoxy resin were placed between the jaws (Figure 2a). Resin was maintained during the entire test enabled vertebrae to be set in the vertical axis of the compression device. A velocity of 5 mm/min was imposed during the test. The failure load (in kN) and displacement were measured. During the mechanical test, the mean failure load was (2.4kN). The initial crack always occurred in the middle of the vertebra .Experimental ultimate load, final displacement, BMD for all vertebrae are reported in Table 1.

The next picture (Figure 2b) taken during the compression of one of the specimens, show clearly that the failure occurs in the middle of the vertebra with the wedge fracture which is close to the group A3 with a reduction of the height of the vertebra.

In Figure 3, the experimental data obtained for the six specimens are reported. The value of the tangent modulus E T used in the second method for young's modulus estimation 2 (Eq.2) are mentioned in the same curves. The maximum value (E T = 39 MPa) is obtained for the specimen 3-L2, whereas the lowest one (E T = 7 MPa) is obtained for the specimen 2-L2. It seems clear that this value depends on the classification of vertebra (healthy or osteoporotic). Indeed, specimen 3, which is reported as healthy, has a BMD of 0.884, whereas, specimen 2, which is reported as osteoporotic, has a BMD of 0.616. The load-displacement curves obtained for three different vertebrae (L2, L3, L4) taken for the same donor (specimen 3), showed the same tendency with a variation in the fracture load magnitude [START_REF] Clouthier | Finite element analysis predicts experimental failure patterns in vertebral bodies loaded viaintervertebral discs up to large deformation[END_REF].Indeed, taken in the same segment, the different lumbar vertebrae will not have the same orientation with respect to the vertical axis of the standing or sitting position. Therefore the same compression forces will act differently on them depending on the position of the vertebra in the spine segment.

In order to reproduce the obtained experimental data as the localization of the crack, the state of damage , the ultimate load failure, a new CDM model coupled with the three young modulus (equations1, 3, 4) is presented in the next section.

Bone constitutive models

Constitutive framework: A quasi brittle damage law

The approach of irreversible thermodynamics with internal variables (Chaboche [START_REF] Chaboche | Continuum Damage Mechanics ", Part I and II[END_REF], Germain [START_REF] Germain | Mécanique des Milieux" continus[END_REF], Krajcinovic [START_REF] Kracinovic | Damage Mechanics[END_REF], Saanouni [START_REF] Saanouni | On the Anelastic flow damage[END_REF], Lemaitre [START_REF] Lemaitre | Formulation Unifiée de lois d'évolution d'endommagement[END_REF], Kachanov [START_REF] Kachanov | Introduction to Continuum Damage Mechanics[END_REF]), is chosen to present a coupled damage elastic model to describe the initiation and the accumulation of the damage in bone structure, more precisely the vertebra. In this work, the damage behavior law describing a quasi-brittle behavior is proposed using an isotropic Continuum Damage Mechanics (CDM) based on Marigo modeling of the damage for brittle and quasi-brittle material behavior [START_REF] Marigo | Formulation of a damage law for an elastic material[END_REF]. This theory is formulated to describe the progressive degradation of material (CDM) models. The new energy based model is described throughout state variable (external and internal).The state variables describing the constitutive equations are represented by the external and the observable variables, respectively the elastic strain tensor and the Cauchy stress tensor . For the sake of simplicity, damage is supposed isotropic described by a couple of scalar internal variables (D, Y) where Y is damage associated variable. The effective variables and including the damage effect which are defined in the framework of the elastic strain equivalence assumption are presented hereafter. The expression of the stored elastic energy density is given by: [START_REF] Marshall | Meta-analysis of how well measures of bone mineral density predict occurrence of osteoporotic fractures[END_REF] Where is the symmetric fourth-rank tensor of elastic properties of the virgin (not affected by damage) material which in the isotropic case can be written in terms of the well-known Lame's constants  and  according to:

A 1 1 2 1       
Where 1 is the second-rank identity (Kronecker) tensor while 1 is a fourth-rank unit tensor.

According to the theory of Marigo, the free energy per unit mass depends only on the two state variables namely the elastic strain tensor and the damage. The state laws and Y are classically derived from the state potential are obtained from the free energy by: (6.1) (6.2)

 1  1  1 (6.3) (7)
The damage criterion (or damage yield function) is described by Y:

(8) , s and m are the material parameters.

The parameters s and m are related to the damage "hardening" of the material. It is here to be noticed that the damage yield function [START_REF] Doblare | Anisotropic bone remodeling model based on a continuum damage-repair theory[END_REF] can describe the initiation of micro-cracks starting from undamaged state (D=0).

In the present model, the dissipation potential is reduced to the yield function f (associative theory): [START_REF] Doblare | On Numerical Modeling of Growth, Differentiation and Damage in Structural Living Tissues[END_REF] The evolution laws derived from the dissipation potential are : [START_REF] Ramtani | Computer simulation of an adaptive damage-bone remodeling law applied to three unit-bone bars structure[END_REF] For this approach, the coupling between damage and elasticity is completed with the following damage evolution law. [START_REF] Garcia | A three-dimensional elastic plastic damage constitutive law for bone tissue[END_REF] With: ee Y : :

    & & (12)
The above presented model emphasizes on the fact that the damage evolution is the influent parameter since it will affect the young modulus degradation as well as the evolution of the Cauchy stress. It will implicitly produce the damage propagation. Solving the nonlinear problem described by equations (6, 7,8,9,10,11 and 12) in order to determine the unknowns of the problem is performed through an approximation of these variables in total time interval . t is the increment between two successive time steps. This approximation is done for every integration point related to every finite element following the robust finite element method.

Thus knowing the initial variables at t n , the discretized problem is solved giving the final solution at the final time t n+1 .The discretization leads to the following expressions of the problem variables at t n+1 =t n +t, the end of the step time:

  (13) ( 14
) (15) (16)
From this last equation the "admissible" value of the damage variable is deduced as:

In this isotropic damage modeling, some remarks can be made:  If the scalar variable describing damage D=0, then the material state is described by the classical isotropic elastic model.  If the fracture condition of the critical value of D=1 is reached. D=0.999 is then assigned to theoretical damage value for numerical consideration and the state of final fracture.

We will call the damage prediction combined with the first method (eq 1) model 1, the one that is coupled with the second method (eq2), model 2. Finally, the damage model associated with the third method (eq 4) will be called model 3. These three damage models and the cracks propagation technique have been implemented into the FE code Abaqus/Standard via the user subroutine VMAT which uses an explicit integration method to solve the nonlinear problem.

Simulations  Boundary and Loading conditions

In general, the fracture of the vertebra occurs when the bone is subjected to a high compression load. A compressive pressure is applied then to the upper surface and the lower surface is constrained, (Figure4). The computations were carried out using the kill element method. This method simulates the propagation direction by setting the stiffness matrix to zero when the critical damage value (D=0.999) is reached inside an element of the mesh. This leads to the distribution of the stress state in the damage area. The pressure loading the vertebra is calculated from the experimental load obtained throughout the previous testing presented in section 1.

Results

As mentioned previously, the purpose of this work is to predict effectively the damage localization as well as the ultimate fracture strength for different specimens tested experimentally and presented in section 1.2. The ultimate strength load value obtained experimentally was applied for the three numerical models presented in section 2.1. The nonlinear analysis showed that the localized damaged zones were found to be different for each model. The results of the fitting of computed and experimental data are described in the next paragraph. The parameters Y 0 , m and s used in equation [START_REF] Garcia | A three-dimensional elastic plastic damage constitutive law for bone tissue[END_REF] have been chosen for every vertebra to best fit the experimental data.

The predicted load-displacement curves based on the variation of the young's modulus models (model1, model 2, model 3) presented in figure 5 have been obtained for the specimen 2-L2 and specimen 3-L2. The two chosen curves have been selected among the six experimental ones from figure 3. They have the highest and the lowest value of ultimate loads. The figure showed a similarity between the experimentally quasi-brittle failure and the numerical obtained results obtained and the best results are found when the quasi-brittle law is combined with for the Hernandez law (model 3) for the healthy vertebra (specimen 3-L2). It showed also that the predicted crack path with different stages of damage localization for the same vertebra. It is obvious that the two curves fit when damage reaches the value of 0.22 which means that the two curves tend to be similar when the damage starts to grow. The ultimate fracture load under compression (3.8±0.21 KN) that occurs for the critical value D=0.52. Indeed, for compression the critical damage can be reduced to the value D= 0.5 [START_REF] Hambli | Multiscale prediction of crack density and crack length accumulation in trabecular bone based on neural networks and finite element simulation[END_REF]. The strong damage accumulation occurs in inner of the vertebra as mentioned in [START_REF] Garcia | A three-dimensional elastic plastic damage constitutive law for bone tissue[END_REF] leading to the crush of the vertebral body as observed experimentally (Figure1b).The numerical curve obtained by model 2 showed a similarity in the value of the ultimate fracture load (4.79 ± 0.91 KN) but with a permanent shift toward the experimental curve. For the osteoporotic vertebra (specimen 2-L2) the comparison between the numerical result predicted by the model 1 and the experimental data (0.682 ±0.265 KN) shows a good agreement. However, the damage localization in that case is on the top of the vertebra (Figure 5). A better similarity is obtained for the same osteoporotic vertebra from the model 2 with crack localization in the middle of the vertebra as it will be detailed in Figure 6.

Table 3 summarizes the obtained computed ultimate loads for the six vertebrae. The standard deviation SD for the models 1, 2 and 3 were respectively included in the intervals [0.265; 2.14], [0.027; 0.91], and [0.048; 0.446]. Model 3 is then the most accurate one since it gives the lowest difference in comparison with the experimental data.

Figure 6 shows that in case we choose model 1, the damage localization is on the top of the vertebra with low value of critical damage [0.177; 0.392]. When we choose model 2 and 3, a good localization is obtained [0.478; 0.57]; i.e. in the middle of the vertebra (as observed experimentally). The choice of model 1 leads to a localization zones located on the top of vertebrae for all computed cases.

From the latest results, it is concluded that the choice of model 1 (apparent density based Young's modulus model) is more relevant for the osteoporotic vertebra; however it gives a localization on the top of the specimen. Numerical model3 expressing exclusively the dependency of Young's modulus on the tissue morphology bone fracture BV/TV is suitable for healthy vertebra (good crack localization). The choice of model 2 (including the value of E T and BV/TV) did give the acceptable results for both vertebra cases (osteoporotic or healthy) with good damage localization.

Discussion

The aim of this work was to develop a simple quasi-brittle model to describe the process of vertebral fracture and to compare it to the experimental load-displacement curves. Constitutive equations were developed using a CDM model to be the best to fit to the experimental data. The results of the numerical investigation showed the capabilities of the proposed FE element model to describe and predict the localization of vertebra failure based on the choice of the right modulus, parameters of the evolution damage law (Y o , s, m) which are hard to obtain experimentally. The accuracy of the results is improved by the use of the BM/TV ratio measured for each vertebra. The quasi brittle isotropic models which are so simplistic can then predict plausibly the compression of bones. The results show that the numerical force-displacement curve using the proposed behavior laws are adequate. This is due to the quasi-brittle isotropic models that made plausible prediction about the results related to the experimental force-displacement curve. The obtained results in this work clearly indicate the dependency of Young's modulus model on the vertebra classification. Adjustment of models can be then applied depending on the BMD of every specific studied case.

It has to be noticed that the numerically predicted load-displacement curves was different from the experimental curves when the value of D was very low. This means that the damage models overestimate the damage at the beginning of the loaddisplacement curves. The presented numerical models based on an isotropic elastic behavior for bone structure give some good quantities and qualitative results. A correlation between micro and macro mechanisms is suitable to better understand the whole physical process of damage in order to enhance the damage modeling. Other behavior models should be then chosen like the ones proposed in [START_REF] Kinzi | Identification of a crushable foam material model and application to strength and damage predication of human femur and vertebral body[END_REF] based on crushable foam plasticity models (an approximation of Tsai-Wu yield function). In addition, for vertebrae highly affected by osteoporosis, 3matic presents some limitation on the quality of the bone structure meshing. Besides, it is the continuum level is largely depending on both elements size and orientation and accordingly the mesh size should be accurately determined for each type of material [43, 44, and 45]. Another limitation of this study was the inclusion of only six samples. This was due to the complexity to generate finite elements meshes. Further studies including more appropriate FE meshing algorithm are indeed necessary to treat the mesh dependency problem. However, despite these limitations, the CDM modeling did show some real potential to predict the right localization of fracture for different types of vertebra (healthy or osteoporosis).

Conclusion

The purpose of this work was to develop and validate a simple FE model based on continuum damage mechanics in order to simulate the complete force-displacement curve of lumbar vertebra failure. The quasi-brittle law prediction bone failure corroborates the experimental response of vertebral bone failure under a quasi-static compression loading. The obtained results show important similarities for both models (model 2 and 3) for healthy vertebra and the significant differences emerged for model 1. For osteoporotic vertebrae, the model 1 seems to be more appropriate. Despite these differences, the potential of the CDM model combined with the right young modulus law is proved through the computed results. In general, concerning FE simulations, the combination of the parameters like geometry and BV/TV (QCT) and the BMD measurement (DXA) did improve the numerical results.

  ) were CT scanned by using a GE Medical Systems scanner available in La Timone University hospital, (Marseille. DICOM images files generated by the scanner are constituted by pixels with different gray intensities). The investigations were approved by research ethics board at the University Hospital.  Step 1: 3D reconstruction has been accomplished using the technique of density segmenting with the research software Mimics 17.0. The generation of the surface and volume meshes was made by the research software 3Matic 9.0.0.231. It was used to create a triangular mesh on the surface of the vertebra, and to generate the volume mesh with linear tetrahedral elements (C3D4, six degrees of freedom per node). The volume mesh was imported again in Mimics, in order to assign the material parameters based on the Hounsfield Units 'HU' or Units Gray values on the scanned images. The density  and the young Modulus E are obtained at the end of this step and the geometry is used in the finite element simulation.

	1.1 Sample preparation
	Six lumbar vertebrae obtained from four donors (female=1, male=3, age = 82 ± 9
	years [range: 63 -91]



Step 2: Through the software Image J combined with the software Bone J, the architectural parameters were extracted: BV/TV (Bone Volume over Total volume; i.e. bone volume fraction), Tb.Th (trabecular thickness), Tb.Sp (Trabecular Spacing) and Tb.N (Trabecular Number

Table 1 . : Failure values for the six lumbar vertebrae.

 1

	Specimen	vertebra	BMD (g/cm2)	Load (KN)	Displacement (mm)
	1	L2(healthy)	1.016 1.615	7
	2	L2(Osteoporotic)	0.616 1.057	9.87
	3	L2(healthy)	0.884 3.501	4.6
	3	L3(healthy)	0.938 3.074	9.14
	3	L4(healthy)	0.881 2.839	5.77
	4	L2 (Osteoporotic) 0.819 2.437	8.4
		SD	0.136 0.924	2.037

Table 2 .1 : Material properties for bone used for the simulation (specimen 2-L2) GR : E = 3500 ρ 2.2

 2

		Density (g/cm3)		Young 's Modulus (MPa)	Poisson ratio
	Cortical bone	0.64-2.09		1330-17859	0.3
	Concellous bone 0.04-0.41		4-517	0.3
	Yang´s Law : E = 1240 * ET * (BV/TV) 1.80
		Density	BV/TV	Young 's Modulus (MPa)	Poisson ratio E T
		(g/cm 3 )
	Vertebra	0.03-1.74	0.3	993.89	0.3	7

Hernandez´s Law : E = 84370*(BV/TV) 2.58

		Density	BV/TV	Young 's Modulus (MPa)	Poisson ratio
		(g/cm 3 )
	Vertebra	0.03-1.74	0.3	3777.12	0.3

Table 2 .2 : Material properties for bone used for the simulation (specimen 3-L2)

 2

	GR : E = 3500 ρ 2.2
		Density (g/cm 3 )		Young's Modulus (MPa)	Poisson ratio
	Cortical bone	0.619-1.89		1219-14299	0.3
	Cancellous bone	0.05-0.414		6-505	0.3
	Yang´s Law : E = 1240 * ET * (BV/TV) 1.80
		Density	BV/TV	Young's Modulus	Poisson ratio	E T
		(g/cm 3 )		(MPa)
	Vertebra	0.05-1.89	0.3	5537.38	0.3	39
	Hernandez´s Law : E = 84370*(BV/TV) 2.58
		Density	BV/TV	Young Modulus (MPa) Poisson ratio
		(g/cm 3 )
	Vertebra	0.05-1.89	0.3	3777.12	0.3
	Resin	1.2		4220	0.4

Table 3 : Ultimate computed load (SD) for the six specimens

 3

		Specimen 1-	Specimen 2-	Specimen 3-	Specimen 3-	Specimen 3-	Specimen 4-
		L2	L2	L2	L3	L4	L2
	Numerical	1.037(0.408) 0.682	0.467(2.145) 1.131(1.373) 1.898(0.665)	0.759(1.17)
	(Model 1)		(0.265)
	Numerical	2.074	0.977(0.056) 4.796(0.915) 3.035(0.027) 3.576(0.5211) 1.478(0.661)
	(Model 2)	(0.324)
	Numerical	1.712	1.274	3.805(0.215) 3.142(0.048) 1.782(0.446)	1.782(0.446)
	(Model 3)	(0.068)	(0.153)


W M Keck

Garreth Ruane

Henry Ngo

Dimitri Mawet

Olivier Absil

Elodie Choquet

Therese Cook

Carlos Gomez Gonzalez

Elsa Huby

Keith Matthews

Tiffany Meshkat

W M Keck Observatory

Élodie Choquet

Maddalena Reggiani

Eugene Serabyn

Nicole Wallack

Reference star differential imaging of close-in companions and circumstellar disks with the NIRC2 vortex coronagraph at

Keywords: techniques: high angular resolution, planets and satellites: detection, protoplanetary disks, stars: imaging, stars: individual: MWC 758, stars: individual: 2MASS J16042165-2130284

INTRODUCTION

Adaptive optics (AO) imaging surveys determine crucial occurrence statistics for giant exoplanet populations with masses 1 M Jup at orbital separations >10 au Corresponding author: Garreth Ruane gruane@caltech.edu Current and future infrared, high-contrast imaging programs seek to (1) bridge the gap between the search completeness of direct imaging and radial velocity surveys, (2) characterize the atmospheres of lower mass planets, and (3) understand the interaction between planets and the circumstellar disks from which they form. All of the above science cases benefit from minimizing the angular separations from the host star at which instruments are sensitive enough to detect and characterize the planets and disks of interest.

Reducing the inner working angle (i.e. the separation at which the transmission is 0.5) of ground-based, high-contrast imaging systems to their theoretical limits (i.e. ∼ λ/D, where λ is the wavelength and D is the telescope diameter; Guyon et al. 2006) opens a critical search space at solar system scales (1-30 au) for stars within a few hundred parsecs. From a technical point of view, it will also pave the way to studies of temperate exoplanets in reflected light with future large aperture telescopes (D ∼30 m), especially in the habitable zones of nearby, late-type stars where the planet-to-star flux ratio is expected to be ∼10 -8 for planet radii of ∼ 1R ⊕ (see e.g. Guyon 2018). Looking further into the future, these developments will also help maximize the number of Earth-like planets orbiting solar-type stars available for study with future exoplanet imaging space missions, such as the Habitable Exoplanet Observatory (HabEx; Gaudi et al. 2018) and Large UV/Optical/IR Surveyor (LUVOIR; The LUVOIR Team 2018) concepts.

The vortex coronagraph mode (Foo et al. 2005;Mawet et al. 2005) of the NIRC2 instrument at W.M. Keck Observatory (Serabyn et al. 2017;Mawet et al. 2017;Ruane et al. 2017;Huby et al. 2017;[START_REF] Guidi | [END_REF][START_REF] Reggiani | [END_REF]Xuan et al. 2018) has the optical throughput needed to search for giant planets down to λ/D, or ∼100 mas, for L and M s bands (λ=3.4-4.8 µm). However, pairing the Keck adaptive optics system with a small inner working angle coronagraph does not immediately provide sensitivity to point sources with ∆mag 3 within 0. 5 of the host star in our experience. Surpassing this level also requires optimized observing strategies and post-processing [START_REF] Mawet | Proc. SPIE[END_REF].

The limitations of ADI

High contrast imaging surveys with Keck/NIRC2 traditionally make use of angular differential imaging (ADI, [START_REF] Marois | [END_REF] in the search for giant planets, such as the four planets discovered orbiting HR 8799 (Marois et al. 2008(Marois et al. , 2010;;Konopacky et al. 2016). In this strategy, the telescope's field rotator is set to vertical angle mode, which keeps the beam fixed in azimuthal orientation with respect to the elevation axis and detector while the field of view rotates over the course of an observation. This assists in differentiating between true astrophysical objects, which revolve in the image, and starlight that leaks through the coronagraph, which typically leaves a residue of speckles in the image.

Although ADI is the most successful observational strategy for detecting giant planets at wide separations to date (see recent review by [START_REF] Chauvin | Proc. SPIE[END_REF], it has some innate limitations. In order to clearly identify companions, each target must be observed for long enough to allow enough parallactic angle (PA) rotation such that a companion would move by a substantial fraction of the width of the point spread function (PSF). For instance, at the inner working angle of a vortex coronagraph (0. 08-0. 1 for NIRC2 L and M s bands), the minimum desired PA rotation is ∼30 • , which moves the PSF of a companion by an arclength of roughly half of its full width at half maximum (FWHM). However, due to practical limitations, observations from ongoing surveys with Keck/NIRC2 vortex coronagraph have a median PA rotation of ∼11 • (Xuan et al. 2018).

To illustrate the problems associated with scheduling ADI observations, Figure 1 shows the PA rotation achieved from Maunakea, Hawaii (latitude 19.8 • N) in a one hour observation versus the target's declination and the timing of the observation. An "on-time" observation (thick, blue, dashed line) in this case starts 30 min before the target crosses the meridian. The thin black lines shows how the PA rotation degrades when observation is initiated too early or too late by multiples of ∆t=15 min.

Figure 1 has several important implications:

• Timing constraints: ADI observations must be carried out during a relatively narrow time window to achieve a PA rotation needed to avoid selfsubtraction effects [START_REF] Marois | [END_REF]. For example, in order for a planet to move 30 , ruling out more than 90% of the sky.

• Inefficient surveys of star forming regions: Many of the most attractive targets for highcontrast imaging in the infrared are in nearby star forming regions (distances 120-150 pc), such as the Taurus and ρ Ophiuchi molecular clouds (see e.g. [START_REF] Bowler | [END_REF]). However, these regions only extend over a few square degrees, allowing for one or two ADI sequences on these targets per night. Therefore, an ADI survey of these regions would need to be carried out over many nights.

• Limited "effective" inner working angle: Unless the star falls within the declination ranges described above and the observations are welltimed, self-subtraction effects at small angular separations [START_REF] Marois | [END_REF] ultimately limit the sensitivity near to the star and inner working angle of the coronagraph.

• Erases rotationally symmetric circumstellar disks: Imaging scattered light from close to pole-on protoplanetary, transitional, and debris disks with ADI may not be possible because rotationally symmetric features self-subtract in postprocessing and complex dust disk distribution become confounded by artifacts (Milli et al. 2012).

The limitations of SDI

An alternate speckle subtraction method known as spectral differential imaging (SDI; Sparks & Ford 2002) is applied on instruments with integral field spectrographs, such as Palomar/P1640 (Hinkley et al. 2011), Gemini/GPI [START_REF] Macintosh | Proc. Natl. Acad. Sci[END_REF], VLT/SPHERE (Vigan et al. 2016), and Subaru/SCExAO/CHARIS (Jovanovic et al. 2015;[START_REF] Groff | Proc. SPIE[END_REF]. SDI decouples speckles from astrophysical objects based on their wavelength dependence rather than azimuthal rotation and therefore alleviates some of the aforementioned limitations of ADI. Between the shortest and longest wavelengths, speckles move radially by α × ∆λ/λ, where α is the angular separation in units of λ/D. For instance, a radial change in the position of > λ/D is only seen at angular separations >5 λ/D using a typical filter with ∆λ/λ ≈ 0.2. Therefore, SDI also suffers from selfsubtraction effects and is not ideal for imaging close-in point sources or pole-on circumstellar disks. Additionally, Keck/NIRC2 does not have an integral field spectrograph for this purpose, making contemporaneous SDI impossible.

Reference star differential imaging

Reference star differential imaging (RDI; Lafrenière et al. 2009;Soummer et al. 2011;[START_REF] Gerard | Proc. SPIE[END_REF]) is an alternative approach that uses images of other stars to build a model of the stellar PSF. RDI is a commonly used observational strategy for Hubble Space Telescope imaging observations of debris disks [START_REF] Golimowski | [END_REF]Schneider et al. 2009;Schneider et al. 2014;[START_REF] Choquet | [END_REF]) and vortex coronagraph observations with the Hale telescope at Palomar Observatory (Mawet et al. 2010;Serabyn et al. 2010;Mawet et al. 2011) where ADI is not possible due to its equatorial mount. In addition, RDI will very likely be applied for highcontrast imaging with future space telescopes, including the James Webb Space Telescope (JWST; [START_REF] Green | Proc. SPIE[END_REF]) and the Wide Field Infrared Survey Telescope (WFIRST; Spergel et al. 2015;[START_REF] Bailey | Proc. SPIE[END_REF] as well as the HabEx (Gaudi et al. 2018) and LUVOIR (The LUVOIR Team 2018) mission concepts. The latter two may push to very small angular separations using vortex or other small inner working angle coronagraphs (Ruane et al. 2018;[START_REF] Pueyo | Proc. SPIE[END_REF].

Here, we demonstrate how RDI with Keck/NIRC2 mitigates practical issues associated with ADI (and potentially SDI on other instruments). We demonstrate the benefits of the RDI observing strategy for groundbased imaging studies, especially the detection of point sources at small angular separations and mapping of scattered light from circumstellar disks. In the following sections, we detail the RDI method and strategy (Sec. 2) using three example observing nights that illustrate the benefits of RDI for point source detection (Sec. 3) and for imaging of circumstellar disks (Sec. 4) with Keck/NIRC2. Section 5 discusses our findings and the limitations of RDI, including error sources and potential artifacts. Section 6 summarizes our conclusions.

RDI METHOD AND STRATEGY

RDI requires the observer to decide which reference frames and algorithms to use in order to model the stellar PSF. Reference frames may be specific reference stars purposely taken on the same night as the primary observation (e.g. Ruane et al. 2017), frames from other stars that were coincidentally observed on the same night in the same observing mode (e.g. Xuan et al. 2018), or from an archive of frames across observing programs (e.g. [START_REF] Choquet | [END_REF]. Example PSF reconstruction algorithms include subtracting a scaled version of the reference image (Schneider et al. 2009), principal components analysis (PCA; Soummer et al. 2012), and nonnegative matrix factorization (NMF; Ren et al. 2018). Xuan et al. (2018) showed that using all frames from a given night (excluding the target of interest) as reference frames and applying PCA improves the sensitivity to point sources at small angular separations with respect to ADI for many of the targets in our database of Keck/NIRC2 vortex coronagraph observations. In the following, we will discuss improvements to this strategy by combining this approach with frame pre-selection.

For each observation, we apply basic pre-processing steps (see Appendix A for details) and PCA to estimate and subtract the starlight from the images using the Vortex Image Processing (VIP) software package (Gomez Gonzalez et al. 2017). Following Soummer et al. (2012), a single frame may be written as X = I + A, where I is the stellar PSF after the coronagraph (a speckle field) and A is the planet PSF. The stellar PSF is

I = K k=1 X, Z (k) Z (k) , (1)
where {Z (k) } k=1,...,N are a basis set of images derived from a series of N reference frames {R (k) } k=1,...,N , X, Z (k) is the projection of the frame onto the kth basis image, and K represents the number of basis vectors used for the PSF model with K ≤ N . Here, the basis modes are computed using PCA. In the case of point sources, we empirically determine the number of basis vectors, K, that provides the best sensitivity to fake companions injected in post-processing. Whereas in ADI-processed images the value of K can drastically affect the sensitivity to point sources (Pueyo 2016) and the occurrence of disk artifacts (Milli et al. 2012), RDI results tend to be relatively insensitive to the choice of K (Soummer et al. 2012). Rather, the choice of which reference frames to use to model the stellar PSF has a dominant effect on the quality of the final image.

We explored a number of different methods for reference frame pre-selection. In each case, we compared the potential reference frames, X i , with the temporal median of the science frames, M , over a speckle noise dominated region of the image (typically within 5-10 λ/D). We assigned a score to each reference frame using three metrics: the mean square error (MSE), the Pearson correlation coefficient (PCC), and the structural similarity index metric (SSIM) as defined in Wang et al. (2004). The MSE of the kth reference frame is given by

MSE (k) = 1 N pix Npix i=1 X (k) i -M i 2 , (2)
where N pix is the number of pixels in the comparison region. The PCC is

PCC (k) = cov(X (k) , M ) std(X (k) )std(M ) , (3)
where cov(.) and std(.) represent the covariance and standard deviation. Specifically, the covariance is given by

cov(P, Q) = 1 N pix -1 Npix i=1 (P i -P )(Q i -Q), (4)
where P and Q are the images being compared, with pixel-wise means P and Q. The standard deviation is

std(P ) = 1 N pix -1 Npix i=1 (P i -P ) 2 .
(5)

The final metric we consider is the mean SSIM:

SSIM (k) = 1 N pix Npix i=1 L (k) i C (k) i S (k) i , (6)
where L, C, and S are the luminance, contrast, and structural terms, which are computed over a FWHM×FWHM window centered on pixel i. The luminance term is the contrast term is

L i = 2 X M + c 1 X2 + M 2 + c 1 , (7)
C i = 2 std(X)std(M ) + c 2 std(X) 2 + std(M ) 2 + c 2 , (8)
and the structural term is

S i = cov(X, M ) + c 3 std(X)std(M ) + c 3 . (9)
c 1 , c 2 , and c 3 are small constants chosen to prevent instability when the denominator is otherwise close to zero. Qualitatively, L i is the relative change in luminance and responds similarly to MSE to differences in the pixel counts. S i is almost identical to the PCC. Thus, SSIM can be thought of as a mixture between the MSE and PCC metrics. In the following section, we investigate the performance gains achieved by modeling the stellar PSF using the principal components (PCs) of a set of reference frames selected based on their MSE, PCC, and SSIM values, as compared to the median science frame.

POINT SOURCE DETECTION

To illustrate the utility of and trade-offs associated with using RDI for point source detection, we revisit Keck/NIRC2 observations taken during commissioning of the vortex coronagraph mode on UT 2016 Apr 13. Table 1 gives the full list of observations, consisting of 11 targets observed over the course of a full observing night for various scientific purposes. The targets span a large range of elevations, spectral types (Spec. type), and magnitudes, which we list in V and Wise W 1 bands (3.4 µm; Wright et al. 2010). The observations also vary in the number of frames (N frames ) and total integration times (t int ). Hinkley et al. (2015) previously identified low-mass stellar companions orbiting at ∼20 au from two of the observed targets, HIP 79124 and HIP 78233, using aperture masking interferometry. HIP 79124 and HIP 78233 are classified as A0V and F0V stars, which reside in the Upper Scorpious subgroup of the Scorpius-Centaurus (USco) association at a distance of ∼120-150 pc (de Zeeuw et al. 1999) and, thus, have an estimated age of 5-20 Myr (Pecaut et al. 2012;Song et al. 2012). HIP 79124 is a triple system with an additional companion, HIP 79124 B, at ∼1 from HIP 79124 A (Lafrenière et al. 2014). RDI allows us to directly image their close-in companions, HIP 79124 C and HIP 78233 B, even though their angular separation from the host star is <0. 2, or 2 λ/D in L band. We use these previously reported astrometric measurements along with data presented in this study (Table 2) to check for common proper motion and confirm that HIP 78233 B, HIP 79124 B, and HIP 79124 C are all gravitationally bound to their host star.

Figure 2 shows the MSE, PCC, and SSIM for each possible reference frame from UT 2016 Apr 13, compared to the median frame in the HIP 79124 and HIP 78233 observations. As predicted, the MSE and PCC respond to different attributes, whereas the SSIM shows similar features to both MSE and PCC. In the case of HIP 79124 and HIP 78233, the brighter stars (HIP 66704, DX Leo, HIP 41152) had considerably worse MSE and SSIM. Also, HIP 79124 and HIP 78233 were mutually poor reference stars for one another according to the PCC and the SSIM despite being observed consecutively. This is likely due to variable observing conditions at that point in the night.

Figures 3 and4 show the ADI and RDI images of HIP 79124 C and HIP 78233 B, respectively. HIP 79124 C was previously imaged by Serabyn et al. (2017) using RDI. However, we have recovered HIP 78233 B for the first time. We show both examples here to demonstrate that, using PCA, the companions are detected at higher significance with RDI compared to our ADI reduction.

In order to determine which reference frames to use, we ordered the reference frames according to the MSE, PCC, and SSIM and varied the number of frames from which we derive the PCs (see Fig. 5). In each case, we compute a forward model of the companion's PSF (see e.g. Lagrange et al. 2009;Soummer et al. 2012). We subtracted a copy of the off-axis PSF at the location of the planet in each pre-processed science frame, varied its location and brightness, and repeated the PCA reduction until the values were minimized in the final images in a 2×FWHM radius about the companions position using a downhill simplex algorithm. The forward model is defined as the difference between the original PCA reduction without injected or subtracted companions and the best-fit residuals.

The signal-to-noise ratio, S/N , is defined as the ratio between estimated flux, F , and the error in the flux fit, σ F , as calculated using the method outlined in Appendix B. We find that the MSE and SSIM metrics lead to higher S/N for these objects compared to PCC, with SSIM offering a slight improvement. The optimal number of frames according to SSIM is 237 (out of 336) and 202 (out of 323) for HIP 79124 C and HIP 78233 B, respectively. We also optimized K in this fashion for each number of reference frames considered, but found that K has much less of an influence on the S/N . At the optimal number of reference frames, we used 11 and 15 PCs for HIP 79124 C and HIP 78233 B, respectively. Figures 3 and4 show the result using the optimal combination of parameters. The S/N has a steep drop when using more than ∼230 frames ordered by MSE and SSIM, which occurs when the reduction includes frames from the brightest stars (HIP 66704, DX Leo, HIP 41152). This is expected since we found previously in Xuan et al. (2018) that the contrast in Keck/NIRC2 images is strongly correlated with stellar magnitude. Table 2 reports the astrometry and photometry values derived from the highest S/N case. There are a few important differences between the values in Fig. 5 and the final fits. First, we used the flux and position we retrieved using all reference frames to compute the fitting error in Fig. 5. In the final fits, we used the actual flux and position retrieved at the optimal number of reference frames and PCs for the injected, fake companion. Second, we applied a correction to the retrieved flux values to account for the throughput of the coronagraph, which was 81% and 68% for HIP 79124 C and HIP 78233 B, respectively, according to our optical model. In the end, we arrived at lower S/N values for each companion than the estimates in Fig. 5. The final S/N was better with RDI than ADI for both stars. Specifically, for HIP 79124 C, the S/N in RDI was 36 whereas the best we could achieve in ADI was 7.9 with K = 2 and a PA rotation of 17 • . Similarly, for HIP 78233 B, the S/N in RDI was 8.6, but the best we achieved in ADI was 5.1 with K = 6 despite having a PA rotation of 32 • .

Comparing to Hinkley et al. (2015), the change in astrometry since 2010 is consistent with orbital motion, leading to the first confirmation that HIP 78233 B is bound. We updated the distances to HIP 79124 and HIP 78233 based on Gaia DR2 (Bailer-Jones et al. 2018), which gave 132±1 pc and 256±4 pc for the respective stars, rather than the distances of 123 pc and 145 pc used in Hinkley et al. (2015). Our measured flux ratio for HIP 79124 C, ∆L = 4.01 ± 0.03, is inconsistent with both Hinkley et al. (2011) and Serabyn et al. (2017) by a small margin, who report ∆L = 4.3 ± 0.1 and ∆L = 4.2 ± 0.1, respectively. The increased flux estimate and updated stellar distance predicts a higher mass for both companions. Assuming the BT-Settl evolutionary model (Allard et al. 2012) and an age of 10 Myr, we estimate minimum masses of 199±4 M Jup and 198±19 M Jup , which are well above the hydrogenburning limit.

For each point source companion, we use previously reported astrometry along with new data presented in this study to constrain the orbits of each companion around their host stars. Although two data points do not provide a well-constrained orbit, the long baseline provides some constraint on the actual semimajor axis and eccentricity. We use the "Orbits of the Impatient" (OFTI; [START_REF] Blunt | [END_REF]) algorithm as implemented by the orbitize software package 1 . OFTI is a Bayesian rejection-sampling method that excels at finding constraints on long period orbits with sparse data points. We determine the median semimajor axis and 68% confidence intervals to be 37 +31 -12 au for HIP 78233 B, 157 +82 -37 au for HIP 79124 B and 26 +20 -8 au for HIP 79124 C. We also found that the data favors lower eccentricity orbits over high eccentricity orbits. Details of these orbit fits and the asymmetrical posteriors are discussed in Appendix C. The posteriors are also available online as supplementary data for use in other studies.

In the two examples above, we have demonstrated that RDI offers significantly improved S/N for companions at angular separations of 2 λ/D and that frame pre-selection is critical for maximizing the performance achieved with RDI. However, previous versions of the data reduction pipeline used by our team for Keck/NIRC2 vortex coronagraph observations (Xuan et al. 2018) used all available reference frames from a given night by default. Although there is often an S/N improvement at small separations using all frames over ADI, we are working towards using SSIM or similar metrics to perform the pre-selection of frames in a robust way and reprocessing archival NIRC2 data to reveal new point sources at close separations.

CIRCUMSTELLAR DISK IMAGING

In addition to improving the sensitivity to point sources at small angular separations, RDI is also advantageous for disk imaging. In the following two examples, we use Keck/NIRC2 observations of circumstellar disks. In each case, reference stars were specifically chosen to be close in elevation when observed and to have similar brightness in R and L bands, where the wavefront sensor and science camera are sensitive. [START_REF] Grady | [END_REF][START_REF] Benisty | [END_REF][START_REF] Reggiani | [END_REF]Ren et al. 2018). ADI images with Keck/NIRC2 from [START_REF] Reggiani | [END_REF] showed evidence of a potential third spiral arm.

Here, we reprocessed one of the epochs presented in Reggiani et al. (2018) using reference stars observed directly before and after the previously published ADI sequence. Table 3 gives the list of observations we used in the reduction. For ADI, we used all of the science frames and K = 6. For RDI, we used similar methods as the previous section keeping the best 70% (50 out of 75) of reference frames according to the SSIM metric and used the projection of 25 out of 50 PCs in the PSF model. However, since the reference stars were a good match to the science targets and conditions were stable during these observations, the variation in MSE, PCC, and SSIM was similar to that of one of the wellcorrelated observations in Fig. 2. Therefore, the RDI image of the disk is visibly unchanged, barring some background variation, whether we use 10% more or less frames or change the number of PCs by ∼20%, while the ADI image is very sensitive to number of PCs.

Figure 6 shows the ADI and RDI images of MWC 758. In addition to the famous North and South spirals (labeled '1' and '2' in Fig. 6b), we recovered the proposed third spiral arm feature (labeled '3') with both ADI and RDI for the first time. The ADI image is dominated by self-subtraction effects, while the RDI image is likely more representative of the true morphology of scattered light features. Although we have confirmed that the proposed third spiral arm was not an ADI artifact, it is unclear from the RDI image whether the feature is truly a third spiral arm or a continuation of the north arm wrapping almost 360 • about the star. Nevertheless, this observation demonstrates that RDI is beneficial for mapping scattered light distributions in circumstellar disks without the sharpening of azimuthal features seen in the ADI images.

J16042165-2130284

We present an observation of the young stellar object 2MASS J16042165-2130284 (from here on 'J1604'). J1604 is a K2 star at a distance of 150 pc (Bailer-Jones et al. 2018) and a likely member of USco (Carpenter et al. 2014) suggesting an age of 5-20 Myr (Pecaut et al. 2012;Song et al. 2012). J1604 is surrounded by a nearly pole-on (7 • inclination; Mathews et al. 2012) pre-transitional disk (Espaillat et al. 2007), the most massive known (∼ 0.1 M Jup ) in USco (Mathews et al. 2012). Several previous works identified a ring-like disk in scattered light (Mayama et al. 2012;Pinilla et al. 2015;Canovas et al. 2017;Pinilla et al. 2018). This disk has received considerable attention owing to its variability (Pinilla et al. 2015(Pinilla et al. , 2018) ) and evidence of planetinduced dust filtration (Rice et al. 2006;Canovas et al. 2017).

The observing list from UT 2017 May 10 (see Table 4) consists of seven observations, including J1604 and a reference star, 2MASS J19121875-2137074 (from here on 'J1912'), specifically selected prior to the observing night because its similar brightness and elevation. The remaining five targets were observed as part of an ongoing survey of M stars (PI: Mawet). The first three M star targets were significantly worse reference stars according to their SSIM values because they were much fainter (W 1 > 9) than J1604 (W 1 = 7.61). The other two M stars had SSIM values close to that of J1912.

Figure 7 shows the ADI and RDI images of J1604. Since the disk is nearly rotationally symmetric, the ADI image is consistent with noise for all values of K. However, the ring of scattered light is clearly visible in the (2018) show azimuthal dips attributed to shadowing by an inner disk. We see similar features labeled '1', '2', and '3' in Fig. 7b. These dips change with time indicating that the inner regions are highly dynamic, which may be evidence of a close-in massive companion. In the southern part of the disk, we estimate the inner and outer edges of the bright scattering feature appear at 440 and 500 mas, respectively. Compared to J-band (1.2 µm) observations taken with VLT/SPHERE by Pinilla et al. (2018) who measured a scattered light peak at 430 mas, the scattering surface appears slightly further from the star in L band (3.8 µm), which may be evidence of spatial segregation of dust particle sizes and/or lower opacities at longer wavelengths. In future work, we will fit a forward model of the disk using radiative transfer and further investigate the wavelength dependence of the scattered light surface and shadowing effects in the disk (Wallack et al., in prep.), all of which would not be possible using ADI.

DISCUSSION

The above sections exemplify that (1) combining the small inner working angle vortex coronagraph and RDI is a sensitive technique for discovering faint companions ∼100 mas from the star and imaging scattering light from circumstellar disks with Keck/NIRC2 and (2) optimizing the reference frame selection process leads to significant improvement in the quality of the stellar PSF subtraction. In this section, we describe the noise characteristics, artifacts, and limitations unique to this observing mode.

Treatment of false positives

Given the typical number of reference frames available per night and the limited diversity they provide, subtraction of the stellar speckle field is imperfect, and the residual speckles in RDI images tend to obey a highly non-Gaussian noise distribution [START_REF] Goodman | Laser Speckle and Related Phenomena[END_REF]Aime & Soummer 2004), especially within a few λ/D of the star. However, the standard detection methods used in high-contrast imaging, including this work, estimate the standard deviation of the speckle noise, σ n , as a function of angular separation from the star and set a fixed detection threshold (typically 5σ n ), which leads to a substantial increase in the number of false positives compared to what is expected from normally distributed noise. In addition, since the spatial scale of speckles is roughly the FWHM of the off-axis PSF, σ n is estimated using only a few independent samples at small angular separations [START_REF] Mawet | [END_REF]. It is therefore especially challenging to differentiate between bright speckles and true companions in a single observation and to set robust upper limits on the brightness of unseen companions (see e.g. Ruffio et al. 2018) within a few λ/D of the star. More work is needed to accurately model the residual speckle noise distribution given a set of science and reference frames.

Typical artifacts

In addition to the non-Gaussian properties of the stellar speckle noise, there are other artifacts that may lead to false positives. For instance, several dust spots on the vortex focal plane mask in Keck/NIRC2 resemble point sources in RDI images because their brightness is not correlated with that of the star and background and therefore they are not well subtracted in post-processing. This requires the observer rule out detections corresponding to the location of the dust. We find that having a large amount of PA rotation helps blur these effects in the de-rotation step.

The RDI approach assumes that stellar speckle noise dominates in both the science and reference frames. Generally speaking, RDI tends not to provide improved sensitivity at small angular separations on targets where spatially resolved circumstellar material dominates the speckle noise. In fact, projecting reference frames of point sources onto the science frames in an attempt to model the stellar PSF modifies the appearance of the circumstellar disk and generates additional speckle noise, which may obscure the true scattered light features and generate point-like false positives. Spatially resolved sources also have a mismatch in the radial size of the main stellar residuals that leaves behind a ringlike residual at the outer edge of the central lobe of the stellar PSF. A similar effect may occur for objects with vastly different colors, though we have not come across a clear example where chromatic effects dominate.

Lastly, there are speckles that appear in NIRC2 images at ∼ 7 λ/D from the star whose brightness depends on the telescope elevation. We attribute these speckles to segment-to-segment phasing errors that depend on the direction of the gravity vector with respect to the primary mirror and the azimuthal position depends on the position of the rotator. For this reason, it may be beneficial to use reference stars at similar declination to the science target and similar instrument settings, including the rotator angle.

Sources of night-to-night variance

The methods outlined above may be generalized to include reference frames from multiple observing nights. However, although the stellar PSF is relatively stable throughout an observing night, it tends to vary considerably from night-to-night for a number of reasons.

The calibration of static wavefront error within NIRC2 uses a phase retrieval algorithm to reconstruct the wavefront from defocused images of the internal source. The solution changes on a nightly basis. Additionally, wavefront errors due to the primary mirror must be sensed and corrected by the Shack-Hartmann wavefront sensor in the Keck AO system (Wizinowich et al. 2000), but segment piston errors (or the "terrace" modes) are not seen by the wavefront sensor and are therefore not corrected. This issue may be mitigated by the pyramid wavefront sensor under development at Keck Observatory [START_REF] Bond | Proc. SPIE[END_REF].

For these reasons, reference frames are generally not well correlated from night-to-night. Nevertheless, there may be correlations with instrument settings (e.g. the rotator angle and wavefront sensor parameters) or observing conditions that can assist in frame pre-selection, which will be investigated in future work.

CONCLUSION

Using the examples of HIP 79124 C and HIP 78233 B, we have demonstrated that RDI offers improved detection of point sources at small angular separation, by up to a factor of 5 with respect to ADI, further unlocking the potential of the vortex coronagraph on Keck/NIRC2. Furthermore, the observations of MWC 758 and 2MASS J16042165-2130284 show that RDI is beneficial for imaging of circumstellar disks because it preserves the morphology, which allows observers to directly interpret the imaged scattered light distribution. We find that frame pre-selection using image comparison metrics, such as the SSIM, significantly improves the performance of RDI, particularly for point source detection, astrometry, and photometry. RDI may help other current and future ground-based instruments achieve better performance at small angular separations, especially when used in conjunction with small inner working angle coronagraphs, and perform more efficient surveys by relaxing PA rotation requirements. RDI may also be the primary strategy for high-contrast imaging with future space missions, including JWST, WFIRST, HabEx, and LUVOIR.

G. Ruane is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1602444. The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration (NASA). The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. In this section, we outline the details of our observations and the image processing steps used in our team's pipeline (Xuan et al. 2018). All observations were carried out using the Keck II telescope, natural guide star adaptive optics, the NIRC2 instrument, and the vortex coronagraph mode in L band (3.7 µm). The angular resolution was 0.08 and the plate scale was 0.01 per pixel (Service et al. 2016). The field rotator was set to vertical angle mode, such that the telescope pupil tracks the elevation axis, to enable ADI. In addition to the science frames (see Tables 1, 3, and4), we took several images of the off-axis PSF with a discrete integration time (DIT) of 0.008 seconds and 100 coadds as well as images of the sky background with integration times matching that of the science and off-axis PSF frames. The sky, off-axis PSF, and background frames were taken approximately every 30 min during the observing sequence. The alignment of the star and the center of the vortex focal plane mask was maintained by the QACITS tip-tilt control algorithm (Huby et al. 2015(Huby et al. , 2017)).

In the pre-processing stage, bad pixels identified in the dark frames and sky flats were replaced by the median of neighboring values in all frames. The sky flat was the median of 10 images of a blank patch of sky with the coronagraph focal plane mask removed (DIT of 0.75 s, 10 coadds each). We then subtracted sky background frames, taken with the coronagraph in place, from each frame individually using a scale factor to account for background variability. The frames were centered based on the position of the optical vortex core in the median of the science frames and each of individual science and reference frames were registered using the peak of the cross correlation with the median science frame. Our pipeline crops the raw 1024×1024 pixel frames to 587×587, which is the largest allowable square frame centered on the vortex focal plane mask.

We applied PCA (Soummer et al. 2012) to estimate and subtract the stellar contribution from each image using the VIP software package (Gomez Gonzalez et al. 2017). In the case of ADI, the PCs of a subset of pixels in the science frames make up the stellar PSF model. Although our standard pipeline uses a set of default annular reduction regions, we generally modify these regions to bound the companion or disk of interest on a case-by-case basis. For RDI, we derive the PCs and the stellar PSF model only from the reference frames. Each frame is de-rotated such that North is in the vertical direction and the ADI/RDI image is taken to be the temporal median of the de-rotated frames. We repeat these steps for all possible number of PCs. Later, an optimal number of PCs may be defined by maximizing the detection significance of a real or injected companion, depending on the scientific goal. For instance, an observer may wish to use a different number of PCs for point source detection and disk imaging.

B. ERROR BARS FOR ASTROMETRY AND PHOTOMETRY OF COMPANIONS

To estimate the errors in the position and flux estimates, with the best-fit companion subtracted from the data, we re-injected the PSF into the pre-processed science frames at the same separation and brightness, but varied the azimuthal angle and retrieved the photometry and astrometry using the same method for each injected companion tracing a full circle about the star. The step size in azimuthal angle was 360 • /(2πr), where r is the radial position in units of the FWHM, to ensure that each fit was performed at the location of an quasi-independent speckle sample. We took the standard deviation of the measured flux and position of the injected companions as the error on each parameter. The signal-to-noise ratio, S/N , is defined as the ratio between estimated flux, F , and the uncertainty in the flux, σ F .

C. ORBIT FITTING OF POINT SOURCE COMPANIONS WITH ORBITIZE AND OFTI

In this section, we provide more details on the orbit fit and posteriors on the point source companions from the systems discussed in Section 3. Table 5 summarizes the astrometry used in our fit. We use stellar distances from Bailer-Jones et al. (2018) of 132±1 pc and 256±4 pc for HIP 78233 and HIP 79124, respectively. We use total system mass estimates from Hinkley et al. (2015), corresponding to 1.7 ± 0.1 M Sun for HIP 78233 and 2.48 ± 0.45 M Sun for HIP 79124. We choose uninformed priors for all parameters. Namely, for the semimajor axis, we choose a prior probability distribution that is uniformly linear in log space (Jeffreys Prior); for inclination, we choose a prior probability distribution corresponding to p(x) ∝ sin(x); and for all other parameters, we choose a linearly uniform prior probability distribution. We fit the companions as point masses and for the HIP 79124 triple system, we treat each companion independently as their large mutual separation would make interactions between the B and C companions Using orbitize's OFTI algorithm, we sample 100,000 orbits for each of the three point source companions and compute posterior probabilities. Figure 8 shows the marginalized posterior probabilities for orbital semimajor axis, eccentricity, and inclination. These figures also include the prior probability to demonstrate how well the data constrains the orbits. The marginalized posterior on semimajor axis is the best constraint on the companion's orbit size when accounting for all available astrometry. For HIP 78233 B and HIP 79124 C, while the eccentricity and inclination posteriors are not well constrained and are not very different from the prior. The eccentricity posterior distribution does show that lower values are more likely than higher values and the inclination posterior distribution disfavors a completely edge-on system (inclination of 90 • ). For HIP 79124 B, the majority of the eccentricity distribution is at values below 0.5 and the inclination distribution strongly favors an inclination less than 90 • , indicating that the object is orbiting counter-clockwise from our viewing angle. These posteriors are available as part of this article's online data for use in further calculations as a table of orbit samples.

Figure 1 .

 1 Figure 1. Parallactic angle rotation during a one hour observation from Maunakea, Hawaii (latitude 19.8• N) as function of target declination δ. The dashed blue line represents an observation initiated exactly 30 min before the object crosses the meridian. The thin black lines represent observation windows shifted by ∆t=15 min. The other thick dashed lines indicate the observation with maximum allowable timing error to achieve a PA rotation that moves a planet by an arclength equal to 1.0×FWHM and 0.5×FWHM at an angular separation of λ/D. Specifically, achieving 60 • and 30 • PA rotations at the most favorable δ requires the timing error to be <36 min and <68 min, respectively.

Figure 2 .

 2 Figure2. The MSE (shown as negative logarithm), PCC, and SSIM for all frames from UT 2016 Apr 13 (grouped by target in an arbitrary order) compared to the median science frame of the science observations: (a) HIP 79124 and (b) HIP 78233. MSE and PCC responds to differences in the image counts and structure, respectively; SSIM is correlated with both.

Figure 3 .

 3 Figure 3. Observation of HIP 79124 C. (a) Image after subtracting stellar PSF using ADI. (b) Best forward model of the point source. (c) Residuals after subtracting the point source from the data. (d)-(f) Same as (a)-(c), but using RDI for the stellar PSF subtraction. North is up and East is left. The companion was retrieved at higher S/N in RDI compared to ADI. (a)-(c) and (d)-(f) are shown on the same scale.

Figure 4 .

 4 Figure 4. Same as Fig. 3, but for HIP 78233 B.

Figure 5 .

 5 Figure5. Estimated photometric S/N = F/σF achieved with RDI versus the number of reference frames for the three frame selection methods. For the sake of finding the best set of reference frames, we injected a fake companion into the RDI residual cube with the best-fit flux and separation obtained using the maximum number of reference frames.

  4.1. MWC 758 MWC 758 (HD 36112) is a young (3.5±2.0 Myr; Meeus et al. 2012) stellar object in the Taurus star-forming region (distance of 160 pc; Bailer-Jones et al. 2018) surrounded by a prominent protoplanetary disk with at least two massive spiral arms seen in scattered light observations

Figure 6 .

 6 Figure 6. The circumstellar disk around MWC 758 after subtracting the stellar PSF using (a) ADI and (b) RDI. Labels '1,' '2,' and '3' indicate the prominent North and South spiral arms and the potential third arm, respectively.RDI reduction. For the RDI image in Fig.7b, we used the best 88 reference frames, according to their SSIM values, out of the 150 available and projected 44 PCs to build the stellar PSF model. Again, the morphology of the disk did not change much when changing the number of reference frames or PCs by ∼ 20% in this case. After removing poorly matching reference frames according to the SSIM, the projection coefficients for the first few PCs dominate.Polarized intensity maps of J1604 fromPinilla et al. (2018) show azimuthal dips attributed to shadowing by an inner disk. We see similar features labeled '1', '2', and '3' in Fig.7b. These dips change with time indicating that the inner regions are highly dynamic, which may be evidence of a close-in massive companion. In the southern part of the disk, we estimate the inner and outer edges of the bright scattering feature appear at 440 and 500 mas, respectively. Compared to J-band (1.2 µm) observations taken with VLT/SPHERE byPinilla et al. (2018) who measured a scattered light peak at 430 mas, the scattering surface appears slightly further from the star in L band (3.8 µm), which may be evidence of spatial segregation of dust particle sizes and/or lower opacities at longer wavelengths. In future work, we will fit a forward model of the disk using radiative transfer

Figure 7 .

 7 Figure 7. The 2MASS J16042165-2130284 disk after subtracting the stellar PSF using (a) ADI and (b) RDI. Labels '1,' '2,' and '3' indicate the the position of potential shadow features in the disk.

Facilities:

  Keck:II (NIRC2) Software: VIP(GomezGonzalezetal.2017),QACITS(Huby et al. 2015, 2017), Astropy (The Astropy Collaboration et al. 2018), orbitize (Blunt et al. 2018) APPENDIX A. OBSERVING AND IMAGE PROCESSING DETAILS

Figure 8 .

 8 Figure 8. Marginalized probability density functions (pdfs) for three orbital parameters from orbit fits of three point source companions. The blue histogram shows the posterior pdf while the orange curve shows the prior pdf. For the semimajor axis pdf, the median and 68% confidence region centered on the median is labelled as a point with horizontal bars. The posterior samples used to generate these pdfs are available as part of the online data.

Table 1 .

 1 Keck/NIRC2 vortex observations in L band on UT 2016 Apr 13

	Star name	Alt. Name	RA	Dec	Spec. Type N frames tint (min) PA rot.	V	W 1 a
	HIP 41152 HD 70313	08 23 48.5 +53 13 11.0	A3V	30	10.0	12.9 •	5.54 5.21
	DX Leo	HD 82443	09 32 43.8 +26 59 18.7	G9V	30	15.0	55.0 •	7.01 5.15
	TW Hya	TWA 1	11 01 51.9 -34 42 17.0	K6V	70	35.0	21.6 •	10.5 7.10
	V1249 Cen TWA 25	12 15 30.7 -39 48 42.6	M0.5	30	15.0	7.91 •	11.2 7.26
	V1252 Cen TWA 10	12 35 04.2 -41 36 38.5	M2V	30	15.0	7.82 •	13.0 8.09
	HIP 66704 HD 119124	13 40 23.2 +50 31 09.9	F7.7V	30	15.0	17.9 •	6.32 4.92
	HIP 78233 HD 142989	15 58 29.3 -21 24 04.0	F0V	36	18.0	32.4 •	9.10 7.63
	HIP 79124 HD 144925	16 09 02.6 -18 59 44.0	A0V	23	11.5	16.9 •	7.78 6.96
	2MASS J16430128-1754274	16 43 01.3 -17 54 27.4	M0.5	30	15.0	11.8 •	12.6 8.44
	2MASS J17520294+5636278 17 52 02.9 +56 36 27.8	M3.5V	30	15.0	11.4 •	13.3 8.20
	V4046 Sgr HD 319139	18 14 10.5 -32 47 34.4	K5+K7	20	10.0	5.63 •	10.7 7.12
	Note-

a Wise W 1 mag (3.4 µm; Wright et al. 2010).

Table 2 .

 2 Best fit photometry and astrometry for the detected companions.Note-a Bailer-Jones et al. (2018). b Lower bound on mass; assumes the BT-Settl model(Allard et al. 2012) and an age of 10 Myr.

	Companion Dist. (pc) a	∆L	Sep (mas)	PA ( • )	F/σF Mass (MJup) b Proj. sep (au)
	HIP 79124 B	132±1	2.49±0.04	971.5±1	100.61±0.03 >31	630±17	128±1
	HIP 79124 C	132±1	4.01±0.03	192±2	246.5±0.8	36	199±4	25.3±0.3
	HIP 78233 B	256±4	4.78±0.12	141±3	6±1	8.6	198±19	36.0±0.9

Table 3 .

 3 Keck/NIRC2 vortex observations in L band on UT 2016 Oct 24.

	Name	Alt. Name	RA	Dec	Type N frames tint (min) PA rot.	V	W 1
	TYC 1845-2048-1 HD 284115 05 02 41.3 +24 07 19.5	K2	25	16.7	23.1 •	8.13 3.98±0.10
	MWC 758	HD 36112	05 30 27.5 +25 19 57.1 A8V	80	53.3	173 •	8.27 4.60±0.08
	TYC 1867-221-1	HD 249769 05 58 47.2 +25 15 28.8	K7	50	33.3	171 •	8.83 4.14±0.09

Table 4 .

 4 Keck/NIRC2 vortex observations in L band on UT 2017 May 10.

	Name	RA	Dec	Type N frames tint (min) PA rot.	V	W 1
	2MASS J11110358-3134591	11 11 03.6 -31 34 59.1	M	20	15.0	6.99 •	14.4 9.35±0.02
	2MASS J11431742+1123126 11 43 17.4 +11 23 12.6	M	20	15.0	10.1 •	12.3 9.12±0.03
	2MASS J13412668-4341522	13 41 26.7 -43 41 52.2 M3.5	20	15.0	6.18 •	14	9.70±0.02
	2MASS J16042165-2130284	16 04 21.7 -21 30 28.5	K2	73	54.8	36.2 •	11.9 7.61±0.03
	2MASS J18580415-2953045	18 58 04.2 -29 53 04.5 M0V	20	15.0	7.25 •	11.8 7.86±0.02
	2MASS J19121875-2137074	19 12 18.8 -21 37 07.4	M	50	37.5	19.8 •	11.3 7.42±0.03
	2MASS J20013718-3313139	20 01 37.2 -33 13 13.9	M1	20	15.0	6.95 •	12.3 8.14±0.02

Table 5 .

 5 Input astrometry for orbit fitting via orbitize The inclination convention is that 0 • ≤ i < 90 • orbits are counter-clockwise and 90 • ≤ i < 180 • are clockwise.

	UT date	Separation (mas) Position Angle ( • )	Reference
	HIP 78233 B astrometry
	2010 Apr 26	133 ± 3	20 ± 1	Hinkley et al. (2015)
	2016 Apr 13	141 ± 64	6 ± 1	this work
	HIP 79124 B astrometry
	2008 May 25	990 ± 1	98.11 ± 0.05	Lafrenière et al. (2014)
	2016 Apr 13	971.5 ± 1	100.61 ± 0.03	this work
	HIP 79124 C astrometry
	2010 Apr 5	177 ± 3	242 ± 1	Hinkley et al. (2015)
	2016 Apr 13	192 ± 2	246.5 ± 0.8	this work
	very small.

https://orbitize.readthedocs.io/