© ESO, 2012

1. Introduction

Rosetta is ESA’s comet orbiter mission, launched in March 2004 and currently en route to Jupiter-family comet 67P/Churyumov-Gerasimenko (hereafter 67P/C-G). The probe will rendez-vous with the comet in 2014 and remain in the vicinity of the nucleus for on-going detailed physical and compositional analysis. While the original target of this mission was studied extensively, we are only now beginning to build up a detailed picture of the bulk physical properties of the nucleus of the new target, as well as its heliocentric lightcurve behaviour. Such knowledge is invaluable with regard to mission planning and science scenarios. Both the size and shape of the nucleus are critical issues that bear directly on the operation of the Rosetta orbiter and on the success of the lander package. Also, such a mission will provide a unique opportunity to test ground based observing techniques. For example, using multiple rotational lightcurves obtained at a range of observational geometries can allow a detailed shape model and pole orientation to be derived. Additionally, using mean magnitude estimates at a wide range of phase angles can allow one to characterize the surface roughness of the nucleus via Hapke modelling (Hapke 2002; Lederer et al. 2005). The only way to definitively test these models is via in-situ analysis.

Comet 67P/C-G orbits the Sun every 6.57 years and is a member of the low inclination Jupiter-family of comets. Prior to its consideration as a potential Rosetta target, very little was known about the nucleus properties. Ground-based observations in April–May 1991 (Mueller 1992) suggested that the mean effective radius is 2.8 ± 0.1 km (assuming an albedo of 0.04) and that the nucleus has an axial ratio a/b > 1.7, where a and b are the axes of the assumed prolate ellipsoid nucleus. However, a complete description of the analysis and results was never published. Subsequent observations using the 1 m JKT only established an upper limit for the radius of 2.9 km (Lowry et al. 1999). In March 2003, Philippe Lamy and colleagues were granted director’s discretionary time to observe 67P/C-G using the planetary camera on the Hubble Space Telescope (HST; Lamy et al. 2006). The comet was at 2.52 AU from the Sun at the time of the observations and was substantially active. By taking advantage of the high angular resolution offered by the HST Lamy et al. were able to detect the rotational signature of the nucleus by applying a coma subtraction method. Applying both phase dispersion minimization and Fourier analysis techniques to their time series R-filter magnitudes yielded a best fit nucleus rotational period of 12.30 ± 0.27 h. The resulting mean effective radius was found to be 1.98 ± 0.02 km, for an assumed albedo of 0.04.

The first direct observations of the nucleus were obtained using ground-based telescopes by our team. Optical imaging observations were performed at the ESO 3.6 m NTT, Chile, in February & April 2004 and May 2005 (also see Lowry et al. 2006), as the comet was passing through aphelion. Observing comets at large heliocentric distances (>3 AU) has proven to be an effective means of directly observing the nucleus due to much decreased or negligible coma activity. Rosetta will encounter the comet at around 4 AU from the Sun. When our observations were taken it was not yet known how active the comet would be in this part of its orbit. Our observations showed that the comet was barely active, if at all, and certainly not one of a significant group of JFCs that are vigorously active over their entire orbits (e.g. Lowry et al. 2010; Fernández et al. 2008; Kelley et al. 2008; Snodgrass et al. 2010). This is rather fortunate for the mission, as there is an opportunity to monitor the progressive increase in activity much closer to an inactive state than may otherwise have been possible with comet 46P/Wirtanen.

An extensive database of nucleus observations has since been acquired, not only from large ground-based telescopes like the ESO VLT (Tubiana et al. 2008, 2011), but also at thermal-IR wavelengths from Spitzer (Kelley et al. 2006, 2009; Lamy et al. 2008). More details on all of these campaigns are discussed later in the paper. Agarwal et al. (2010) presents observations of the dust trail at optical and thermal-IR wavelengths as the comet passed through aphelion from 2004–2006. The ESO 2.2 m and Spitzer were used.

In this paper we present all of our observations from the ESO NTT. Preliminary analysis of part of our complete dataset has been presented elsewhere (Lowry et al. 2006; Lamy et al. 2007). Our aim here is to present a complete and finalized analysis of our data, with the ultimate aim, and the long-term aim of our programme, of linking all of the current and published nucleus photometry that may be available into a single shape and spin model, which can also be used to perform a detailed thermophysical analysis of its surface.

In Sect. 2 we discuss our observations in detail, as well as the photometric analysis methods that were applied to the data to extract the nucleus rotational lightcurves. We then go on to discuss other published optical datasets that are to be used in our shape and spin state modelling, discussed in Sect. 3. Using this shape model, we then applied a new thermophysical modelling method to characterize the thermal properties of the nucleus surface, such as albedo, thermal inertia, surface roughness, as well as nucleus physical dimensions (Sect. 4). All available Spitzer data on the nucleus was incorporated into this analysis. No thermal-IR data exists from the ground due to the intrinsic faintness of the inactive nucleus of this comet. An overall summary of the main results obtained is presented in Sect. 5.

2. Optical observations and data analysis

We conducted a series of observing runs at the 3.6 m New Technology Telescope (NTT) at the European Southern Observatory, La Silla, Chile. Our main aim was to obtain time-series optical-wavelength imaging in order to detect and extract the comet’s nucleus rotational lightcurve at as wide a range of observational geometries as possible, for the purpose of shape and spin state modelling. Our observing runs span a range of dates from February 2004 to September 2007 when the comet passed through aphelion. See Table 1 for the range of observational geometries attained.

The main bulk of our data comes from two separate runs conducted in May 2005 and July 2007, but other smaller datasets are discussed. Three nights of optical imaging data were taken on May 10-12-14, 2005 (R_h = 5.60 AU [aphelion = 5.71 AU]). We used the “red arm” of the ESO Multi-Mode Instrument (EMMI) for all but one of our runs. EMMI uses two 2k  ×  4k MIT/LL CCDs with a pixel scale of 0.33′′/pixel (in 2 × 2 binning mode) and a full effective FOV of  ~9.1′ × 9.1′ on the red arm of the detector. The “blue arm” was not used during this programme. At the time of this set of observations Lamy and colleagues had already observed the comet while in its active stage with the HST in 2003 (Lamy et al. 2006) and noted a rotation period of just over 12 h. This information was used to plan our observing strategy for this run. We observed the comet for 3 nights over a 5 night period (with 1-night intervals) to allow 100% of the rotational phase to be covered. We cycled through the filter sequence RVRIR each night, and full rotational lightcurves in all three passbands were obtained, allowing detailed inspection of surface colours at all rotational phases.

Extraction of the R-filter photometric lightcurves involved first dividing the full R-filter image list into groups of 6, or more in a few cases if the signal level from the comet was relatively low. These groups were then co-added first according to the rates of motion of the comet during the time-range of that group, and then a second co-added image was created, by registering and stacking this same set of images on the stars. As all images were sidereally tracked, with the exposure times set low enough that the motion across the chip within a given 100 s exposure produced negligible trailing, we were then able to reliably reconstruct the background PSF as if we had a single continuous 600 s exposure differentially tracked at the comet’s rates. This in turn allowed us to apply aperture corrections much more effectively. For a given co-added image we then set the cometary photometric aperture radius at some multiple of the seeing for that combined image, and repeated this in a consistent manner for each group. A selection of stars was identified and used as comparison stars. Once relative R-filter lightcurves were successfully extracted we set about calibrating them using standard stars from the Landolt catalogue (Landolt 1992). Our best night was on night 2, which displayed photometric conditions for the entire night. We took the opportunity to obtain calibration images at the comet’s position on the preceding and succeeding nights, in all three filters, for calibration purposes. Analysis of the data was successful, and the rotational signature of the nucleus was clearly detected at all passbands. For the coma search we stacked all R-filter images for a given night on the comet and compared the brightness profile with the stellar profile from the same image set stacked on the stars (see Fig. 1). The pre-coaddition method for the lightcurve extraction was very necessary to detect the rotation lightcurve of the nucleus. Simply extracting the relative magnitudes directly from each 100 s exposure and taking a 6-point average afterwards resulted in low S/N random scatter. The extracted magnitudes are listed in Table 2. The dust trail was not detected on any of our NTT images including the stacked R-filter images. This is not surprising given our surface brightness lower limit for detection of  ~26.4 mag/arcsec², for the co-added R-filter image for night 2. Tubiana et al. (2008) measure an R-filter surface brightness for the trail of  ~28 mag/arcsec² using the larger VLT facility, when the comet was at smaller heliocentric and geocentric distances of 4.89 AU and 4.29 AU, respectively, and phase angle of 10.3°. Scale this value according to the observational geometry of our NTT observations, and the trail is  ~1.7 mag/arcsec² too faint for detection.

Fig. 1 Composite R-filter images obtained at the NTT in May 2005. Each panel consists of 71 × 100 s, 123 × 100 s, and 82 × 100 s, images stacked according to the comet’s projected rate of motion. In each case, as with all our data, the comet appears as a sharp point source on our images. Celestial North and East are marked along with the projected solar direction and heliocentric velocity vectors.

Selected key preliminary results from these data were presented in Lowry et al. (2006) and Lamy et al. (2007), and included a Fourier analysis of the R-filter lightcurve which gave a best-fit synodic rotation period of 12.72 ± 0.05 h. The observed brightness range of 0.38 ± 0.04 mag implied a projected nucleus axial ratio of 1.42 ± 0.05. The mean R-filter magnitude was 22.41    ±    0.03, which corresponds to mean nucleus effective radius of 2.26 ± 0.03 km, for an assumed albedo and linear phase coefficient of 0.04 and 0.04 mag/deg, respectively. Assuming the same albedo gives projected semi-axial dimensions of around 2.94 ± 0.15 and 2.07 ± 0.04 km, although this should be revised downwards as the albedo is somewhat larger as we show later (also see Kelley et al. 2009; Lamy et al. 2008). These observations were the first time-series imaging to be taken of the comet at large heliocentric distances, and showed that the nucleus was likely to be inactive or close to inactive during the initial approach phase of Rosetta. The comet was observed at opposition with an average phase angle of just 0.46 degrees. Our three nights of VRI-band imaging span a phase angle range of 0.07–0.84 degrees, allowing the opposition surge to be sampled.

The 2nd of our main data sets include 3 nights of observations taken in July 16-19-22, 2007 (R_h = 4.63 AU [inbound]). As with the May 2005 run, we used EMMI and incorporated an interval between each night of observation to ensure 100% coverage of the rotational phase, while making optimal use of the telescope. On this occasion a 2-night interval was used due to NTT scheduling constraints. The rotational signature of the nucleus was again clearly detected at all passbands and there was no indication of resolved coma activity. These July ’07 observations are particularly important as they were taken around the same orbital position that the comet is expected to have during the initial approach phase of Rosetta. Extraction of the relative photometric lightcurves was performed using the same approach as in the May 2005 data set. However, in this case the comet was passing through a dense stellar region. Care was taken to identify those frames that may have some flux excess at the comet’s position due to contamination from a nearby star, and subtract from each affected image a star-only image of that field (obtained through a median combination of comet frames taken before and after the image). This worked well, but the quality of the resulting photometry was not quite as good as the May 2005 case. Unfortunately the data from July 22 was deemed too unreliable due to the dense stellar field, and excluded from the analysis. Photometric calibration was performed in the normal manner using Landolt standards. A mean R-filter absolute magnitude of 15.35  ±  0.04 was observed (“Reduced” magnitude m_R(1,1,α)    =    15.80), along with a colour index of (V − R) = 0.548 ± 0.050.

Fig. 2 Upper panels: left – periodogram produced via least-squares fitting of a Fourier series to the raw lightcurve data. The most likely period solutions are 12.40971 ± 0.00005 h and 12.76236 ± 0.00006 h. Middle and right – folded lightcurves of the uncorrected data points for each of these two best-fit rotation periods. The fit in the 12.41 h case is unconvincing with various datasets not following the phase of the folded curve. The fit for the 12.76 h case is better, but can be further improved by correcting for the differing observational geometry between the various datasets. Lower left: periodogram produced via least-squares fitting of a Fourier series to the lightcurve data which has been corrected for the movement of the PAB. In this case the most likely solution was 12.76137 ± 0.00006 h. Lower right: folded lightcurve using the data corrected for the movement of the PAB vector. The scatter and phasing of the datapoints is improved, leading to a more reliable period measurement.

A further 1/2-night of NTT time was allocated to this programme and scheduled for Sept. 13, 2007 to obtain VRI-band data (R_h = 4.40 AU [inbound], Δ = 4.25 AU, α = 13.21°). The aim here was to acquire a partial lightcurve at VR-passbands to measure the mean brightness at a relatively large phase angle. 3 × 210 s, 2 × 400 s, 3 × 600 s images were taken in each of the V and R filters spanning a continuous 6 h period, to measure the full apparent magnitude brightness range of the nucleus, as opposed to a full lightcurve that could be used in any shape modelling analysis. Unfortunately, seeing was very poor throughout the run and a detection was not possible even by co-adding the data. Fortunately we were able to use unpublished brightness measurements obtained at the VLT based on an actual detection, and taken on Aug. and Oct. 2007 at similar phase angles to our Sep. 2007 NTT data set. These provided important substitute brightness measurements at phase angles of  ~10.1–12.5° (Snodgrass et al., in prep.). However, the two datasets include just very short image sequences at various rotational phases. Therefore before we can successfully incorporate these additional unpublished VLT measurements into any updated analysis of the phase darkening curve, we must first scale the observed apparent magnitudes to the lightcurve midpoints. For this we need an adequate shape and spin state model, which is presented later in Sect. 3. For this reason we delay presenting an updated phase darkening analysis until Sect. 3.2.

Finally, a small data set was acquired using the SuSI2 imaging system on the NTT, on 26th February 2004. On this occasion the phase angle was at  ~11.5°, giving us another important datapoint at a larger phase angle than our other NTT datasets where we have full lightcurve data. As with the VLT data from Snodgrass et al. (in prep.), a magnitude correction is needed based on the rotation phase and shape model orientation at the time the data was taken. SuSI2 uses two 2k  ×  4k CCDs, and has a FOV of  ~5.5′ × 5.5′ and pixel scale of 0.166′′/pixel. Four 300 s images were taken in each of the BVR filters, while cycling through the filters. Conditions were photometric and the comet was clearly detected in each filter. Each set was stacked on the comet and apparent magnitudes were measured using aperture corrections on several relatively bright stars nearby. Apparent magnitudes in each filter were measured as m_B = 23.53, m_V = 22.71 and m_R = 22.29 (± 0.06), corresponding to (B − V) and (V − R) colour indices of 0.81 and 0.43 (± 0.09), respectively.

2.1. Other optical data included in the analysis

In addition to the unpublished short image sequences from the VLT (see Sect. 2), several nucleus photometry datasets are available in the literature, which we can incorporate into our analysis, especially the shape and spin state modelling component. The first of these are from Lamy et al. (2006). Observations were performed with HST’s Wide Field and Planetary Camera 2 between 11.4–12.3 March, 2003 (21 h time interval). The comet was at heliocentric and geocentric distances of  ~2.51 and 1.53 AU, respectively, with the comet still in its active phase. Lamy et al. applied a detailed coma-subtraction method to successfully extract the nucleus flux to reveal a double-peaked lightcurve consistent with a rotation period of 12.41    ±    0.41 h. A total of 61 data points were used, taken with the HST F555W and F675W optical-wavelength filters, while focussing mainly on the latter. This is a particularly important dataset given the close proximity to Earth and thus the very different observational geometry attained as compared to observations from Earth when the comet was at large heliocentric distances. Expanding the range of observational geometries over which nucleus photometry is obtained is crucial for the shape modelling work.

Another excellent data set was acquired at the ESO Very Large Telescope by Cecilia Tubiana and colleagues (Tubiana 2008; also see Tubiana et al. 2008; 2011), from April 2004–July 2007, roughly the same time period as our data. As with our campaign, the comet was targeted when it was at large heliocentric distances. The heliocentric distance and phase angle ranged from  ~4.62–5.62 AU, and 0.5–10.5°, respectively. The observations focused on time-series R-filter imaging with the FORS instrument. Due to the much larger aperture used no co-addition of multiple images was required for the extraction of a useful lightcurve data point. A total of 10 nightly datasets were used in our analysis, comprising 295 individual lightcurve magnitudes. A full listing of datasets used can be seen in Table 1.

3. Nucleus shape and spin state modelling

Typically, a large number of lightcurves spanning a wide range of observing geometries are required to derive a reliable shape and spin state model. The observing geometry of 67P/C-G does not change significantly during the course of the ground-based observations presented here, due to the low orbital inclination of the comet and the necessity of performing the observations when the comet was at large heliocentric distances. At these distances we believe the comet was inactive, thus negating the effects of coma contribution to the nucleus brightness. This is due to the comet’s heliocentric lightcurve behaviour throughout the entire ground-based datasets considered here, which was consistent with an inert body. The observer-centered Ecliptic longitude varied by approximately 100 degrees while the observer-centered Ecliptic latitude varied by only 15 degrees. This makes accurate shape modelling more difficult since we are receiving shape information from a reduced proportion of the nucleus surface.

Fig. 3 Contour plot showing the expected location of the rotation pole produced by inverting the lightcurves at each potential pole location. Darker colours here indicate lower chi-squared values. Ecliptic latitude and longitude are denoted by β and λ, respectively. The plot indicates that the pole is most likely located at high latitude, with two prominent and well-defined minima at λ = 78° and β =  + 58° and a 2nd unrelated solution at λ = 278° and β =  + 58°. (A colour version of this figure is available online.)

Reconstruction of a shape model from lightcurve inversion techniques also naturally relies upon the accurate measurement of the sidereal rotation period. As an initial step in this process we performed a least squares Fourier fit of all of the lightcurve data to determine an estimate for the rotation period. The resulting periodogram is plotted in Fig. 2. The best fit solution was found to be approximately 12.41 h with another potential solution at 12.76 h. The resulting folded lightcurves for both rotation periods are also plotted in Fig. 2. Although the fit appears reasonable in the case of the 12.41 h solution, on closer inspection we found that the phase of individual lightcurves did not appear to correspond to the expected phase of the folded lightcurve. We therefore rejected this period solution in favour of the solution at 12.76 h, which folds quite well. As an additional check we isolated some of the more densely-sampled and well-covered lightcurves, and again applied least squares fitting to a Fourier series. This resulted in period solutions consistently in the range 12.6 to 12.9 h. Of course, this method does not give the sidereal rotation period, but rather an average of multiple synodic rotation periods.

To determine an initial estimate for the sidereal rotation period, we used the methods of Kaasalainen et al. (2001,b) and our customized version of the software described by Ďurech et al. (2010). This works by fitting a sidereal rotation period and shape model simultaneously to the lightcurve data for six random poles spread throughout the celestial sphere. At each pole location the rotation period and shape model are optimised using a Levenberg-Marquardt algorithm as described by Press et al. (1992). Essentially, model lightcurves are generated and compared with the observational data with the fitted parameters adjusted simultaneously to minimize fit residuals. The spin rate of the model is the sidereal rotation period, as the geometry vectors between the Earth, comet and Sun are computed beforehand. The resulting artificial lightcurves are then equivalent to the observed “synodic” lightcurves, automatically negating the need to apply corrections to the JDs of the data points due to varying geometry between the datasets via some approximation method like the “PAB” method (see below). We set our initial sidereal period scanning range to 12.4–12.9 h based on our initial estimate above. The best fitting period was found to be at 12.761 h. Since this method relies on an optimisation method and is not a direct measurement per se, no formal error was established, nor needed at this stage. The rotation period corresponded fairly well to previous estimates (see below).

To determine the pole solution we performed the same operation over the entire pole space. To do this we constructed a 5 × 5 degree grid of potential pole locations and used lightcurve inversion to optimise a shape model at each trial pole, using the estimated sidereal rotation period of 12.761 h. The lightcurves were treated as relative throughout our analysis so scattering parameters were not optimised in the process. This does not affect the shape of the comet nucleus, only its size scale. By using calibrated lightcurves the scattering function could also be estimated, but this was not critical for our analysis. Again a Levenberg-Marquardt algorithm is used to determine the optimal values for the various parameters which determine the shape of the nucleus at each pole location. The parameters varied during this process were the coefficients used in the Laplace series to generate the facet areas. The period and pole were not optimised during this process and 50 iterations were used to determine the best potential shape for the given pole and rotation period. A chi-squared value was returned for each pole location indicating the level of agreement between the potential model and data. The lowest chi-squared value and thus best pole solution was located at λ = 65°β = 65°.

With a confident pole solution at hand we can now apply a semi-independent method to constrain the sidereal rotation period, and iterate this back into the pole search for a more refined pole measurement, and thus more reliable shape. The main difficulty in successfully measuring the sidereal rotation period of the nucleus is due to the changing geometry of the various observations. As already noted the observer-centred Ecliptic latitude does not vary strongly due to the low inclination of the comet, but the observer-centred Ecliptic longitude varied by up to 100 degrees. The result of this variation is that lightcurve maxima or minima might be observed earlier or later than expected. Upon phasing the lightcurves from various observational geometries we can arrive at an erroneous rotation period solution. This can be avoided by correcting the times of the lightcurve datapoints to the same observing geometry. This is achieved by determining the movement of the phase angle bisector (PAB) between the various observations (e.g. Magnusson 1989). We then determine the small shifts in the times of the data points. However, knowledge of the spin axis orientation is required before this can be achieved, since it is the movement of the PAB vector projected onto the equatorial plane of the body which is important.

Fig. 4 Comparison of our rotation pole with others from the literature. For consistency with other authors, mainly Lamy et al. (2007), we compare the various measurements in RA-Dec space. We first converted all previous values to Ecliptic coordinates if not done so already, then took the regions defined by these values and their associated uncertainties and mapped each region onto the celestial sphere. Key: S06 – Schleicher (2006), C04 – Chesley (2004), V12 – Vincent et al. (2012), L07 – Lamy et al. (2007), DG05 – Davidsson & Gutierrez (2005). Our primary solution at λ = 78° and β =  + 58° (or RA = 51°, Dec = +79° – Solution A), is in excellent agreement with those values determined from completely independent methods. We have a 2nd unrelated solution at λ = 278° and β =  + 58° (or RA = 275°, Dec = +35° – Solution B), which is well outside the range of all the other solutions. With the constraints of the other studies in mind, we discount Solution B. An inherent property of lightcurve inversions is the inevitable existence of “conjugate” solutions, z-axis “mirrors” of the best fit shape model(s), for which alternate best-fit pole solutions are possible (i.e. at λ ± 180°, Kaasalainen & Lamberg 2006). For completeness, we investigated these conjugates by generating the mirrored shape models for each of the two best fit models (at Solutions A+B in Fig. 3), to see where their preferred pole solutions lie (marked as Ac and Bc on the plot). We find that these solutions result in much poorer fits to the data, especially solution Ac, which is barely defined at all. Note that the Lamy et al. (2007) value is also based on lightcurve inversion, but on just a small subset of the data used here. They also find a 2nd solution and its conjugate, marked as “L07 (B)” and “L07 (Bc)”, respectively.

The best pole solution was then used to determine the movement of the PAB vector with respect to the comet’s equatorial plane for each dataset. The time of each lightcurve data point was then shifted by the appropriate amount. The sidereal rotational period was measured to be 12.76137 ± 0.00006 h via least-squares fitting of a Fourier series, and is somewhat larger than a previous measurement of 12.7047 ± 0.0011 h (Tubiana et al. 2008). The periodogram and folded lightcurve are shown in Fig. 2. The match between the individual lightcurves and the fitted lightcurve was much improved over previous attempts. To refine our pole measurement the pole grid was searched again, this time over a 1 × 1 degree grid using the new period measurement, which was of course held fixed. The best pole was found to be at Ecliptic coordinates: λ = 78° and β =  + 58°. Formal error bars are  ± 10°, the determination of which will be discussed in Sect. 3.1. New time-shifts were calculated for the lightcurve data points using this updated pole solution. The new PAB-based sidereal period solution did not change significantly and further iterations were deemed unnecessary. Figure 3 shows the variation in chi-squared with pole location. A second unrelated, yet equally good, solution is uncovered at around λ = 278° and β =  + 58°. However we believe we can exclude this second solution based on previous observational constraints (see Fig. 4). Observation and modelling of jet activity on the comet indicates that the pole should lie at either and β = 70 ± 10° or and β = −70 ± 10° (Lamy et al. 2007, and references therein). A new independent measurement, based on modelling of dust jets from the 2009 apparition, produces pole solutions at RA = 40°, Dec = 70° (prograde, λ = 66° and β =  +51°), and RA = 220°, Dec = –70° (retrograde, λ = 246° and β = −51°) (Vincent et al. 2012). Since our potential pole solutions are both located at positive latitudes implying prograde rotation, at odds with previous observational constraints, we opt to discount the potential solution at longitude  ~280°. Our final shape solution was constructed using the programme conjgrad (Ďurech et al. 2010) with the measured period of 12.76137 ± 0.00006 h and corresponding pole at λ = 78° and β =  + 58°. The program was allowed to use 100 iterations to arrive at the best fit shape. The final convex shape model is displayed in Fig. 5 with views along the x, y, and z axes. Also see Fig. 7, which compares each observational lightcurve with its corresponding model lightcurve, generated using the shape and spin-state model.

Fig. 5 Best-fit convex triangular-facet shape model produced via lightcurve inversion using the measured sidereal rotation period of 12.76137 h and the most likely pole orientation, λ = 78° and β =  + 58°. Axial ratios are well constrained at b/a = 1.239 and c/a = 0.819.

3.1. Shape model variants and pole uncertainties

While the overall combined dataset is one of the most extensive datasets obtained for a presumably inactive nucleus, the data is of relatively low S/N compared to asteroid datasets, for which this kind of shape modelling is normally applied. It is therefore particularly important to investigate how reliable our best-fit shape model is, or at least its overall relative axial lengths. We therefore set about quantifying; a) how flexible we can be in our best-fit shape while still producing a statistically meaningful fit to the data; b) the formal statistical uncertainties in our rotation pole measurement, given this measure of flexibility; c) and as a final step we generated several shape model variants to assess the stability of our shape model as the pole is varied across its formal 3-σ uncertainties.

For step “a”, we took our best-fit shape model and transformed the relative sizes of the axial dimensions of the model by stretching the model along each principal axis at random increments to include all realistic axial ratios. We then used all the new shape model permutations to generate new sets of model lightcurves for comparison with the data via chi-squared fitting. Our best-fit sidereal rotation period and pole solution remained fixed at each iteration. Note that we decided not to add small-scale features as a degree of freedom in this analysis, given the low S/N of the data. Our aim was to generate a 3D grid of chi-squared values at all realistic stretching factors along the three principal axes. In order to display the results from this analysis in a convenient manner, we plotted the 2D chi-squared grids in each of the a − b, a − c, and b − c planes (in each case the third axis is kept constant at its best-fit value, i.e. the scaling factor is kept at a value of 1). This is illustrated in Fig. 6. We display the results this way because if we let all three axes vary simultaneously we would get an infinite number of best fit values when all three axes are increased or decreased by the same fraction, as opposed to a 3D error ellipsoid. The best-fit axial lengths are located at the origin – the abc correspond to the xyz axes in Fig. 5. We find that a, b, and c can vary by as much as 7%, and still result in a statistically meaningful fit, at the 3-σ level. As we are not considering smaller scale features, this percentage may be slightly larger.

Fig. 6 To investigate how “flexible” we can be with our best-fit shape while still producing a statistically meaningful fit to the data, we used our best-fit shape model and transformed the relative sizes of the axial dimensions of the model by stretching the model at each principal axis at random increments to include all realistic axial ratios. We then took all the new shape model permutations to generate new sets of model lightcurves for comparison with the data via chi-squared fitting. We find that a, b, and c can vary by as much as 7%, and still result in a statistically meaningful fit, at the 3-σ level. Lower right panel: given this measure of flexibility in the shape, we find that the uncertainties in the pole coordinates are revised upwards to  ~±10°, indicated by the outer curve. The grey band indicates how the outer curve can vary when the sidereal rotation period is allowed to vary within its associated uncertainty. See Sect. 3.1 for full details.

To investigate how flexible the pole can be in a more formal manner, we first took our best-fit shape and assigned a series of pole orientations to it, covering a 1 × 1 degree grid around our best fit solution. Again, the sidereal rotation period always remains fixed throughout. For each pole location, we generated new model lightcurves for chi-squared fitting with the data. This returned a formal 3-σ error ellipse, illustrated in Fig. 6. We determined how this 3-σ region grows when we factor in the full 3-σ range of the shape-model axial lengths. As expected, this enlarges the 3-σ confidence region considerably. The formal uncertainties are  ± 8° in Ecliptic longitude and  ± 10° in Ecliptic latitude, which we round off to  ± 10° for both parameters, given the approximate nature of this approach. We regard these uncertainties as good estimates but not necessarily definitive values.

Finally, we generated a series of shape model variants for various extreme values of the pole location, based on these formal 3-σ uncertainties. While there are some minor differences, the overall shape remains satisfactorily stable across each variant. While lightcurve inversion techniques cannot provide a detailed picture of the nucleus, we find that this technique, with these data, do provide strong constraints on the overall relative dimensions of the nucleus, specifically the axial ratios (which can be appropriately size-scaled with sufficient supporting thermal-IR data). Some additional detail in the surface features is required to match the lightcurves more closely. In some cases a perfect match between data and model was not possible, possibly indicating broad concave features that would otherwise not show up in our modelling. In any case we can safely assume that when Rosetta eventually images the nucleus we will see a more rounded body like 9P/Tempel 1 (A’Hearn et al. 2005), as opposed to a highly elongated and possibly “bi-lobed” nucleus similar to comets 103P/Hartley (A’Hearn et al. 2011) or 8P/Tuttle (Harmon et al. 2010).

Fig. 7 Observed relative lightcurves used in our shape modelling analysis, compared with best-fit model lightcurves. The two models represent two slightly different scattering functions. They are the Lommel-Seeliger (solid line), and a combination of Lommel-Seeliger and Lambertian scattering laws (dashed line). The overall fit is quite good with the exception of a few lightcurve segments. These cases may be due to several factors including: the presence of one or more large concavities on the surface, which this modelling procedure is not sensitive to; non-uniform surface albedo; and imperfect scaling of the lightcurves to the data, accentuating differences between the two.

3.2. Phase darkening behaviour

We used our new data to determine an independent value for the phase darkening coefficient. Our model allowed us to reliably incorporate our snapshot observations taken at the NTT in February 2004 (α = 11.47°) and unpublished magnitudes from the VLT, taken at phase angles of 10.1–12.5° and thus greatly extending our phase angle coverage. The problem with these datasets is that the rotational phase is not known, which therefore could lead to greatly increased uncertainties in the measured phase darkening coefficient if such data are included (it’s the lightcurve midpoints that we need to determine). Our shape and spin-state model provides a solution to this issue. For each nightly dataset we generated a model lightcurve, that allowed us to ascertain what shift factor needed to be applied to the average measured magnitude of each dataset to ascertain the lightcurve midpoint magnitude for each dataset. We repeated the process for all other datasets even though lightcurves were obtained. In these cases we generated the model lightcurve for each nightly data set, shifted each magnitude datapoint accordingly, then took the nightly averages of the corrected magnitudes. The uncertainty is the standard deviation of the resulting magnitudes.

Fig. 8 We used just the HG system (Bowell et al. 1989) to characterize the phase darkening behaviour of the nucleus and its absolute magnitude, HR(1,1,0). Our best-fit values for both parameters are: G = 0.11 ± 0.12 and HR(1,1,0) = 15.31 ± 0.07 which are determined from chi-squared fitting, results of which are shown in the left panel along with the 1–3σ co-variance ellipses. This fit is very satisfactory with very little scatter, as illustrated on the right panel. We point out that our NTT data from 10th May, 2005 (α = 0.1°) has been excluded from the analysis above due to a suspected calibration error of unknown source leading to magnitudes 0.35 mag fainter than might be expected. This has no bearing on our shape modelling analysis as we use relative lightcurves. While not included in our fit, we overplot the HST data (Lamy et al. 2006) and the VLT data (Tubiana 2008) for comparison. See Sect. 3.2 for more discussion.

We used the HG system (Bowell et al. 1989) to characterize the phase darkening behaviour of the nucleus and its absolute magnitude, H_R(1,1,0). Our best-fit values for both parameters are: G = 0.11 ± 0.12 and H_R(1,1,0) = 15.31 ± 0.07 which are determined from chi-squared fitting, the results of which are shown in Fig. 8 along with the 1–3σ co-variance ellipses. This fit is very satisfactory with very little scatter. We also fit a linear phase function, in the form m_R(1,1,α) ∝ βα, for comparison with previous results. A least squares fit gives a phase function slope of β = 0.059 ± 0.006 and H_R(1,1,0) = 15.43 ± 0.04. The fit is reasonably acceptable when the extremes of the errors are taken into account, which implies that an opposition surge may not be required to explain the phase behaviour. The derived slope is shallower than that found by Tubiana et al. (2008) (β = 0.076 ± 0.003), but closer to the average value for Jupiter Family comets of β = 0.053 (Snodgrass et al. 2011). While the formal uncertainty in our β value is larger than Tubiana et al., we note that the error bar on their result is likely underestimated, as it is derived from a global search to solve simultaneously for phase function and rotation period, rather than the formal uncertainty on the linear fit to the points shown in their Fig. 5.

While not included in our fit, we overplot the HST data (Lamy et al. 2006) and the VLT data (Tubiana 2008) for comparison. The scaling method outlined above was also applied to these data for consistency. The scatter is increased somewhat, perhaps due to imperfections in the calibration process, or over-removal of coma in the case of the HST data. An independent treatment of the data for consistency with our reduction methods is beyond the scope of this work, so we opt to quote best-fit HG parameters using just our NTT data and the unpublished magnitudes from the VLT. These parameters are important for the thermophysical modelling discussed in Sect. 4 below.

3.3. Additional discussion

Our shape model with a constant sidereal rotation period fits all data very well and now allows accurate predictions of rotational phase, lightcurve shape, and thermal-IR flux levels throughout the observing period of 2004–2007 and beyond (also see Sect. 4 below). The quoted formal 1-σ uncertainty for the rotation period will lead to a degradation of rotational phase prediction of just  ~1.2° per year. New observational lightcurve data can easily be incorporated into our model for more accurate phase and flux predictions at encounter by Rosetta. Furthermore, a change in spin rate between the two successive aphelion passages may even be detected. Such a change in rotation period brought about by localised outgassing would not be surprising, and has been detected on several comets to date. These are comets 9P/Tempel 1 (Belton et al. 2011), 103P/Hartley 2 (Samarasinha et al. 2011), and 10P/Tempel 1 (Knight et al. 2011).

New lightcurve extraction procedures are constantly being developed and refined, which are particularly powerful at extracting rotational signatures from moving objects in very crowded stellar fields (Snodgrass et al., in prep.). We encountered this issue during the July 2007 data set, as did Tubiana et al. (2008) during their near-simultaneous observations at the VLT. While each group made every effort to negate this issue, there is scope to improve the quality of the lightcurves obtained and to rerun the modelling procedure.

4. Thermophysical modelling

The lightcurve-derived shape model and spin state can be combined with thermal-IR flux measurements to derive the size, albedo, thermal inertia, and surface roughness of the nucleus of 67P/C-G by using a suitable thermophysical model. Rotation-resolved thermal-IR flux measurements of 67P/C-G have been obtained by the Spitzer Space Telescope (SST) in February 2004 at 24 μm (Lamy et al. 2008), and in May 2007 at 8 and 24 μm (Kelley et al. 2009). Both data sets were taken post-aphelion when 67P/C-G was at heliocentric distances of 4.5 and 4.8 AU respectively (see Table 1). A dust trail was still visible in the thermal-IR images of both data sets, which required subtraction of a modelled dust trail flux from the total flux in order to arrive at the flux emitted by the nucleus. Using an earlier version of the 67P/C-G shape model and a thermophysical model, Lamy et al. (2008) derive a mean effective radius of 1.93–2.03 km, and an R-band albedo of 0.039–0.043 using an R-band absolute magnitude of 15.46. These quantities were derived by making the assumptions that the level of thermal inertia was low (i.e.  <20 J m^-2 K^-1 s^−1/2) and that the degree of surface roughness was high (i.e. a beaming parameter of 0.85), as these properties could not be inferred from the single wavelength observations. They also note that their model cannot be fitted to the entire set of data points. In contrast, Kelley et al. (2009) use the simpler NEATM model of Harris (1998) to derive a mean effective radius of 2.04 ± 0.11 km, a beaming parameter of 0.68 ± 0.06, and an R-band albedo of 0.054  ±  0.006 using an R-band absolute magnitude of 15.35 ± 0.04. The low value of the beaming parameter strongly suggests that the surface has a low thermal inertia and a high degree of surface roughness. They also attempt to phase their thermal-IR observations with those of Lamy et al. (2008) and with R-band optical lightcurves, and notice several disagreements between the different data sets. It is speculated that the disagreements may indicate deviations in scattering and/or thermal behaviours of the surface, such as variations in albedo, thermal inertia, or surface roughness. In this section, we use these thermal-IR observations with the new shape model and spin state to provide the first measurements of the level of thermal inertia and degree of surface roughness of the 67P/C-G nucleus by using the Advanced Thermophysical Model (Rozitis & Green 2011). Updated values of the mean effective radius and R-band albedo are also derived during the model fitting.

4.1. Thermal-infrared flux fitting

The Advanced Thermophysical Model (ATPM) uses a detailed shape model, described by a mesh consisting of 1000 s of triangular facets, to predict the thermal-IR emission for a range of surface and nucleus properties. For each facet of the shape model, unresolved surface roughness is represented by a fractional coverage, f_R, of hemispherical craters (consisting of 100 facets), with the remaining fraction, (1 − f_R), representing a smooth flat surface. The hemispherical crater is a simple way to accurately reproduce the thermal-IR beaming effect produced by a range of surface roughness morphologies and spatial scales, and has been verified by application to lunar data (Rozitis & Green 2011). Other approaches may be applicable to comet surfaces but we do not consider them here. For specified surface and nucleus properties, the ATPM computes the surface temperature variation for each shape and roughness facet during a rotation by solving the 1D heat conduction equation with a surface boundary condition that includes direct and multiple scattered solar radiation, shadowing, and re-absorbed thermal-IR radiation from interfacing facets. Sublimation and outgassing of volatile ices is not included; however, this is a valid approximation for interpreting the 67P/C-G Spitzer observations since it is expected to be very minimal and limited to very small areas of the surface at the heliocentric distance 67P/C-G was at during the time of the observations. To determine the thermal-IR emission, a Planck function is applied to the derived temperatures and summed across visible facets to give the emitted flux as a function of observation wavelength, rotation phase, and surface and nucleus properties. The model thermal-IR fluxes can then be compared with observations to infer what the likely properties are.

For 67P/C-G, the free parameters to be constrained by fits to the thermal-IR observations are the effective radius, albedo, thermal inertia, and surface roughness. The effective radius and albedo are related to one another and can be considered a single free parameter. The geometric albedo, p, is related to the effective radius in km, R, by: (1)where H is the absolute magnitude of the nucleus, and M is the apparent magnitude of the Sun (Russell 1916; Jewitt 1991). Here the R-band magnitudes are used, which are H_R = 15.31 ± 0.07 for 67P/C-G from above and M_R = −27.14 for the Sun. The geometric albedo is related to the effective Bond albedo of the nucleus surface, , by: (2)where q is the phase integral that measures the angular dependence of scattered solar radiation. Lamy et al. (2008) assume a value of q = 0.27 that was measured for comet 19P/Borrrelly by Buratti et al. (2004). However, we choose to use our measured value of G to derive q based on the following relation q = 0.290 + 0.684G (Bowell et al. 1989). Since the degree of surface roughness influences the effective Bond albedo it is related to the Bond albedo of a smooth flat surface, A_B, by; (3)where f_R is the roughness fraction as before. Therefore, each effective radius and roughness fraction combination leads to a unique Bond albedo value to be used in the ATPM. However, it is computationally expensive to run ATPM for each value separately. Like Wolters et al. (2011), we run ATPM for just one Bond albedo value and perform a pseudo-correction to the predicted flux for different values. This flux correction factor, FCF, is given by; (4)where A_B is the smooth surface Bond albedo determined by inversion of Eq. (3), and is the model Bond albedo used in the ATPM. A value of is assumed, and as the albedo of 67P/C-G is very low, most of the flux correction factors used in the fitting process were very near unity. Like Lamy et al. (2008), the thermal-IR emissivity, ε, is assumed to be 0.95. Separate thermophysical models were run for thermal inertia values, Γ, ranging from 0 to 100 J m^-2 K^-1 s^−1/2 using the best-fit shape model and spin state derived from the lightcurve inversion. Using 400 time steps and 40 depth steps to solve the 1D heat conduction equation, the ATPM typically takes between 10 and 100 revolutions to reach convergence.

The model flux predictions, F_MOD(λ_n, ϕ_n, Γ, R, f_R), were then compared with the observations, F_OBS(λ_n, ϕ_n), and observational errors, σ_OBS(λ_n, ϕ_n), by varying the effective radius, thermal inertia, and roughness fraction to give the minimum-chi-squared fit; (5)for a set of N observations with wavelength λ_n and rotation phase ϕ_n. The total absolute errors given in Lamy et al. (2008) and Kelley et al. (2009) were used during the fitting process. As the nucleus flux measurements required subtraction of a modelled dust trail flux from the total flux to arrive at them, the nucleus flux measurement errors may not be normally distributed, which would result in elevated χ² values. Therefore, we chose a parameter region bounded by a constant Δχ² value at the 3-σ confidence level to define the range of possible parameters. To ensure that no outlier data points and/or underestimated error bars were placing heavy constraints on this parameter region, data points that were more than 3-σ error bars away from the best-fit models were excluded during the fitting process. Also, to take into account the absolute magnitude and phase parameter uncertainties in the albedo calculations, a range of possible H and G values to use were randomly selected from normal distributions with means and widths equal to their nominal values and uncertainties respectively. As the Lamy et al. (2008) and Kelley et al. (2009) data sets used slightly different methods/techniques to obtain the nucleus fluxes we do not combine them in the fitting process. Instead, chi-square fitting is performed on them separately to obtain consistency checks on the derived parameter values. As the Lamy et al. (2008) data set is only at a single wavelength, it cannot be used to derive the level of thermal inertia and, instead, the mean value derived from the Kelley et al. (2009) data set is assumed for it.

4.2. Results and discussion

Figure 9 shows example model fits to the thermal-IR data as a function of rotation phase, Figs. 10 and 11 show the parameter distributions for the possible models, and Table 3 summarises the derived properties for 67P/C-G. During the model fitting process, only 1 out of 16 data points of the Lamy et al. (2008) data set was found to be more than 3-σ away from the best-fit model and was therefore excluded. All 22 data points of the Kelley et al. (2009) data set were found to be within 3-σ error bars of the best-fit model and none were therefore excluded. The best-fit model to the Kelley et al. (2009) data set has zero thermal inertia and a fully rough surface (i.e. f_R = 1). This indicates that the surface is essentially in instantaneous equilibrium between absorbed sunlight and emitted thermal-IR radiation. However, within the 3-σ confidence region, thermal inertia values up to 15 J m^-2 K^-1 s^−1/2 and roughness fractions as low as 0.65 are possible but with decreasing probability. The high level of the 8 μm flux (i.e. before the blackbody peak) relative to the 24 μm flux (i.e. after the blackbody peak) places strong constraints on the possible values of thermal inertia and roughness fraction, which forces them to be low and high respectively. As shown by the different line styles in Fig. 9, increasing the thermal inertia slightly decreases the overall level of the 8 μm flux but does not adjust the 24 μm flux level. By weighting each possible model equally, mean values of the thermal inertia and roughness fraction of 4 J m^-2 K^-1 s^−1/2 and 0.9 are obtained. The effective radii and albedo are almost normally distributed and have values of 1.97 ± 0.04 km and 0.060 ± 0.005 respectively. By assuming a thermal inertia value of 4 J m^-2 K^-1 s^−1/2, chi-square fitting to the Lamy et al. (2008) data set determines values of 2.01  ±  0.05 km, 0.058  ±  0.005, and 0.7  ±  0.2 for the effective radius, albedo, and roughness fraction respectively. These values are very consistent with those derived from the Kelley et al. (2009) data set.

Fig. 9 ATPM fits to the 67P/C-G a) Kelley et al. (2009) data set; and b) Lamy et al. (2008) data set. The data points are shown by the markers where their error bars represent the total absolute 1-σ uncertainty, and the model fits are shown by the lines. The triangle represents a data point that is more than a 3-σ error bar away from the best model fit. The long-dashed, solid, short-dashed, and long-dashed-dotted lines represent thermal inertia values of 0, 4, 10, and 15 J m-2 K-1 s−1/2 respectively. The greyed out regions represent the 3-σ range of possible thermal-IR fluxes caused by non-uniform surface roughness for a thermal inertia of 4 J m-2 K-1 s−1/2.

As shown in Fig. 9, model fits to the Kelley et al. (2009) data set recreate the shape and phase of its variations only roughly, whilst the model fits to the Lamy et al. (2008) data set recreate those variations perfectly. The fit to the Lamy et al. (2008) data is much better than the model presented in their paper in which they were only able to fit to the main peak. Although all of the data points in the Kelley et al. (2009) data set fell within 3-σ error bars of the best-fit model, the flux measurements could still be contaminated slightly by dust tail and/or coma emissions which might explain why a perfect fit to the 2007 data was not obtained. The non-perfect fit could also be caused by non-uniform surface properties as already suggested by Kelley et al. (2009). To see how much variation in the emitted thermal-IR flux could be caused by non-uniform surface properties, we allowed the roughness fraction to vary across the surface. Following the method that was first introduced in Rozitis & Green (2012), we randomly divided the surface into 10 patches, and then assigned each patch a unique roughness fraction. These roughness fractions were randomly selected from a normal distribution with a mean of 0.9 and a 1-σ width of 0.1. However, where a randomly selected roughness fraction was less than zero it was assumed to be zero, and where a randomly selected roughness fraction was greater than one it was assumed to be one. 1000 independent realisations of the possible surface roughness distributions were then made to evaluate the degree of scatter in the emitted thermal-IR flux. The greyed out regions in Fig. 9 show the 3-σ range of the possible thermal-IR fluxes for a thermal inertia of 4 J m^-2 K^-1 s^−1/2 as a function of rotation phase. The degree of scatter is not that large and we are unable to tell if there are any variations in surface roughness from it; however, it may explain one or two slightly anomalous data points.

67P/C-G’s derived low thermal inertia of  <15 J m^-2 K^-1 s^−1/2 and high roughness fraction of  >0.65 are consistent with its low beaming parameter of 0.68  ±  0.06 determined by Kelley et al. (2009). To accurately compare this thermal inertia value to that of other planetary bodies, it must be scaled to a common heliocentric distance to take into account the temperature dependence of thermal conductivity (Delbo et al. 2007). Assuming a common heliocentric distance of 1 AU, the scaled thermal inertia value for 67P/C-G is  <50 J m^-2 K^-1 s^−1/2, which is comparable with that measured for other comets, e.g.  ~5 J m^-2 K^-1 s^−1/2 for Comet 2P/Encke (Lamy et al. 2003),  <70 J m^-2 K^-1 s^−1/2 for Comet 9P/Tempel 1 (Lisse et al. 2005; Groussin et al. 2007),  <100 J m^-2 K^-1 s^−1/2 for Comet 22P/Kopff (Groussin et al. 2009), and  <10 J m^-2 K^-1 s^−1/2 for Comet 8P/Tuttle (Boissier et al. 2011). The scaled 67P/C-G value is less than the thermal inertia of  ~50 J m^-2 K^-1 s^−1/2 measured for lunar regolith (Rozitis & Green 2011, and references therein), which suggests its surface consists of finer grained material. It is also much less than the average thermal inertia of 200 ± 40 J m^-2 K^-1 s^−1/2 determined for km-sized near-Earth asteroids (Delbo et al. 2007). Finally, the roughness fraction of  >0.65 produces an rms slope of  >40°, which is greater than the  ~32° measured for the Moon (Rozitis & Green 2011). This is to be expected as the thermal skin depth of  ~1 mm for 67P/C-G is much shorter than  ~1 cm for the Moon, and therefore allows thermal-IR beaming to be induced by a wider range of surface roughness spatial scales.

Fig. 10 67P/C-G parameter distributions for the possible models derived from ATPM fitting to the Kelley et al. (2009) data set. The black, white, and grey columns correspond to parameter distributions derived at the 1-, 2-, and 3-σ confidence levels respectively.

Fig. 11 67P/C-G parameter distributions for the possible models derived from ATPM fitting to the Lamy et al. (2008) data set. The black, white, and grey columns correspond to parameter distributions derived at the 1-, 2-, and 3-σ confidence levels respectively. The thermal inertia is kept fixed at 4 J m-2 K-1 s−1/2.

5. Summary of main conclusions

We report on optical observations that were acquired of comet 67P/C-G at the 3.6 m ESO New Technology Telescope from February 2004 to September 2007, when the comet was at heliocentric distances of 4.4–5.6 AU. We performed a detailed photometric analysis of the data to extract the nucleus lightcurves, for linkage with data available in the literature as well as some unpublished snapshot observations, for shape-model and thermophysical characterisation of the nucleus. The extent of the data that is now available on this comet presents a unique opportunity for an analysis of this kind. Our main results are as follows:

• 1.

  Our initial attempt at defining the side-real rotation period simultaneously with thepole orientation and shape yielded a value ofP_sid = 12.761 h. Using our best fit pole solution we refined the sidereal rotation period measurement and associated uncertainty considerably to P_sid = 12.76137 ± 0.00006 h. This was achieved by applying the Phase Angle Bisector (PAB) approximation method, which entailed shifting each datapoint in time according to its PAB orientation at the time of the observation and some initial measurement of its rotation period, then applying a straight forward Fourier analysis on the time-shifted data. This was iterated until convergence. The corresponding pole solution had the following Ecliptic coordinates: λ = 78°(± 10°), β = 58°(± 10°) (J2000). Uncertainties are 3-σ. A second prograde pole solution was evident but was not consistent with separate measurements determined by completely independent methods, and was therefore discounted.

• 2.

  A best-fit convex triangular-facet shape model was generated with axial ratios b/a = 1.239 and c/a = 0.819. These values can vary by as much as 7% in each axis and still result in a statistically significant fit to the observational data, at the 3-σ level.

• 3.

  The nucleus phase darkening behaviour was measured and best characterized using the IAU HG system. Best fit parameters are: G = 0.11 ± 0.12 and H_R(1,1,0) = 15.31 ± 0.07. A linear phase function is also acceptable within the measurement uncertainties. A least squares fit gives a phase function slope of β = 0.059 ± 0.006 and H_R(1,1,0) = 15.43 ± 0.04. The phase coefficient is close to the average value for the Jupiter-Family comet population.

• 4.

  We performed a detailed surface thermal analysis, by incorporating our shape and spin state model and optical photometry into the new Advanced Thermophysical Model developed by Rozitis & Green (2011, 2012), along with all available Spitzer thermal-IR data, in this case 8–24 μm flux measurements from Lamy et al. (2008) and Kelley et al. (2009). Our shape model combined with the ATPM can satisfactorily reconcile all optical and thermal-IR data, with the fit to the Lamy et al. data being exceptionally good. We derive a range of mutually-consistent physical parameters for each thermal-IR data set, including effective radius, geometric albedo, surface thermal inertia and roughness fraction. These values are summarized in Table 3.

• 5.

  67P/C-G’s derived low thermal inertia of  <15 J m^-2 K^-1 s^−1/2 and high roughness fraction of  >0.65 are consistent with its low beaming parameter of 0.68  ±  0.06 determined by Kelley et al. (2009). This value of thermal inertia is comparable with that measured for other comets once scaled to a common heliocentric distance, and implies a surface regolith finer than lunar surface material.

• 6.

  Our shape model with a constant sidereal rotation period fits all data very well and now allows accurate predictions of rotational phase, lightcurve shape, and thermal-IR flux levels throughout the observing period of 2004–2007 and beyond. The quoted formal 1-σ uncertainty for the rotation period will lead to a degradation of rotational phase prediction of  ~1.2° per year.

• 7.

  The comet’s dust trail was too faint for detection during our observations at the NTT.

Acknowledgments

We thank the anonymous referee for their very helpful comments. All image reduction and processing were performed using the Image Reduction and Analysis Facility (IRAF) (Tody 1986, 1993). IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. JPL Horizons. Support is acknowledged from the UK Science Technology Facilities Council (STFC) (S. Lowry, S. Duddy, S. Green, and B. Rozitis), and the Southeast Physics Network (SEPnet) (S. Lowry). The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 268421 (C. Snodgrass).

References

All Tables

All Figures

Fig. 1 Composite R-filter images obtained at the NTT in May 2005. Each panel consists of 71 × 100 s, 123 × 100 s, and 82 × 100 s, images stacked according to the comet’s projected rate of motion. In each case, as with all our data, the comet appears as a sharp point source on our images. Celestial North and East are marked along with the projected solar direction and heliocentric velocity vectors.

In the text

Fig. 2 Upper panels: left – periodogram produced via least-squares fitting of a Fourier series to the raw lightcurve data. The most likely period solutions are 12.40971 ± 0.00005 h and 12.76236 ± 0.00006 h. Middle and right – folded lightcurves of the uncorrected data points for each of these two best-fit rotation periods. The fit in the 12.41 h case is unconvincing with various datasets not following the phase of the folded curve. The fit for the 12.76 h case is better, but can be further improved by correcting for the differing observational geometry between the various datasets. Lower left: periodogram produced via least-squares fitting of a Fourier series to the lightcurve data which has been corrected for the movement of the PAB. In this case the most likely solution was 12.76137 ± 0.00006 h. Lower right: folded lightcurve using the data corrected for the movement of the PAB vector. The scatter and phasing of the datapoints is improved, leading to a more reliable period measurement.

In the text

Fig. 3 Contour plot showing the expected location of the rotation pole produced by inverting the lightcurves at each potential pole location. Darker colours here indicate lower chi-squared values. Ecliptic latitude and longitude are denoted by β and λ, respectively. The plot indicates that the pole is most likely located at high latitude, with two prominent and well-defined minima at λ = 78° and β =  + 58° and a 2nd unrelated solution at λ = 278° and β =  + 58°. (A colour version of this figure is available online.)

In the text

Fig. 4 Comparison of our rotation pole with others from the literature. For consistency with other authors, mainly Lamy et al. (2007), we compare the various measurements in RA-Dec space. We first converted all previous values to Ecliptic coordinates if not done so already, then took the regions defined by these values and their associated uncertainties and mapped each region onto the celestial sphere. Key: S06 – Schleicher (2006), C04 – Chesley (2004), V12 – Vincent et al. (2012), L07 – Lamy et al. (2007), DG05 – Davidsson & Gutierrez (2005). Our primary solution at λ = 78° and β =  + 58° (or RA = 51°, Dec = +79° – Solution A), is in excellent agreement with those values determined from completely independent methods. We have a 2nd unrelated solution at λ = 278° and β =  + 58° (or RA = 275°, Dec = +35° – Solution B), which is well outside the range of all the other solutions. With the constraints of the other studies in mind, we discount Solution B. An inherent property of lightcurve inversions is the inevitable existence of “conjugate” solutions, z-axis “mirrors” of the best fit shape model(s), for which alternate best-fit pole solutions are possible (i.e. at λ ± 180°, Kaasalainen & Lamberg 2006). For completeness, we investigated these conjugates by generating the mirrored shape models for each of the two best fit models (at Solutions A+B in Fig. 3), to see where their preferred pole solutions lie (marked as Ac and Bc on the plot). We find that these solutions result in much poorer fits to the data, especially solution Ac, which is barely defined at all. Note that the Lamy et al. (2007) value is also based on lightcurve inversion, but on just a small subset of the data used here. They also find a 2nd solution and its conjugate, marked as “L07 (B)” and “L07 (Bc)”, respectively.

In the text

	Fig. 5 Best-fit convex triangular-facet shape model produced via lightcurve inversion using the measured sidereal rotation period of 12.76137 h and the most likely pole orientation, λ = 78° and β =  + 58°. Axial ratios are well constrained at b/a = 1.239 and c/a = 0.819.
In the text

Fig. 6 To investigate how “flexible” we can be with our best-fit shape while still producing a statistically meaningful fit to the data, we used our best-fit shape model and transformed the relative sizes of the axial dimensions of the model by stretching the model at each principal axis at random increments to include all realistic axial ratios. We then took all the new shape model permutations to generate new sets of model lightcurves for comparison with the data via chi-squared fitting. We find that a, b, and c can vary by as much as 7%, and still result in a statistically meaningful fit, at the 3-σ level. Lower right panel: given this measure of flexibility in the shape, we find that the uncertainties in the pole coordinates are revised upwards to  ~±10°, indicated by the outer curve. The grey band indicates how the outer curve can vary when the sidereal rotation period is allowed to vary within its associated uncertainty. See Sect. 3.1 for full details.

In the text

Fig. 7 Observed relative lightcurves used in our shape modelling analysis, compared with best-fit model lightcurves. The two models represent two slightly different scattering functions. They are the Lommel-Seeliger (solid line), and a combination of Lommel-Seeliger and Lambertian scattering laws (dashed line). The overall fit is quite good with the exception of a few lightcurve segments. These cases may be due to several factors including: the presence of one or more large concavities on the surface, which this modelling procedure is not sensitive to; non-uniform surface albedo; and imperfect scaling of the lightcurves to the data, accentuating differences between the two.

In the text

Fig. 8 We used just the HG system (Bowell et al. 1989) to characterize the phase darkening behaviour of the nucleus and its absolute magnitude, HR(1,1,0). Our best-fit values for both parameters are: G = 0.11 ± 0.12 and HR(1,1,0) = 15.31 ± 0.07 which are determined from chi-squared fitting, results of which are shown in the left panel along with the 1–3σ co-variance ellipses. This fit is very satisfactory with very little scatter, as illustrated on the right panel. We point out that our NTT data from 10th May, 2005 (α = 0.1°) has been excluded from the analysis above due to a suspected calibration error of unknown source leading to magnitudes 0.35 mag fainter than might be expected. This has no bearing on our shape modelling analysis as we use relative lightcurves. While not included in our fit, we overplot the HST data (Lamy et al. 2006) and the VLT data (Tubiana 2008) for comparison. See Sect. 3.2 for more discussion.

In the text

Fig. 9 ATPM fits to the 67P/C-G a) Kelley et al. (2009) data set; and b) Lamy et al. (2008) data set. The data points are shown by the markers where their error bars represent the total absolute 1-σ uncertainty, and the model fits are shown by the lines. The triangle represents a data point that is more than a 3-σ error bar away from the best model fit. The long-dashed, solid, short-dashed, and long-dashed-dotted lines represent thermal inertia values of 0, 4, 10, and 15 J m-2 K-1 s−1/2 respectively. The greyed out regions represent the 3-σ range of possible thermal-IR fluxes caused by non-uniform surface roughness for a thermal inertia of 4 J m-2 K-1 s−1/2.

In the text

	Fig. 10 67P/C-G parameter distributions for the possible models derived from ATPM fitting to the Kelley et al. (2009) data set. The black, white, and grey columns correspond to parameter distributions derived at the 1-, 2-, and 3-σ confidence levels respectively.
In the text

	Fig. 11 67P/C-G parameter distributions for the possible models derived from ATPM fitting to the Lamy et al. (2008) data set. The black, white, and grey columns correspond to parameter distributions derived at the 1-, 2-, and 3-σ confidence levels respectively. The thermal inertia is kept fixed at 4 J m-2 K-1 s−1/2.
In the text