Giant impact on early Ganymede and its subsequent reorientation

The degree of isostatic compensation and the value of \({{\varvec{g}}}_{20}\)

Following Melosh¹⁹, we assumed the following conditions for a crater basin: (i) the ejecta blanket of excavated material has a uniform thickness, \({t}_{e}\), and a simple annulus shape with inner and outer radii of \(R\phi\) and \(R\theta\), respectively, where \(\phi\) and \(\theta\) are the angles from the center of the object, respectively, and the outer radius is twice the inner radius, \(\theta =2\phi\); and (ii) the basin floor has a uniform depth, \(h\), and a simple circular shape with a radius of \(R\phi\). In this case, a simple geometric solution gives the volume of the ejecta blanket and hole as follows:

In general, both volumes are equal because of the conservation of mass¹⁹, even if the very shallow basin was created after the collapse of a transient crater; subsequently, we can use the simple relation of \(h={t}_{e}(\text{cos}\phi -\text{cos}\theta )/(1-\text{cos}\phi )\). We considered the variations in the density of a nearly spherical object expressed by variations over a series of concentric spherical shells, and the density perturbations were azimuthally symmetric. The density perturbation in this model of the ejecta blanket model is expressed as follows:

where \({\rho }_{0}\left(r\right)\) is the mean density for radius \(r\), \({\rho }_{c}\) is the density of the crust, and \(\varphi\) is the colatitude from the crater center. The degree-2 order-0 coefficient of the density perturbation is expressed as follows⁶²:

We can obtain.

The degree-\(l\) order-\(m\) spherical harmonic coefficients of the gravitational potential at an arbitrary distance, \(s (s\ge R)\), are given by⁶²:

where \(s\) is the distance from the center. Therefore, if this density perturbation is supported by an infinitely rigid lithosphere, the degree-2 order-0 gravitational potential coefficient at the object surface can be described as follows:

If we assume that \(h={t}_{e}(\text{cos}\phi -\text{cos}\theta )/(1-\text{cos}\phi )\)), this gravitational coefficient is always negative, \({g}_{20}, which gives \(Q.

The actual shapes of the eject blankets in Fig. 3 are not simple annulus shapes. However, this simple annulus provides a good approximation of the ejecta distribution because the wavelength of the degree-2 gravitational potential is sufficiently longer than that of the basin radius and mostly averages inhomogeneity or asymmetry in the ejecta distribution. As a simple example, if we define \(\theta =1.3\phi\) or \(\theta =2.5\phi\) instead of \(\theta =2\phi\), the value of \(Q\) changes by only approximately 20–30%. Inhomogeneity in the tangential direction should have little effect on \({g}_{20}\) if the total mass of ejecta as a function from the center is equal to the simple annulus shape because Eq. (6) does not depend on the tangential direction.

A mass excess/deficit that is compensated for by isostasy does not create a free-air gravity anomaly, whereas the ejecta mass supported by the bending and membrane stresses of the lithosphere creates a gravity anomaly. Therefore, the gravity anomaly decreases by a factor of \((1-{C}_{n})\), where the degree of isostatic compensation, \({C}_{n}\), is described as follows⁶³:

where \(N=n\left(n+1\right)\), \(n\) is the degree of spherical harmonics of the horizontal width of the load (\(n=2\pi R/w\), where \(w\) is the horizontal width of the load),\(\overline{\rho }\) is the mean density of the satellite, and \(\tau\) and \(\sigma\) are nondimensional parameters that are defined as follows:

where \(t\) is the elastic thickness;\({\rho }_{m}\) is the density of the mantle underlying the lithospheric plate; \(E\) and \(\nu\) are the Young’s modulus and Poisson’s ratio of the lithosphere, respectively; and \(g\) is the surface gravity. Here, \({C}_{n}=0\) indicates that the load is fully supported by the lithosphere, and \({C}_{n}=1\) indicates that the load is fully compensated by isostasy. We assumed that \(E\) = 9 GPa, \(\nu\) =0.32 (ref. 64), \({\rho }_{m}\) = 1000 kg m⁻³, and \({\rho }_{c}\) = 930 kg m⁻³ (ref. 65). We assumed that \(w=\) 1000 km because the ejecta blanket spreads across several basin radii. Note that the value of \(w\) hardly affects \({C}_{n}\): \({C}_{n}=0.871\) (\(w=\) 3000 km),\({C}_{n}=0.875\) (\(w=\) 1000 km), and \({C}_{n}=0.868\) (\(w=\) 500 km) for \(t=\) 10 km, and \({C}_{n}=0.992\) (\(w=\) 3000 km),\({C}_{n}=0.993\) (\(w=\) 1000 km), and \({C}_{n}=0.993\) (\(w=\) 500 km) for \(t=\) 0.5 km. Consequently, we assume \({C}_{n}=0.87\) for \(t=\) 10 km and \({C}_{n}=0.99\) for \(t=\) 0.5 km.

If the value of \({C}_{n}\) is homogeneous inside and outside the basin, the gravity anomaly is always negative (\({g}_{20}); thus, reorientation that orients the basin center to the tidal axis does not occur. At least \({C}_{n}\) outside the basin must be smaller than \({C}_{n}\) inside the basin to obtain \({g}_{20}>0\). Nimmo et al.¹² assumed that isostasy is archived inside the Sputnik Planitia basin. Similarly, in the case of Ganymede, the gravity anomaly becomes positive (\({g}_{20}>0\)) if (i) isostasy is archived inside the basin, which has the same attribute as \(h=0\), and (ii) it is not archived outside the basin (i.e., the load of the ejecta is partially supported by the lithosphere). This assumption is natural for the following reasons. First, numerical simulations²⁶ or impact experiments⁶⁶ have proposed that the impactor can penetrate Europa’s ice shell and create conduits to the underlying ocean. Similarly, if Ganymede has an ocean and its ice shell is sufficiently thin, isostasy is easily achieved within the basin floor by the inflow of water. The furrow-forming impact alone may not create a large ocean because the size of a local melt pool created by the impact is up to 6 times the impactor radii⁶⁷, which is less than 10% of the total volume of Ganymede. However, isostasy would be achieved within the basin floor even if Ganymede’s interior is warm convective ice without an ocean. Numerical simulations of warm convective ice without the ocean^11,27 have shown that cold lithospheric material is replaced by warm subsurface material and that the postimpact lithospheric thickness becomes quite thin within the transient crater radius when the transient crater size is sufficiently larger than the lithospheric thickness. The lithospheric thickness recovered by cooling from the surface is roughly determined by the characteristic thermal diffusion depth, \(\sqrt{\kappa {t}_{d}}\), where \(\kappa\) is the thermal diffusivity and where \({t}_{d}\) is the timeframe⁶⁸. The relaxation time of a load with a horizontal width of 100 km is 25 years (see below). If we use \(\kappa =\) 2.9 \(\times\) 10^–5 for ice^69,70 and \({t}_{d}=\) 25 years, we obtain \(\sqrt{\kappa {t}_{d}}\)= 150 m. An elastic thickness of \(t=\) 150 m does not support a load above the lithosphere at all \(({C}_{n}=1)\). In other words, the mass anomaly within the basin floor quickly disappears before the lithospheric thickness recovers. Therefore, isostasy can be achieved within the basin floor of the furrow system, regardless of the two surface models used. Second, numerical simulations^11,27 have shown that the postimpact lithospheric thickness and temperature do not become thin and warm outside the basin. Numerical simulations for the Sputnik Planitia⁷¹ also show that hot preimpact thermal structures or cold thermal structures with an ocean greater than 150 km thick produce basins that are broadly isostatically compensated inside the basin and that are not fully isostatically compensated outside the basin. Therefore, it is likely that isostasy is not fully archived outside the basin. Third, furrow formation is unlikely to cause a significant change in the rheology of the lithosphere overall before and after the furrow-forming impact because individual furrows occur in the brittle part of the lithosphere^6,7, and the brittle part above the brittle‒ductile transition layer should already include many cracks and faults due to impacts or tectonics. Notably, if the lithosphere is fully damaged by impact and isostasy is fully and broadly archived even outside the basin, any gravity anomalies will not be created by ejecta. If so, we should conclude that a coincidence or another hypothesis, such as the impactor-origin silicate mascon beneath the basin or the variation in the lithospheric thickness, is responsible for the current position of the center of the furrow system coinciding with the tidal axis.

Based on the above assumption (\(h=0\)), the dimensionless parameter becomes.

or

The radius of the transient crater of the Gilgamesh basin is 135 km (ref. 46), which indicates a load of \(Q\)=0.027 (\(\Delta {k}_{2}\)=0.6) or 0.54 (\(\Delta {k}_{2}\)=0.03). The Gilgamesh basin does not cause significant reorientation of Ganymede, although a reorientation caused by the Gilgamesh basin has previously been proposed²⁸. As another example, the mound at the sub-Jovian point, discovered by the Juno flyby and legacy data reanalysis, with a height of 3 km and an oval of 450 km × 750 km (refs. 49,50), has a value of \(Q\)=0.018 (\(\Delta {k}_{2}\)=0.6) or 0.36 (\(\Delta {k}_{2}\)=0.03), assuming \(h\)=-3 km, \({t}_{e}\)= 0, \(R\phi\)=300 km, and \({C}_{n}=0\).87 (t = 10 km). Following Eq. (41) in Matsuyama and Nimmo¹⁶, the mound leads to very small reorientation of less than 2.6° if its latitude is 45°. We obtain a value of \(Q\)=0.17 for the mound if we assume that \(h\)=-3 km, \({t}_{e}\)= 0, \(R\phi\)=300 km, \({C}_{n}=0.\) 99 (t = 0.5 km), and \(\Delta {k}_{2}\)=0.005.

The degree 2 Love number of Ganymede

The value of \(\Delta {k}_{2}={k}_{2}^{T*}-{k}_{2}^{T}\) is determined by the internal structure of Ganymede and has a significant effect on the value of \(Q\). Many studies have investigated the tidal response of Ganymede, assuming various interior structure models^72,73,74. For example, Moore and Schubert⁷² reported that \({k}_{2}^{T}\) is approximately 0.01 if the present-day Ganymede has no ocean and its mantle has a high viscosity of \(\eta >\) 10¹⁴ Pa s or approximately 0.4 to 0.6 if Ganymede has an ocean or its mantle has a low viscosity of \(\eta 10¹² Pa s. Although it is unclear whether Ganymede’s interior is differentiated or undifferentiated, we estimate the Love number on the basis of a four-layer Ganymede model. Here, we assume that Ganymede consists of (1) a 10 km-thick elastic icy shell, (2) a viscoelastic icy mantle, (3) a viscoelastic rocky mantle, and (4) a liquid metallic core, where the rigidities of the icy shell and rocky mantle are 10 and 100 GPa, respectively; the densities of the icy shell, icy mantle, rocky mantle, and metallic core are 1050, 1050, 1050, 3100, and 5150 kg/m³, respectively; the outer radii of the icy shell, icy mantle, rocky mantle, and metallic core are 2638, 2628, 1745, and 710 km, respectively; and the viscosity of the rocky mantle is 10²⁰ Pa s. The value of \({k}_{2}^{T*}\) is defined as the value when the elastic icy shell is assumed to be a liquid ocean layer. Those rigidities, densities, outer radii, and viscosities are from Kamata et al.⁷³. We use a code released by Kamata⁷⁵ to calculate the periodic spheroidal deformation of a planetary body to obtain the Love number. As a result, when the viscosity of the icy mantle is 10¹¹, 10¹³, or 10¹⁵ Pa s, we obtain (\({k}_{2}^{T},{k}_{2}^{T*})\) = (0.579,0.604), (0.335,0.605), or (0.018,0.605), respectively. If the subsurface ocean layer (with a density of 1050 kg/m³) is present instead of the viscoelastic icy mantle, we can obtain (\({k}_{2}^{T},{k}_{2}^{T*})\) =(0.579,0.604). Considering the large uncertainties, we studied the two cases of \(\Delta {k}_{2}\)=0.03 and \(\Delta {k}_{2}\)=0.6 for t = 10 km. If we assume a 0.5 km-thick elastic icy shell instead of a 10 km-thick elastic icy shell, when the viscosity of the icy mantle is 10¹¹, 10¹³, or 10¹⁵ Pa s, we obtain (\({k}_{2}^{T},{k}_{2}^{T*})\) =(0.599,0.604), (0.446, 0.605), or (0.018,0.605); therefore, we study the two cases of \(\Delta {k}_{2}\)=0.005 and \(\Delta {k}_{2}\)=0.6 for t = 0.5 km. This combination of \(\Delta {k}_{2}\)=0.6 and t = 0.5 km may be unlikely because the very thin lithosphere and cold interior are not consistent.

Tidal decay of rotation

Following a large impact event, tidal friction quickly dampens the motions of the satellite’s orientation, such as librational oscillations and nonsynchronous rotation, and the satellite attains its new orientation on a short timescale¹⁹. The timescales of the tidal decay of librational rotation and nonsynchronous rotation have been estimated in previous studies^19,76,77, and all three studies were equivalent. Following Noyelles et al.⁷⁷, the timescale is expressed as follows:

where \(C\) is the maximum moment of inertia about the center of mass; \({n}_{o}\) is the mean orbital motion; and \({Q}_{d}\) is the specific dissipation factor of the satellite. Assuming a homogenous interior (\(C=0.4 M{R}^{2}\)) and reasonable values for warm icy satellites,\({k}_{2}^{T}\)=0.5 and \({Q}_{d}\)=100 (ref. 78), we can obtain \({T}_{damp}\)=870 years. If Ganymede was cold and rigid (\({Q}_{d}/{k}^T_{2}\) ~ 10⁵), we can obtain \({T}_{damp}\)=4.4 million years.

Viscous relaxation of the load and bulge

The mass anomaly not supported by the lithosphere disappears owing to viscous relaxation on a timescale of \(\eta /{\rho }_{m}gw\) (ref. 68), where \(w\) is the horizontal width of the load. The icy materials underlying the elastic layer are convecting, and their reference viscosity is as low as \(\eta\)= 10¹² to 10¹⁷ Pa s (ref. 79). For example, a 100 km-wide load not supported by the lithosphere would vanish in less than 25 years. Similarly, the density perturbations below the lithosphere created by heating and redistribution of the impact crater would vanish on a short timescale. If the relaxation time of the remnant tidal/rotational bulge is given by \(\eta /{\rho }_{m}gR\) (ref. 17), the bulge would vanish in 1 year. However, this estimate is suitable only for an isoviscous case, whereas the internal structure of Ganymede is much more complicated and, in particular, may have a rigid outer shell that can maintain a remnant bulge over geological time. In short, the actual bulge relaxation timescale is extremely uncertain. It is estimated that old craters on the dark terrain of Ganymede have relaxed viscously for 1 billion years^55,80. The timescale of relaxation of the mass anomaly of the ejecta blankets and the tidal/rotational bulge may also be approximately 1 billion years.

The Valhalla basin

The original crater radius of the Valhalla basin was estimated to be 500 km on the basis of the mapping of ejecta and secondaries⁸¹. The thickness of the lithosphere was estimated to range from 15 to 20 km during the formation of the Valhalla basin⁶. Assuming \(R\phi\)=500 km, \(t=\) 20 km,\(h=\) 0, \(w\)= 500 km, \(n\)= 15, \({C}_{n}\)=0.66, and physical parameters for Callisto, we can obtain \(Q\)=199.8 (\(\Delta {k}_{2}\)=0.03) or \(Q\)=10.0 (\(\Delta {k}_{2}\)=0.6) for the Valhalla basin, which should lead to a significant reorientation of Callisto. Possible explanations for this include Callisto’s rotational period, incomplete isostatic compensation within the basin floor, and other surface loads. If Callisto’s rotational period was 50 h (currently 400 h), we could obtain \(Q\)=0.16 (\(\Delta {k}_{2}\)=0.6). The Valhalla basin is one of the oldest surface features on Callisto^21,82, and Callisto at the time of the formation of the Valhalla basin may not have yet reached a synchronous rotation state. However, this scenario may be unlikely because Callisto attains a synchronous rotation state in 2 million years if we use \({Q}_{d}/{k}_{2}\) ~ 20 (ref. 83). If Callisto was cold and rigid (\({Q}_{d}/{k}_{2}\) ~ 10⁵), Callisto obtained a synchronous rotation state 1 billion years after its formation. Note that because Ganymede attains a synchronous rotation state at 0.66 million years for \({Q}_{d}\)=100 after its formation⁷⁶, Ganymede’s rotational period at the time of furrow formation should not differ much from the present period. If isostasy is not achieved within the basin floor of the Valhalla basin (i.e., \(h>0\)), incomplete isostatic compensation of the negative topographic profile within the basin would decrease \(Q\). Finally, as many areas of Callisto have not yet been imaged with sufficient resolution, we cannot rule out the possibility that large undiscovered loads exist on Callisto’s surface. Additionally, Callisto may have had a large remnant bulge, such as the equatorial bulge of Iapetus, because of Callisto’s undifferentiated and inert interior.

Note for the calculation of the ejecta blanket

Many studies have focused on the shape of the ejecta blanket, which shows that the mass distribution of ejecta is strongly sensitive to the impact incidence angle and highly asymmetric^{84,85,86,87,88}. Therefore, most naturally, the impact incidence angle is responsible for the asymmetric load to explain the location of the center of the furrow system. However, there could be various alternative explanations, as we described above.

It is not unnatural that a thick ejecta blanket would be superimposed on where the furrows lie because this is also the case in the Valhalla basin. The Valhalla basin can be divided into three distinct zones outward from the center: a central smooth zone (r_b r_b is the distance from the basin center), an inner ridge and trough zone (360 km r_b r_b 3,21. If troughs and grabens were created first, they would be damaged by ejecta particles. Therefore, troughs and grabens were likely created after the accumulation of ejecta. Similarly, no structures older than the furrows have been found on the Ganymede surface, and no craters cut by the furrows have been found³. Therefore, all craters on Ganymede were likely created after furrow formation (perhaps, all ancient craters created before furrow formation were removed by ejecta). Compared with ejecta particles flying in a ballistic trajectory, the timing of furrow formation formed by the asthenospheric flow would be delayed.

The thickness of the ejecta blankets at the rim is 100 km in our calculation; however, this is unrealistic. This is mainly because our ejecta model does not consider any subsequent movement after the landing of ejecta particles, which creates a point of discontinuity in the thickness along the basin rim. In reality, some dispersion induced by mass wasting or re-ejection should occur. The ejecta volume at a point of discontinuity is not large; thus, we expect that it hardly affects our results. Nevertheless, we recalculated the ejecta distribution in the vicinity of the basin rim, assuming a realistic surface slope. Specifically, we removed ejecta steeper than a given slope (\(\uplambda\)) around the basin rim, calculated the moment of inertia tensor of the mass distribution of the ejecta blanket whose steep part was removed, and remapped the ejecta blanket after the most stable reorientation (Supplementary Fig. 1). Here, we assume that the removed ejecta do not create a mass anomaly. Although this calculation does not assume a particular geologic process, it is known that a portion of the ejecta rim moves into the basin floor due to lateral movement, such as mass wasting. Thus, the ejecta mass within the basin does not create a free-air gravity anomaly because we assume that isostasy is archived inside the basin. The 15 examples in this figure correspond to the 15 examples shown on the right side of Figs. 3, 4, and 5. As a result, among the 15 examples, only #3 for \(\uplambda \le\) 20 degrees is different from the original example. For the rest, we find that almost the same principal axis rotation state is obtained. It is unclear how much gradient is actually allowed on Ganymede’s surface; however, this figure shows that the discontinuities around the basin rim hardly affect our results. This is mainly because the unrealistically thick ejecta rim is created by the accumulation of ejecta in a very small area, and the removed ejecta mass was approximately less than 10% of the total mass when \(\uplambda \ge\) 20 degrees and less than 20% when \(\uplambda =\) 5 degrees.

We performed similar numerical simulations for Pluto. The simulations are almost the same as those for Ganymede, other than the trajectories of the ejecta particles. The trajectories were solved via the equation of motion for the circular restricted three-body motion (Pluto, Charon, and an ejecta particle)⁴⁸. Supplementary Fig. 2 shows a similar probability as a function of the transient crater radius, similar to that in Fig. 7. The mean radius of the Sputnik Planitia is 550 km (ref. 13), which corresponds to a transient crater radius of 420 km. For this transient crater radius, the probability that the impact center is within 5° from the point that deviates poleward by 20° from the tidal axis is approximately 20–30%. Therefore, the center of the Sputnik Planitia of Pluto, which, together with the furrow system, deviates poleward by 20° from the tidal axis, can be explained by a similar model.