The impact crater at the origin of the Julia family detected with VLT/SPHERE?

Context. The vast majority of the geophysical and geological constraints (e.g., internal structure, cratering history) for main-belt asteroids have so far been obtained via dedicated interplanetary missions (e.g., ESA Rosetta, NASA Dawn). The high angular resolution of SPHERE/ZIMPOL, the new-generation visible adaptive-optics camera at ESO VLT, implies that these science objectives can now be investigated from the ground for a large fraction of D ≥ 100 km main-belt asteroids. The sharp images acquired by this instrument can be used to accurately constrain the shape and thus volume of these bodies (hence density when combined with mass estimates) and to characterize the distribution and topography of D ≥ 30 km craters across their surfaces. Aims. Here, via several complementary approaches, we evaluated the recently proposed hypothesis that the S-type asteroid (89) Julia is the parent body of a small compact asteroid family that formed via a cratering collisional event. Methods. We observed (89) Julia with VLT/SPHERE/ZIMPOL throughout its rotation, derived its 3D shape, and performed a reconnaissance and characterization of the largest craters. We also performed numerical simulations to ﬁrst conﬁrm the existence of the Julia family and to determine its age and the size of the impact crater at its origin. Finally, we utilized the images/3D shape in an attempt to identify the origin location of the small collisional family. Results. On the one hand, our VLT/SPHERE observations reveal the presence of a large crater ( D ∼ 75 km) in Julia’s southern hemisphere. On the other hand, our numerical simulations suggest that (89) Julia was impacted 30–120 Myrs ago by a D ∼ 8 km asteroid, thereby creating a D ≥ 60 km impact crater at the surface of Julia. Given the small size of the impactor, the obliquity of Julia and the particular orientation of the family in the (a,i) space, the imaged impact crater is likely to be the origin of the family. Conclusions. New doors into ground-based asteroid exploration, namely, geophysics and geology, are being opened thanks to the unique capabilities of VLT/SPHERE. Also, the present work may represent the beginning of a new era of asteroid-family studies. In the ﬁelds of geophysics, geology, and asteroid family studies, the future will only get brighter with the forthcoming arrival of 30–40 m class telescopes like ELT, TMT, and GMT.


Introduction
Our understanding of the surface composition of asteroids and its distribution across the asteroid belt has improved enormously over the last two decades (see Burbine 2016 andBeck 2017 for recent reviews); however, this is not the case for the surface topography and shape of asteroids. This stems from the fact that the visible and near-infrared disk-integrated spectral properties of asteroids (and to a lesser extent the mid-infrared properties) could be accurately measured via Earth-based (either ground-based or space-based) telescopic observations for a large number of bodies, whereas disk-resolved observations of asteroids could only be obtained with sufficient spatial resolution for a handful of bodies via (i) rare in situ space missions (Galileo, NEAR, Hayabusa, Rosetta, Dawn; Thomas et al. 1994Thomas et al. , 1999Thomas et al. , 2002Veverka et al. 1994Veverka et al. , 1997Veverka et al. , 1999Sullivan et al. 1996;Demura et al. 2006;Saito et al. 2006;Keller et al. 2010;Sierks et al. 2011;Jorda et al. 2012;Jaumann et al. 2012;Marchi et al. 2012;Schenk et al. 2012;Russell et al. 2013;Buczkowski et al. 2016;Hiesinger et al. 2016;Russell et al. 2016); (ii) high angular resolution imaging with the HST or large gound-based telescopes for a few of the largest asteroids, (1) Ceres and (4) Vesta in particular (e.g., Binzel et al. 1997; Thomas et al. 1997aThomas et al. ,b, 2005Zellner et al. 1997;Carry et al. 2008Carry et al. , 2010bMerline et al. 2013;Drummond et al. 2014;Marsset et al. 2017); or (iii) radar echos (e.g., Ostro et al. 2000;Shepard et al. 2017Shepard et al. , 2018. The high angular resolution of SPHERE/ZIMPOL (∼20 mas at lambda = 600 nm; Schmid et al. 2017), the new-generation visible adaptive-optics at ESO/VLT, implies that such a science objective can now be investigated from the ground for a large fraction of D ≥ 100 km main-belt asteroids because most of these bodies display an angular diameter around opposition larger than 100 mas. The sharp images acquired with this instrument can be used to accurately constrain the shape and thus the volume of these bodies, and hence their density (when combining the volume with existing mass estimates) and to characterize the distribution and topography of D ≥ 30 km craters across their surfaces.
Density is the physical property that constrains best the internal compositional structure of large (D ≥ 100 km) asteroids. Most D ≤ 100 km asteroids are seen as collisionally evolved objects (Morbidelli et al. 2009) whose internal structure can be largely occupied by voids (called macroporosity, reaching up to ∼50% in some cases; Fujiwara et al. 2006;Carry 2012;Scheeres et al. 2015), thus limiting our capability to interpret meaningfully their bulk density in terms of composition(s). On the contrary, most large bodies (D ≥ 100 km) are seen as primordial remnants of the early solar system (Morbidelli et al. 2009); in other words, their internal structure has likely remained intact since their formation, and they can be considered the smallest protoplanets. For these objects, the macroporosity is expected to be minimal (≤10%) and their bulk density is therefore an excellent tracer of their initial bulk composition and thus of their formation location. This last characteristic determines the nature of their building blocks: rock in the case of the inner solar system and a mixture of ice and rock in the case of the outer solar system. The topography of the largest craters offers complementary information to the bulk density in the sense that it allows us to characterize the response of the outer shell (of the first 20-40 km below the surface) to large impacts, the latter response being directly related to the nature and composition of the outer shell.
To make substantial progress in our understanding of the shape, internal compositional structure (i.e., density), and surface topography of large main-belt asteroids, we are carrying out an imaging survey (via an ESO Large program entirely performed in service mode with seeing constraints ≤0.8 ; PI: P. Vernazza; ID: 199.C-0074) of a statistically significant fraction of all D ≥ 100 km main-belt asteroids (∼35 out of ∼200 asteroids; our survey covers the major compositional classes) at high angular-resolution with VLT/SPHERE throughout their rotation (typically six epochs per target).
Here, we present the first results of our survey. We focus on asteroid (89) Julia (hereafter simply Julia), a large (D ≥ 100 km; Tedesco et al. 2004) S-type main-belt asteroid that was discovered in 1866 by French astronomer Edouard Stephan. Whereas Julia had originally received significantly less attention than other large S-type asteroids such as (6) Hebe (see Marsset et al. 2017 and references therein) or (8) Flora (e.g., Florczak et al. 1998;Nesvorný et al. 2002;Vernazza et al. 2008Vernazza et al. , 2009Dykhuis et al. 2014;Vokrouhlický et al. 2017), it came under scrutiny in the early 2000s when it became a potential flyby target for the ESA Rosetta mission (Birlan et al. 2004). Ultimately, however, the main-belt asteroids (2867) Steins and (21) Lutetia were selected as flyby targets (Barucci et al. 2005). Later on, spectroscopic observations of Julia at visible and near-infrared wavelengths (0.4-2.5 µm) have shown that its surface composition is dominated by olivine and low calcium pyroxene and that LL chondrites are its closest meteoritic analog (Vernazza et al. 2014). Recently, it has been proposed that Julia is the parent body of a small compact asteroid family that consists of 33 known members (Nesvorný et al. 2015) and that formed via a cratering collisional event.
Here, we evaluate this hypothesis via several complementary approaches. We observed Julia with VLT/SPHERE/ZIMPOL throughout its rotation and derived its 3D shape, as well as the topography of the largest craters. We also performed numerical simulations to first confirm the existence of the Julia family and to determine its age and the size of the impact crater at its origin. Finally, we utilized the images/3D shape in an attempt to identify an origin location consistent with the proposed small collisional family.
Our paper is organized as follows. In Sect. 2, we present the ZIMPOL data and our survey methods in terms of data reduction and image deconvolution. In Sect. 3, we present the archive data that was used in addition to the ZIMPOL data for the shape reconstruction. It comprises Keck/NIRC2 disk-resolved data and optical disk-integrated data. The shape reconstruction is described in Sect. 4.1. We discuss Julia's physical properties in Sects. 4.2 and 4.3. In Sect. 5, we present the results of the dynamical simulations that were performed to constrain the age of the Julia family and the size of the impact crater. Finally, we conclude our work in Sect. 6.

Observations, data reduction, and deconvolution
Julia was observed with the SPHERE instrument (ESO/VLT; Beuzit et al. 2008) around its opposition at eight different epochs (see Table B.1 for a complete list of the observations). We used ZIMPOL (Thalmann et al. 2008) in narrowband imaging mode (N_R filter; filter central wavelength = 645.9 nm, width = 56.7 nm). Each observational sequence consisted of a series of five images, where each image corresponded to a series of detector integration times (DITs) of 10 s, during which Julia was used as a natural guide star for adaptive optics (AO) corrections. Observations were performed under good seeing conditions For the epochs with the highest contrast (typically epochs 2 and 3), we observe a variation in the intensity between the brightest and darkest areas on the order of 50%. We interpret these rather circular, low-intensity areas as craters (see Fig. 3).
(≤0.8 ) with an airmass usually below 1.6. After every asteroid observation, we observed a nearby star for deconvolution purposes to estimate the instrument point spread function (PSF). Finally, standard calibrations, which include detector flat-fields and darks, were acquired in the morning as part of the instrument calibration plan.
Data reduction of the asteroid and PSF observations were performed with the ESO pipeline Esorex and included the following steps: (1) pre-processing of the images, which consists in combining the two detector segments into a single trimmed image and splitting the images into odd and even sub-frames. In the imaging mode, all the information is contained in the odd subframe, the rest being masked; (2) background subtraction; (3) flat-fielding; (4) bad-pixel correction; and (5) stacking of the cubes of DITs. For a given epoch, the final product of the reduction procedure consists of five images. We note that for each epoch we also produced an average image to increase the signal-to-noise ratio and perform a companion search. Finally, the optimal angular resolution of each image was restored with MISTRAL, a myopic deconvolution algorithm optimized for images with sharp boundaries (Fusco et al. 2003), using the stellar PSF acquired on the same night as our asteroid data. The deconvolved images of asteroid Julia are shown in Figs. 1 and B.1.
The images acquired for a given epoch are very consistent with each other with the exception of one case (image 5 of epoch 4; the difference here may be due to a sudden change in seeing conditions; see Fig. B.1). The minimal difference between the images for a given epoch can be explained by the relatively slow rotation period of Julia (∼11.4 h) and the short time interval between the first and last image of a given epoch (∼14 min), implying that there is only a difference of ∼7 • in rotation phase between the first and last image of a given epoch. Overall, the fact that fine details are repeatedly visible in the different images of each series gives us confidence that they are genuine features and not artifacts.

Disk-resolved images
In addition to our VLT/SPHERE data, we used one Keck image of Julia (see Hanuš et al. 2017 and Table B.1 for additional information). Given the difference in spatial resolution between the Keck image and our VLT/SPHERE images (a factor of ∼2 in one spatial direction), we used the Keck image only as a consistency check (see Fig. B.2).

Shape models and optical lightcurves
The Database of Asteroid Models from Inversion Techniques (DAMIT 1 ;Ďurech et al. 2010) contains the most recent shape models of asteroid Julia and the photometric data used to produce those models: 22 lightcurves sampling three different apparitions, namely, those in 1968 (Vesely & Taylor 1985), 1972 (Schober & Lustig 1975), and 2009 (Hanuš et al. 2013). In the present work, we used these 22 lightcurves along with 15 new lightcurves recorded in 2017. Our dataset of disk-integrated optical photometry therefore contains 37 lightcurves from four different apparitions spanning a time interval of 49 yr. The aspect data and the references of the lightcurves are listed in Table B.2.

Stellar occultations
We also considered two stellar occultations recorded in 2005 and 2006 and already published in Hanuš et al. (2017). However, because of their rather low quality, these measurements did not bring significant constraints on the shape. As in the case of the single Keck image (see above), the occultations were mostly used as a consistency check (see Fig. B.2). A&A 618, A154 (2018) Fig. 2. Comparison between the deconvolved images of Julia (bottom panels) and the corresponding shape model projections (top panels). The AO epochs were ordered 1-8-2-5-6-7-3-4 and correspond to an increasing rotation phase. Julia's spin axis (red) is also shown.

Reconstruction of Julia's 3D shape with the ADAM algorithm
The All-Data Asteroid Modeling (ADAM) algorithm (Viikinkoski et al. 2015;Viikinkoski 2016) is a multi-data inversion technique that allows the 3D shape of an asteroid to be reconstructed by combining both disk-integrated and diskresolved observations. While the disk-resolved data contribute the most to fine-scale details, optical lightcurves are essential for shape reconstruction since they stabilize the optimization process and provide information on the regions of the asteroid that are not covered by disk-resolved observations. In ADAM, a shape model is represented as a polyhedron with triangular facets, and the vertex locations are given by a parametrization. Given initial estimates for the shape parameters and rotation state, the optimization process searches for parameters minimizing the objective function where χ 2 lc and χ 2 ao measure the model fit to lightcurves and adaptive optics images, respectively. The last term corresponds to regularization functions γ i (penalizing fine-scale concavities) and their weights λ i . These functions are described in detail in Viikinkoski et al. (2015).
In general, the reconstructed shape is not unique. The use of different shape supports (i.e., subdivision surfaces and octanoids; see Viikinkoski et al. 2015) and regularization functions allows peculiarities caused by parametric representations to be distinguished from model features supported by the data. In particular, real features should be visible in all the shape models with identical χ 2 fits.
We applied the ADAM algorithm to our dataset (40 VLT-SPHERE/ZIMPOL images from 8 different epochs, 37 optical lightcurves, 2 occultations) and produced various shape models of Julia by considering different combinations of (i) shape support (i.e., octanoids and subdivision), (ii) AO data types (deconvolved or non-deconvolved images) and occultation measurements (with or without), and (iii) shape model resolutions (see Fig. B.3). In all cases, we obtained qualitatively similar shape model solutions and volumes. The range of values for the physical parameters such as spin axis orientation, sidereal rotation period, volume-equivalent diameter, and dimensions along the major axes obtained for the various shape models allowed us to estimate their average and their uncertainty
( Table 1). The comparison between the eight AO epochs and the corresponding shape model projections is shown in Fig. 2; we show the shape model based on the deconvolved AO data, the octanoid shape support, and a high-resolution shape model.

Volume and bulk density
We combined the derived volume with current mass estimates to constrain the density of Julia. Because of the few mass estimates available and given their disparity (see Table 2), the average mass estimate for Julia has a significant error bar [(4.30 ± 3.59) × 10 18 kg] leading to a density estimate with a large uncertainty (ρ = 3.0 ± 2.6 g cm −3 ). It follows that no reliable conclusions can be reached regarding the density of Julia given that it is consistent with the density of all classes of meteorites, except perhaps iron meteorites if we assume reasonable (≤20%) macro-porosity values.

Reconnaissance and analysis of impact craters
We first performed a thorough visual inspection of the images to identify likely impact craters present at Julia's surface and used the 3D shape model to determine the corresponding location of each crater. We only considered the impact craters present in the front and side views and not the ones visible in the image contours as the latter could be confused with depressions due to Julia's irregular shape. We found one highly probable large impact basin clearly visible at rotation phases 0.18 and 0.49 both in the images and in the shape model. Two putative craters A and B are visible at rotational phase 0.67; we consider the craters   Parameter Value putative because we only detected them in one epoch (see Fig. 3). We named the large impact basin after the village Nonza where Saint Julia of Corsica was born ((89) Julia was named after her). We note that Nonza is not clearly visible at rotation phase 0.53, possibly because of degraded observing conditions or an unfavorable viewing angle.
Second, we fitted Julia's 3D shape model with an ellipsoid in order to produce an elevation map. The map reveals three major depressions, one of them being Nonza (located at zero longitude and −32 • latitude). We searched for the origin of the two large northern hemisphere depressions (located at longitudes Fig. 4. Elevation map (in km) of the shape model of Julia obtained by best-fitting the shape of Julia with an ellipsoid. The gray area (∼20% of Julia's surface) highlights the latitudes that were not covered by our ZIMPOL images. Given the fact that Julia's shape deviates significantly from ellipsoid, the amplitude of the relief is significantly enhanced with respect to the local topographic relief. Nonetheless, the main surface features are clearly visible (e.g., Nonza at −32 • in latitude). We defined the zero longitude with respect to the center of the Nonza crater. ∼70 • and ∼200 • ) by comparing the 3D model to the images. We found that they are most likely due to the overall shape of Julia (we can see them clearly in the contours of epoch 1) rather than to unambiguous impact features, although we cannot firmly rule out this origin given the limited spatial resolution (e.g., epoch 1) of some of our images.
Finally, we used Julia's 3D shape model to constrain the diameter, depth, and volume of Nonza (summarized in Table 3). It is the only impact crater whose size is compatible with the one at the origin of the Julia family that we could identify in our images (see Sect. 5.3). The diameter (D N ) was estimated from the 3D shape model as the mean value of the maximum distances between two rim points. The depth (d N ) was defined as the maximum distance between a crater point and the plane best-fitting the rim disk. Its rather large uncertainty reflects the complexity of the surface due to the large crater diameter with respect to Julia's diameter. Finally, we used two approaches to constrain the excavated volume. First, we computed the volume from the derived diameter D N and the depth d N assuming two simple crater shapes, namely, a spherical bowl shape and a parabolic shape. The third method consisted in the numerical integration of the spatial coordinates of the crater points. Before proceeding with the integration, we performed a correction of the local slope of the crater (Nonza's rim disk is inclined by approximately −19 • with respect to the local tangent plane). The volumes derived with the three methods fall in the 9000-11 000 km 3 range; we list the mean value in Table 3 (the large uncertainty on the volume is due to the uncertainty on both the diameter and crater depth). We note that this value is a lower limit because the shape of Julia did not allow us to easily take into account the local curvature of the terrain (i.e., because of the high ratio of Nonza's diameter to Julia's diameter).

Constraining the age of the Julia family and the size of the impact crater
The Julia family was first reported in Nesvorný et al. (2015). It is located in the middle main asteroid belt, at high inclination where the number density of asteroids is already low (see Fig. 5). Such a favorable location (i.e., low-density region) clearly facilitated its identification. Nesvorný et al. (2015) identified 33 members and suggested that it must have formed via a recent collisional event. There is therefore the intriguing possibility that one of the craters on the surface of Julia is directly related to the Julia collisional family. In the following sections, we first constrain the basic properties of the family and actually confirm its existence. We then use an N-body orbital model, a Monte Carlo collisional model, and a smoothed particle hydrodynamics (SPH) collisional model to place constraints on the family age and the resulting crater size.

Basic properties of the family
As a first step, we identified the members of the family using the standard hierarchical clustering method (HCM; Zappalà et al. 1995) applied on a recent catalog (version June 2017) of asteroid synthetic proper elements (Knežević & Milani 2003). The cutoff velocity v cut was treated as a free parameter as it should correspond to local conditions, in our case a low background. Because the number of family members N fam depends on v cut and since this dependence turns out to be flat around v cut = 80 m s −1 , we adopted this cutoff velocity and identified 66 family members (see Table B.2 for a complete list) with diameters falling in the 1-2.5 km size range (see Figs. 5 and 6).
No additional criteria were used to remove interlopers, even though we checked the WISE catalog (Masiero et al. 2011) for visual geometric albedos p V , the Sloan Digital Sky Survey (SDSS) Moving Object Catalog (MOC; Ivezić et al. 2001) for color indices a , i−z, and the semimajor axis versus absolute magnitude (a p , H) plot for any outliers. Because of their faintness, we found physical data for only a handful of these bodies; in particular, colors are only available for nine members, which indicate an S-type nature for these objects in agreement with the taxonomic class of Julia. Independently, using the population at slightly higher inclination (sin I p = 0.30-0.32) as a representative background, we do expect at most one interloper among family members.
The family in the space of proper semimajor axis a p , proper eccentricity e p , and proper inclination sin I p is shown in Fig. 6. proper inclination sin I p plot. All asteroids with eccentricities e p = 0.11-0.14 are shown. Symbol sizes are proportional to the logarithm of diameters D. Colors correspond to the geometric albedos: p V < 0.07 is blue, p V > 0.15 (brown and yellow, gray if unknown). The family is close to the 3/1 mean-motion resonance with Jupiter and there are several secular resonances between the J3/1 separatrix and the family (2.51 and 2.55 au). The libration center of the 3/1 mean-motion resonance with Jupiter and its approximate width is denoted by the vertical strip. Approximate positions of the secular resonances z 1 ≡ g − g 6 + s − s 6 , z 2 , and z 3 computed for the mean eccentricity e = 0.125 are plotted as dotted curves. The black dashed rectangle indicates the region where the Julia family is located and which is analyzed in detail in the present work. The region delimited by the gray rectangle (above the black one) was selected as a suitable background population. The large family at sin I p < 0.27 is Maria. The largest asteroid in the neighborhood is (13) Egeria. Julia family members identified at the cutoff velocity v cut = 80 m s −1 are shown in orange. The family is mainly composed of small asteroids (D ∈ (1; 2.5) km) and is well separated from background asteroids.
In the (a p , e p ) plane, it looks like an inclined structure, somewhat similar to an ellipse, which is actually expected for an isotropic ejection if the true anomaly f at the time of impact was close to 180 • . The left-hand part is cut at 2.54 au, most likely due to the proximity of the J3/1 resonance. The right-hand part seems more scattered in e p . Moreover, its overall size is comparable to the escape velocity from Julia, v esc = √ 2GM/R 100 m s −1 if we assume a bulk density ρ = 3300 kg m −3 , considering LL chondrites as a suitable meteorite analog (Vernazza et al. 2014). We note that there is a notable offset in inclination by ∆I = 0.002 rad, likely arising from a cratering event and subsequent ejection of all fragments into a half-space.

Constraining the family age
We constructed an orbital-evolution model in order to constrain the age of the family. It is based on a symplectic N-body integrator from the Swift package (Levison & Duncan 1994), modified according to Laskar & Robutel (2001). This integration scheme makes it possible to use a time step of ∆t = 91 days and a time span that reaches up to 4 Gyr. Our dynamical model contains perturbations by (i) the four giant planets (with a barycentric correction applied to initial conditions) and (13)   which induce a systematic drift in a, but interactions with resonances may also induce drifts in e or I (Brož & Vokrouhlický 2008); and (iii) the YORP effect with thermal torques from Capek & Vokrouhlický (2004), corresponding secular evolution of spins, random collisional reorientations, and a mass shedding at a critical spin rate, as summarized in Brož et al. (2011).
To compute the mean elements, we performed sampling of osculating elements at 1 yr and used the convolution filters A, A, A, B from Quinn et al. (1991), with decimation factors 10, 10, 5, 3. Proper elements were subsequently computed by the frequency-modified Fourier transform (Šidlichovský & Nesvorný 1996) from 512 samples. Planetary frequencies (forced oscillations) were removed to obtain the required amplitudes of free oscillations. Finally, we used a running-window filter 10 Myr wide to suppress remaining oscillations due to longperiod secular resonances. This averaging involves a certain amount of data and delays the initial output (by about 5 kyr for the mean-element filter, 0.8 Myr for the proper-element filter, and 10 Myr for the running window).
Initial conditions of the planets were taken from the DE405 ephemeris, and osculating elements of asteroids from the Astorb catalog (version Dec 2017). All massive bodies and (89) Julia were integrated from the epoch of osculation until the true anomaly reached f = 180 • and ω + f = 80 • (ω: argument of perihelion). We then generated 660 bodies, i.e., ten times the observed number of family members, in order to get a large sample of orbits, so that these could be easily resampled off-line. For simplicity, we assumed an isotropic velocity field, and the velocity distribution of Farinella et al. (1994) with a slope α = 1.25 and a maximum velocity v max = 500 m s −1 (to also obtain some outliers). This distribution peaks at about the escape velocity v esc . We also assumed initially isotropic spins, and a uniform distribution of spin periods, P ∈ (2; 10) h. Thermal parameters of synthetic bodies were selected as follows: bulk density ρ = 3300 kg m −3 , surface density (regolith) ρ surf = 1500 kg m −3 , heat capacity C = 680 J kg −1 K −1 , thermal conductivity K = 10 −3 W m −1 K −1 , Bond albedo A = 0.10, and infrared emissivity = 0.9. Diameters were taken from the observed size frequency distribution (SFD), with ten clones for each.
The orbital evolution of the synthetic family is shown in Fig. 7 (left column). Qualitatively, it is clear that in about 100 Myr the synthetic family becomes too dispersed compared to the observations. However, we cannot directly (quantitatively) compare the outcome of our N-body simulations with the observations as there are three limitations at least (Brož 2016): (i) observed asteroids have to be carefully selected, not only the family identified by the HCM, but also its surroundings where the bodies can be scattered to (cf. Brož & Morbidelli 2013); (ii) the synthetic SFD is different and even changes over the course of the simulation; and (iii) there is an inevitable contribution of background asteroids, and this background can be variable.
To overcome these limitations, we post-processed the output of our N-body simulations. As a preparation, we had to select all observed asteroids which encompass the family, i.e., with a p ∈ (2.52; 2.58) au, e p ∈ (0.11; 0.14), and sin I p ∈ (0.28; 0.30). A suitable background population was chosen with nearly the same parameters, except sin I p ∈ (0.30; 0.32) (see the rectangles in Fig. 5).
We then proceeded with the method described in detail in Brož (2018). It accounts for the surroundings, includes the background population, and matches the SFDs in every single time step, using a random selection of synthetic asteroids for all given size bins (D, D + dD). Finally, we computed the numbers of both synthetic and observed asteroids in four boxes, defined by ∆a = 0.03 au, ∆e = 0.015 (see Fig. 7 middle column), and the metric where the uncertainties were assumed Poisson-like (σ = √ N).
According to Fig. 7 (right column), it is clear that systematic good fits are only possible for ages t 120 Myr, our optimal solution being for t = 30 Myr. The lower limit t 10 Myr stems from the fact that the left-hand part of the family inside 2.54 au needs to be dispersed. We note that we obtained similar results with an entirely different approach (see Appendix A), which shows that our estimated age is robust.

Constraining the size of the impact crater and its location
As a next step we performed numerical simulations of the collision that led to the currently observed Julia family using an SPH model. These simulations allowed us to place constraints on the total ejected mass and the crater size.
As a preparatory task, we employed a set of Durda et al. (2007) size-frequency distributions and a simple scaling law (from D = 100 km up to D = 150 km) discussed therein to estimate preliminary impact parameters from the observed SFD. Given the size/mass difference between the largest remnant, i.e., (89) Julia, and the second largest fragment within the Julia family, it is not surprising that the best-fit SFD corresponds to a small cratering event that occurred with an oblique angle and imparted only a limited amount of kinetic energy on to the target. Optimal A154, page 7 of 16 A&A 618, A154 (2018) P. Vernazza et al.: Asteroid (89) Fig. 7. Proper semimajor axis a p vs. the proper eccentricity e p for the synthetic family (left column). It is an output of the N-body simulation, which initially included 660 bodies (i.e., 10 times more than observed) and later rescaled to the same size-frequency distribution as the observed one; a random background population was also included (middle column). Colors correspond the numbers of bodies N box in the respective boxes. Finally, there is the observed family, and the corresponding χ 2 metric (right column). The "initial" conditions are shown at t = 0 Myr (top row), the best fit at 30 Myr (middle row), and a very poor fit at 520 Myr (bottom row) for comparison. Dotted lines on the χ 2 (t) plot correspond to the 1-σ and 3-σ levels. A reasonable match is only possible for young ages t 100 Myr. values for the expected projectile size, the impact velocity, and the impact angle were found to be d = 8 km, v imp = 6 km s −1 , and ϑ imp = 75 • , respectively.
We then used the following tools for computations: the SPH5 code for the fragmentation phase (Benz & Asphaug 1994) and PKDGRAV for the reaccumulation phase (Richardson et al. 2000;Stadel 2001). For SPH5, we assumed the Tillotson (1962) equation of state, the von Mises (1913) yielding criterion, and the Grady & Kipp (1980) fracture model. Material parameters in our simulations correspond to those of basalt 2 , except the zero-pressure density ρ 0 = 3.3 g cm −3 . The time step ∆t was controlled by the Courant criterion with a factor of 1.0. The time span (100 s) was set to a higher value than the back and forth travel time of the shock wave (2D/v imp ), yet to a lower value than the ejection timescale (v ej /g). Artificial viscosity parameters were set to slightly higher values than the standard ones (α av = 4.0, β av = 8.0). We used a modification of the scalar damage D as explained in Ševeček et al. (2017, Appendix C). The number of SPH particles was set to N part 7.0 × 10 5 .
The handoff, or conversion of the SPH particles to solid spheres, was performed with the relation R i = [3m i /(4πρ i )] 1/3 , which allows for some expansion below ρ 0 . During the reaccumulation phase, the self-gravity was computed approximately, with gravitational moments up to hexadecapole order and a tree opening angle θ = 0.5 rad. We assumed a perfect merging so that no shape information was preserved during this stage. The time step was set to 10 −6 (in G = 1 units) and the time span to 50 000 ∆t so that the reaccumulation phase was definitively over.
The results of our simulations are as follows: -The diameter of the resulting crater (see Fig. 8) is at least 60 km. We note that measuring the diameter of the resulting crater is not a straightforward procedure as we only have a transient crater (which could become larger due to later landslides) at the end of the fragmentation phase. -The size of the second largest family member matches the observed one quite well (3 km versus 2 km). -The excavated volume at the end of the fragmentation phase is V ex 7600 km 3 0.005 V Julia , and the ejected volume computed after the gravitational reaccumulation is V ej 176 km 3 , i.e., only 10 −4 V Julia (the uncertainties are at least on the order of 10 %). Our results thus show in a clear manner that V ex V ej , which implies the formation of a large crater and a small total volume for the family members (excluding Julia), in agreement with the currently observed value. Also, the excavated volume is in excellent agreement with that measured for Nonza (see Table 3). The crater also sets constraints for the material strength and other properties (ρ, µ, Y, k, m), but this would require a large number of simulations, and a better physical model with gravity in the fragmentation phase, to describe a collapse of the crater and eventual surface features. This is postponed to a future work.
As a final step, we analyzed the velocity field resulting from the SPH model. As expected for a cratering event, ejecta can only fly to a half-space. The ejection velocity with respect to the barycenter (without outliers like projectile fragments) is on the order of v ej 100 m s −1 . Using the Gauss equation (e.g., Murray & Dermott 2000; n being the mean motion and ∆v W the out-of-plane velocity component), this translates to for a suitable orientation of the ejection half-space, i.e., just above the orbital plane. This is indeed comparable to the ∆I value observed for the Julia family (as discussed in Sect. 5.1).
Since the orbital angular momentum L orb of the impactor is only about 1/10 of the rotational angular momentum L rot of Julia, it is reasonable to assume that even a highly oblique impact could not have significantly modified Julia's rotation axis. Given the orientation of Julia's spin axis and the orientation of its orbit, the obliquity of Julia is γ = −17 • . Because the latitude φ = γ is the one at which ejecta can fly the most above the orbital plane (in W direction), it appears that Nonza (located around φ −32 • ) is an ideal candidate for being the impact crater at the origin of the Julia family. Specifically, our observations covered almost completely (99.4%) the surface that lies between latitudes 33 • and −67 • (which is the region between γ + 50 • and γ − 50 • ) and Nonza is the only D ≥ 60 km crater that we could recognize in this region of interest. We cannot exclude that we missed a similarly sized crater in this region of interest as it was partly covered by lower resolution images (e.g., rotational phase: 0.0). We also recall here that Nonza is not clearly visible at rotational phase 0.53, illustrating that the reconnaissance of craters is highly dependent on the observing conditions (e.g., seeing) and/or the viewing geometry in the case of groundbased AO observations. Independently, the current main-belt impact rate suggests that only a few large craters (D ≥ 60 km) should have formed at the surface of Julia over the course of the solar system's history. In the case of Lutetia (D ∼ 100 km), it was suggested that an impact such as the one that formed the ∼55 km wide Massilia impact basin should occur every 9 billion years (Sierks et al. 2011), implying that D ≥ 60 impact basins should be rare at the surfaces of D ∼100-150 km asteroids (see Appendix A for more details). The low frequency of such collisional events gives additional support to the fact that Nonza may indeed be the impact crater at the origin of the ∼100 Myr old Julia family.

Conclusions
Here, we evaluated via several complementary approaches the recently proposed hypothesis that the S-type asteroid (89) Julia is the parent body of a small compact asteroid family that formed via a cratering collisional event. We observed Julia with VLT/SPHERE/ZIMPOL throughout its rotation during eight apparitions between July 8 and October 7, 2017, and derived its 3D shape as well as the topography of the largest craters. We also performed numerical simulations to first confirm the existence of the Julia family and to determine its age as well as the size of the impact crater at its origin. Finally, we utilized the images/3D shape in an attempt to identify an origin location consistent with the proposed small collisional family.
Our numerical simulations suggest that (89) Julia was impacted 30-120 Myr ago by a D ∼ 8 km asteroid, thereby creating a D ≥ 60 km impact crater at its surface. Given the small size of the impactor, the obliquity of Julia, and the orientation of the family in the (a,i) space, the impact most likely occurred in the southern hemisphere of Julia, close to latitude φ = −17 • . Our VLT/SPHERE observations reveal that Julia's southern hemisphere hosts one such large crater (D ∼ 75 km; centered around φ = −32 • ) that therefore appears to be an ideal candidate for being at the origin of the ∼100 Myr old Julia family.
The same type of science investigation that could be performed 20 yr ago with HST in the case of (4) Vesta and the discovery of its south pole impact basin at the origin of the Vesta family (Thomas et al. 1997a), can now be performed for many D ≥ 100 km main-belt asteroids with VLT/SPHERE. New opportunities for ground-based asteroid exploration, namely, geophysics and geology, are becoming available thanks to VLT/SPHERE's unique capabilities. Also, the present work may represent the beginning of a new era of asteroid-family studies. In these fields (geophysics, geology, and asteroid-family studies), the future will only get brighter with the arrival of 30-40 m class telescopes (ELT, TMT, GMT). Here, we used a method independent from the one presented in the main text to constrain the age of the family. To perform this research, we used a Monte Carlo code called Boulder (Morbidelli et al. 2009), which computes the evolution in time of size-frequency distributions due to fragmentation (and reaccumulation). We performed several modifications to the original code, the most relevant being that we took the size-dependent dynamical decay due to the Yarkovsky effect and neighboring resonances into account, including a handling of fractional probability (as in Cibulková et al. 2014). We also realized that the parametric relation describing the mass of the largest fragment M lf (Q/Q D ) and its dependence on the specific impact energy Q, divided by the strength Q D , does not fit our high-resolution SPH simulations from Sect. 5.3 in the region of interest Q/Q D 0.002. Indeed, the original fit was performed in linear space and its extrapolation towards low Q values was offset, predicting M lf too large by a factor of 10, or D lf by a factor of 10 1/3 2. In the Monte Carlo model, this offset would push up the number of Julia family forming events. We thus performed the fit again in the log-log space; we incorporated our low-Q SPH simulations, which provided upper limits for the value of M lf , because the largest fragment was barely resolved, and obtained the following relation (see also Fig. A.1): In order to account for the dynamical decay, we computed the relative numbers of asteroids N(t)/N 0 located within the respective regions. They were approximated by exponentials (N/N 0 = exp(−t/τ)); the respective timescales τ for the five different bin sizes (four sizes in the D ∈ (0.5; 2.5) km range and 150 km which corresponds to the size of Julia) were as follows: 226,434,595,424,and 42 200 Myr. In the next step, we computed the intrinsic collisional probabilities of the Julia region and of the main belt, using the algorithm developed by Bottke & Greenberg (1993): P i = 2.99 × 10 −18 km −2 yr −1 within the main belt, 3.84 × 10 −18 km −2 yr −1 for mutual collisions, and 6.10 × 10 −18 km −2 yr −1 within the family (even though these collisions are clearly negligible). The complementary median impact velocities are v imp = 5.04, 6.62, and 6.57 km s −1 . The main belt decay was taken from Bottke et al. (2005) and Julia's decay from the current paper.
Physical properties of the populations were described by the scaling law Q D (r) = Q 0 r a + Bρr b , where we adopted the parameters from Benz & Asphaug (1999) for basalt at 5 km s −1 , except for the Julia family, for which we used a slightly higher density (ρ = 3.3 g cm −3 ) in agreement with its LL-like spectral properties (Vernazza et al. 2014).
The initial SFDs were rather similar to the observed values except the tails, which were prolonged down to D min = 0.005 km or eventually bent down below a slope q = −3 to prevent a divergence in mass. The time step was set to ∆t = 10 Myr (or smaller if necessary) and the time span up to 4 Gyr. The model had to be run multiple times because of the fractional probabilities.
The results are shown in Fig. A.2 (left column). It turns out that Julia's neighborhood is replenished continually with small asteroids as collisions with (89) Julia are relatively frequent. On the other hand, if we focus on a single event (right column), without (89) Julia itself which continually creates fragments, this transient population decays within a few 100 Myr, partly because of the dynamical decay and partly due to collisional grinding. This young age is in perfect agreement with that derived via our orbital model.
We also computed an overall statistics of such cratering events, producing largest fragments D lr ≥ 2.6 km, and we expected up to 10 events over 4 Gyr (from 100 MC runs). Even though the total ejected mass is still relatively small, and the bulk of the 150 km parent body is preserved, it may indicate a substantial resurfacing and create an overall irregular shape if the number of events is 1. The main belt SFD is plotted in blue and the family SFD in red; the respective initial conditions are plotted in cyan and yellow, and the observations in gray; the bend at around D 1 km is mostly due to the observational incompleteness. The situation at time t = 30 Myr (top row), and t = 520 Myr (bottom row) is shown. There are always 10 runs with a different random seed in order to see lower probability events. The first model demonstrates that the Julia family is by no means an exception; the second model shows that small fragments (from a single breakup) decay quickly within a few 100 Myr at most.