Predictions for Sparse Photometry of Jupiter-Family Comet Nuclei in the LSST Era

The Legacy Survey of Space and Time (LSST) at Vera C. Rubin Observatory will deliver high-quality, temporally-sparse observations of millions of Solar System objects on an unprecedented scale. Such datasets will likely enable the precise estimation of small body properties on a population-wide basis. In this work, we consider the possible applications of photometric data points from the LSST to the characterisation of Jupiter-family comet (JFC) nuclei. We simulate sparse-in-time lightcurve points with an LSST-like cadence for the orbit of a JFC between 2024-2033. Convex lightcurve inversion is used to assess whether the simulation input parameters can be accurately reproduced for a sample of nucleus rotation periods, pole orientations, activity onsets, shapes and sizes. We find that the rotation period and pole direction can be reliably constrained across all nucleus variants tested, and that the convex shape models, while limited in their ability to describe complex or bilobed nuclei, are effective for correcting sparse photometry for rotational modulation to improve estimates of nucleus phase functions. Based on this analysis, we anticipate that LSST photometry will significantly enhance our present understanding of the spin-state and phase function distributions of JFC nuclei.


INTRODUCTION
Much of our present understanding of comet physical properties has been derived from observations of their nuclei in broad-band optical filters.Jupiter-family comets (JFCs) are an intriguing and informative population of comets to study in this way.The source populations of these objects lie beyond the orbit of Neptune in dynamically unstable regions of the Kuiper Belt before being propelled to low-inclination, short-period (< 20 years) orbits following gravitational interactions with the giant planets (Levison & Duncan 1997).Studying the physical characteristics of JFCs provides insights into the processes they experience during their time in the inner Solar System, which is essential for contextualising their evolution from their reservoirs of origin.
As comet nuclei are km-sized, these objects typically remain observable by 1-4m class observing facilities over most of their orbits.Near the Sun, the temperature is high enough that their surfaces are shrouded to the Corresponding author: A. Donaldson a.donaldson@ed.ac.uk ground-based observer by a diffuse coma of gas and dust driven by the sublimation of volatile species in the nucleus.For most JFCs, sublimative activity diminishes with increased heliocentric distance, at which point the dominant signal from the comet becomes sunlight reflected from the nucleus surface.When this is the case, the ground-based observer can estimate nucleus properties depending on the amount and quality of available data.As a starting point, multi-band snapshot photometry allows for the absolute magnitude of the nucleus to be measured in filters of interest, placing constraints on the colour and effective size of the nucleus (e.g.Licandro et al. 2000;Lowry et al. 2003).This technique is somewhat limited in that it captures nucleus information at only a single instance of its rotational phase.With sufficiently long sequences of dense-in-time imaging, typically in a single filter, the rotation period of the nucleus about its spin axis can be measured.This estimate improves in reliability with lightcurves collected across successive nights, probing larger portions of the object's rotational phase without notably altering the observing geometry of the nucleus with respect to the observer.
From lightcurves collected across multiple epochs spanning a range of solar phase angles, the dependence of nucleus brightness on phase angle can be estimated (e.g.Lamy et al. 2004;Snodgrass et al. 2011;Kokotanekova et al. 2017).Combined with reliable nucleus size measurements, the resulting phase function allows for the geometric albedo (p V ) of the nucleus to be estimated.With a large enough number of observational epochs, it is possible to fit a convex model to the collection of lightcurves to identify a likely shape and pole orientation for the nucleus (Kaasalainen & Torppa 2001;Kaasalainen et al. 2001).This was recently demonstrated by Donaldson et al. (2023) for a JFC, using absolutely-calibrated dense-in-time lightcurves of 162P/Siding Spring spanning more than fifteen years of observations.The resulting shape model indicated an elongated shape for the nucleus of 162P and was used to correct the phase function for rotational variation.A similar approach was used by Lowry et al. (2012) and Mottola et al. (2014) to accurately predict the spin state of JFC 67P/Churyumov-Gerasimenko prior to the arrival of the Rosetta spacecraft.Using lightcurve inversion to model the shape of a comet nucleus provides the most detailed information possible from the ground on its shape and spin state.Such comprehensive information for a greater number of JFCs may provide insights into the distribution of these properties for a more representative sample of the population and offer new perspectives into the mechanisms driving their activity behaviours and evolutionary processes.
Given the limited number of comets for which large amounts of reliable, dense-in-time lightcurve data are presently available, we consider the role that large-scale optical surveys might play in enhancing the determination of JFC nucleus properties.Optical surveys can obtain intermittent photometry of small bodies contained within their field of view -often fortuitously -with the regularity of data points for each object observed dependent on survey cadence and orbital parameters.In recent decades, survey telescopes have become increasingly powerful, reaching fainter limiting magnitudes and covering greater areas of sky, resulting in the discovery of thousands more small bodies.Kaasalainen (2004) showed it was theoretically possible to obtain reliable asteroid shape and spin estimates from solely sparse-intime photometry.Data from surveys such as the Asteroid Terrestrial-impact Last Alert System (ATLAS) and Zwicky Transient Facility (ZTF) have produced reliable shape and period estimates of asteroid populations solely using sparse photometry over long temporal baselines (e.g.Ďurech et al. 2020;Ďurech & Hanus 2023).
Such observations of comet nuclei are typically more challenging.To avoid contamination from activity in photometry apertures, we are restricted to observing them at large heliocentric distances (typically > 3 au) where they are faintest.Small survey telescopes are therefore unlikely to detect inactive comet nuclei at sufficient signal-to-noise for reliable photometry (for a recent summary of the contributions of telescopic surveys to comet studies, see Bauer et al. 2022).
The Legacy Survey of Space and Time (LSST) at the Vera C. Rubin Observatory on Cerro Pachón, Chile, will survey 18,000 deg 2 of the southern sky in six ugrizy filters over its 10-year period of operation (Ivezić et al. 2019;Bianco, F. & SCOC 2022).At the time of writing, the survey is expected to commence in mid-2025.The primary observing mode of the LSST is a 'widefast-deep' (WFD) optical survey covering the entire observable sky in approximately three-day intervals.Additional time will be allocated to several micro-surveys probing smaller sky regions to greater depth.The survey will be conducted with the Simonyi Survey Telescope, which has a primary mirror diameter of 8.4-metres.The design is such that the single exposure limiting object brightness in the r-band for the WFD mode is around 24th magnitude.It is currently predicted that the survey will discover more than five million new solar system objects, and will provide unprecedented numbers of calibrated photometric data points for the moving objects it detects.The large survey area coupled with the faint limiting magnitude have profound implications for the study of the rotational lightcurves of faint small bodies.LSST will likely be the first survey capable of reliably acquiring sparse photometry of inactive short period comets.Given the challenges presented by observing comet nuclei at a wide range of observing geometries with targeted observations, will nucleus properties such as pole orientation and shape be accessible with LSST?
The LSST science priorities encompass many astronomical disciplines in addition to inventorying the Solar System.Finding an observing strategy that will optimise the scientific returns in each of these areas is not straightforward.The intricacies of the survey cadence are therefore not yet finalised, and will likely continue to be tweaked throughout the early years of the LSST.This work does not aim to contribute to the discussion on survey cadence specifics, and the methods and results discussed here by no means reflect the exact final output from the LSST.Our goal is to examine the potential of sparse-in-time LSST photometry for the characterisation of JFC nuclei.To date, no survey has been capable of obtaining high-quality nucleus photometry on such vast scales, and our focus is the application of this photometry to measuring nucleus properties on an individual object basis.
In this work we simulate observations of typical JFC nuclei with an LSST-like cadence.Section 2 describes the process for generating synthetic sparse-in-time nucleus lightcurves using existing survey simulator tools.Section 3 details the lightcurve analysis methods we employ by application to a noiseless example case.By individually varying the input parameters of the nucleus model, we examine how intrinsic nucleus properties affect the outcomes of lightcurve analysis in Section 4. In Section 5 we incorporate the effects of photometric uncertainty on the synthetic lightcurve points and analyse these using a Monte Carlo approach.We also evaluate combining the sparse survey data points with simulated dense-in-time lightcurves.Finally in Section 6 we assess our findings in a realistic context beyond the simplifying assumptions made to generate the lightcurve datasets, and discuss their implications for how reliably we can expect to measure JFC properties using LSST photometry.

Simulating lightcurve points with an LSST-like cadence
We first consider the temporal spacing of the synthetic lightcurve points.The timestamps i.e. when the comet is contained within the LSST field of view (FOV), depend on both the observing cadence of the survey and the orbital parameters of the comet.To create these timestamps, we use the objectsInField (OIF) module within the Asteroid Survey Simulator package (Cornwall et al. 2020).For an input population of small bodies, OIF identifies field pointings in which the specified objects are detected.We chose the orbit of comet 67P/Churyumov-Gerasimenko (hereafter 67P) as it is representative of a 'typical' JFC orbit and is not expected to undergo any significant orbital modifications during the LSST operation period.We used a polygonal survey FOV with an effective area of 10.5 deg 2 to simulate the array of pixel sensors on the LSST Camera (Ivezić et al. 2019, Fig. 12).The survey cadence used was baseline 3.0, the most current database of simulated cadence information available at the time.The database contains simulated observational properties of each pointing such as filter, visit start time, and the limiting magnitude at that pointing.It should be noted that baseline 3.0 is based on the recommendations of the Survey Cadence Optimization Committee in December 2022 and does not reflect the exact final observing strategy for LSST -rather, it incorporates updates from pre-vious versions of the observing strategy based on feedback (Bianco, F. & SCOC 2022).
The resulting dataset from which we construct the nucleus lightcurves consists of 613 individual field pointings between January 2024 and September 2033.The median effective seeing FWHM is 1.03 arcseconds for these pointings, and throughout this work we assume that the object's proper motion is sufficiently slow that it would not be trailed in the ∼30-second LSST exposures.In the baseline 3.0 cadence, most of the grizy pointings are comprised of two 15-second exposures in each filter, separated by a shutter-close time of 1 second and a readout time of 2 seconds.To define the lightcurve timestamps, we take the midpoint of these pointings to be 16.5 seconds after the start time of the field pointing.Some of the r-and i-band exposures are single 15-second exposures, for which we take the timestamp to be 7.5 seconds after the pointing start time.The uband pointings in the baseline 3.0 cadence are entirely single 30-second visits, chosen to minimise read noise in this filter.We take the midpoint as being 15s from the start time of the pointing.We correct each timestamp for light-travel time using the geocentric distance of the comet so that each timestamp depicts the time at which the light left the nucleus surface at the exposure midpoint.

Nucleus magnitudes
To associate a nucleus brightness with each timestamp, we must account for the nucleus spin state and photometric properties which vary from comet to comet.In the following we describe our procedure for generating magnitudes applied to what we refer to as the nominal case, the properties of which were chosen to represent a 'typical' comet nucleus.The values chosen for the parameters of the nominal case are those of 67P, with the exception of rotation period.This we took as the mean JFC rotation period from the compilation of periods by Knight et al. (2023), P i = 14.8978 h, rather than the true value of ∼12h.This was simply to rule out the possibility of confusing real periodic signals identified in the lightcurves with alias signals due to the 24-hour day in later analysis.The values for all the initial parameters are given in Table 1, and those defining the nominal case are indicated in bold.
We generated a value for the reduced magnitude H 0 (1,1,α) as a function of solar phase angle α at each light-time corrected timestamp using Equation 1.
The absolute magnitude H(1, 1, 0) is defined as the theoretical magnitude an object would have at heliocen-   tric distance (r h ) and geocentric distance (∆) 1 au and phase angle α = 0 • , where the object brightness scales linearly with phase angle by phase function β.We denote the reduced magnitude with subscript 0 to indicate that the effects of rotation have not yet been incorporated.We used the linear phase function parameters for 67P determined by Lowry et al. (2012) H R,0 (1,1,0) = 15.43 and β L 1 = 0.059 mag/deg.The LSST grizy filter bandpasses are similar to those of Pan-STARRS1 (PS1), therefore we convert all magnitudes in this work to the PS1 filter system.We assumed a nucleus colour (B − V ) = 0.87, the mean colour for JFC nuclei determined by Jewitt (2015), and used the conversions described in Tonry et al. (2012) to derive an absolute magnitude of H r,0 (1, 1, 0) = 15.64 in the PS1 r-band.
To these initial static reduced magnitudes we added the shape-driven effects of nucleus rotation.For the nominal case, we synthesised lightcurves using 67P shape model SHAP8, produced by spectrophotoclinometry of images taken in situ of the entire comet nucleus acquired with the Optical, Spectrocopic and Infrared Remote Imaging System (OSIRIS) instrument onboard Rosetta (Gaskell et al. 2008;Jorda et al. 2016).Various resolutions of this model are available -to minimise computation time, we used the lowest resolution shape model with 3000 facets.
1 subscript L to differentiate between model linear phase function, β L and initial ecliptic latitude, β i .
We assigned to the shape model a constant rotation period P i , and a fixed pole orientation (λ i , β i ) -values are given in Table 1.We used the Earth-comet and Suncomet vectors (Giorgini et al. 1996) at each timestamp along with the pole orientation to determine the viewing geometry of the nucleus at each epoch.The light scattering properties of the nucleus surface were modelled using a Lommel-Seeliger (LS) function 2 .Starting with an arbitrary initial rotational phase, we determined the total LS flux at each timestamp.We accounted for self-shadowing due to the non-convex shape of 67P by including the flux contributions of only facets that were illuminated, visible to the observer and unobscured by any other facets.Converting the total flux at each timestamp to a relative magnitude centred on zero gave the brightness shift to be added to the static reduced magnitudes to produce the rotational lightcurve.These reduced magnitudes were converted to apparent magnitudes by accounting for the comet's heliocentric and geocentric distances.
To associate an apparent magnitude uncertainty with each lightcurve point, we returned to the OIF output.This includes a fiveSigmaDepth value, providing an estimate of the limiting magnitude in the specified filter at which an object would be detected with a signal-tonoise ratio (SNR) of 5.At every timestamp we compared this limiting magnitude to the corresponding apparent 2 The flux contribution from each model facet is µµ 0 µ+µ 0 , where µ 0 and µ are the cosines of the angles of incidence and reflection respectively.magnitude, first converting from PS1 r-band into the filter associated with that pointing.We assume that the nucleus photometric properties are constant across all wavelengths.For the grizy filters we converted the average JFC nucleus colours from Jewitt (2015) to PS1 g − r, r − i, r − z and r − y colours using linear combinations of the polynomial colour conversions described in Tonry et al. (2012).As PS1 has no direct counterpart to the LSST u-band, to estimate the colour in this band we took the average JFC (V − R) = 0.51 from the compilation by Knight et al. (2023), noting that this is consistent with the value given in Jewitt (2015) but based on a greater sample of 49 objects.We determined the spectral slope S ′ corresponding to this colour following the method given in Luu & Jewitt (1990), yielding S ′ = 11.6 ± 0.1%/1000 Å.Assuming a constant spectral slope from the centre of the u passband to the centre of the r passband, we use this value for S ′ to determine a value for u − r.This value is shown with the other colour indices in Table 1.By approximating magnitude uncertainties as 1/SNR, we extrapolated magnitude error values from the 5-σ limiting magnitude of each pointing.If the uncertainty of a given point was larger than 0.2, it was assumed that the object was too faint at that instance to be detected in an LSST exposure and so it was discarded from further analysis.Note that we converted the synthetic apparent magnitudes to the filter corresponding to their timestamp solely to associate a representative uncertainty value with each lightcurve point.The final lightcurve is comprised of only r-band magnitudes and incorporates all points with sufficient signal-to-noise to be detected by LSST over the entire ten-year time frame.
Thus far we have assumed at all timestamps that the brightness is due entirely to reflected light from the nucleus surface, which is not the case for most JFCs.For targeted observations of comet nuclei, the standard strategy to minimise contamination from activityrelated brightening is to plan and undertake observations only when the comet is further than 3 au from the Sun where temperatures are insufficient for the significant sublimation of water ice.This strategy is not without flaw, and the PSF of all targeted observations should be carefully checked to ensure that the comet's radial profile is stellar-like.To simulate this strategy, for the nominal case lightcurve we implement a simple heliocentric distance cut-off, discarding lightcurve points where the comet is at r h < 3.0 au.This gives a final nucleus lightcurve as might be acquired by LSST, shown in the right panel of Figure 1.In the final dataset, the median number of lightcurve points per night is 2 for the 167 individual epochs in which the nucleus is both visible and detected.The comet spans a range of 183 • to 298 • in ecliptic longitude, −4 • to 2 • in latitude and a phase angle range 0.2 • to 16.7 • .

Simple phase function fit and period search
With the synthetic sparse-in-time lightcurves for the nominal case in hand, our aim was to establish whether the known input parameters could be reproduced by standard methods of lightcurve analysis.One of the most exciting prospects for nucleus studies with LSST is the unprecedented observational coverage across entire JFC orbits.In particular, measuring a nucleus phase function has historically required years of targeted data collection at as large as possible a range of solar phase angles.However, with LSST it will be possible to gather calibrated, temporally-sparse observations spanning a large α range for many objects.
The phase angle variation covered by the nominal lightcurve is shown in Figure 2, spanning 0.2 • − 16.7 • at r h > 3 au.Note that there is no noise in this dataset (this is considered in Section 5) -scatter in the points is due solely to nucleus rotation.We fit a phase function by linear regression to the reduced magnitudes, yielding a phase function coefficient β = 0.0565 ± 0.0014 mag/ • , slightly shallower than the input value β L = 0.059 mag/ • .We corrected the points for this phase function by Equation 1, and searched for the best rotation period by performing a Lomb-Scargle (LS) period search (Lomb 1976;Scargle 1982).LS is typically used to determine synodic rotation periods for observational datasets in which the object geometry is unchanging.Performing a period search over the long time-base lightcurve yields an estimate for the 'average' synodic period.We tested a single-term Fourier fit over a frequency space corresponding rotation periods between ∼1 h and 45 h.The resulting periodogram is shown in Figure 3.The highest periodogram peak corresponds to a rotation period P = 14.8981 h.This is close in value to the model input sidereal period P i = 14.8978 h, and phasing the lightcurves to the most likely Lomb-Scargle period results in a similar profile to that shown in the right panel of Figure 1.With the a priori knowledge of the model inputs, this is consistent with the input sidereal period since the difference between synodic and sidereal periods is typically small for such observations (Harris et al. 1984).For observational datasets, due to the large range of nucleus observing geometries that occur over the long   temporal baseline, it is challenging to select the best period from the possible LS peaks by visual inspection of the lightcurve morphology, particularly regarding the choice of a single or double-peaked lightcurve.

Convex lightcurve inversion
To test the benefits of the availability of a large range of observing geometries, we next assessed the synthetic lightcurves using convex lightcurve inversion (CLI).The rotation period, shape and aspect of an object directly impact the characteristics of its lightcurve measured by the ground-based observer, and thus as the observing geometry of the object changes across its orbit, so too does the appearance of its lightcurve (as can be seen in Figure 1).CLI (Kaasalainen & Torppa 2001;Kaasalainen et al. 2001) uses a Levenberg-Marquardt minimisation algorithm to iteratively fit the shape and spin state that most closely match the input lightcurves.This method is most effective when a wide range of object viewing geometries is covered by the input lightcurves, a diversity that is provided naturally by intermittent detections with survey telescopes.Convex inversion of solely sparse-in-time lightcurves is capable of producing extremely reliable estimates of asteroid spin states and good approximations of their overall shapes (Kaasalainen 2004).The primary requirement is that each data point across the entire baseline is calibrated to the same absolute brightness scale.The sparse points can then be combined and treated as a single lightcurve with a long temporal baseline.
For this work, we used the implementation of the convexinv software package for sparse lightcurves (see e.g.Hanuš et al. 2011).The lightcurve inversion software parameterizes an object by its rotation period (P ), pole orientation (ecliptic longitude λ and latitude β), light scattering properties and shape attributes.In the optimization process the shape is represented by the expansion of a spherical harmonic series.The parameters to be optimised are the coefficients of the series with degree l and order m.To strike a balance between reconstructing a detailed shape and avoiding over-fitting a model to a relatively small number of lightcurve points, we used l = m = 6, giving 49 shape parameters to be optimised.The light scattering function is described by an empirical three-coefficient linearexponential phase function (Kaasalainen et al. 2001).Following the method of Ďurech & Hanus (2023) the exponential coefficients of the phase function were held fixed, fitting solely the linear coefficient of the phase function.
We searched initially for the best-fitting sidereal rotation period.The periodsearch program within convexinv steps through a series of user-defined trial periods and iteratively fits a shape at each period from six starting pole orientations.The quality of the fit to the lightcurves at each trial period is quantified by a χ 2 value.The best period value is taken as the trial period corresponding to the minimum value of χ 2 .Figure 4 shows the resulting periodogram from trial periods in the range 5 − 45 h.The best period was found to be consistent with P i , with a value of P = 14.8978 ± 0.0001 h, where we express the uncertainty as the period values within 10 per cent of the χ 2 minimum (Rożek et al. 2022).The uncertainties found by this method are gen-erally very small as a result of the small period step size (described by e.g.Kaasalainen 2004).
We next searched for the optimal orientation of the object's rotation pole, with the starting value of P set to the above value.At each trial pole orientation on a 5 × 5 deg 2 grid of ecliptic coordinates, the shape was optimised over 50 iterations.The χ 2 evaluation of the fit quality for each possible pole location is shown in the right panel of Figure 4.There are two regions of equally-likely pole orientations, separated by approximately 180 deg in ecliptic longitude, as is expected for an object on a low-inclination orbit by the ambiguity theorem (Kaasalainen & Lamberg 2006).We defined unique pole solutions as those with a χ 2 value at least 10 per cent lower than the surrounding solutions, and which were separated from other minima by at least 30 deg.This resulted in two possible pole orientations corresponding to the regions of darkest blue shown in Fig- ure 4. Solution 1 at λ = 84 • , β = 40 • is consistent with the input coordinates used for the pole orientation, and Solution 2 at λ = 274 • , β = 35 • gives the mirror solution.These values are permitted to vary as part of the final shape optimisation procedure described below, but testing showed that the pole orientation at which the final model converges is close to the coordinates of the minima presented here.Therefore, to associate a formal uncertainty with these pole orientations, we consider the uncertainty region as that with χ 2 values within ten per cent of the χ 2 minimum value.This yields final pole solutions of λ = (84 ± 6) • , β = (40 ± 5) • for Solution 1, and λ = (274 ± 7) • , β = (35 ± 5) • for Solution 2.
Solution 1 can be directly compared to the results of other published convex inversions of 67P lightcurves: the first acquired using dense-in-time photometry by Lowry et al. (2012) who identified a pole orientation (λ, β) = (78 • ± 10 • , 58 • ± 10 • ); and with the addition of lightcurves from the approach phase of the Rosetta spacecraft, Mottola et al. (2014) obtained a most likely pole orientation (λ, β) = (65 • ± 15 • , 59 • ± 15 • ).These solutions are consistent with one another within their specified uncertainty ranges, and with the true orientation of 67P's rotation pole.
For each of these pole orientations, we developed a convex shape model.We used 50 iterations of the Levenberg-Marquardt loop to fit each model, unless the χ 2 values had not satisfactorily converged within this number of iterations.In most cases it was found that iterating beyond 50 when the χ 2 values had already converged (varied by < 10 −3 with each iteration) made minimal difference to the axis ratios and fit to the lightcurves of the final shape models, and only served to make the models 'blockier' and less physically realistic in appear- The axis ratios of the convex models are substantially less elongated than those of the input 67P model (see Table 1).This is to be expected for an object as strongly non-convex as 67P -the lightcurve inversion procedure cannot recreate surface concavities and instead tends to mask over them, which has the effect of artificially lengthening the rotation axis of the convex model in this case.We test the impact of different nucleus shapes later in Section 4.4.It is therefore more illustrative to compare these shapes to the previously published convex models of 67P obtained by Lowry et al. (2012) and Mottola et al. (2014).The axial ratios based on the inertia tensors of these models are a/b = 1.19, b/c = 1.25 and a/b = 1.14, b/c = 1.22 respectively.The sparse-intime model axis ratios derived in this work are remarkably similar to the previously published models.This is extremely promising for characterising a greater number of JFCs with LSST photometry: the diversity of observing geometries covered over the survey's 10-year period of operation is capable of providing shape and pole information that is equivalent to years of targeted monitoring.
While the convex models can provide valuable insights into the object's global shape properties, these insights are somewhat limited when the original shape is intrinsi-  cally non-convex, such as that of 67P.One can convince oneself that Model 1 is a gift-wrapped rendering of the shape of the bilobed nucleus, close to a faithful reconstruction of its convex hull.In contrast, since the model was optimised to fit each input lightcurve point, we can use the model to correct the lightcurve points for the effects of brightness variation due to the rotation of the irregular shape (analogous to the reference phase function of Kaasalainen et al. 2001).Following the method of Donaldson et al. (2023), all lightcurve points within a single rotation were binned and shifted by the mean correction.The shifted magnitudes were plotted as a function of phase angle α, shown in Figure 6.While the phase function of the rotationally-scattered points is slightly shallower than the input phase function β L , following the rotational-correction with the convex model β L is exactly recovered.

TESTING A RANGE OF COMETARY PROPERTIES
We have thus far tested LSST-like photometry for a single combination of possible comet physical properties.In reality, the properties of the well-studied JFCs span a wide range of values.We selected five parameters to vary that we considered likely to have the largest impact on the simulated lightcurves: rotation period; pole obliquity; nucleus shape; nucleus size; and heliocentric distance cutoff (below which it is assumed that bare nucleus observations are not possible due to the dominance of sublimation-driven activity in photometry apertures).This section aims to test whether the intrinsic nucleus properties affect the likelihood that the 'correct' results, i.e. the input parameters, are reproduced by lightcurve analysis.
For each test case, we regenerate a nucleus lightcurve by changing only one model parameter at a time.The parameter variant tests are described below, and the final lightcurve analysis results for all cases can be found in Table 2.The starting rotation of the model is kept as the same arbitrary time T 0 as was used for the nominal case.In every test case, we found that the pole search returned two equally-likely unique solutions for the pole orientation, corresponding to the best solution and its mirror as described in Section 3.2.For brevity, in the following we report only the pole coordinates and shape parameters of the solution which corresponds to the 'correct' orientation of the rotation pole, and we discuss the implications of the pole solutions further in Section 6.1.3.

Activity distance
To examine how cometary activity impacts the nucleus lightcurves, we implemented a simple threshold to dictate the minimum r h value at which lightcurve points could be included in the final dataset.For the nominal case, this threshold was set at 3 au.We tested two additional cutoffs: one at 4.3 au, the heliocentric distance of 67P below which which its activity could be detected by ground-based facilities (Snodgrass et al. 2013); and one at 0 au i.e. the comet was inactive for the duration of its orbit.The latter case could be representative of a devolatised JFC at the end of its dynamical lifetime (Jewitt 2004).
Generating synthetic lightcurves for r h > 0 au and r h > 4.3 au as outlined in Section 2.2, the number of data points remaining in the final lightcurves was 423 and 221 respectively.Including points across the entire orbit meant a maximum phase angle α max = 26 • was covered by the lightcurves, and the initial phase function fit to these scattered lightcurves was a close match to the input value of β L .Discarding all points below 4.3 au reduced the maximum observed phase angle to α max = 12.6 • , yielding a slope value of β scat = (0.054 ± 0.002) mag/deg, 7.8 per cent shallower than the input value of β L .
For both cases, the period search component of the convex lightcurve inversion identified P i as the period corresponding to the global χ 2 minimum.The pole searches successfully identified coordinates for the pole positions that were consistent with the simulation inputs.For the case r h > 4.3 au, it was found that the uncertainty in the pole coordinates was larger than for the nominal case -examination of the distribution of χ 2 values over the entire ecliptic sphere showed a large range of coordinates had a χ 2 value within 1.1× the minimum.This is likely a result of the reduced number of data points in the lightcurve covering a narrower range of nucleus aspects -in the dataset, both the nominal case and r h > 0 au cover an aspect angle range of ∼ 46 • , while imposing an r h > 4.3 au cutoff covers only 30 • in aspect.
The axis ratios of both final shape models are similar to those found for the nominal case, implying that the inversion procedure robustly identifies the optimal shape when lightcurve points are both added and removed.This indicates that a sufficient number of observing geometries are covered by the lightcurves in the r h > 4.3 au test case, that the final shape and spin solution is not noticeably improved with the addition of lightcurves closer to the Sun.The rotationally-corrected the phase functions β corr for both cases deviated from  β L by only 1 per cent, and the uncertainty in the phase function fit was reduced.

Size
We next considered the effect of nucleus size on our ability to extract nucleus properties from the LSST photometry.The effective nucleus radius can be expressed in terms of its geometric albedo p r and mean absolute magnitude Hr by Equation 2: The constant m ⊙ refers to the brightness of the Sun in the filter of measurement, and k = 1.496 × 10 8 km per au.For an effective radius of approximately 1 km, we use a value of Hr = 17.1 mag, and for 10 km Hr = 12.2 mag.These radius values were chosen to order-ofmagnitude approximate the extremes of the size range of the currently-known JFCs (see Table 1A of Knight et al. 2023).
The lower limit of 1 km was imposed by our method of lightcurve generation -we initially trialled 0.5 km as the small size case, but found that this resulted in the disposal of 97 per cent of the initial 613 timestamps due to the low signal-to-noise of the lightcurve points at large r h .This left fewer observed points than model parameters to fit in the lightcurve inversion process, and so we opted to scale the small radius up to 1 km.The final synthetic lightcurves consisted of 91 and 515 lightcurve points for the r n ∼ 1 km and r n ∼ 10 km models respectively.
For both cases, the initial phase function fit to the scattered lightcurve points was 3-4 per cent shallower than the initial β L value.The period and pole searches correctly identified the input period and pole orientation for both size models -this is particularly noteworthy for the small value of r n with only 91 data points total.The shape models were very similar.In both cases, using the convex model to correct the phase curve for rotational scatter allowed the input phase function β L to be precisely recovered, within the uncertainties on the slope.

Rotation period
Varying the rotation period allows us to examine the effects of spin rate on the quality of the recovered nucleus properties.To sample a representative range of JFC rotation periods, we generated lightcurves using nucleus models rotating with the minimum (P i = 2.8 h), maximum (P i = 41.3 h) and median (P i = 11.1 h) rotation periods for all known JFCs (from the compilation in Knight et al. 2023, Table 1A).We also tested a nucleus with a spin rate identical to Earth's sidereal period (P i = 23.9344696h) as an extreme case of sampling only a fraction of the object's lightcurve profile.
As was the case for the previous tests, the lightcurve inversion successfully identified the input period value as the most likely rotation period for each trial lightcurve, and the pole searches returned the input pole orienta- tion.This was true also for challenging case of P i ∼ 24 h, although the distribution of χ 2 values for all pole coordinates tested (Figure A1) suggests four regions of similar probability as a result of the ambiguous shape of the lightcurve.We have accounted for these additional possibilities in the uncertainties quoted for the pole coordinates in Table 2.The initial phase function estimate for the P = 41.3hlightcurves deviated significantly from β L by 5.8 per cent and the pole search revealed a relatively large uncertainty in polar longitude.The latter finding can be attributed to the sparse observations sampling comparatively less of the rotational phase of the slow rotator, leading to increased uncertainty in the estimate of the body meridian.

Shape
The choice of 67P as the typified JFC nucleus with which to simulate lightcurves allowed us to probe whether the non-convex shape model still allowed us to constrain physical parameters such as period and pole.The meagre sample of five JFCs with resolved nucleus imaging (9P, 19P, 67P, 81P and 103P) exhibit a wide diversity of shapes.As such, 67P was a natural choice of shape model as a bilobed object that is not particularly elongated.
To examine the impact of nucleus shapes on the lightcurve analysis outcomes, we trialled three additional shape models: JFCs 9P/Tempel 1 and 103P/Hartley 2; and trans-Neptunian object (486958) Arrokoth.Comets 9P and 103P represent the two extremes of the known JFC shapes, with 9P being irregularly-shaped but primarily rounded (Farnham & Thomas 2013a), while 103P is extremely elongated and dogbone-like in shape (Farnham & Thomas 2013b).The shape model of Arrokoth is based on images returned by the NASA New Horizons mission (Stern et al. 2019), which revealed a bilobed object with an extremely large a/c axis ratio.To date there have been no in situ detections of a JFC with such a shape.For 9P and Arrokoth, the initial phase function fit to the scattered lightcurve points already matched β L well within the slope uncertainties.In both cases, P i was successfully identified as the best period and the input pole orientation was recovered by the lightcurve inversion, however the uncertainties associated with the pole longitude are sizeable.This is likely a result of the small range of lightcurve amplitudes covered: for 9P, the lightcurve ∆m ranged from 0.12−0.24,and for Arrokoth 0.05 − 0.29.The convex model for 9P was slightly less elongated than the true shape (axis ratios given in Table 1) but is visually similar to the true shape (see Figure A2).In contrast, the convex model for Arrokoth was an extremely poor match to the input shape model.The lightcurve inversion favoured a close to spherical shape to fit the input lightcurves (see Figure A2).This is likely due to the extremely low ∆m values of the lightcurves.The phase function fits following correction by the con-vex model were comparable to the initial scattered fit within their uncertainties.
The initial phase function fit to the 103P lightcurves resulted in a shallow β scat = 0.0542 ± 0.0018 mag/ • , more than 8 per cent shallower than β L .The lightcurve inversion successfully determined the period and constrained the pole orientation to within a very small region of uncertainty.For the shape optimisation, we found it was necessary to iterate 400 times before the χ 2 fits to the lightcurves began to converge.The inversion identified an elongated shape for 103P (Figure A2), though the axis ratios were substantially less elongated than the true shape by nature of the convex model artificially elongating the rotation axis.The difference between β corr and β L was reduced to 2.5 per cent following rotational correction.

Pole
Spin axis orientation is known for only a small number of JFCs and is typically constrained by tracking evolving coma morphology.For low-activity comets, pole orientation is coupled with shape information, and thus simultaneous shape and pole optimisation by lightcurve inversion, such as demonstrated here, is one of the few ways of estimating this parameter for such objects.
For the above tests we used the real pole orientation of 67P (Preusker et al. 2015) which corresponds to an obliquity 3 I ∼ 52 • .For a given orbit, obliquity has a large impact on the measured lightcurve, as the region of the nucleus surface reflecting sunlight towards the observer changes as the polar aspect varies across the orbit.To test this effect, we synthesised lightcurves for nuclei in stable rotation states with obliquities of I ∼ 0 • , I ∼ 180 • and I ∼ 90 • , i.e. the polar direction aligned with the orbital normal vector in both the positive and negative directions, and pole perpendicular to the orbit normal vector.The ecliptic longitude and latitude values used to generate these lightcurves for 67P's orbit are given in Table 1.
For all three pole orientations, the β scat fit was notably shallower than β L .The best period was again consistent with P i in every case.For each test the pole search was able to converge to two unique solutions.At the maximum obliquity angle of I = 90 • , the pole orientation is tightly constrained and the convex shape properties are comparable with those from the previous test cases.For the I ∼ 0 • and I ∼ 180 • tests, the unique pole solutions were close to the ecliptic poles where longitude values are more closely spaced, shown in Figure 7.The 3 Angle subtended by the orbital momentum vector (perpendicular to orbit plane) and the vector in the direction of the spin axis.regions of uncertainty on the pole span 20 − 25 • degrees in ecliptic latitude, centred on the poles.While there is seemingly more ambiguity in the pole coordinates when considering the uncertainty regions, there is very little physical difference between different longitude values at ecliptic latitudes close to ±90 • .The axis ratios of the best-fit shape for the I = 0 • test deviate most from the nominal case than any other varied parameter.The same cannot be said for the op-posing obliquity at I = 180, but visual inspection of the best-fit model reveals a physically improbable shape with more than 100 • offset between the rotation pole and the shortest axis of inertia of the shape's equivalentvolume ellipsoid.These effects are likely due to the limited variation in object geometry caused by the low obliquity.Across these datasets the pole aspect varies by approximately 4 • , resulting in only small portions of their surface ever being illuminated to the observer.It is unlikely that any physically-meaningful shape information can therefore be extracted from lightcurves of low-inclination objects with low obliquity.The shape models were still effective for correcting the lightcurves for rotation, and in both cases brought the β corr fit much closer to the true value.

Sparse lightcurves
We demonstrated above that convex lightcurve inversion reproduces rotation periods and pole orientations extremely effectively when the lightcurves precisely describe the input model parameters -i.e. the lightcurve points include no photon noise.The small sizes of JFCs mean they are faint (20-24 mag) at the large heliocentric distances typically required to observe them with minimal activity, thus real nucleus brightness measurements will include significant photon noise.In the following we attempt to evaluate the effect of photon noise on our ability to extract nucleus parameters from LSST JFC observations using a Monte Carlo (MC) approach.Noisy versions of the nominal lightcurve are generated by replacing each lightcurve point with a randomly-generated magnitude clone taken from a normal distribution with mean and standard deviation respectively equal to the apparent magnitude and uncertainty value determined in Section 2.2.We repeat this process creating n = 100 lightcurves to sample a large range of the noisy parameter space.
We first tested the simple case of fitting a phase function and rotation period to the sparse lightcurves, as in Section 3.1, for n = 100 lightcurves.For each lightcurve we transformed the randomised apparent magnitudes back into reduced magnitudes H r (1, 1, α), and used weighted linear regression to fit the phase function β scat as a function of α.The resulting distribution is shown in the left panel of Figure 8.A normal distribution fit to this histogram yields a mean value βMC = 0.0566 ± 0.0012 mag/ • .Each randomised lightcurve was converted to absolute magnitudes by its respective β scat value, and we performed a Lomb-Scargle period search as described in Section 3.1.The distribution of best synodic period values is shown in the right panel of Figure 8.The period most often identified by the Lomb-Scargle algorithm (60 per cent of tests) was consistent with P i , and a period of ∼ 21 h was preferred for approximately 40 per cent of the lightcurves.We found that across all 100 lightcurves, the range of possible values for the best period identified by the algorithm was very small -within each histogram bin, the period values were similar to within 2 decimal places.This was true even for period searches with a large number of frequency samples per peak, and is likely a result of the sparse rotational sampling of the data over an extensive temporal baseline.
Applying convex inversion to the random lightcurves, we repeated the procedure outlined in Section 3.2 with 100 newly-generated lightcurve clones.We found that for every lightcurve, a period was identified that was close to P i , with the values ranging from 14.8976 − 14.9018 h.For 36 of the 100 lightcurves, a period multiple in the range 29.7953 − 29.8033 h was identified as equally likely, highlighting the challenge of correctly identifying the period from noisy data.We performed the pole search for each random lightcurve with a starting period value of P i , and allowed this value to be optimised along with the pole over 300 Levenberg-Marquardt iterations.The distribution of the top two unique pole solutions for all lightcurves is given in Figure 9. Considering only the models consistent with the input pole orientation (λ i = 78 • , β i = 41 • ) i.e. not the mirror solutions, the mean pole was located at (λ = 91 • , β = 40 • ) with interquartile ranges for the longitude and latitude values of 10.5 • and 4 • respectively.This implies that it is highly likely that after 10-years of the LSST, the rotation periods and pole orientations for a large number of JFCs will be known to a high degree of certainty.We optimised the convex shapes at the top two pole solutions for every lightcurves.Distributions of the model axes ratios are shown in Figure 10.Most models are similar in elongation to the noiseless models, and lie within 1.05 < a/b < 1.25.
For every random lightcurve, we selected the convex model corresponding to the pole direction closest to the true orientation, to correct for rotation and fit the phase function.The β scat and β corr distributions are shown in Figure 11.We found that the phase function values shifted positively to be centred on β L following the correction, with a mean and standard deviation of βcorr,s = 0.0591 ± 0.0015 mag/ • .While the correction improves the likelihood of correctly estimating the phase function, the shape model fit to the noisy data points resulted in an increased spread in the final distribution of β corr values.

Addition of dense-in-time lightcurves
As a final test, we repeated the MC convex inversion with the addition of dense-in-time lightcurves.Lightcurves covering large segments of the rotational phase help to constrain the rotational parameters at that geometry, and the inclusion of a single dense lightcurve has been demonstrated to measurably improve inversion outcomes for sparse asteroid photometry (e.g.Santana-Ros et al. 2015).
We created synthetic lightcurves across three dark nights in April 2029 when the comet would be observable with a 2-metre class telescope in the northern hemisphere for around half of the night.In these 'observa- tions', the target is at a heliocentric distance of 3.6 au, spanning phase angles 4.6−5.2• with a pole aspect ∼80 • .We estimated the magnitude uncertainties as 0.04 mag, and created 100 random instances of these lightcurves to be combined with the randomised LSST points.As the uncertainties in the dense lightcurves were approximately half as large as the average uncertainty in the sparse lightcurve, we weighted them to contribute twice as much as the sparse lightcurve to the fit at each stage of the optimisation process.Across all 100 of the combined dense and sparse lightcurves, the period search identified values consistent with P i as the best-fit rotation period solutions.Unlike in the case of the solely sparse lightcurves, no alias periods were returned as equally-likely, indicating that the inclusion of dense lightcurves substantially reduced the ambiguity in the period estimation from the noisy sparse datasets.The distribution of pole orientations resulting from the pole searches is shown in Figure 12.Considering the points closest to the true pole orientation, the mean ecliptic coordinates are λ = 84 • , β = 41 • .The distribution of possible pole orientations are demonstrably more tightly constrained when the lightcurve inversion has access to dense-in-time photometry: compared with the sparse-only lightcurves the inter-quartile range (IQR) of the model polar latitude is unchanged at 4 • , while the longitude IQR is reduced to 6 • and there is less scatter across the entire set of models.The rotationallycorrected phase function for the combined sparse and dense lightcurves is shown in Figure 11, with a mean and standard deviation βcorr,sd = 0.0591 ± 0.0012 mag/deg.The distribution mean is unchanged from the phase function fits to the sparse-only lightcurves, but there is a reduction in the distribution standard deviation.

Predictions for the LSST
The broad conclusion of this work is that accurate characterisations of nucleus properties will be possible from LSST photometry.The range of nucleus variants for which the model parameters were successfully reproduced suggests that these characterisations will enhance our understanding of JFCs on an unprecedented scale.Throughout this analysis we have made several simplifying assumptions in order to generate the synthetic lightcurves.Below we contextualise our findings outside of these assumptions, and predict what can realistically be expected from the output of LSST.

When will reliable results be attainable?
This analysis was based on lightcurves comprised of data points acquired across the complete operational period of LSST.This prompts the question of how much data is required before we can accurately measure the different JFC rotational properties, and at which point of the survey this amount of data will be available.The answer to this question clearly depends upon a host of factors including the survey start date, observational cadence, and orbital properties on an object-by-object basis.While we cannot succinctly answer this within the scope of this analysis, we tested our ability to determine the input parameters from the lightcurves obtained for our nominal test case after one and three years of survey commencement.This emulated the potential data availability after the first and third annual data releases.
Comet 67P has an orbital period of 6.4 years and its perihelion occurs approximately halfway through the LSST dataset in 2028.The first three years therefore cover the comet around aphelion and inbound to ∼ 4.6 au.The Year 1 dataset was comprised of 51 lightcurve points spanning ecliptic longitudes of 224 • to 231 • and latitudes −1 • to −0.5 • .The lightcurve amplitudes (∆m) ranged from 0.36 to 0.39.For the Year 3 dataset there were 120 lightcurve points spanning 224 • to 267 • in ecliptic longitude and −3 • to −0.5 • in latitude.Lightcurve amplitudes were in the range 0.30 to 0.39.
For both cases with noiseless lightcurves, the mostlikely rotation period was found to be consistent with P i .For the Year 1 test the pole solution did not converge to a unique solution, but for Year 3 two unique poles were identified at λ 1 = 280 −15 .The latter solution is close to the input pole orientation λ i = 78 • , β i = 41 •  We also ran an MC period and pole search on the Year 3 lightcurve.Across 100 random lightcurves, there were an average of 2.25 equally-likely period solutions with a maximum of 10.Every lightcurve found a period value within 0.005 h of P i as a possible solution.For 70 per cent of the lightcurves there were one or more equally-likely rotation periods in the range 5.1 − 39.5 h.Analysing the pole solutions as in Section 5, using an initial period of P i and disregarding the mirror solutions, the mean pole orientation was λ = 104 • , β = 26 • with an IQR of 15 • for longitude and 24 • for latitude respectively.
While there is significantly more uncertainty in the Year 3 dataset for period and pole compared with the results for the 'complete' dataset described in Section 5, these findings are encouraging and point to the possibility of measuring the nucleus rotational parameters within a single orbit.Our method therefore indicates that, assuming observations are obtained at large heliocentric distance, it is possible for the lightcurve inversion to identify the correct period and approximate pole orientation from three years of LSST observations (or > 100 detections with sufficient SNR).

Rotation period
For every nucleus parameter variant tested in this work, the input model rotation period could be accurately recovered solely from the sparse lightcurve points.Incorporating photometric noise led to difficulty distinguishing between the true period and period multiples, Figure 12.Same MC analysis as in Figure 9 but with the addition of a dense lightcurve to the sparse survey data points.
which was alleviated with the addition of dense-in-time lightcurves.For real nucleus photometry with no prior knowledge of the object's rotation period, it may be necessary to perform lightcurve inversion for all possible periods and rule out values by visual inspection.Nevertheless this is extremely promising for augmenting the number of JFCs with well-constrained rotation periods, without the need for extensive targeted monitoring campaigns.
The simplifying assumption made in this analysis was that the model nucleus was in a stable spin state for the survey duration, rotating with a constant period and pole orientation.It is well established that JFCs can experience period changes on measurable timescales (e.g.Knight et al. 2023, their Table 2 and references therein).Such spin changes are likely to occur close to perihelion where sublimative activity is maximised, and empirical models have indicated that the spin rates of smaller, more active nuclei are more strongly impacted than larger nuclei (Samarasinha & Mueller 2013;Jewitt 2021).
In the event that significant spin changes take place, employing sparse long-time baseline observations such as those simulated here may aid in constraining the extent of the variation.We have shown that changing nucleus geometry can be satisfactorily accounted for by lightcurve inversion, and that the rotation period can be extracted from a single orbit (Section 6.1.1).Therefore, if a rotation period is measured by lightcurve points obtained over a time frame including a perihelion passage, and it produces a poor fit to lightcurve points, this may indicate a change in spin rate.This approach was used to detect the Yarkovsky-O'Reefe-Radzievskii-Paddack (YORP) effect in near-Earth asteroids from dense-intime photometry (Kaasalainen et al. 2007).YORP induces linear-in-time variation in the spin states of these objects in contrast with the step function around perihelion that is exhibited by JFCs.
Given that the magnitude of spin change is believed to scale inversely with nucleus size, we expect that if spin changes do take place throughout the operational period of LSST, the nuclei for which the effect is strongest will be too faint for lightcurve inversion by LSST (see discussion below).Most JFCs will make at most one perihelion passage over the ten-year survey duration.Based on our findings we anticipate that accurate rotation periods will be accessible for many km-sized JFCs, provided they have been observed at large heliocentric distance for a significant portion of the survey.If the observations encompass a perihelion passage, we do not rule out the possibility of period changes being detectable by lightcurve inversion.

Shape and pole
The strength of lightcurve inversion is the simultaneous optimisation of shape and pole in addition to rotation period.We found that in the test cases with a low obliquity angle, the pole orientation was correctly constrained to its respective ecliptic hemisphere.For every test case with an intermediate pole obliquity, the best-fit pole orientation returned the input direction as well as the λ+180 • mirror solution.Where there was a large region of uncertainty around the best-fit pole coordinates, either the lightcurves exhibited low amplitude values or, in the case of the comet rotating with P i ∼ 24 h, there was fundamental uncertainty in the lightcurve profile.
These findings suggest that for most low-inclination JFCs, unless they have obliquities close to 0 • or 180 • , lightcurve inversion of LSST photometry will result in dual equally-likely pole orientations.However, this ambiguity is only in polar longitude as the pole latitude is similar (within ∼ 10 • ) for both the true and mirror pole orientations.As discussed by Ďurech et al. (2009), for low-inclination orbits the obliquity can be estimated directly from the polar latitude.
A larger sample of JFCs with known obliquities will provide insights into seasonal activity variations from ground-based observations.Constraining the rotation poles for a more representative sample of the JFC population may also allow further development of models of comet evolution.Samarasinha (1997) and Neishtadt et al. (2002) predicted that the rotation poles of short period comets would tend towards a point of rotational stability oriented along the direction of maximum outgassing i.e. towards the perihelion or aphelion directions.
Regarding nucleus shapes, we had variable success in reproducing the overall shape profiles of the input shape models by convex inversion.The best-fit convex models for the lightcurves generated with a 67P shape model had comparable axis ratios to the literature convex models, with the exception of the low obliquity tests which varied in aspect by only ∼ 4 • .We would caution against relying on shape information from sparse data points for such objects with low orbital inclinations as a result.
The lightcurve inversion favoured more rounded shapes even when the input shape model was elongated in the cases of 103P and Arrokoth.For every case the convex model axis ratio was less elongated than the true a/b value.We attribute this to the convex regularisation having the effect of masking over potential concavities and artificially extending the shorter model axes.We expect that elongated shapes like 103P and Arrokoth are likely to require substantial dense-in-time lightcurve coverage to give the lightcurve inversion adequate constraints for a more representative shape model.Interpretation of how well convex shape models generated by sparse photometry represent non-convex comet nuclei is beyond the scope of this work.

Activity
With the inclusion of a dense lightcurve in Section 5 we weighted the lightcurves' contribution to the final shape and spin model according to their uncertainties, which stem from the apparent brightness of the data points.However, for real active objects with poorlyunderstood activity patterns, brighter objects tend to be closer to the Sun and more likely to be contaminated by activity.Treatment of uncertainties in brightness measurements of potentially active objects should always be carefully considered alongside an assessment of their activity levels.
The simplifying assumption typically applied to JFCs (and used for the majority of the tests in this work) is that they are inactive beyond 3 au from the Sun where temperatures are insufficient for water ice to sublimate.However, a significant number of JFCs have exhibited activity at heliocentric distances beyond this (e.g.Mazzotta Epifani et al. 2007, 2008;Kelley et al. 2013).We demonstrated that only incorporating data points acquired at r h > 4.3 au still produced accurate results in the noiseless test case, but the drivers of activity at large distances from the Sun are poorly understood and implementing a simple distance cutoff does not guarantee data points will be free from activity contamination, for example in the case of outbursts.Interpreting photometric observations of comet nuclei therefore always requires a careful assessment of the comets' activity levels.
There are various challenges in assessing the activity of faint sources, but it is an important consideration as undetected activity in photometry apertures may lead to erroneous interpretation of lightcurves.If unaccounted for, significant coma flux could suppress the rotational signature of the nucleus and yield anomalously bright lightcurve points, which may be mistaken for shape effects if incorporated into datasets for lightcurve inversion.This scenario will be difficult to recognise solely from sparse survey data points.A deep activity search combining frames acquired at similar geometries may reveal unresolved coma.It is likely that the frequent detections and signal-to-noise capabilities of LSST for many JFCs at large r h will significantly enhance our understanding of activity at such distances -for example, calibrated brightness measurements across large portions of the orbit may aid in the future identification of anomalous brightening.

Nucleus size limitations
Using our method of synthetic lightcurve generation, we found it was necessary to use a minimum effective radius of 1 km.Setting the low end of the nucleus size range to 0.5 km, which more closely reflects the size of the smallest observed JFC nuclei, left only ∼ 20 data points for lightcurve analysis.This meant fewer data points than model parameters to fit by the chosen method of lightcurve inversion, and thus we opted not to test it in this work.
If the true photometric uncertainties in the survey are larger than our simple SNR approximation, it is possible that only those JFCs with r n > 1 km will be detected by the survey at the large heliocentric distances typically necessary for nucleus characterisation.Using the Bauer et al. (2017) debiased JFC size distribution, and accounting crudely for the fraction of inactive JFCs as 75 per cent those passing aphelion during the LSST that do not exhibit activity at large heliocentric distance (Kelley et al. 2013), we estimate that this analysis will be possible for a few hundred JFC nuclei.For objects with r n ∼ 0.5 km, 20 lightcurve points at a range of geometries may be sufficient to deduce a first approximation to the rotation period and pole orientation if a model with fewer fit parameters is employed.These could be refined with targeted follow-up observations by a designated observing program to characterise small nuclei discovered by LSST.

Phase functions
Varying the input parameters for the nucleus model in Section 4, it was found that the phase function fit β scat to the reduced magnitudes was often not consistent with the input phase function value β L .The average difference (∆β scat ) between the phase function fit to the scattered lightcurves and β L across all varied parameters was 4 per cent.This deviation can be attributed simply to rotational scatter: fitting a linear function to data with y-axis scatter yields greater uncertainty in the slope (and intercept) values.Even with an MC approach to sample the distribution of phase function fits to the noisy lightcurves as in Figure 8, recovering the true slope value was unlikely.
In most cases, the difference between β L and the measured β scat was insignificant, but for several cases the deviation was of order 10 per cent and the phase function measured by this method may lead to erroneous interpretations of the nucleus surface properties.The large deviation in the 103P slope value is likely a result of the high lightcurve amplitudes (∆m) which, when sampled by the survey simulator, resulted in a larger spread of possible magnitude values for the data points.A similar argument can be applied to the phase function fit to the I = 0 • lightcurves: such an obliquity for an object with a low-inclination orbit would result in a mostly equatorial view of the rotating nucleus.For 67P, this implies viewing the bilobed nucleus face-on, resulting in long periods at the flattened lightcurve maxima and hence more scatter in y-axis and greater slope uncertainty.The same should also be true for the opposite pole orientation at I = 180 • , which has a comparably-sized slope uncertainty.The slope value for the lightcurves with r h > 4.3 au also deviated substantially from β L .Notably, this was the lightcurve that covered the lowest range of phase angles (0.2 − 12.6 • ) which likely impacted the straight line fit.
This highlights the challenge of fitting phase functions to photometric time series without lightcurve correction.When the lightcurve scatter is reduced -i.e. when we correctly account for the rotational brightness modulation of each lightcurve point and remove its corresponding amplitude from the reduced lightcurve -the slope value is consistently close to the true value.This was demonstrated in Section 4, in which the mean difference between β L and the corrected β corr across all noiseless lightcurves was reduced to 0.9 per cent.As shown in Figure 11, this is true on average even when the rotational variation is partially obscured by photon noise.
To improve the reliability of literature phase function values, where possible these values should be corrected for rotation.For dense lightcurves, this is typically in the form of rotational averaging, or correction using a convex model.It was recently demonstrated that removing rotational modulation had a significant impact on the phase function fit to JFC 162P/Siding Spring, which increased by 9 per cent from 0.0468 ± 0.0001 mag/ • to 0.051 ± 0.002 mag/ • (Donaldson et al. 2023).For populations with large amounts of available sparse photometry such as asteroids, methods are still being developed to fit reliable phase curves (e.g.Martikainen et al. 2021;Wilawer et al. 2022).As identified by Dobson et al. (2023) and Robinson et al. (2024), there are non-negligible differences between phase function fits to sparse lightcurves based on whether any attempt to account for the rotation is included.Currently there are only around 20 JFCs for which it has been possible to estimate the nucleus phase function (Knight et al. 2023, and references therein).Our analysis demonstrates that sparse LSST photometry alone is sufficient to accurately measure the rotationally-corrected phase function.This should be possible for many objects, leading to a vast improvement in our population statistics and a meaningful way to compare JFCs with other populations, to trace potential origins, evolution and dynamical links (e.g.Kokotanekova et al. 2018).

SUMMARY
We generated synthetic JFC nucleus lightcurves for objects on a 67P orbit observed with an LSST-like cadence.Each synthetic lightcurve tested a different variant of nucleus characteristics, including rotation period, size, shape, obliquity and activity distance.In every case it was assumed that the comet was inactive during the observations, and that the nucleus was in a stable configuration with a constant rotation period.We processed the sparse-in-time lightcurve points using convex lightcurve inversion to determine whether the input model properties could be reproduced.The key findings and conclusions are summarised as follows: 1. LSST will provide unprecedented coverage of the orbits of observable JFCs, detecting km-sized nuclei at large heliocentric distances with sufficient quality for lightcurve analysis.Due to signal-tonoise limitations, we anticipate that such analysis will be possible for a few hundred objects with effective radii ∼ 1 km and larger.
2. For every lightcurve tested, the best sidereal period identified by lightcurve inversion was consistent with the model input period.The input period was identifiable even with the addition of photon noise.This is extremely promising for accurately determining rotation periods for many JFCs following the ten-year operational period of LSST.
3. The range of observing geometries covered by LSST observations for a typical low-inclination JFC is sufficient to constrain the nucleus pole orientation to within two possible directions.This will allow us to measure pole obliquities for more JFCs than is currently feasible.
4. The properties of the convex shapes determined from lightcurve inversion of sparse survey data are comparable with those found by dense-in-time targeted monitoring campaigns.
5. It is possible to identify the correct rotation period and approximate pole orientations from lightcurve points spanning a three-year period.
6. Throughout the survey duration, many JFCs will be observed at a larger range of phase angles than has been realistically attainable using groundbased facilities to date.In all test cases it was possible to rotationally-correct the lightcurves using the best-fit convex model, and recover a phase function value close to the model input value.We suggest that this becomes the standard practise for quoting JFC phase functions where possible to reduce uncertainty in the slope estimates from inherent lightcurve scatter.This will allow for more valid comparisons to be made between JFCs and other populations.
7. Incorporating dense-in-time nucleus lightcurves with the sparse survey photometry measurably improves the accuracy of the period, pole and phase function estimates.
8. Cometary activity is unpredictable and JFCs can exhibit activity at large heliocentric distances.Lightcurve points should carefully checked for signatures of sublimation-induced brightening, ideally automatically, prior to inclusion in any rotational lightcurve analysis.

Figure 1 .
Figure 1.Left: filter (LSST ugrizy) of every pointing covered by the simulated nominal dataset.Lightcurve points are discarded from the final dataset first by SNR, then heliocentric distance r h .Right: final lightcurve for the nominal case of this study, with 325 individual points remaining within imposed 5σ and r h limits.Marker colour represents time from first observation.Points were corrected for solar phase by linear phase function slope value βL, and phased to period Pi = 14.8978 h.The misalignment of lightcurve profiles is due to synodic effects (i.e. the object viewing geometry varying with respect to the observer over the 10-year lightcurve baseline).

Figure 2 .
Figure 2. Reduced magnitude as a function of phase angle α for the nominal lightcurve.Solid line gives best fit to lightcurve points from linear regression with a slope value β = (0.0565 ± 0.0014) mag/ • .Dashed line is phase function used as input (slope βL) to generate lightcurves.Scatter in points is due solely to lightcurve amplitude, no photometric noise incorporated at this stage and error bars are not used to weight the phase function fit.

Figure 3 .
Figure 3. Lomb-Scargle periodogram for single-peaked rotation period search on nominal lightcurve.The highest peak corresponds to a most-likely synodic rotation period P = 14.898088 h.

Figure 4 .
Figure 4. Left: periodogram for period search using the nominal lightcurve.The χ 2 minimum is consistent with Pi and the next strongest peaks occur at Pi multiples and sub-multiples.Other peaks are due to sampling aliases.Right: distribution of rotation pole solutions.Darker regions indicate a lower value of χ 2 and therefore a better fit to the input lightcurves.Yellow dot depicts the ecliptic coordinates corresponding to the pole orientation with best fit to the lightcurves, yellow circles show regions within which χ 2 values are within 10% of this minimum value.ance, likely due to the small number of lightcurve data points (typically 1-2) to fit to at each geometry.The two final shape models are shown in Figure 5.The axis ratios of Model 1 are a/b = 1.15 and b/c = 1.15, and for Model 2 a/b = 1.13 and b/c = 1.16.The axis ratios of the convex models are substantially less elongated than those of the input 67P model (see Table1).This is to be expected for an object as strongly non-convex as 67P -the lightcurve inversion procedure cannot recreate surface concavities and instead tends to mask over them, which has the effect of artificially lengthening the rotation axis of the convex model in this case.We test the impact of different nucleus shapes later in Section 4.4.It is therefore more illustrative to compare these shapes to the previously published convex models of 67P obtained byLowry et al. (2012) andMottola et al. (2014).The axial ratios based on the inertia tensors of these models are a/b = 1.19, b/c = 1.25 and a/b = 1.14, b/c = 1.22 respectively.The sparse-intime model axis ratios derived in this work are remarkably similar to the previously published models.This is extremely promising for characterising a greater number of JFCs with LSST photometry: the diversity of observing geometries covered over the survey's 10-year period of operation is capable of providing shape and pole information that is equivalent to years of targeted monitoring.While the convex models can provide valuable insights into the object's global shape properties, these insights are somewhat limited when the original shape is intrinsi-

Figure 5 .
Figure 5.A comparison of input 67P shape model (top) with the final convex models corresponding to pole Solution 1 (centre) and the mirror at Solution 2 (bottom).Three orthogonal views are given for each model, with the viewing angle oriented along the axis indicated in bold above the model.The rotation axis is aligned with the Z-axis.

Figure 6 .
Figure 6.Nucleus phase function for the nominal lightcurves.Blue-coloured points show the synthetic magnitudes corrected for the effects of rotation by the Solution 1 shape model, coloured according to heliocentric distance r h .The translucent grey points depict the lightcurve points before the rotational correction.The solid line shows the best linear fit to the corrected points (slope value βcorr = 0.0588 ± 0.0002 mag/deg, intercept H(1, 1, 0) = 15.634 ± 0.002), and the dashed line the fit to the original, scattered points (slope value βscat = 0.0565 ± 0.0014 mag/deg, intercept H(1, 1, 0) = 15.665 ± 0.012).

Figure 8 .
Figure 8. Histograms from the Monte Carlo (a) phase function and (b) period analysis for the nominal case.The true phase function slope βL is given by the dashed vertical line in the left panel.The period values shown are twice the most-likely period value identified by the Lomb-Scargle analysis.

Figure 9 .
Figure 9. Pole orientations corresponding to the top two unique pole solutions from lightcurve inversion of 100 sparsein-time lightcurves.Colour corresponds to density of points where lighter regions have a larger number of points.Dashed lines mark the input pole coordinates (λi, βi).

Figure 10 .
Figure10.Distribution of principal axis ratios a/b and b/c for the top two shape models from 100 MC randomised lightcurve trials.Colour scale depicts pole longitude: darker coloured points correspond to pole orientations closest to the simulation input pole, lighter colours are shapes at the mirror pole.Also shown are axis ratios for the noiseless models developed in Section 3, and previously published models byLowry et al. (2012, L12)  andMottola et al. (2014, M14).

Figure 11 .
Figure 11.Distribution of phase function fits to 100 MC lightcurves for three different cases.Dashed vertical line shows the input phase function value βL.Orange line (βscat) demonstrates weighted linear regression on reduced magnitudes with lightcurve amplitude scatter.The mean of this distribution is 0.0567 ± 0.0005 mag/ • .Red line (βcorr,s) shows distribution of phase functions for lightcurves corrected for rotation by a convex model constructed from purely sparse lightcurves with mean 0.0591 ± 0.0015 mag/ • .Purple line (labelled β corr,sd ) results from the phase function fit using a convex model constructed from combination of sparse and dense lightcurves, with mean 0.0591 ± 0.0012 mag/ • .

Table 1 .
Values of physical properties used to generate synthetic nucleus lightcurves.Note-Parameter values shown in boldtype are those which were used to generate the lightcurve for the nominal case.* Parameters for which the values were kept fixed across all artificial lightcurves.* * Alternative model parameters which are varied for testing in Section 4. † Heliocentric distance r h below which it is assumed the comet is active.For case 0 au, the comet is inactive throughout the entire orbit i.e. points at all values of r h are included in the lightcurve dataset.

Table 2 .
Resulting values for parameter fits when model inputs were varied.Range of model lightcurve amplitudes i.e. ∆mmax − ∆mmin. 2 Absolute percentage difference between measured phase function and input phase function β L .The large uncertainties in pole orientation for cases Pi = 23.93.. h, I = 0 • and I = 180 • correspond to error ellipses with equivalent radii 24 • , 23 • and 25 • respectively.