Exocomet Orbit Fitting

Exocomet Orbit Fitting: Accelerating Coma Absorption During Transits of Pictoris


Comets are a remarkable feature in our night sky, visible on their passage through the inner Solar system as the Sun’s energy sublimates ices and liberates surface material, generating beautiful comae, dust, and ion tails. Comets are also thought to orbit other stars, and are the most promising interpretation of sporadic absorption features (i.e. transits) seen in spectra of stars such as  Pictoris and 49 Ceti. These “exocomets” are thought to form and evolve in the same way as in the Solar system, and as in the Solar system we may gain insight into their origins by deriving their orbits. In the case of Pictoris, orbits have been estimated indirectly, using the radial velocity of the absorption features coupled with a physical evaporation model to estimate the stellocentric distance at transit . Here, we note that the inferred imply that some absorption signatures should accelerate over several hours, and show that this acceleration is indeed seen in HARPS spectra. This new constraint means that orbital characteristics can be obtained directly, and the pericentre distance and longitude constrained when parabolic orbits are assumed. The results from fitting orbits to 12 accelerating features, and a handful of non-accelerating ones, are in broad agreement with previous estimates based on an evaporation model, thereby providing some validation of the exocomet hypothesis. A prediction of the evaporation model, that coma absorption is deeper for more distant transits, is also seen here.

comets: general — planets and satellites: detection — planetary systems — planet-disc interactions — circumstellar matter — stars: individual: Pictoris

1 Introduction

Comets are a well-known and important component of our Solar system, and while a scientific focus for astronomers, some are visible to the naked eye and capture the imagination of millions. These icy bodies become visibly more active as they near the Sun, providing a window into the primordial conditions in the outer Solar system (e.g. Blum et al., 2017). Comets formed beyond the Asteroid belt, and are only brought within a few astronomical units through interaction with the giant planets; in the case of shorter period comets (e.g. Jupiter Family Comets) these interactions are recent and ongoing (e.g. Duncan & Levison, 1997), while in the case of long period comets, these interactions occurred billions of years ago and resulted in the formation of the Oort cloud, the source of long period comets (e.g. Duncan et al., 1987).

While our detailed knowledge of planetary systems around other stars is relatively limited in comparison, it is clear that processes analogous to those postulated for the Solar system’s formation and evolution do occur. Most obviously, planets exist in spectacular abundance (e.g. Petigura et al., 2013), and the cometary source regions known as “debris discs” are seen around at least 30% of other stars (e.g. Eiroa et al., 2013; Sibthorpe et al., 2018). Of particular relevance here is the fact that the first clear evidence of another planetary system – the image of the edge-on disc orbiting Pictoris (Smith & Terrile, 1984) – prompted a series of papers that may mark the discovery of extrasolar comets (e.g. Kondo & Bruhweiler, 1985; Lagrange et al., 1987; Ferlet et al., 1987; Lagrange-Henri et al., 1988). These “exocomets” manifested as transient red-shifted ionic absorption lines in high-resolution ultraviolet and optical spectra; the signature of a large coma transiting the face of the star. The exocomet interpretation was (and perhaps still is) controversial, so these events have also become known by the more prosaic term “falling evaporating bodies”, or FEBs.

These discoveries were backed up by the development of a physical model, whereby material is ejected from a cometary nucleus into a coma dominated by neutral hydrogen. Observed elements such as calcium and magnesium are ionised when visible to stellar radiation, and are subsequently blown into a tenuous ion tail by radiation pressure, thus forming a parabolic front of metallic ions around the nucleus (Beust et al., 1989, 1990). The effect of radiation pressure was found to be critical to explain the observations; the closer the coma to the star, the stronger the radiation pressure and the smaller the front of ions, and therefore the shallower the absorption. While subject to systematic uncertainties from model assumptions, this property means that the distance of the coma to the star for a given transit could be estimated from the depth of the absorption line. The key modelling conclusion was that the velocity of the observed absorption corresponds to the radial velocity of the coma (Beust et al., 1990), meaning that by measuring these velocities, and estimating the distance to the star at transit, one may derive constraints on the comet orbit.

A common inference from application of the model to the data for Pictoris, and indeed which could be inferred from the data themselves, is that the tendency of the FEBs to be redshifted implies that the cometary orbits are not randomly distributed. The relative lack of blue shifted events (Lagrange-Henri et al., 1992; Crawford et al., 1998) means that most comets are approaching pericentre, with values for the longitude of pericentre (measured from the line of sight to the star), being between -50 and 130 (Beust et al., 1990; Kiefer et al., 2014).

This conclusion led to proposals for three main theories for the origins of Pictoris’ transiting comets. The first is that the events represent a single family of comets on common orbits, the aftermath of the breakup of a single larger object (Beust & Lissauer, 1994). The main issue with this theory is that we are unlikely to witness the aftermath of such an event, suggesting that mechanisms that work on longer timescales would provide more satisfactory explanations.

One such theory appeals to the Solar system, noting that the eccentricities of Asteroids whose pericentres precess at the same rate as Saturn (i.e. those in the secular resonance) increase in a way that is correlated with Saturn’s pericentre. Thus, these bodies make their closest approaches to the Sun with a preferred pericentre direction, and the same explanation was offered for Pictoris’ FEBs (Levison et al., 1994). The criticism levelled at this theory is that it is very specific, and has not been shown to generically reproduce the FEB phenomenon (see Beust & Morbidelli, 1996).

The final theory, and the one that has received the most attention, is that the 4:1 and/or 3:1 resonance with an outer planet lies within a source belt (Beust & Morbidelli, 1996; Quillen & Holman, 2000; Beust & Morbidelli, 2000; Thébault & Beust, 2001; Pichierri et al., 2017). The attraction of this theory is that objects in these resonances can be excited to star- grazing orbits with a very specific range of pericentre longitudes. A potential weakness of this theory (which also applies to secular resonances) is that the timescale to reach very high eccentricity orbits may be long relative to the time taken to destroy the comet by evaporation, limiting the ability of the resonance to produce the wide range of pericentre distances that is inferred from the models described above. Various solutions to this issue exist. One is that some comets are larger than about 30 km in radius, so can survive long enough to reach small pericentre distances (Beust & Morbidelli, 2000). Another invokes the existence of additional planets that either perturb the first planet, or perturb the comets directly, thus causing larger orbit to orbit changes in the comet’s pericentre distances (Beust & Morbidelli, 1996, 2000). One of the ways that this theory has been tested and developed is by aiming to match the derived properties of FEBs, for example their distributions of periastron distances, pericentre angles, and stellocentric distances at transit (Beust & Morbidelli, 2000).

2 Motivation

Thousands of spectra of Pictoris have been taken in pursuit of FEBs, but one aspect that has received little attention is the prospect for directly deriving the cometary orbits. While orbital characteristics have been derived based on the evaporation theory (Beust et al., 1990; Kiefer et al., 2014), these are model-dependent. That is, while the evaporation theory clearly reproduces the observations very well, it necessarily includes many assumptions for quantities that may vary from comet to comet, and it would be far preferable to derive the orbits in a less model-dependent way. Such an ability would provide further tests of the exocomet hypothesis, and more certain tests of comet origin models. Orbital constraints mean that the evaporation model could be refined, leading to a better physical understanding of the comets and their comae.

To derive an orbit from the absorption lines, enough information to constrain the orbital elements is needed (here we use pericenter distance , eccentricity , inclination , line of nodes , longitude of pericentre , and true anomaly ). All six are not necessary however, because we may assume an edge-on orbit, thus eliminating and . We may further eliminate by defining to be the angle from the line of sight to the star at transit centre (i.e. ). An assumption, that the orbits are parabolic () may be made based on the expectation that comets cannot survive more than a few periastron passages near Pictoris before evaporating. That is, the comets are not spiraling in on near-circular orbits, but are appearing at just a few stellar radii via relatively large orbit-to-orbit changes in (Beust & Morbidelli, 1996). At this point, two further independent pieces of information are needed to constrain the orbit. In previous work this has been the observed radial velocity of the absorption line, and the distance to the star at transit estimated from the evaporation model.

Given that FEBs are estimated to lie within a few tens of stellar radii when they are observed, one may ask whether acceleration could be detected in the course of a night’s observations. For an object at from Pictoris (using and ), the acceleration due to gravity is 2 m s, or 7 km s h. The magnitude of this acceleration is within easy reach of high-resolution spectrographs such as the High Accuracy Radial velocity Planet Searcher (HARPS), which has a spectral resolution equivalent to a few km s. Therefore, given that absorption features inferred to be caused by absorption of material at less than ten stellar radii are regularly seen, the expected accelerations should be present and readily detectable in sequences of spectra spanning more than about one hour. Non-detection of acceleration where sufficient temporal coverage exists would be a serious issue for the exocomet hypothesis.

Aside from the expectation of acceleration given the likely orbits, hints that it should be detectable can be seen in previous works; Ferlet et al. (1987) show a series of spectra taken over several hours, in which a redshifted line clearly accelerates by about 5 km s, and temporal variation was also seen by Lagrange et al. (1996) and Petterson & Tobin (1999) With the detection of such an acceleration, the distance derived from the evaporation model would no longer be needed to constrain a parabolic orbit (or could be used to refine the evaporation model).

3 Detection of Exocomet Acceleration

As a nearby (19pc) naked eye A6V star (), Pictoris is observable at high signal to noise ratios with a variety of instruments. One particular focus has been high resolution () spectra, which in addition to FEBs, reveals a relatively stable circumstellar calcium absorption line (e.g. Hobbs et al., 1985) and stellar pulsations (Koen et al., 2003), and yields constraints on the possible orbits of unseen planets (Galland et al., 2006) and on the mass of Pictoris b (Lagrange et al., 2012).

Specifically, the latter study obtained about 1000 HARPS spectra between 2004 and 2011. In the time since, observations have continued, and in late 2017 all 2642 public HARPS spectra of Pictoris were downloaded from the ESO archive for use in this study. These spectra have been processed by the HARPS Data Reduction Software, and have been shown by Kiefer et al. (2014) to be of sufficient quality for work related to FEBs; the spectra have a stable wavelength calibration that far exceeds our needs here, and the relative flux calibration across the spectral lines of interest and over the years since 2004 is high enough that variations at the 1% level are easily detected (see for example their Extended Data Fig. 2).

Figure 1: Examples of five HARPS spectra, showing the calcium K and H lines centred at the radial velocity of Pictoris. The dashed line shows the stellar reference spectrum. The absorption at zero km s is thought to originate in a stable circumstellar disc, while the varying components are attributed to exocomet comae. Features at both low and high velocity are seen, with those at low velocity tending to be deeper.

Here, we focus on the calcium K & H lines at 3933.66 and 3968.47 Å. The only post-processing applied to the spectra here is to extract the region within 500 km s of the calcium lines, and apply a multiplicative normalisation in the wings of the line so that all spectra are on the same (arbitrary) flux scale. At this stage a “stellar” reference spectrum near each line is computed; the method follows very closely that described by Kiefer et al. (2014), which used the highest flux values minus an estimated noise level across all spectra to derive a spectrum free of absorption. An additional step added here is interpolation across the circumstellar line, which is always present but varies in width/velocity slightly, so that reference divided (or subtracted) spectra do not show spurious variations that might be interpreted as FEB activity near zero velocity.

An example of five “randomly” chosen spectra is shown in Figure 1, where the stellar systemic velocity has been set to 21 km s.1 The obvious feature of all spectra is that the calcium lines are broadened by a few hundred km s by the fast stellar rotation, and in addition every spectrum shows the stable circumstellar line and several other narrow (few to a few tens of km s) and/or broad (few tens to hundreds of km s) absorption features; the FEBs. The stable circumstellar line is present in all spectra, and is seen to move and/or change in width at the km s level. Spectra generally appear different night-to-night, because the transit time for comets is of order hours. Absorption features are generally present in both Ca lines, though are deeper and more apparent in K than H unless the absorption is optically thick (i.e. the optical depth can be estimated from the line ratios, e.g. Lagrange-Henri et al., 1992; Petterson & Tobin, 1999).

Figure 2: Example of an accelerating Ca K line absorption feature, seen at 60-80 km s. The upper panel shows the level of absorption on a log stretch, and the lower panel shows 1d spectra that have been divided by the reference and offset vertically for clarity. The 1d spectra have been binned by a factor of two in time, and four in velocity, relative to the image in the upper panel. Two strong lines at and 25 km s, plus a weaker one moving from 60 to 80 km s, are visible, and the latter is the accelerating absorption feature. In the image this moving feature is seen as a faint stripe moving up and to the right over time. In the series of spectra the feature appears as a series of small absorptions that move up and to the right over time. Both panels illustrate clearly that this absorption feature is moving over the observation sequence.

To detect accelerations we start with the graphical technique of plotting temporal sequences of spectra as images, an example of which is shown in the upper panel of Figure 2. The spectral dimension is simply that output by the HARPS reduction pipeline, but spectra are binned temporally; each vertical pixel is 1/250th of a day (about 6 minutes), so rows where multiple spectra fall are averaged. On nights with continuous monitoring, bins typically contain five spectra. A weak temporal variation is present on most nights, with a period of approximately 0.6 hours; this has been removed by subtracting a smoothed image with the absorption lines removed (see section 4). A final feature present is pulsation, whose primary signal is diagonal stripes with a gradient of approximately 120 km s h; it is strongest outside the deepest part of the lines (200 km s at about the 1% level) so not visible in Figure 2, and not corrected for. The pulsation signal is much weaker than the accelerating features considered here, and has a greater velocity gradient than all fitted features, so is unlikely to affect the results.

An alternative way to visualise the results is simply to plot the spectra, as is done in the lower panel of Figure 2. This panel shows the same spectra as in the upper panel, but binned temporally by a factor of two, and spectrally by a factor of four. These spectra can highlight temporal variations of absorption lines by creating the illusion of a surface; static absorption lines result in a series of features that are directly above each other (i.e. a vertical ‘valley’), while those that accelerate show features that move in velocity with each successive spectrum (i.e. a diagonal valley).

Figure 2 shows three absorption features, two at 0 and 25 km s that do not accelerate, and one at 60-80 km s that does. The accelerating feature is seen as either a stripe in the upper panel, or a rightward-moving dip in the lower panel.

Figure 3: Accelerating features identified for modelling, indicated by transparent dashed lines. The modified Julian date at the start of each observing sequence is given. Each panel is the same as the upper panel in Figure 2 and has the same scale, where lighter colours indicate stronger absorption.
Figure 4: Spectra of accelerating features identified for modelling. The modified Julian date at the start of each observing sequence is given. Each panel is the same as the lower panel in Figure 2. As discussed in the text, weak accelerating features are not always easily discernible.

A set of accelerating features were identified in spectra by eye using images similar to Figure 2, from 35 nights’ data with more than about an hour of continuous observation, and these are shown in Figures 3 and 4. The dates on which these sequences started are shown in each panel, and given in Table 1. In most cases the accelerating feature is apparent, though the relatively low contrast relative to the deeper circumstellar line makes their visualisation difficult in some cases. Also, in a few cases the acceleration is not large enough to be clearly visible (e.g. the 40 km s feature in second panel from top is accelerating slightly, but significantly). While not all features are clearly visible in the figures, the identified features are verified to i) have average absorption significantly greater than zero (i.e. exist), and ii) be accelerating significantly, by the fitting method outlined in section 4.

All accelerating features found are red shifted, and are accelerating towards the star. No automated detection methods were attempted, though we found that some features could also be discovered by fitting Gaussians to absorption features in individual spectra, and plotting the temporal evolution of their radial velocities (see Petterson & Tobin, 1999, for such plots). In addition, distinct unblended absorption features that did not show apparent acceleration were also selected for further analysis, with the expectation that useful orbital constraints could still arise (e.g. a minimum distance at transit). For these, only red shifted events were selected because the aim is to illustrate orbital fitting of accelerating features, and to compare the orbital constraints for both accelerating and non-accelerating features. There is a degree of subjectivity in the subset of spectra used below, and more accelerating features might be discovered with a more systematic approach. We are not however attempting to make any statistical inference from the subset of features analysed, so any biases introduced by our approach do not influence the results that follow.

4 Orbit Fitting

Given an accelerating absorption feature it is relatively simple to construct and fit a model of an orbit. The primary assumption here is that the orbit crosses the centre of the stellar disc, leaving three elements to be constrained (eccentricity , longitude of pericenter , and pericentre distance ). Thus, an orbital fit will yield constraints, but with a strong degeneracy between these within the allowed parameter space. However, as discussed above the reasonable assumption of a parabolic orbit (), means that an accelerating feature contains enough information to constrain the orbit fully. Here we illustrate this degeneracy by leaving to vary, but also restrict the results to parabolic orbits.

Radial velocity is the variable being fitted here, with respect to Earth once Pictoris’ systemic velocity has been subtracted (i.e. is only the same as the stellocentric radial velocity at transit center). This velocity can be derived using the motion of an orbit in a frame , where the pericentre lies along the positive coordinate, and Pictoris is at the origin. The system is viewed from the negative direction, where the , coordinates are rotated clockwise by from , (e.g. Murray & Dermott, 2000). That is, the orbit in the frame is rotated such that the pericentre direction is at an angle from (the reader may wish to refer ahead to Fig. 6.) In the frame of the orbit


where , is the true anomaly, , and for parabolic orbits and for elliptical orbits. The rotation yields


We fit absorption features directly (i.e. instead of fitting accelerations and then deriving orbital constraints), so deriving an orbit from the images requires five or six free parameters, the orbital elements , , and (for elliptical orbits), the time of mid-transit relative to the start of the observation sequence, and the (constant) depth and velocity width of the absorption line, which is assumed to be Gaussian. The radial velocity at a given time is found by first finding the true anomaly (using Barker’s equation in the case of parabolic orbits), and then using equation (4). Thus, a model of an absorption line analogous to the accelerating feature in Figure 2 is created by computing at the centre of each temporal bin, and the image comprises a vertical sequence of (negative) Gaussians of depth (where positive values mean absorption) and width centred at . At each time, an absorption line is only assumed to be present if the centre of the body is in front of the star (i.e. ). Additional parameters that can be derived from the models2 are the radial velocity at mid-transit , the distance to the star at mid-transit (), and the acceleration .

As can be seen from Figures 3 and 4, most spectra contain fixed absorption lines that are much deeper than the accelerating ones. Fitting the model and obtaining reasonable values from residual images that are also easily interpreted visually therefore requires that these are removed. First, we find the median spectrum over the nightly sequence. The noise level in each temporal bin is then computed as the standard deviation of the residuals when this median is subtracted from the initial image. Next, a second order polynomial is fitted to the median in the region near the line of interest (typically 20 km s on either side, though in cases with other nearby lines (e.g. for lines with low ) this region is necessarily smaller). The outer edges of this region have twice the weight relative to the centre in the fit, and this polynomial is taken as an estimate of the local time-independent background near the accelerating absorption feature. This polynomial fit to the median was found to be much better than simply using the median, as the absorption line of interest was commonly partially subtracted with the latter method. During the fitting we subtract this polynomial and the proposed model from the image to obtain a residual image.

Fitting orbits to the spectra is done in two steps, where the goodness of fit metric is the sum of squares of the residuals divided by the noise. First, a model is fitted to the original image with the polynomial subtracted to obtain a preliminary fit. A preliminary residual image is created by subtracting the polynomial and this model, and the residuals are smoothed with a small temporal extent of 8.5 minutes, but a wide spectral extent of 100 km s. These smoothed residuals (which vary at the percent level relative to the absorption features) are then subtracted from the original image before fitting.

As noted above, very strong degeneracies between fitted parameters exist when the eccentricity is not fixed, so the second step is to use a Markov Chain Monte Carlo (MCMC) method to explore the parameter space. Specifically, we use the ensemble sampler emcee (Goodman & Weare, 2010; Foreman-Mackey et al., 2013) with 256 parallel chains (“walkers”). The chains are run for 100,000 steps to ensure the parameter space is fully explored, as in some cases the autocorrelation lengths were of order 10,000 steps. The only restrictions on the fitted parameters are that orbits are bound or parabolic (i.e. ), do not hit the star at pericentre (i.e. ), and that the depth and width are positive.

After these runs, the walkers from the final step were used to generate distributions of orbital elements and other derived quantities of interest. This fitting was run twice; once for elliptical orbits, and then again for parabolic orbits, thus obtaining two sets of orbital elements from each modelled absorption line. The parabolic orbits are of course a subset of the elliptical orbits, but fitting them is more satisfactory than simply applying an arbitrary cut in eccentricity to the elliptical orbit results (at , say).

Figure 5: Fitting results, showing the same data as in Figure 2, but with the best-fit model subtracted. The upper panel shows the level of absorption after model subtraction, where brighter colours indicate deeper absorption. The lower panel shows the binned spectra, with the black lines showing the original data, and the grey lines the model subtracted data.

An example of a model fit is shown in Figure 5, for the same feature in Figure 2 (and the top row of Figure 3). While the accelerating feature is not subtracted perfectly, the model adequately reproduces the accelerating absorption line, and while the depth of the absorption can appear to vary (Beust & Lissauer, 1994), any variation is not at a sufficiently strong level that the addition of another model parameter is warranted for the purposes of orbital fitting. This procedure was applied to all eleven sets of spectra shown in Figure 3, plus others where non-accelerating lines were well-separated from others, thus obtaining sets of elliptical and parabolic orbital elements from a set of 27 absorption features. Of these, twelve are accelerating features (there are two in the second panel from the top in Figure 3).

Figure 6: Distribution of possible elliptical orbits for the example in Figure 2. Each line shows one of the 256 walkers at the final step in the MCMC chains, and each has a corresponding circle at pericenter. The axes here are , , so pericentre is at an angle anticlockwise from . The thick black circle in the centre is the star, and the system is viewed from below (i.e. ).

5 Results

The results of the orbital fits may be visualised in various ways; Figure 6 shows a top-down (, ) view of the final set of orbits for the example above. From this figure the degeneracies are clear. A wide range of orbits is possible, and those with smaller pericentre distances reach pericentre at a greater longitude from the line of sight to the star. The degeneracy with eccentricity is more subtle and harder to discern here, as the most eccentric orbits are in fact those with the greatest pericentre distances and smallest pericentre longitudes.

Figure 7: Accelerations for fitted absorption features, plotted against the distance to the star at mid-transit. Each cluster of dots shows the final step of the 256 walkers in the MCMC chain for the fit to a given absorption feature. The dot colours are used consistently across Figs. 7-10. The dashed line shows the acceleration expected for bodies in free-fall from infinity. Those without detected acceleration extend down to the free-fall line, and their symbols are transparent.

These orbits may also be viewed in parameter space plots, each of which shows clusters of 256 points from the final step in the MCMC chains. These clusters therefore show the distribution of possible parameters for a given absorption feature. The first example is Figure 7, which illustrates the distinction between absorption features that do and do not show measurable accelerations. While features with significant acceleration yield a direct measurement of the distance to the star, features with undetectable accelerations are consistent with orbits that pass in front of the star at arbitrarily large distances. However, the restriction that fitted orbits are bound or parabolic means that for a given velocity at mid-transit the smallest possible acceleration coresponds to a object in free-fall from infinity (which is ). Thus, non- accelerating features are easily identified in Figure 7 by their lack of clustering at a given , and the extension of their distributions down to the expected limit based on the fitting method.

Figure 8: Distributions of possible orbital elements for all fitted features. Each cluster of dots shows the final step of the 256 walkers in the MCMC chain for the fit to a given absorption feature. Those without detected acceleration are transparent and dot colours are used consistently across Figs. 7-10. The left panel shows results for elliptical orbits, and the right panel shows results for parabolic orbits. The range in the right panel is restricted (to the grey box in left panel), as the poorly constrained parabolic orbits look very similar to the elliptical ones. For elliptical orbits the orbital element degeneracies are strong, but the elements are well constrained for parabolic orbits when acceleration is seen.

In the left panel of Figure 8, all fitted orbits are shown, and again the distinction is made between orbits with and without detected acceleration. It is clear that degeneracies apply in all cases. The same trends visible in Figure 6 can be identified, but it is now apparent that each fit reaches a minimum eccentricity that is consistent with the data. This can be understood from equation (4), because at transit , and therefore and the maximum radial velocity at transit for a given orbit occurs when is near 90. Thus, the lowest eccentricity orbit consistent with the data should also be near 90. The eccentricities rise to larger pericentre longitudes, but stop short of because of the restriction that the pericentre distance be larger than . In cases with detected acceleration and are strongly correlated, but this correlation is not present for those without and the distributions form a broad cloud of points centered near .

Restricting the orbits to be parabolic yields tighter constraints, as shown in the right panel of Figure 8. The plotted range is smaller here, because non-accelerating features yield very similar constraints to those for elliptical orbits. For accelerating features however, the constraints localise the orbits to relatively small regions in space. A summary of the parameters derived for parabolic orbits is given in Table 1.

Figure 9: Distributions of velocity versus stellocentric distance at mid-transit, for elliptical (left panel) and parabolic (right panel) orbits. Each cluster of dots shows the final step of the 256 walkers in the MCMC chain for the fit to a given absorption feature. Those without detected acceleration are transparent and dot colours are used consistently across Figs. 7-10. The grey boxes show ‘high’ and ‘low velocity features’, following Beust & Morbidelli (2000). The constraints in this space are only moderately improved when parabolic orbits are assumed.

A final way to view the orbit distributions is shown in Figure 9, primarily for comparison with the models presented in Beust & Morbidelli (2000) (e.g. see their Figure 4). Here the radial velocity at the time of mid-transit is used; while Figure 3 shows that in general the velocity of an accelerating feature at any given time is well constrained, this is not necessarily true at the time of mid-transit (which may occur during, before, or after the observation sequence). Thus, as shown in the left panel, for elliptical orbits the distribution of can be rather wide. When the orbits are restricted to be parabolic, the right panel shows that the constraints become tighter, though not so much that this assumption yields significantly more information than when eccentricity is a free parameter. That is, while parabolic orbits must be assumed to obtain tight constraints on orbital elements, this is less true in the parameter space of Figure 9.

Non-accelerating features result in poor constraints in Figure 9, most of which extend beyond the right side of the plot. The distances at transit are only constrained to be larger than some value that depends on the strength of the signal (i.e. the upper limit on the level of acceleration that was not detected). However, the fact that the absorption features are seen at all means that is no more than a few tens of stellar radii for evaporation to be occuring (Beust et al., 1996, derive 35), so the true of these orbits is likely below this value, and consistent with the orbits with stronger constraints.

In Figure 9 the general trend is that orbits that transit at smaller stellocentric distances have a greater radial velocity. The grey boxes show the estimated regions in which previously identified populations of FEBs lie, where the evaporation model was used to estimate the distance at transit, and the distinction between “high” and “low velocity features” (HVFs and LVFs) is largely historic.3 The first conclusion is therefore that the orbits derived here in general have similar properties to those derived using the evaporation model.

While our orbits tend to lie within the HVF box, they are on average above the LVFs. However, the lack of orbits falling within the LVF box is probably the result of a bias, as features that are the result of transits at greater distances more commonly occur at low velocity, and these are more common and more likely to be blended with other features. Thus, the main conclusions from placing the orbits derived here on this plot is that there is not strong evidence for a lack of orbits within the HVF and LVF boxes, but that orbits certainly exist outside them, in particular at greater velocities for a given distance at transit when . The existence of such orbits is indeed predicted by the resonance model (Beust & Morbidelli, 2000), and the orbits derived from this method will provide useful new constraints for this model.

Figure 10: Correlation between distance at transit and absorption fraction for elliptical orbits. Each cluster of dots shows the final step of one of the 256 walkers in the MCMC chain for the fit to a given absorption feature. Those without detected acceleration are transparent and dot colours are used consistently across Figs. 7-10. This link is predicted by the evaporation model because exocomet comae are more strongly ablated closer to the star, and therefore cover a smaller fraction of the star. The constraints are not improved with the assumption of parabolic orbits (so that plot is omitted here).

The last view of the fitting results includes the absorption depth derived in each case, shown against the distance at transit in Figure 10. The values for are in units of normalised flux (i.e. as in Figure 1), but can be converted to fractional peak absorption using the reference spectrum with , where is the reference (stellar) flux at . The relevance of this plot is that the evaporation model described above predicts that comet absorption is deeper for more distant transits, because the effect of stellar radiation pressure on the calcium ions is less. Thus the comet comae are larger and cover a greater fraction of the projected stellar surface. This prediction is borne out by Figure 10, which shows a clear tendency for the closest transits to have shallower absorption. As with Figure 9, the likely distances for poorly constrained orbits are towards the low ends of their distributions, so the true correlation may be as tight as suggested by the clouds of points for orbits derived from accelerating features.

6 Discussion and Conclusions

This work is motivated by the realisation that previously estimated orbital elements for Pictoris’ comets imply acceleration of the absorption features should be visible in spectra taken over several hours. Figures 2, 3 and 4 show that these accelerations are indeed seen, and using a simple model these accelerating features can be used to derive orbital elements for individual comets. The constraints on the pericentre distance and longitude , and eccentricity are degenerate, but constraints on the distance to the star at the time of transit are less so. Thus, useful constraints that can be directly compared to models of the comets’ origins can be derived in a model-independent way.

The constraints can be further improved with the reasonable assumption that the orbits are parabolic, because comets on low eccentricity orbits spend too much time near the star and may not survive long enough to be observed. With this assumption, the pericentre distance and longitude are well constrained, and the latter lies in the range 25 to 75 for the orbits modelled here (Fig. 8). That is, the tendency for the pericentre longitudes to lie in a restricted range suggested in the past based on orbits estimated using an evaporation model (e.g. Beust et al., 1990) is also found when the orbits are derived directly.

A further test of the evaporation model is of the prediction that comet absorption is deeper for more distant transits, because comae of closer comets are more strongly ablated and therefore cover less of the stellar surface. This correlation is also borne out by our analysis (Fig. 10).

Our technique of exocomet orbit fitting is analogous to other techniques that measure the radial velocity of planets directly, by using cross-correlation with models of CO and HO absorption to show that their spectra are double-lined (e.g. Snellen et al., 2010). Because they can measure the planet velocity with a precision of a few percent, rather than the reflex motion of the host star, these techniques yield stringent orbital constraints, for example on the eccentricity and inclination of non-transiting planets (Brogi et al., 2012). These observations have detected both absorption in the thermal emission spectrum of the planet (e.g. Birkby et al., 2017), and absorption in transmission spectra (i.e. during transit) (e.g. Snellen et al., 2010; Brogi et al., 2018). The latter detections are analogous to the results presented here, and future transmission spectroscopy of transiting planets may also yield independent and useful constraints on the orbital elements.

Our results are reasonably straightforward because the approach is simple. However, a few caveats and inbuilt assumptions should be noted. In fitting orbits, the stellar radius and mass have been assumed fixed, so there remains a small systematic uncertainty in the results. We further assumed that the transits cross the centre of the stellar disc, but there are no cases where both ingress and egress are seen so this assumption should not influence the results significantly. In some cases the absorption lines are probably blended, meaning that a single model may have been applied where two might have been more appropriate; the most likely examples are the 80 km s feature in the second panel and the two broad features between 50 and 150 km s in the fifth panels in Figures 3 and 4. Again, the results are unlikely to be significantly changed by use of a more complex model, but the uncertainties on the derived orbital parameters may increase. The most important issue for further work to address is probably a systematic and objective identification of accelerating features. Other issues include i) simultaneous fitting of H and K lines where possible, ii) simultaneous fitting of possibly bended lines, and iii) modelling variable absorption during transit.

This work has various possible implications. Most obviously, the accelerating absorption features provide strong evidence that the FEB phenomenon indeed originates in bodies that pass in front of Pictoris at close distances (for which comets is the leading hypothesis). By fitting models to these features the orbital properties are found to be consistent with previous estimations that were based on an evaporation model, and therefore these results largely validate that model as a reasonable physical explanation of the comet comae. Further work could attempt to refine the evaporation model by attempting to again fit the absorption features, but now using orbital element constraints derived when acceleration is detected.


GMK is supported by the Royal Society as a Royal Society University Research Fellow. I am grateful to Matteo Brogi and Daniel Bayliss for discussions during the preparation of this work, and to the referee for a valuable review.

The code used in this research is available on github at https://github.com/drgmk/feb-accel, including the final steps of the chains used in Figures 6 to 10.

d km s h km s rad %
54542 9 0.4 68 2.0 8 0.2 0.6 0.02 9 0.2 0.03 0.001
54829 2 0.3 88 1.8 12 0.5 1.2 0.04 18 1 0.04 0.0009
54829 2 0.3 39 1.1 20 2 0.6 0.03 22 2 0.05 0.001
54913 38 0.7 90 0.6 4 0.04 0.6 0.004 4 0.04 0.03 0.0005
55170 1 0.1 73 1.2 18 0.5 1.3 0.03 28 1 0.2 0.002
55567 22 1.1 68 0.5 5 0.1 0.5 0.008 6 0.1 0.02 0.0008
55597 65 6.0 117 3.3 3 0.1 0.7 0.02 3 0.2 0.02 0.001
56685 3 0.4 86 2.8 11 0.5 1.1 0.05 15 1 0.04 0.002
56694 1 0.3 52 1.4 20 2 0.8 0.06 23 3 0.06 0.002
56982 5 1.0 73 4.1 10 1 0.8 0.06 12 1 0.06 0.003
56988 13 1.3 125 7.3 5 0.4 1.1 0.06 7 0.4 0.02 0.002
57344 21 3.4 89 9.7 5 0.5 0.7 0.07 6 0.5 0.02 0.002
Table 1: Median parameters derived for accelerating features assuming parabolic orbits, and their standard deviations. MJD refers to the modified Julian date on which the spectral observations began, and the table is listed in the same order as Figure 3.


  1. i.e. the spectra shown were not chosen completely randomly, but only a few attempts need to be made to obtain a similarly varied set of spectra where both deep and blueshifted lines are present.
  2. We do not distinguish between and because the distance to the star changes little during a transit. We use instead of for the same reason.
  3. Here the ‘very low velocity features’ (VLVFs) are not considered, as none in this class were found to be accelerating.


  1. Beust, H., Lagrange, A.-M., Plazy, F., & Mouillet, D. 1996, A&A, 310, 181
  2. Beust, H., Lagrange-Henri, A. M., Vidal-Madjar, A., & Ferlet, R. 1989, A&A, 223, 304
  3. Beust, H. & Lissauer, J. J. 1994, A&A, 282, 804
  4. Beust, H. & Morbidelli, A. 1996, Icarus, 120, 358
  5. —. 2000, Icarus, 143, 170
  6. Beust, H., Vidal-Madjar, A., Ferlet, R., & Lagrange-Henri, A. M. 1990, A&A, 236, 202
  7. Birkby, J. L., de Kok, R. J., Brogi, M., Schwarz, H., & Snellen, I. A. G. 2017, AJ, 153, 138
  8. Blum, J., Gundlach, B., Krause, M., Fulle, M., Johansen, A., Agarwal, J., von Borstel, I., Shi, X., Hu, X., Bentley, M. S., Capaccioni, F., Colangeli, L., Della Corte, V., Fougere, N., Green, S. F., Ivanovski, S., Mannel, T., Merouane, S., Migliorini, A., Rotundi, A., Schmied, R., & Snodgrass, C. 2017, MNRAS, 469, S755
  9. Brogi, M., Giacobbe, P., Guilluy, G., de Kok, R. J., Sozzetti, A., Mancini, L., & Bonomo, A. S. 2018, ArXiv e-prints
  10. Brogi, M., Snellen, I. A. G., de Kok, R. J., Albrecht, S., Birkby, J., & de Mooij, E. J. W. 2012, Nature, 486, 502
  11. Crawford, I. A., Beust, H., & Lagrange, A.-M. 1998, MNRAS, 294, L31
  12. Duncan, M., Quinn, T., & Tremaine, S. 1987, AJ, 94, 1330
  13. Duncan, M. J. & Levison, H. F. 1997, Science, 276, 1670
  14. Eiroa, C., Marshall, J. P., Mora, A., Montesinos, B., Absil, O., Augereau, J. C., Bayo, A., Bryden, G., Danchi, W., del Burgo, C., Ertel, S., Fridlund, M., Heras, A. M., Krivov, A. V., Launhardt, R., Liseau, R., Löhne, T., Maldonado, J., Pilbratt, G. L., Roberge, A., Rodmann, J., Sanz-Forcada, J., Solano, E., Stapelfeldt, K., Thébault, P., Wolf, S., Ardila, D., Arévalo, M., Beichmann, C., Faramaz, V., González-García, B. M., Gutiérrez, R., Lebreton, J., Martínez-Arnáiz, R., Meeus, G., Montes, D., Olofsson, G., Su, K. Y. L., White, G. J., Barrado, D., Fukagawa, M., Grün, E., Kamp, I., Lorente, R., Morbidelli, A., Müller, S., Mutschke, H., Nakagawa, T., Ribas, I., & Walker, H. 2013, A&A, 555, A11
  15. Ferlet, R., Vidal-Madjar, A., & Hobbs, L. M. 1987, A&A, 185, 267
  16. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  17. Galland, F., Lagrange, A.-M., Udry, S., Chelli, A., Pepe, F., Beuzit, J.-L., & Mayor, M. 2006, A&A, 447, 355
  18. Goodman, J. & Weare, J. 2010, Comm. App. Math and Comp. Sci, 5, 65
  19. Hobbs, L. M., Vidal-Madjar, A., Ferlet, R., Albert, C. E., & Gry, C. 1985, ApJ, 293, L29
  20. Kiefer, F., Lecavelier des Etangs, A., Boissier, J., Vidal-Madjar, A., Beust, H., Lagrange, A.-M., Hébrard, G., & Ferlet, R. 2014, Nature, 514, 462
  21. Koen, C., Balona, L. A., Khadaroo, K., Lane, I., Prinsloo, A., Smith, B., & Laney, C. D. 2003, MNRAS, 344, 1250
  22. Kondo, Y. & Bruhweiler, F. C. 1985, ApJ, 291, L1
  23. Lagrange, A.-M., De Bondt, K., Meunier, N., Sterzik, M., Beust, H., & Galland, F. 2012, A&A, 542, A18
  24. Lagrange, A. M., Ferlet, R., & Vidal-Madjar, A. 1987, A&A, 173, 289
  25. Lagrange, A.-M., Plazy, F., Beust, H., Mouillet, D., Deleuil, M., Ferlet, R., Spyromilio, J., Vidal-Madjar, A., Tobin, W., Hearnshaw, J. B., Clark, M., & Thomas, K. W. 1996, A&A, 310, 547
  26. Lagrange-Henri, A. M., Gosset, E., Beust, H., Ferlet, R., & Vidal-Madjar, A. 1992, A&A, 264, 637
  27. Lagrange-Henri, A. M., Vidal-Madjar, A., & Ferlet, R. 1988, A&A, 190, 275
  28. Levison, H. F., Duncan, M. J., & Wetherill, G. W. 1994, Nature, 372, 441
  29. Murray, C. D. & Dermott, S. F. 2000, Solar System Dynamics
  30. Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Science, 110, 19273
  31. Petterson, O. K. L. & Tobin, W. 1999, MNRAS, 304, 733
  32. Pichierri, G., Morbidelli, A., & Lai, D. 2017, A&A, 605, A23
  33. Quillen, A. C. & Holman, M. 2000, AJ, 119, 397
  34. Sibthorpe, B., Kennedy, G. M., Wyatt, M. C., Lestrade, J.-F., Greaves, J. S., Matthews, B. C., & Duchêne, G. 2018, MNRAS, 475, 3046
  35. Smith, B. A. & Terrile, R. J. 1984, Science, 226, 1421
  36. Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
  37. Thébault, P. & Beust, H. 2001, A&A, 376, 621
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description