# Erosion and accretion by cratering impacts on rocky and icy bodies

Ryuki Hyodo ISAS, JAXA, Sagamihara, Japan  Hidenori Genda Earth-Life Science
Institute, Tokyo Institute of Technology, Tokyo, 152-8550, Japan

###### Abstract

During planet formation, numerous small impacting bodies result in cratering
impacts on large target bodies. A fraction of the target surface is eroded,
while a fraction of the impactor material accretes onto the surface. These
fractions depend upon the impact velocities, the impact angles, and the escape
velocities of the target. This study uses smoothed particle hydrodynamics
simulations to model cratering impacts onto a planar icy target for which
gravity is the dominant force and material strength is neglected. By
evaluating numerical results, scaling laws are derived for the escape mass of
the target material and the accretion mass of the impactor material onto the
target surface. Together with recently derived results for rocky bodies in a
companion study, a conclusion is formulated that typical cratering impacts on
terrestrial planets, except for those on Mercury, led to a net accretion,
while those on the moons of giant planets, e.g., Rhea and Europa, led to a net
erosion. Our newly derived scaling laws would be useful for predicting the
erosion of the target body and the accretion of the impactor for a variety of
cratering impacts that would occur on large rocky and icy planetary bodies
during planet formation and collisional evolution from ancient times to today.

planets and satellites: formation – planets and satellites: dynamical
evolution and stability – methods: numerical

††journal: ApJ

##  1 Introduction

Cratering impacts – collisions of small impacting bodies on large planetary
bodies – erode a portion of the large planetary surface by ejecting materials
with a speed greater than the escape velocity for the large planetary body.
Simultaneously a fraction of impactor’s material is accreted onto the target
surface (e.g., Housen et al., 1983; Svetsov, 2011; Housen & Holsapple, 2011;
Hyodo & Genda, 2020, and references therein). The erosion of the target’s
surface materials and the accretion of the impactor’s materials are key
processes that dictate the mass gain/loss of planetary bodies and that
characterize the surface composition (e.g., Melosh, 1984, 1989).

Using a point-source model, the crater size is assumed to be much larger than
that of the impacting body, an analytical scaling law that predicts the amount
of the ejecta $m_{\rm eje}(>v_{\rm eje})$ whose ejection velocity is larger
than a given value $v_{\rm eje}$ is derived (e.g., Housen & Holsapple, 2011;
Holsapple & Housen, 2012). The amount of the target material that is eroded by
escaping from its own gravity upon cratering impacts (hereafter, the escape
mass $m_{\rm esc,tar}$) increases with increasing impact velocity $v_{\rm
imp}$ compared to the escape velocity of the target $v_{\rm esc}$.

In our companion study (Hyodo & Genda, 2020, hereafter Paper 1), the fraction
of escape material from the target body $m_{\rm esc,tar}$ was studied for a
variety of cratering impacts for rocky bodies. Using smoothed particle
hydrodynamics (SPH) simulations, Paper 1 demonstrated that the point-source
scaling law agreed with the numerically derived escape mass of the target when
$v_{\rm imp}\gtrsim 12v_{\rm esc}$. Otherwise, it was found to overestimate
the escape mass when the ejection velocity of the target material is only
moderately larger than the escape velocity of the target.

The escape mass as well as the applicability of the point-source scaling law
have not been studied for the case of icy bodies. In the outer solar system,
icy bodies are ubiquitous. Following the same approach as Paper 1, an
extensive number of cratering impact simulations onto a planar icy target were
performed. These simulations consider a case where gravity dominates and the
material strength may be neglected. We independently study the high-speed
escaping ejecta (i.e., $v_{\rm eje}>v_{\rm esc}$) that originate from the icy
target or the icy impactor. A goal of this study is to derive scaling laws for
the escape mass of the target material $m_{\rm tar,esc}$ and that of impactor
$m_{\rm imp,esc}$. Considering the mass balance, the mass originating from the
impactor that are deposited onto the target surface upon cratering is given as
$m_{\rm imp,acc}$ as $m_{\rm imp,acc}=m_{\rm imp}-m_{\rm imp,esc}$ ($m_{\rm
imp}$ is the impactor mass).

The motivation of this work is to provide scaling laws that would be useful
for the planetary community. These scaling laws would better predict the
fraction of the eroded target body materials which escape the gravity and the
fraction of impactor materials deposited on the target surface during a
cratering impact in the gravity-dominated regime. The scaling laws are
provided as a function of the impact velocity $v_{\rm imp}$, the impact angle
$\theta$, and the escape velocity of the target $v_{\rm esc}$ (see also Paper
1 for rocky bodies).

The numerical methods are presented in section 2. In section 3, the numerical
results of the ejecta that originate from the target (the target material) are
presented. Section 3 also presents the derivations of new scaling laws for the
escape mass of the target material. In section 4, the numerical results of the
ejecta originating from the impactor (the impactor material) are presented. A
new scaling law for the accretion mass of the impactor material onto the
target surface is also presented. In Section 5, a discussion of applications
using the newly derived scaling laws for planet formation is presented.
Finally, in section 6 a summary is presented.

##  2 Numerical methods

The fate of high-speed ejecta with exceeding the escape velocity of the icy
target bodies for a given impact velocity and impact angle were evaluated.
This was performed using 3D SPH simulations (Lucy, 1977; Monaghan, 1992) for
cratering impacts (see also Paper 1 for rocky bodies). The relevant equations
of state (EOSs) was a 5-Phase EOS for icy materials (Senft & Stewart, 2008).
This 5-Phase EOS included the effects of melting and vaporization for H2O. The
numerical code was identical as that used to study cratering impacts (Kurosawa
et al., 2019; Hyodo et al., 2019) and as that used to study a catastrophic
disruption of icy bodies (Hyodo & Charnoz, 2017). In this work, high-speed
ejecta whose ejection speed exceeds the escape velocity was considered.
Material strength were neglected in the simulations.

The icy impactor was idealized as a spherical projectile with a radius of
$R_{\rm imp}=10$ km. The icy target was idealized as a flat surface contained
in a half-sphere with a radius of typically greater than that of the impactor
(projectile) by a factor of 10-20. For the numerical resolution, the total
number of equal-mass SPH particles (for the target and for the impactor) used
in our simulations was $\sim 1.7\times 10^{7}$. As a resolution test, the
total number of particles was kept constant but kept changing the target size
and the number of particles to resolve impactors. $\sim 3.4\times 10^{4}$ and
$\sim 4.8\times 10^{3}$ SPH particles were used for the impactors which
correspond to $20$ and $10$ SPH particles per projectile radius (PPPR),
respectively. Convergence of the numerical results was confirmed by comparing
10 and 20 PPPR cases for the total mass of the high-speed ejecta of $v_{\rm
eje}>v_{\rm esc}$ (with $v_{\rm esc}=1-10$ km s-1).

Different sources of the escape mass, the target material (Section 3) and
impactor material (Section 4), were independently investigated. Ejecta
particles were defined as those existing above the surface of the target for
the analysis. Snapshots were utilized when the shock reached the boundary of
the target for the analysis.

Various impact velocities typically ranging from $6-62$ km s-1 (7 km s-1
interval) and at various impact angles from $15-90$ degrees ($15$ degrees
interval), where a 90-degree impact is a head-on/vertical impact, were
analyzed. In dimensionless form, model fidelity increases as the ratio of the
large planetary body size to that of the cratering impacts increases such that
the curvature of the target is of negligible importance (Genda et al., 2017).
The results can be related to any impactor size because all hydrodynamic
equations can be rewritten in a dimensionless form without gravity and
strength (Johnson & Melosh, 2012).

Figure 1: Schematic illustration of cratering impacts (left) and velocity
distribution of the ejecta (right). High- and low-speed ejecta are launched
closer to and farther from the impact point, respectively; that is, the escape
mass is launched from the vicinity of the impact point (left). Irregular-
shaped objects in the left panel indicate ejecta of target and impactor
materials (the total amount and their fractions depend on the impact
condition). The point-source scaling law is characterized by a unique slope
$-3\mu_{\rm HH11}$ in the velocity distribution (the black line in the right
panel; Equation 1), while direct SPH simulations of cratering impacts of icy
bodies demonstrated that high-speed ejecta would have a steeper slope for the
velocity distribution (the gray line in the right panel; see also Paper 1).
This discrepancy leads to an overestimation in the escape mass when using the
point-source scaling law, depending on the considered escape velocity of the
target $v_{\rm esc}$ (the same was reported for rocky bodies; Paper 1).

##  3 Escape of target material

In this section, the escape mass that originates from the target (the target
material) is discussed. First, the point-source scaling law (Section 3.1) is
explained. In Section 3.2, the ejection velocity distributions for different
impact angles are investigated. In Section 3.3, a new scaling law for the
escape mass of the target materials is derived. Finally, the newly derived
scaling law is compared to the classical point-source scaling law (Section
3.4).

###  3.1 The point-source scaling law

Cratering impacts lead to an ejection of the target material from the surface.
Pioneering studies (Housen et al., 1983; Holsapple & Schmidt, 1987; Holsapple,
1993; Holsapple & Housen, 2007; Housen & Holsapple, 2011) predicted $m(>v_{\rm
eje,tar})$ as a function of the impact velocity $v_{\rm imp}$ and impact angle
$\theta$ assuming that the size of the impact crater is assumed to be much
larger than that of the impactor (hereafter $M_{\rm HH11,esc,tar}$)111We note
that the point-source scaling law implicitly assumed negligible contribution
from the impactor material.. The point-source scaling law can potentially
predict the escape mass by using the relation, $v_{\rm eje,tar}=v_{\rm esc}$,
where $v_{\rm esc}$ is the escape velocity of the target, as follows
(hereafter HH11):

$\frac{M_{\rm HH11,esc,tar}(>v_{\rm esc})}{m_{\rm imp}}=C_{\rm
HH11}\left(\frac{v_{\rm esc}}{v_{\rm imp}\sin\theta}\right)^{-3\mu_{\rm
HH11}},$ (1)

where $C_{\rm HH11}=(3k/4\pi)C_{0}^{3\mu_{\rm HH11}}$ and $\mu_{\rm HH11}$ are
constants. From the experiments, $C_{\rm 0}=1.5$ and $\mu_{\rm HH11}=0.55$ for
nonporous icy materials and rocky materials, respectively (e.g., Holsapple &
Housen, 2007). $k=0.2$ for ice and $k=0.3$ for rock, respectively (Gault et
al., 1963), respectively. The densities of the impactor and target are assumed
to be identical.

###  3.2 Velocity distribution of the ejecta

Figure 2: Cumulative ejecta mass ($>v_{\rm eje}$; target material) that
exceeds a given ejection velocity $v_{\rm eje}$ for a different impact angle
$\theta$. The black, blue, green, and red lines represent the results of SPH
simulations for $v_{\rm imp}=6,13,34$, and $62$ km s-1, respectively. The
black line represents the point-source scaling law ($-3\mu_{\rm HH11}$;
Equation 1). Numerical results do not converge for small velocities where the
point-source scaling assumption is valid, but this does not affect the
conclusion of the current study (since our primary interest is the high-speed
ejecta of $v_{\rm eje}>v_{\rm esc}$).

The impact shock propagates through the interior of the target after the
impactor makes contact with the target surface. The crater forms as the target
material is sheared, moving upward and outward along the bowl-shaped crater
edge (Figure 1; see also Housen & Holsapple, 2011). The speed of the ejecta
$v_{\rm eje}$ depends on its launch position (e.g., Piekutowski, 1980), that
is, high- and low-speed ejecta are launched closer to and farther from the
impact point, respectively (Figure 1).

Figure 2 shows the cumulative mass of ejecta originating from the target as a
function of ejection velocity $v_{\rm eje}$ for different impact velocities
and impact angles. The ejection velocity is scaled by $v_{\rm imp}\sin\theta$.
Note that the data for the ejection velocity with $v_{\rm eje}\lesssim
0.1v_{\rm imp}\sin\theta$ failed to converge during the performed simulation
time (i.e., more numerical time is needed for the convergence of smaller
$v_{\rm eje}$ because the lower-speed ejecta is launched farther from the
impact point222It was noted that the point-source scaling law would not be
appropriate for a very small ejection velocity because of the effects of
material strength and gravity (e.g., Housen & Holsapple, 2011), although such
a regime does not pertain to this study.). The point-source scaling law
(Equation 1) is plotted as a black line. The results indicate that the
distribution is uniquely scaled by $v_{\rm imp}\sin\theta$, which is in
accordance with the basic principles evident in the point-source assumption
(Holsapple & Housen, 2007).

The slope of the ejection velocity distribution is characterized by
$-3\mu_{\rm HH11}$, a unique constant value, in the case of the point-source
scaling law (Equation 1). This agrees with the results obtained from direct
impact simulations only for the low-speed ejecta (Figure 2). This indicates
that the point-source scaling law is valid only for a limited range of
ejection velocity distributions.

The ejecta velocity distribution has a steeper slope for high-speed ejecta
than that predicted by the point-source scaling law (see also Figure 1). The
ejection velocity distribution for such high-speed ejecta has a unique power-
law dependence depending on different impact angles (Figure 2). This result
indicates that the point-source scaling law overestimates the mass of the
high-speed ejecta with a launch point near the impact point. This was also
reported in the case of rocky bodies (Paper 1). In the next subsection,
updated scaling laws are derived for the high-speed ejecta where the point-
source scaling law is not applicable.

###  3.3 Updated scaling law for high-speed ejecta

Figure 3: The exponent $\mu_{\rm HG21,tar}(\theta)$ (left) and the coefficient
$C_{\rm HG21,tar}(\theta)$ (right) for the new scaling law of the escape mass
of the target material (Equations 3.3 and 3.3) as a function of the impact
angle. The points represent the results of SPH simulations and solid curves
are defined by fitted quadratic and cubic functions of the impact angle for
$\mu_{\rm HG21,tar}(\theta)$ and $C_{\rm HG21,tar}(\theta)$, respectively.
Figure 4: Escape mass of target material (in units of $m_{\rm imp}$) as a
function of $v_{\rm imp}$ for different $\theta$. Points are the results of
the SPH simulations. Solid and dashed lines represent the updated scaling law
(Equation 3.3; $\min\left\\{M^{*}_{\rm HG21,esc,tar},M_{\rm
HH11,esc,tar}\right\\}$). Solid lines represent the cases where the same
function as in Equation 3.3 ($M^{*}_{\rm HG21,esc,tar}$) are used. Dashed
lines depict cases in which the function is the same as in Equation 1 ($M_{\rm
HH11,esc,tar}$). The black, blue, green, and red lines represent cases where
$v_{\rm esc}=1,3,6$, and $10$ km s-1, respectively.

Following the methodology of Paper 1 and from the arguments in the previous
subsection, the escape mass of the target material $M^{*}_{\rm HG21,esc,tar}$
for the high-speed ejecta regime ($\gtrsim 0.1v_{\rm imp}\sin\theta$) is
uniquely expressed by a power-law as a function of impact velocity and impact
angle as follows:

$\displaystyle\frac{M^{*}_{\rm HG21,esc,tar}(>v_{\rm esc})}{m_{\rm imp}}=$
$\displaystyle C_{\rm HG21,tar}(\theta)$ $\displaystyle\left(\frac{v_{\rm
esc}}{v_{\rm imp}\sin\theta}\right)^{-3\mu_{\rm HG21,tar}(\theta)},$ (2)

where $C_{\rm HG21,tar}(\theta)$ and $\mu_{\rm HG21,tar}(\theta)$ are a new
coefficient and a new exponent that depends on the impact angle, respectively.

One may obtain the exact $C_{\rm HG21,tar}(\theta)$ and $\mu_{\rm
HG21,tar}(\theta)$ at varying impact angles using the numerical results
(indicated by blue points in Figure 3). Subsequently, the approximated
$\mu_{\rm HG21,tar}(\theta)$ and $C_{\rm HG21,tar}(\theta)$ (the black lines
in Figure 3) were derived using the quadratic and cubic functions of the
impact angle (the black lines in Figure 3) as

$\displaystyle\mu_{\rm HG21,tar}(\theta)=a_{\rm tar}\theta^{2}+b_{\rm
tar}\theta+c_{\rm tar}$ (3) $\displaystyle C_{\rm
HG21,tar}(\theta)=\exp\left(d_{\rm tar}\theta^{3}+e_{\rm tar}\theta^{2}+f_{\rm
tar}\theta+g_{\rm tar}\right),$ (4)

where $a_{\rm tar}$, $b_{\rm tar}$, $c_{\rm tar}$, $d_{\rm tar}$, $e_{\rm
tar}$, $f_{\rm tar}$, and $g_{\rm tar}$ are the fitted parameters,
respectively (Table 1).

Figure 4 shows the escape mass of the target material (target mass with an
ejection velocity greater than the escape velocity) as a function of the
impact velocity. Points were obtained from the numerical simulations.

As Paper 1 established that for rocky bodies, the original point-source
scaling law for icy bodies (HH11; Equation 1) is valid only at sufficiently
large distances from the impact point (i.e., for a sufficiently small ejection
velocity). The coefficient and exponent deviate from those of HH11 (compare
the solid black line and the other lines in Figure 2) for large ejection
velocities (i.e., small distances to the impact point).

As observed in Figure 2 (see also Figure 1), HH11 for icy bodies overestimates
the escape mass of the target material, especially for a significant ejection
velocity. Therefore, the new scaling law for the escape mass of the target
material that combines with HH11 is given by $\min\left\\{M^{*}_{\rm
HG21,esc,tar},M_{\rm HH11,esc,tar}\right\\}$ and is written as follows
(hereafter, HG21):

$\displaystyle\frac{M_{\rm HG21,esc,tar}(>v_{\rm esc})}{m_{\rm imp}}=$
$\displaystyle\min\left\\{C_{\rm HG21,tar}(\theta)\left(\frac{v_{\rm
esc}}{v_{\rm imp}\sin(\theta)}\right)^{-3\mu_{\rm HG21,tar}(\theta)},\right.$
$\displaystyle\hskip 71.13188pt\left.C_{\rm HH11}\left(\frac{v_{\rm
esc}}{v_{\rm imp}\sin(\theta)}\right)^{-3\mu_{\rm HH11}}\right\\}.$ (5)

The new scaling law (Equation 3.3) is plotted by solid and dashed lines in
Figure 4. A close match between the updated scaling law and impact simulations
(points) for $v_{\rm imp}=6-62$ km s-1 and $\theta=15-90$ degrees for $v_{\rm
esc}=1-10$ km s-1 was observed.

$a_{\rm tar}$ | $b_{\rm tar}$ | $c_{\rm tar}$ |
---|---|---|---
$-3.76\times 10^{-5}$ | $1.45\times 10^{-2}$ | $3.27\times 10^{-1}$ |
$d_{\rm tar}$ | $e_{\rm tar}$ | $f_{\rm tar}$ | $g_{\rm tar}$
$4.13\times 10^{-5}$ | $-6.85\times 10^{-3}$ | $2.53\times 10^{-1}$ | $-4.55$
$a_{\rm imp}$ | $b_{\rm imp}$ | $c_{\rm imp}$ |
$-3.61\times 10^{-4}$ | $4.48\times 10^{-2}$ | $-5.99\times 10^{-1}$ |
$d_{\rm imp}$ | $e_{\rm imp}$ | $f_{\rm imp}$ | $g_{\rm imp}$
$3.97\times 10^{-5}$ | $-4.74\times 10^{-3}$ | $4.99\times 10^{-2}$ | $2.53\times 10^{-1}$
Table 1: Parameters of the fitted polynomial and exponential functions for the
escape masses of the target material (Equation 3.3) and the impactor material
(Equation 4), respectively. Those for rocky bodies are shown in Paper 1.

###  3.4 Comparison to the point-source scaling law

The point-source scaling law (Equation 1) has been widely used. As
demonstrated in the case of the rocky bodies (Paper 1), the point-source
scaling law for the escape mass of the target material of icy bodies also
agrees with the numerical results when $v_{\rm imp}\gg v_{\rm esc}$, whereas
it overestimates the escape mass when $v_{\rm imp}\gtrsim v_{\rm esc}$.

Figure 5 shows the overestimation of the point-source scaling law (HH11;
Equation 1) when compared to the results of the direct impact simulations
(HG21; Equation 3.3) as a function of the impact velocity. Except for impact
angles of $15$ and $30$ degrees, HH11 overestimated the escape mass for
$v_{\rm imp}\lesssim 9v_{\rm esc}$, which exponentially increased as the
impact velocity decreased. As found in the case of the rocky bodies (Paper 1),
the difference became more significant as the impact angle tended toward a
head-on collision. HH11 overestimated the escape mass of the target material
by approximately $\sim 100$ times when $v_{\rm imp}/v_{\rm esc}\sim 1$ at
vertical impact ($\theta=90$ degrees). For the $\theta$-averaged escape mass
over the $\sin 2\theta$ distribution (Equations 10–11), HH11 overestimated by
a factor of $\sim 4$ larger than HG21 when $v_{\rm imp}/v_{\rm esc}\sim 1$.

Figure 5: The ratio of the point-source scaling law (HH11; Equation 1) to the
newly derived scaling law (HG21; Equation 3.3) for the escape target mass of
icy bodies. The $x$-axis represents the impact velocity in a unit of escape
velocity. The red, orange, brown, green, cyan, and blue lines are the cases of
the impact angles of $\theta=$90, 75, 60, 45, 30, and 15 degrees,
respectively. The solid black line represents the $\theta$-averaged case over
$\sin 2\theta$. Those for rocky bodies are shown in Paper 1.

##  4 Accretion of impactor material

In this section, a discussion regarding the accretion mass of the impactor
material onto the target surface upon cratering impacts is presented. As a
consequence of a cratering impact, a fraction of the impactor material (mass
of $m_{\rm imp,esc}$) escapes as its speed exceeds the escape velocity of the
target, and the remaining material accretes onto the target surface.

The accretion mass originating from the impactor, the accretion mass of the
impactor material $m_{\rm imp,acc}$, is given by the mass balance as $m_{\rm
imp,acc}=m_{\rm imp}-m_{\rm imp,esc}$. By evaluating $m_{\rm imp,esc}$ using
SPH simulations and by considering the mass balance, one obtains $m_{\rm
imp,acc}$. A scaling law of the accretion mass of the impactor material onto
the target surface is derived, as shown below.

Figure 6: The exponent $\mu_{\rm HG21,imp}(\theta)$ (left) and the coefficient
$C_{\rm HG21,imp}(\theta)$ (right) for the new scaling law (Equation 4) that
predicts the escape mass of the impactor material as a function of the impact
angle. The points are the results of SPH simulations and solid curves are the
fitted quadratic and cubic functions of the impact angle for $\mu_{\rm
HG21,imp}(\theta)$ and $C_{\rm HG21,imp}(\theta)$, respectively. Note that
$\mu_{\rm HG21,imp}(\theta)=0$ and $C_{\rm HG21,imp}(\theta)=1$ for
$\theta<15$ degrees. Figure 7: Escape mass of the impactor material (in units
of $m_{\rm imp}$) as a function of $v_{\rm imp}$ for different $\theta$. The
points are the results of the SPH simulations. Triangles are used to indicate
cases for which numerical simulations do not converge. Circles are employed
when numerical simulations converge. Solid lines represent the new scaling law
(Equation 4) derived in this study. The black, blue, green, and red lines
represent the cases of $v_{\rm esc}=1,3,6$, and $10$ km s-1, respectively.

The derived escape mass of the impactor material is defined as $M_{\rm
HG21,esc,imp}(>v_{\rm esc})$. Using the same arguments as Paper 1, we assume
that $M_{\rm HG21,esc,imp}(>v_{\rm esc})$ follows a power-law function as
follows:

$\displaystyle\frac{M_{\rm HG21,esc,imp}(>v_{\rm esc})}{m_{\rm imp}}=$
$\displaystyle C_{\rm HG21,imp}(\theta)$ $\displaystyle\left(\frac{v_{\rm
esc}}{v_{\rm imp}\sin(\theta)}\right)^{-3\mu_{\rm HG21,imp}(\theta)},$ (6)

where $C_{\rm HG21,imp}(\theta)$ and $\mu_{\rm HG21,imp}(\theta)$ are
coefficients and exponents that depend on the impact angle. Equation 4 is fit
to the converged numerical results mainly at $v_{\rm esc}=10$ km s-1 to obtain
the coefficient and exponent at varying impact angles.

Figure 6 shows the coefficients and exponents obtained from the numerical
simulations. As in the case of the rocky bodies (Paper 1), we derived
$\mu_{\rm HG21,imp}(\theta)$ and $C_{\rm HG21,imp}(\theta)$ (lines in Figure
6), respectively, as:

$\displaystyle\mu_{\rm HG21,imp}(\theta)=a_{\rm imp}\theta^{2}+b_{\rm
imp}\theta+c_{\rm imp}$ (7) $\displaystyle C_{\rm
HG21,imp}(\theta)=\exp\left(d_{\rm imp}\theta^{3}+e_{\rm imp}\theta^{2}+f_{\rm
imp}\theta+g_{\rm imp}\right),$ (8)

where $a_{\rm imp}$, $b_{\rm imp}$, $c_{\rm imp}$, $d_{\rm imp}$, $e_{\rm
imp}$, $f_{\rm imp}$, and $g_{\rm imp}$ are the fitted parameters,
respectively (Table 1). Note that $\mu_{\rm HG21,imp}(\theta)=0$ and $C_{\rm
HG21,imp}(\theta)=1$ for $\theta<15$ degrees.

Figure 7 shows the escape mass of an impactor material with a speed greater
than $v_{\rm esc}$. Points represent data obtained from the SPH simulations.
Some of the simulations for high impact velocities and impact angles did not
converge within a reasonable computational time (plotted by triangles), while
converged parameters are plotted by circles.

The newly derived scaling law (Equation 4) is plotted in Figure 7 by solid
lines. The law generally agrees with the numerical results, especially for
$\theta=30-75$ degrees. Note that impacts with $\theta=0$ and $90$ degrees do
not statistically occur in planet formation, because the impact angle
distribution is given by $\sin 2\theta$ (Shoemaker, 1962).

The accretion mass of the impactor material onto the target surface upon
cratering impacts is given as:

$\displaystyle\frac{M_{\rm HG21,acc,imp}(<v_{\rm esc})}{m_{\rm imp}}=$
$\displaystyle 1-C_{\rm HG21,imp}(\theta)$ $\displaystyle\left(\frac{v_{\rm
esc}}{v_{\rm imp}\sin\theta}\right)^{-3\mu_{\rm HG21,imp}(\theta)}.$ (9)

Discarding factors of the order of unity difference, we found that our
numerical results are compatible with Pierazzo & Chyba (2002) (case of icy
impactors on Europa’s icy surface with small porosity) and Yue et al. (2013)
(case of rocky impactors on the rocky surface of the Moon). Comparing the
results of Ong et al. (2010) (case of icy impactors on the rocky surface of
the Moon) with our results, we found a notable difference in the impactor’s
accretion mass. Our numerical simulations are limited to the target and
impactor having the same composition. We believe further research needs to be
carried out to study the impact of various factors, such as materials of
target and impactor, numerical methods used on the numerical results (e.g.,
2D/3D, different EOSs, Eulerian/Lagrangian methods).

##  5 Applications for planet formation

###  5.1 $\theta$-averaged escape and accretion masses

Cratering impacts of small bodies on a large planetary body occur much more
frequently than collisions between similar-sized bodies. The statistical
distribution of the impact angle follows $\sin 2\theta$ with a peak of $45$
degrees (Shoemaker, 1962).

Using the newly derived scaling laws for the escape mass of icy target
material (Equation 3.3) and the accretion mass of icy impactor material
(Equation 4), one can derive the $\theta$-averaged scaling laws weighted over
the $\sin 2\theta$ distribution (see Paper 1 for rocky bodies). These are as
follows.

The $\theta$-averaged escape mass originating from the icy target as a
function of impact velocity was derived as

$\displaystyle\left<\frac{M_{\rm ice,esc,tar}(>v_{\rm esc})}{m_{\rm
imp}}\right>_{\theta}=$ $\displaystyle 0.013\left(\frac{v_{\rm imp}}{v_{\rm
esc}}\right)^{2.29}$ (10) $\displaystyle\left[\mbox{Icy bodies;\ }v_{\rm
imp}\lesssim 9v_{\rm esc}\right]$ $\displaystyle\left<\frac{M_{\rm
ice,esc,tar}(>v_{\rm esc})}{m_{\rm imp}}\right>_{\theta}=$ $\displaystyle
0.051\left(\frac{v_{\rm imp}}{v_{\rm esc}}\right)^{1.65}.$ (11)
$\displaystyle\left[\mbox{Icy bodies;\ }v_{\rm imp}\gtrsim 9v_{\rm
esc}\right]$

That for rocky bodies was (derived in Paper 1):

$\displaystyle\left<\frac{M_{\rm rock,esc,tar}(>v_{\rm esc})}{m_{\rm
imp}}\right>_{\theta}$ $\displaystyle=0.02\left(\frac{v_{\rm imp}}{v_{\rm
esc}}\right)^{2.2}$ (12) $\displaystyle\left[\mbox{Rocky bodies;\ }v_{\rm
imp}\lesssim 12v_{\rm esc}\right]$ $\displaystyle\left<\frac{M_{\rm
rock,esc,tar}(>v_{\rm esc})}{m_{\rm imp}}\right>_{\theta}$
$\displaystyle=0.076\left(\frac{v_{\rm imp}}{v_{\rm esc}}\right)^{1.65}.$ (13)
$\displaystyle\left[\mbox{Rocky bodies;\ }v_{\rm imp}\gtrsim 12v_{\rm
esc}\right]$

In the above equations, one uses the smaller one of the two (Equations 10 – 11
or Equations 12 – 13). Equations 11 and 13 are obtained by averaging the
point-source scaling laws333From the point-source assumptions (Equation 1),
the $\theta$-averaged power-law dependences are $\propto(v_{\rm imp}/v_{\rm
esc})^{3\mu_{\rm HH11}}=(v_{\rm imp}/v_{\rm esc})^{1.65}$ for both rocky and
icy bodies. The $\theta$-averaged coefficients (see Equation 1) are $<C_{\rm
HH11}>_{\rm\theta}=C_{\rm HH11}\int_{0}^{\pi/2}(\sin\theta)^{3\mu_{\rm
HH11}}\sin 2\theta d\theta\simeq 0.051$ for icy bodies and $\simeq 0.076$ for
rocky bodies, respectively..

The $\theta$-averaged accretion mass originating from the icy impactor as a
function of impact velocity is

$\displaystyle\left<\frac{M_{\rm ice,acc,imp}(<v_{\rm esc})}{m_{\rm
imp}}\right>_{\theta}=0.88-$ $\displaystyle 0.024\left(\frac{v_{\rm
imp}}{v_{\rm esc}}\right)^{1.66}.$ $\displaystyle\left[\mbox{Icy
bodies}\right]$ (14)

That for rocky bodies was (derived in Paper 1):

$\displaystyle\left<\frac{M_{\rm rock,acc,imp}(<v_{\rm esc})}{m_{\rm
imp}}\right>_{\theta}=0.85$ $\displaystyle-0.071\left(\frac{v_{\rm
imp}}{v_{\rm esc}}\right)^{0.88}.$ $\displaystyle\left[\mbox{Rocky
bodies}\right]$ (15)

These $\theta$-averaged scaling laws are shown in Figure 8. In Section 5.3,
the typical outcome of a cratering impact at a variety of planetary bodies in
the solar system (see Table 2) is discussed.

###  5.2 Critical target size where strength needs to be included

Small impactors are attracted by the gravity of large target bodies (i.e.,
cratering impacts) and impact velocity becomes greater than the escape
velocity of the target ($v_{\rm imp}\gtrsim v_{\rm esc}$). Our numerical
simulations neglected material strength. In this subsection, using a simple
but a well-defined analytical argument, we estimate the critical target size
below which material strength affects the determination of ejecta having
velocity greater than the escape velocity of the target (i.e., ejecta of
$v_{\rm eje}\geq v_{\rm esc}$).

Using the Rankine-Hugoniot relations, the shock pressure $P_{\rm shock}$
during the cratering impact is given by

$\displaystyle P_{\rm shock}=\rho_{\rm 0}\left(C_{\rm 0}+su_{\rm
p}\right)u_{\rm p},$ (16)

where $\rho_{\rm 0}$ is material density before the shock compression, $C_{\rm
0}$ is bulk sound velocity ($C_{\rm 0}\simeq 1000-3000$ km s-1 for rocky and
icy materials, respectively; Melosh, 1989), $u_{\rm p}$ is a particle velocity
($u_{\rm p}$ is related to shock velocity as $u_{\rm s}=C_{\rm 0}+su_{\rm
p}$), and $s\equiv du_{\rm s}/du_{\rm p}$ is a constant ($s\simeq 1-2$;
Melosh, 1989).

For $P_{\rm shock}>P_{\rm HEL}$ where $P_{\rm HEL}$ is the Hugoniot elastic
limit (HEL), a critical pressure above which fluid-like response is expected,
the ejection process is expected to be unaffected by material strength. In
this paper, the minimum ejection velocity of our primary interest is equal to
escape velocity (i.e., $v_{\rm eje}=v_{\rm esc}$ and a larger ejection
velocity accompanies a larger shock pressure). $u_{\rm eje}$ should be larger
than $u_{\rm p}$ due to acceleration during adiabatic expansion. Approximately
$u_{\rm eje}\sim 2u_{\rm p}$ (Melosh, 1989), which leads to $u_{\rm p}=v_{\rm
esc}/2=\sqrt{2G\pi\rho_{\rm 0}/3}R_{\rm tar}$ where $R_{\rm tar}$ is the
target radius.

For $C_{\rm 0}>su_{\rm p}$ (where $u_{\rm p}=v_{\rm esc}/2)$, we can
approximate $C_{\rm 0}+su_{\rm p}\simeq C_{\rm 0}$ to determine the minimum
shock pressure experienced by ejecta of $v_{\rm eje}\geq v_{\rm esc}$ during a
cratering impact. This condition is rewritten as

$\displaystyle R_{\rm tar}$ $\displaystyle<C_{\rm 0}/(s\sqrt{2G\pi\rho_{\rm
0}/3})$ (17) $\displaystyle<2300{\rm\,km}\times$ (18) $\displaystyle\hskip
30.00005pt\left(\frac{C_{\rm
0}}{3000{\rm\,m\,s^{-1}}}\right)\left(\frac{s}{2}\right)^{-1}\left(\frac{\rho_{\rm
0}}{3000{\rm\,kg\,m^{-3}}}\right)^{-1/2}.$ (19)

For a small sized target (Equation 17), we can determine the shock pressure
for above condition by

$\displaystyle P_{\rm shock}$ $\displaystyle=\rho_{\rm 0}\left(C_{\rm
0}+su_{\rm p}\right)u_{\rm p}\simeq C_{\rm 0}\rho_{\rm 0}u_{\rm p}$ (20)
$\displaystyle=\sqrt{2\pi G/3}C_{\rm 0}\rho_{\rm 0}^{3/2}R_{\rm tar}$ (21)
$\displaystyle\simeq 600{\rm\,MPa}\times$ (22) $\displaystyle\hskip
15.00002pt\left(\frac{C_{\rm
0}}{3000{\rm\,m\,s^{-1}}}\right)\left(\frac{\rho_{\rm
0}}{3000{\rm\,kg\,m^{-3}}}\right)^{3/2}\left(\frac{R_{\rm
tar}}{100{\rm\,km}}\right).$ (23)

Typically $P_{\rm HEL}\simeq 2-5{\rm\,GPa}$ for rocky bodies (Sekine et al.,
2008) and $P_{\rm HEL}\simeq 150-300{\rm\,MPa}$ for icy bodies (Stewart &
Ahrens, 2003). The critical size for $P_{\rm shock}=P_{\rm HEL}$ is readily
obtained from Equation 20. For icy bodies, $P_{\rm shock}>P_{\rm HEL}$ (i.e.,
strength is negligible) when $R_{\rm tar}\gtrsim 600$ km (with $P_{\rm
HEL}\sim 300$ MPa, $C_{\rm 0}\sim 1300$ m s-1, and $\rho_{\rm 0}\sim 1000$ kg
m-3). For rocky bodies, $P_{\rm shock}>P_{\rm HEL}$ when $R_{\rm tar}\gtrsim
1000$ km (with $P_{\rm HEL}\sim 5$ GPa, $C_{\rm 0}\sim 2600$ m s-1, and
$\rho_{\rm 0}\sim 2900$ kg m-3).

We have ignored the effects of material strength in our study. However,
further research needs to be carried out on the effects of material strength
on the erosion mass and accretion mass.

###  5.3 Cratering impacts on rocky and icy bodies

Name |  $v_{\rm imp}$ [km s-1] | Description
---|---|---
(rocky/icy) | ($v_{\rm imp}$ [$v_{\rm esc}$]) | (Reference)
Mercury |  $\sim 36$ km s-1 | Instability phase
(rocky) | ($\sim 8.5v_{\rm esc}$; $v_{\rm esc}\simeq 4.25$ km s-1) | (e.g., Hyodo et al., 2021)
Venus |  $\sim 26$ km s-1 | Instability phase
(rocky) | ($\sim 2.5v_{\rm esc}$; $v_{\rm esc}\simeq 10.4$ km s-1) | (e.g., Mojzsis et al., 2019)
Earth |  $\sim 22$ km s-1 | Instability phase
(rocky) | ($\sim 2.0v_{\rm esc}$; $v_{\rm esc}\simeq 11.2$ km s-1) | (e.g., Mojzsis et al., 2019)
Moon |  $\sim 19$ km s-1 | Instability phase
(rocky) | ($\sim 8.0v_{\rm esc}$; $v_{\rm esc}\simeq 2.4$ km s-1) | (e.g., Brasser et al., 2020)
Mars |  $\sim 13$ km s-1 | Instability phase
(rocky) | ($\sim 2.6v_{\rm esc}$; $v_{\rm esc}\simeq 5.0$ km s-1) | (e.g., Mojzsis et al., 2019)
Ceres |  $\sim 5$ km s-1 | Today’s population
(rocky) | ($\sim 9.8v_{\rm esc}$; $v_{\rm esc}\simeq 0.5$ km s-1) | (e.g., Bottke et al., 1994)
Vesta |  $\sim 5$ km s-1 | Today’s population
(rocky) | ($\sim 13.9v_{\rm esc}$; $v_{\rm esc}\simeq 0.35$ km s-1) | (e.g., Bottke et al., 1994)
Rhea |  $\sim 15$ km s-1 | Instability phase
(icy) | ($\sim 23.6v_{\rm esc}$; $v_{\rm esc}\simeq 0.64$ km s-1) | (e.g., Wong et al., 2020)
Europa |  $\sim 24$ km s-1 | Instability phase
(icy) | ($\sim 12.0v_{\rm esc}$; $v_{\rm esc}\simeq 2.0$ km s-1) | (e.g., Wong et al., 2020)
Pluto |  $\sim 2$ km s-1 | Today’s population
(icy) | ($\sim 1.7v_{\rm esc}$; $v_{\rm esc}\simeq 1.2$ km s-1) | (e.g., Dell’Oro et al., 2013)
Table 2: Examples of typical impact velocities for different rocky and icy
bodies in different contexts and epochs. The strong gravity of giant planets
attracts small bodies as high-speed impactors onto their moons. During the
instability phase, the solar system likely experienced putative cataclysm
leading to a rise in impact flux that characterized the surface morphology and
composition of solar system bodies. Figure 8: $\theta$-averaged scaling laws
for rocky (red lines) and icy (blue lines) bodies as a function of $v_{\rm
imp}/v_{\rm esc}$. Left: the escape mass of the target material $M_{\rm
esc,tar}/m_{\rm imp}$ (i.e., erosion mass of target; Equations 10–13). Middle:
the accretion mass of the impactor material onto the target surface $M_{\rm
acc,imp}/m_{\rm imp}$ (Equations 5.1–5.1). Right: total escape mass ($M_{\rm
esc,tot}=M_{\rm esc,tar}+M_{\rm esc,imp}$ in a unit of $m_{\rm imp}$; $M_{\rm
esc,tot}>1$ indicates a net erosion and vice versa for a net accretion).
Examples of typical outcomes of cratering impacts on different solar system
bodies under different contexts are indicated by the labels in each panels
(see Table 2 for the assumed impact conditions).

In our solar system, a variety of orbits exist for planetary bodies including
small bodies (asteroids and comets). For example, the orbital semi-major axis
ranges from $a\sim 0.4$ to $40$ au for Mercury and Pluto, respectively.
Cratering impacts inevitably occur in different contexts and epochs for a
typical impact velocity range of $v_{\rm imp}\sim 1-100v_{\rm esc}$ (Table 2).
High-speed cratering impacts would occur during the instability phase (Mojzsis
et al., 2019; Brasser et al., 2020); these include the ”Nice model” (e.g.,
Gomes et al., 2005; Tsiganis et al., 2005), the ”Grand-tack” hypothesis (e.g.,
Walsh et al., 2011), and the ”early instability” scenario (e.g., Clement et
al., 2018). Giant impacts, such as those that formed the Moon (Bottke et al.,
2015) and/or the Martian moons (Hyodo & Genda, 2018), would distribute impact
debris throughout the inner solar system, and high-speed collisions between
the debris and asteroids (and/or planets) may correspondingly occur ($v_{\rm
imp}\gg 5$ km s-1).

Secular collisional evolution over billions of years and/or a rise of impact
flux during the instability phase are expected to characterize the surface
morphology and composition of solar system bodies. Figure 8 illustrates
examples of typical outcomes of cratering impacts that would potentially occur
on different planetary bodies.

Comparing rocky and icy bodies, rocky target bodies are, in general, more
eroded than icy target bodies (left panel in Figure 8), while less icy
impactor’s material is accreted on icy target surface than the case of rocky
bodies for $v_{\rm imp}\gtrsim 4v_{\rm esc}$ (middle panel in Figure 8).

Cratering impacts on rocky terrestrial planets, except Mercury, generally lead
to a net mass increase as their $v_{\rm imp}/v_{\rm esc}$ ratios are typically
$v_{\rm imp}/v_{\rm esc}\lesssim 5$ even in an instability phase. Mercury is
the innermost planet in the solar system, leading to the highest orbital
Keplerian velocity. Together with its small size, Mercury, cratering impacts
on average lead to a net erosion (see also Hyodo et al., 2021). The typical
impact velocity on the Moon in an instability phase is similar to that on
Earth, whereas its mass is $\sim 0.01$ Earth mass. This leads to $v_{\rm
imp}/v_{\rm esc}\sim 8$, and cratering impacts are generally erosive for the
Moon.

The largest asteroids are Ceres ($\sim 500$ km in radius) and Vesta ($\sim
260$ km in radius). The current average impact velocity among asteroids is
$\sim 5$ km s-1 (e.g., Bottke et al., 1994). Therefore, $v_{\rm imp}/v_{\rm
esc}\sim 10$ and $v_{\rm imp}/v_{\rm esc}\sim 14$ for Ceres and Vesta,
respectively. Such a high-velocity impact is expected to be erosive. However,
we note that our scaling laws might not be applicable to these asteroids
because we have ignored material strength in our study (Section 5.2).

Moons around giant planets and trans-Neptunian objects, such as Rhea, Europa,
and Pluto, are icy bodies. The strong gravity of giant planets attracts small
bodies as high-speed impactors (e.g., $v_{\rm imp}/v_{\rm esc}\sim 23.6$ for
the largest Saturn’s regular moon, Rhea; Wong et al., 2020) and erosive
cratering impacts on the moons occur. The typical impact velocity on Pluto is
$v_{\rm imp}/v_{\rm esc}\sim 1.7$ (Dell’Oro et al., 2013), which leads to a
net accretion of cratering impacts.

##  6 Summary

Cratering impacts inevitably occur during planet formation in different
contexts and epochs. Cratering impacts erode a portion of the target surface
by ejecting materials with a speed larger than its escape velocity while
implanting a fraction of the impactor’s material onto the target surface
(e.g., Housen et al., 1983; Housen & Holsapple, 2011; Hyodo & Genda, 2020, and
references therein).

As was done in the companion study for rocky bodies (Paper 1), SPH simulations
of cratering impacts on icy bodies were performed. The main findings are as
follows:

  * –

Updated scaling laws were derived that predict the escape mass of the icy
target material whose speed exceeds the escape velocity of a target (i.e., the
erosion mass of the target; Equation 3.3). It was demonstrated that a widely
known classical point-source scaling law largely overestimated the escape mass
of icy bodies when $v_{\rm imp}\lesssim 9v_{\rm esc}$ (the same occurs when
$v_{\rm imp}\lesssim 12v_{\rm esc}$ for rocky bodies; Paper 1).

  * –

Correspondingly, previous studies that have adopted the point-source scaling
laws for rocky and icy bodies in an inappropriate $v_{\rm imp}/v_{\rm esc}$
regime should be redone. The conclusions derived in such studies would
overestimate the erosion mass upon cratering impacts.

  * –

We also derived a scaling law that predicts the accretion mass of an
impactor’s material upon cratering impacts of icy bodies (Equation 4; see
Paper 1 for rocky bodies). The $\theta$-averaged scaling law weighted over
$\sin 2\theta$ (Equation 5.1) indicates that no impactor material accretes on
the target surface for $v_{\rm imp}/v_{\rm esc}\gtrsim 9$.

  * –

Statistically speaking (the $\theta$-averaged over $\sin 2\theta$), rocky
bodies are generally more easily eroded than icy bodies (Equations 10–13). Net
erosion for icy and rocky bodies takes place for $v_{\rm imp}/v_{\rm
esc}\gtrsim 5$ (Figure 8).

Diverse cratering impacts have critical roles in the mass gain/loss and in
characterizing the surface composition of planets and asteroids (Figure 8).
Further studies on cratering impacts are warranted for diverse planetary
bodies.

R.H. acknowledges the financial support of JSPS Grants-in-Aid (JP17J01269,
18K13600). H.G. acknowledges the financial support of MEXT KAKENHI Grant
(JP17H06457), and JSPS Kakenhi Grant (19H00726). We would like to thank
Editage (www.editage.com) for English language editing.

## References

  * Bottke et al. (1994) Bottke, W., Nolan, M., Greenberg, R., & Kolvoord, R. 1994, Icarus, 107, 255, doi: 10.1006/icar.1994.1021
  * Bottke et al. (2015) Bottke, W., Vokrouhlický, D., Marchi, S., et al. 2015, Science, 348, 321, doi: 10.1126/science.aaa0602
  * Brasser et al. (2020) Brasser, R., Werner, S., & Mojzsis, S. 2020, Icarus, 338, 113514, doi: 10.1016/j.icarus.2019.113514
  * Clement et al. (2018) Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icarus, 311, 340, doi: 10.1016/j.icarus.2018.04.008
  * Dell’Oro et al. (2013) Dell’Oro, A., Campo Bagatin, A., Benavidez, P. G., & Alemañ, R. A. 2013, A&A, 558, A95, doi: 10.1051/0004-6361/201321461
  * Gault et al. (1963) Gault, D. E., Shoemaker, E. M., Moore, H. J., et al. 1963, Spray Ejected from the Lunar Surface by Meteoroid Impact, NASA technical note (National Aeronautics and Space Administration).  https://books.google.co.jp/books?id=dHyZxkQ2XIsC
  * Genda et al. (2017) Genda, H., Fujita, T., Kobayashi, H., et al. 2017, Icarus, 294, 234, doi: 10.1016/j.icarus.2017.03.009
  * Gomes et al. (2005) Gomes, R., Levison, H., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466, doi: 10.1038/nature03676
  * Holsapple (1993) Holsapple, K. 1993, Annual Review of Earth and Planetary Sciences, 21, 333, doi: 10.1146/annurev.ea.21.050193.002001
  * Holsapple & Housen (2007) Holsapple, K. A., & Housen, K. R. 2007, Icarus, 187, 345, doi: 10.1016/j.icarus.2006.08.029
  * Holsapple & Housen (2012) —. 2012, Icarus, 221, 875, doi: 10.1016/j.icarus.2012.09.022
  * Holsapple & Schmidt (1987) Holsapple, K. A., & Schmidt, R. M. 1987, Journal of Geophysical Research, 92, 6350, doi: 10.1029/JB092iB07p06350
  * Housen & Holsapple (2011) Housen, K. R., & Holsapple, K. A. 2011, Icarus, 211, 856, doi: 10.1016/j.icarus.2010.09.017
  * Housen et al. (1983) Housen, K. R., Schmidt, R. M., & Holsapple, K. A. 1983, Journal of Geophysical Research, 88, 2485, doi: 10.1029/JB088iB03p02485
  * Hyodo & Charnoz (2017) Hyodo, R., & Charnoz, S. 2017, AJ, 154, 34, doi: 10.3847/1538-3881/aa74c9
  * Hyodo & Genda (2018) Hyodo, R., & Genda, H. 2018, ApJ, 856, L36, doi: 10.3847/2041-8213/aab7f0
  * Hyodo & Genda (2020) —. 2020, ApJ, 898, 30, doi: 10.3847/1538-4357/ab9897
  * Hyodo et al. (2021) Hyodo, R., Genda, H., & Brasser, R. 2021, Icarus, 354, 114064, doi: 10.1016/j.icarus.2020.114064
  * Hyodo et al. (2019) Hyodo, R., Kurosawa, K., Genda, H., Usui, T., & Fujita, K. 2019, Scientific Reports, 9, doi: 10.1038/s41598-019-56139-x
  * Johnson & Melosh (2012) Johnson, B. C., & Melosh, H. J. 2012, Icarus, 217, 416, doi: 10.1016/j.icarus.2011.11.020
  * Kurosawa et al. (2019) Kurosawa, K., Genda, H., Hyodo, R., et al. 2019, Life Sciences in Space Research, 23, doi: 10.1016/j.lssr.2019.07.006
  * Lucy (1977) Lucy, L. 1977, AJ, 82, 1013, doi: 10.1086/112164
  * Melosh (1989) Melosh, H. 1989, Impact cratering : a geologic process
  * Melosh (1984) Melosh, H. J. 1984, Icarus, 59, 234, doi: 10.1016/0019-1035(84)90026-5
  * Mojzsis et al. (2019) Mojzsis, S. J., Brasser, R., Kelly, N. M., Abramov, O., & Werner, S. C. 2019, ApJ, 881, 44, doi: 10.3847/1538-4357/ab2c03
  * Monaghan (1992) Monaghan, J. 1992, ARA&A, 30, 543, doi: 10.1146/annurev.aa.30.090192.002551
  * Ong et al. (2010) Ong, L., Asphaug, E. I., Korycansky, D., & Coker, R. F. 2010, Icarus, 207, 578, doi: 10.1016/j.icarus.2009.12.012
  * Piekutowski (1980) Piekutowski, A. 1980, Lunar and Planetary Science Conference Proceedings, 3, 2129
  * Pierazzo & Chyba (2002) Pierazzo, E., & Chyba, C. F. 2002, Icarus, 157, 120, doi: 10.1006/icar.2001.6812
  * Sekine et al. (2008) Sekine, T., Kobayashi, T., Nishio, M., & Takahashi, E. 2008, Earth, Planets, and Space, 60, 999, doi: 10.1186/BF03352857
  * Senft & Stewart (2008) Senft, L. E., & Stewart, S. T. 2008, Meteoritics and Planetary Science, 43, 1993, doi: 10.1111/j.1945-5100.2008.tb00657.x
  * Shoemaker (1962) Shoemaker, E. M. 1962, Physics and Astronomy of the Moon, 283.  https://ci.nii.ac.jp/naid/10009520628/
  * Stewart & Ahrens (2003) Stewart, S. T., & Ahrens, T. J. 2003, Geophys. Res. Lett., 30, 1332, doi: 10.1029/2002GL016789
  * Svetsov (2011) Svetsov, V. 2011, Icarus, 214, 316, doi: 10.1016/j.icarus.2011.04.026
  * Tsiganis et al. (2005) Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. 2005, Nature, 435, 459, doi: 10.1038/nature03539
  * Walsh et al. (2011) Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, Nature, 475, 206, doi: 10.1038/nature10201
  * Wong et al. (2020) Wong, E. W., Brasser, R., & Werner, S. C. 2020, Icarus, 114184, doi: https://doi.org/10.1016/j.icarus.2020.114184
  * Yue et al. (2013) Yue, Z., Johnson, B. C., Minton, D. A., et al. 2013, Nature Geoscience, 6, 435, doi: 10.1038/ngeo1828


# Tiling the plane with hexagons: improved separations for $k$-colourings

Aubrey D.N.J. de Grey; Jaan Parts<EMAIL_ADDRESS>
(Mountain View, California, USA; Kazan, Tatarstan, Russia)

###### Abstract

It has been common knowledge since 1950 that seven colours can be assigned to
tiles of an infinite honeycomb with cells of unit diameter such that no two
tiles of the same colour are closer than $d(7)=\frac{\sqrt{7}}{2}$ apart.
Various authors have described tilings using $k>7$ colours, giving
corresponding values for $d(k)$, but it is generally unknown whether these are
the largest possible for a given $k$. Here, for many $k$, we describe tilings
with larger values of $d(k)$ than previously reported.

##  1 Background

In 1950, Isbell observed [7] that a tiling of the plane using regular hexagons
of unit diameter can be $7$-coloured such that tiles of the same colour are at
least $d(7)=\frac{\sqrt{7}}{2}$ apart (see Fig. 1). If we allow more colours
$k$, this minimum separation $d(k)$ increases. In extension of earlier work,
Chybowska-Sokół et al. recently reported [1] values of $d(k)$ for various
$k>7$. Their analysis, like all prior work of which we are aware (with a few
rare exceptions), was basically restricted to regular hexagons. With such a
tiling, if and only if $k$ is a so-called Löschian number (that is,
$k=a^{2}+ab+b^{2}$ for some integers $a>0,a\geq b\geq 0$), colours can be
assigned such that the tiles of any given colour lie on a (regular) hexagonal
sublattice. Löschian numbers are quite common: the first few are
$k\in\\{1,3,4,7,9,12,13,16,19,21,\dots\\}$. In such cases, $d(k)$ is
$\frac{\sqrt{3}}{2}(a-1)$ when $b=0$, otherwise
$\sqrt{(3a+3b-4)^{2}+3(a-b)^{2}}/4$.

It is left to the reader to verify that, in any tiling with regular hexagons,
the smallest $k^{\prime}>k$ that can have all same-coloured pairs of tiles
more than $d(k)$ apart is always also a Löschian number. When multiple choices
of $a,b$ give the same $k\in\\{49,91,133,147,169,\dots\\}$, the largest $d(k)$
arises from the largest value of $b$.

Figure 1: The classical $7$-tiling of the plane with regular hexagons.

However, regular hexagons are not the only ones that can tile the plane, even
if we restrict ourselves to the case where all tiles are identical (or, more
correctly, congruent). We wondered whether tilings with non-regular hexagons
might achieve improved bounds, i.e. larger values for some $d(k)$. Here we
report that they indeed do. In fact, we now have $d(k)>d(k-3)$ for all $k\leq
175$.

##  2 Preliminaries

### Problem statement and main restrictions.

Without loss of generality, we restrict ourselves to tiles of a unit
diameter111 Sometimes it is more correct to use the term unit width, referring
to the existence of figures with a constant width, such as the Reuleaux
triangle and pentagon, which do not fit into a circle of unit diameter. In
fact, we require that the distance between any pair of points on a tile be
less than one. . Our task is to obtain the maximum possible distance $d(k)$
between tiles of the same colour for each $k$.

We restrict the area of study to the case of periodic tiling of so-called
lattice-sublattice scheme, where i) the plane is partitioned into $k$
congruent sets of tiles, obtained one from the other by translation and
reflection (so for analysis it is enough to consider tiles of one chosen
colour), ii) all tiles have the same shape.

Until now, we have implicitly assumed that $d(k)$ is a non-decreasing function
of the number of colours $k$. The restrictions i) and ii) introduced have a
slightly unexpected effect: some $k$ can give a local decrease in $d(k)$.

### A “hierarchy of irregularity” of hexagons that tile the plane.

From the point of view of maximising $d(k)$, convex hexagons and pentagons are
of main interest. A full taxonomy of tilings of the plane with identical
(congruent) convex hexagons has been described by Gardner [2]. Also, Rao
recently showed [6] that the known classification of tilings using identical
pentagons is complete.

Here we restrict ourselves to a narrow subset of hexagons, namely those with
all three diagonals being the same (unit) length and opposite edges being
parallel. We term these rectilinear hexagons. We explored some other classes,
but did not find a hexagonal tiling that gives a greater $d(k)$ than the best
rectilinear one for any $k$ that we examined. However, our search was in no
way exhaustive, so it remains very possible that such tilings exist.

Among this class of hexagon, there is a subset that we term semi-regular:
these have four edges of equal length, so they can be oriented with two
vertices at $(0,\pm\frac{1}{2})$ and the others at $(\pm x,\pm y)$ with
$x^{2}+y^{2}=\frac{1}{4}$. We decided to explore whether there exist values of
$k$ for which semi-regular hexagons improve on regular hexagons, and similarly
whether rectilinear hexagons ever beat semi-regular ones (see Fig. 2). It
turns out that both such scenarios often occur.

Figure 2: Tilings of the plane with regular, semi-regular or rectilinear
hexagons (respectively from left to right).

### Tiling parameters.

An example of tiling is shown in Fig. 3. Tiles of some (base) colour are
highlighted in grey. An oblique coordinate system $(i,j)$, where
$i,j\in\mathbb{Z}$, is introduced for indexing tiles. For ease of
visualization, one of the axes is oriented horizontally.

In our case, the specific tiling is determined by the position of two tiles of
the base colour relative to the base tile, labeled $(0,0)$. For calculations,
it is convenient to use the tiles $(g,0)$ and $(h,k/g)$. (We will define the
parameters $g$ and $h$ below.) To present the results, it is more convenient
to use the tiles $(i_{1},j_{1})$ and $(i_{2},j_{2})$, which together with the
tile $(0,0)$ form the triple with the smallest distances
$\\{d_{01},d_{02},d_{12}\\}$ for the given tiling222 We used a similar
construction in [4] for the case of colouring the plane in several layers. But
in that work we stipulated equality of distances in each triple of the nearest
tiles, $d_{01}=d_{02}=d_{12}$. Perhaps some results obtained there can be
improved by relaxing this restriction., $\min(d_{01},d_{02},d_{12})=d(k)$. The
relation is fulfilled: $k=i_{1}j_{2}-i_{2}j_{1}$.

The shape of a tile is described by two parameters, for example, the length of
two adjacent edges $\\{r,s\\}$, of which we are interested in the length $r$
of the vertical edge that completely defines the shape of a semi-regular
hexagon.

Figure 3: Notation for the tilings discussed in this study. In the example
shown, $k=22$, $g=11$ and $h=6$.

### Finding the shortest distance between hexagons.

The colouring options for a given $k$ and a given class of hexagon can be
defined in terms of only two parameters, $g$ and $h$, with $g$ a divisor of
$k$ and $0\leq h<g$. Each horizontal row is coloured in a repeating sequence
of $g$ colours, so for any $n$, row $k/g+n$ uses the same colours as row $n$,
offset by $h$ hexagons.

We are thus interested in the minimum distance between a point in hexagon
$(0,0)$ and a point in hexagon $(ah+bg,ak/g)$, $a,b\in\mathbb{Z}$. It suffices
to consider all tiles with coordinates $|\,i\,|+|\,j\,|<2\sqrt{k},\;j\geq 0$.
We wish to identify the shape of hexagon (within the class being considered)
and values of $g$ and $h$ that maximise that distance. For the classes we are
considering, all of a hexagon’s angles must be obtuse, and each pair of
opposite edges forms a rectangle. Thus, as with regular hexagons, $d(k)$ is
either a multiple of the length of the relevant rectangle or a distance
between corners of the two tiles.

### Computation.

Our work was performed using Mathematica, version 13. Briefly, we used
numerical optimisation (the function NMaxValue or FindMaximum) to determine
the best parameter values for a given $k$, and then, if the result improved
upon any for a more restricted class of hexagon with the same or a smaller $k$
or for any class with a strictly smaller $k$, we derived exact expressions for
$d(k)$ and the parameters using Maximize.

##  3 Results

### Main results.

Table 1 compiles our findings, for $3\leq k\leq 175$ (below we will assume
this range of $k$ by default). Note that we show $d^{2}(k)$, not $d(k)$,
because $d^{2}(k)$ is usually rational whereas $d(k)$ rarely is. A blank cell
indicates that the best result obtained is no better than that for a more
restricted class of hexagon with the same $k$ or for some class with a smaller
$k$. However, for each $k$ we show the best result even when it is beaten by a
smaller $k$; cases that do beat all smaller $k$ are denoted by bold-face in
the last column.

The coordinates $(i_{1},j_{1})$, $(i_{2},j_{2})$ correspond to the tiling that
provides the best obtained value of $d(k)$ for a given $k$. For semi-regular
hexagons we also give a parameter $r$, which is the length of the vertical
edges. Where possible, we impose additional restrictions on the parameters:
tile 1 lies between the positive directions of the $(i-j)$ and $i$ axes, tile
2 lies between the $i$ and $j$ axes, and also the condition $d_{01}=d_{02}\leq
d_{12}$ is satisfied. In the case of $k=77$, these conditions contradict each
other, and the distance $d_{01}(77)$ is the largest in the triplet.

The values for rectilinear $k=11,23,45$ and $187$ are roots of quartic
equations. See the next subsection for more details, including the definition
of the function $f(k)$ given as $d^{2}(k)$ for certain $k$.

Table 1: Optimal $d(k)$ for regular, semi-regular and rectilinear hexagons. $k$ | $i_{1}$ | $j_{1}$ | $i_{2}$ | $j_{2}$ |  $d^{2}(k)$ [; $r$] | $d(k)$
---|---|---|---|---|---|---
|  |  |  |  |  regular |  semi-regular |  rectilinear | approx.
3 | 2 | -1 | 1 | 1 |  $(1/2)^{2}$ |  |  | 0.500000
4 | 2 | 0 | 0 | 2 |  3/4 |  |  | 0.866025
5 | 3 | -1 | 2 | 1 |  |  $(5/6)^{2}$; 5/6 |  | 0.833333
6 | 3 | 0 | 0 | 2 |  |  $\frac{17+\sqrt{33}}{24}$; $\frac{\sqrt{33}-3}{12}$ |  $f(6)$ | 0.992076
7 | 3 | -1 | 1 | 2 |  7/4 |  |  | 1.322876
8 | 3 | -2 | 1 | 2 |  |  $(7/5)^{2}$; 1/5 |  | 1.400000
9 | 3 | 0 | 0 | 3 |  3 |  |  | 1.732051
10 | 4 | -1 | 2 | 2 |  |  11/4; 3/4 |  | 1.658312
11 | 3 | -2 | 1 | 3 |  |  11/4; 3/4 |  $1.674155^{2}$ | 1.674155
12 | 4 | -2 | 2 | 2 |  $2^{2}$ |  |  | 2.000000
13 | 4 | -1 | 1 | 3 |  19/4 |  |  | 2.179449
14 | 4 | -2 | 1 | 3 |  |  19/4; 1/4 |  | 2.179449
15 | 4 | -3 | 1 | 3 |  |  153/32; 3/8 |  | 2.186607
16 | 4 | 0 | 0 | 4 |  27/4 |  |  | 2.598076
17 | 5 | -1 | 2 | 3 |  |  627/100; 7/10 |  709371/111724 | 2.519785
18 | 4 | -3 | 2 | 3 |  |  6; 2/3 |  $20-\sqrt{192}$ | 2.478627
19 | 5 | -2 | 2 | 3 |  31/4 |  |  | 2.783882
20 | 5 | 0 | 0 | 4 |  |  $\frac{5931+405\sqrt{129}}{1352}$; $\frac{3\sqrt{129}-15}{52}$ |  $f(20)$ | 2.862821
21 | 5 | -1 | 1 | 4 |  37/4 |  |  | 3.041381
22 | 5 | -2 | 1 | 4 |  |  $3^{2}$; 1/3 |  193496/21275 | 3.015791
23 | 5 | -3 | 1 | 4 |  |  103/12; 1/6 |  $2.942985^{2}$ | 2.942985
24 | 6 | -2 | 3 | 3 |  |  1216/121; 29/44 |  869/86 | 3.178781
25 | 5 | 0 | 0 | 5 |  12 |  |  | 3.464102
26 | 6 | -1 | 2 | 4 |  |  203/18; 2/3 |  2150/187 | 3.390771
27 | 6 | -3 | 3 | 3 |  $(7/2)^{2}$ |  |  | 3.500000
28 | 6 | -2 | 2 | 4 |  13 |  |  | 3.605551
29 | 7 | -2 | 4 | 3 |  |  44719/3844; 47/62 |  481078403/41028908 | 3.424230
30 | 6 | 0 | 0 | 5 |  |  $\frac{1072+128\sqrt{51}}{147}$; $\frac{2\sqrt{51}-6}{21}$ |  $f(30)$ | 3.762292
31 | 6 | -1 | 1 | 5 |  61/4 |  |  | 3.905125
32 | 6 | -1 | 2 | 5 |  |  59/4; 3/8 |  2307733/153928 | 3.871988
33 | 7 | -3 | 4 | 3 |  |  $(227/58)^{2}$; 37/58 |  | 3.913793
34 | 7 | -2 | 3 | 4 |  |  1921/121; 7/11 |  170126/10611 | 4.004121
35 | 7 | 0 | 1 | 5 | $\frac{17559+135\sqrt{12257}}{2048};\frac{\sqrt{12257}-15}{128}$ |  | 3.983917
36 | 6 | 0 | 0 | 6 |  75/4 |  |  | 4.330127
37 | 7 | -3 | 3 | 4 |  73/4 |  |  | 4.272002
38 | 6 | -1 | 2 | 6 |  |  17; 3/5 |  343577/19520 | 4.195388
39 | 7 | -2 | 2 | 5 |  79/4 |  |  | 4.444097
40 | 7 | -4 | 3 | 4 |  |  $(137/31)^{2}$; 11/31 |  | 4.419355
41 | 7 | -2 | 3 | 5 |  |  637/32; 3/8 |  19168435/946444 | 4.500346
42 | 7 | 0 | 0 | 6 |  $\frac{84075+13125\sqrt{33}}{7688}$; $\frac{15\sqrt{33}-35}{124}$ |  $f(42)$ | 4.651075
43 | 7 | -1 | 1 | 6 |  91/4 |  |  | 4.769696
44 | 7 | -1 | 2 | 6 |  |  22; 2/5 |  4700717/209612 | 4.735589
45 | 7 | -3 | 1 | 6 |  |  423/20; 3/10 |  $4.607942^{2}$ | 4.607942

Continuation of Table 1.

$k$ | $i_{1}$ | $j_{1}$ | $i_{2}$ | $j_{2}$ |  $d^{2}(k)$ [; $r$] | $d(k)$
---|---|---|---|---|---|---
|  |  |  |  |  regular |  semi-regular |  rectilinear | approx.
46 | 8 | -2 | 3 | 5 |  |  49051/2116; 57/92 |  73151/3116 | 4.845197
47 | 7 | -3 | 4 | 5 |  |  3224/147; 13/21 |  5560412771/243356228 | 4.780048
48 | 8 | -4 | 4 | 4 |  $5^{2}$ |  |  | 5.000000
49 | 7 | 0 | 0 | 7 |  27 |  |  | 5.196152
50 | 8 | -1 | 2 | 6 |  |  825/32; 5/8 |  2143125/81224 | 5.136669
51 | 7 | -1 | 2 | 7 |  |  6016/243; 19/27 |  32914863743/1281311828 | 5.068374
52 | 8 | -2 | 2 | 6 |  28 |  |  | 5.291503
53 | 9 | -2 | 4 | 5 |  |  358639/13924; 83/118 |  21713005483/826893628 | 5.124307
54 | 8 | -2 | 3 | 6 |  |  335/12; 5/12 |  118349/4160 | 5.333787
55 | 9 | -1 | 1 | 6 |  |  $(449/86)^{2}$; 21/86 |  15510042323/562475260 | 5.251154
56 | 8 | 0 | 0 | 7 |  |  $(277/51)^{2}$; 31/51 |  $f(56)$ | 5.533452
57 | 8 | -1 | 1 | 7 |  127/4 |  |  | 5.634714
58 | 8 | -1 | 2 | 7 |  |  123/4; 5/12 |  7191956/229175 | 5.601959
59 | 8 | -1 | 3 | 7 |  |  89/3; 1/3 |  742918955/24763644 | 5.477261
60 | 9 | -2 | 3 | 6 |  |  118969/3721; 37/61 |  4775023/147232 | 5.694907
61 | 9 | -4 | 4 | 5 |  133/4 |  |  | 5.766281
62 | 8 | -2 | 3 | 7 |  |  49987/1600; 47/80 |  2563175/79001 | 5.696037
63 | 9 | -3 | 3 | 6 |  139/4 |  |  | 5.894913
64 | 8 | 0 | 0 | 8 |  147/4 |  |  | 6.062178
65 | 9 | -1 | 2 | 7 |  |  11431/324; 11/18 |  42769993/1184612 | 6.008714
66 | 9 | -3 | 4 | 6 |  |  23987/675; 17/45 |  1193079661/33188989 | 5.995670
67 | 9 | -2 | 2 | 7 |  151/4 |  |  | 6.144103
68 | 10 | -2 | 4 | 6 |  |  12649/361; 13/19 |  3193957/89125 | 5.986387
69 | 9 | -2 | 3 | 7 |  |  262/7; 3/7 |  6272261/164372 | 6.177292
70 | 10 | -4 | 5 | 5 |  |  44339/1156; 81/136 |  414848/10783 | 6.202613
71 | 10 | -1 | 1 | 7 |  |  1035/28; 5/14 |  3274003629/86168876 | 6.164025
72 | 9 | 0 | 0 | 8 |  $\frac{177233+7203\sqrt{537}}{8664}$; $\frac{7\sqrt{537}-63}{228}$ |  $f(72)$ | 6.411744
73 | 9 | -1 | 1 | 8 |  $(13/2)^{2}$ |  |  | 6.500000
74 | 9 | -1 | 2 | 8 |  |  41; 3/7 |  9026779/215683 | 6.469317
75 | 10 | -5 | 5 | 5 |  $(13/2)^{2}$ |  |  | 6.500000
76 | 10 | -4 | 4 | 6 |  43 |  |  | 6.557439
77 | 11 | 0 | 2 | 7 |  |  1215/32; 5/8 |  | 6.161879
78 | 9 | -2 | 3 | 8 |  |  21844/529; 53/92 |  4833828787/112826683 | 6.545453
79 | 10 | -3 | 3 | 7 |  181/4 |  |  | 6.726812
80 | 10 | -4 | 5 | 6 |  $\frac{63261+972\sqrt{3770}}{2738}$; $\frac{\sqrt{3770}-12}{74}$ |  241351/5371 | 6.703428
81 | 9 | 0 | 0 | 9 |  48 |  |  | 6.928203
82 | 10 | -1 | 2 | 8 |  |  1157/25; 3/5 |  1184696/25029 | 6.879893
83 | 9 | -1 | 2 | 9 |  |  957559/21316; 47/146 |  975386934409/20985431116 | 6.817569
84 | 10 | -2 | 2 | 8 |  $7^{2}$ |  |  | 7.000000
85 | 11 | -5 | 6 | 5 |  |  $(1097/158)^{2}$; 93/158 |  | 6.943038
86 | 10 | -2 | 3 | 8 |  |  40953/841; 17/29 |  30121/610 | 7.026997
87 | 9 | -2 | 3 | 9 |  |  638661/13924; 79/118 |  1080643/23188 | 6.826679
88 | 10 | -2 | 4 | 8 |  |  191/4; 3/8 |  84065/1729 | 6.972847
89 | 11 | -3 | 4 | 7 |  |  194709/3844; 109/186 |  511258004631/9991926364 | 7.153119

Continuation of Table 1.

$k$ | $i_{1}$ | $j_{1}$ | $i_{2}$ | $j_{2}$ |  $d^{2}(k)$ [; $r$] | $d(k)$
---|---|---|---|---|---|---
|  |  |  |  |  regular |  semi-regular |  rectilinear | approx.
90 | 10 | 0 | 0 | 9 |  $\frac{140352+30720\sqrt{19}}{5329}$; $\frac{12\sqrt{19}-20}{73}$ |  $f(90)$ | 7.287290
91 | 10 | -1 | 1 | 9 |  217/4 |  |  | 7.365460
92 | 10 | -1 | 2 | 9 |  |  211/4; 7/16 |  4639657267/86189275 | 7.336963
93 | 11 | -4 | 4 | 7 |  217/4 |  |  | 7.365460
94 | 11 | -2 | 3 | 8 |  |  $(713/97)^{2}$; 57/97 |  666148999/12139951 | 7.407595
95 | 12 | -5 | 7 | 5 |  |  $(277/38)^{2}$; 25/38 |  | 7.289474
96 | 11 | -6 | 5 | 6 |  |  $(67/9)^{2}$; 11/27 |  | 7.444444
97 | 11 | -3 | 3 | 8 |  229/4 |  |  | 7.566373
98 | 12 | -1 | 2 | 8 |  |  115625/2178; 8/33 |  613512713/11338873 | 7.355748
99 | 11 | 0 | 1 | 9 |  $\frac{498775+1575\sqrt{91553}}{16928}$; $\frac{\sqrt{91553}-63}{368}$ |  | 7.590561
100 | 10 | 0 | 0 | 10 |  243/4 |  |  | 7.794229
101 | 11 | -1 | 2 | 9 |  |  28449/484; 13/22 |  12663270675/210818572 | 7.750300
102 | 12 | -5 | 6 | 6 |  |  29651/500; 29/50 |  208565/3509 | 7.709550
103 | 11 | -2 | 2 | 9 |  247/4 |  |  | 7.858117
104 | 12 | -4 | 5 | 7 |  |  176485/2916; 125/216 |  4142419871/67995631 | 7.805245
105 | 11 | -2 | 3 | 9 |  |  1097/18; 4/9 |  447216991/7200892 | 7.880722
106 | 10 | -2 | 3 | 10 |  |  $(84/11)^{2}$; 29/44 |  33151587/560347 | 7.691723
107 | 11 | -5 | 6 | 7 |  |  2161/36; 7/18 |  $\frac{12486512240719}{202886322964}$ | 7.845022
108 | 12 | -6 | 6 | 6 |  $8^{2}$ |  |  | 8.000000
109 | 12 | -5 | 5 | 7 |  259/4 |  |  | 8.046738
110 | 11 | 0 | 0 | 10 |  $\frac{2184003+72171\sqrt{849}}{66248}$; $\frac{9\sqrt{849}-99}{364}$ |  $f(110)$ | 8.160909
111 | 11 | -1 | 1 | 10 |  271/4 |  |  | 8.231039
112 | 11 | -1 | 2 | 10 |  67 |  |  1034089141/15361771 | 8.204618
113 | 11 | -1 | 3 | 10 |  |  451319/6962; 38/59 |  387649933983/5906484172 | 8.101312
114 | 12 | -2 | 3 | 9 |  |  937099/13924; 137/236 |  283698767/4150532 | 8.267550
115 | 13 | -4 | 6 | 7 |  |  4244941/64516; 163/254 |  $\frac{12272144911631}{184524424676}$ | 8.155175
116 | 12 | -4 | 5 | 8 |  |  8723/128; 13/32 |  2045764331/29695775 | 8.300045
117 | 12 | -3 | 3 | 9 |  283/4 |  |  | 8.411302
118 | 11 | -2 | 4 | 10 |  |  451319/6962; 38/59 |  22520512223/336563063 | 8.180049
119 | 13 | -3 | 5 | 8 |  |  $(1643/198)^{2}$; 65/198 |  $\frac{17514883072583}{252372028148}$ | 8.330729
120 | 12 | -3 | 4 | 9 |  |  $\frac{28677+3630\sqrt{58}}{784}$; $\frac{3\sqrt{58}-5}{28}$ |  102121645769/1420956281 | 8.477515
121 | 11 | 0 | 0 | 11 |  75 |  |  | 8.660254
122 | 12 | -1 | 2 | 10 |  |  655/9; 7/12 |  10143575/136512 | 8.620056
123 | 11 | -1 | 2 | 11 |  |  392863/5476; 79/222 |  129728247871/1768473940 | 8.564814
124 | 12 | -2 | 2 | 10 |  76 |  |  | 8.717798
125 | 13 | -2 | 4 | 9 |  |  5551687/77284; 179/278 |  976273547119/13237256716 | 8.587896
126 | 12 | -2 | 3 | 10 |  |  1499/20; 9/20 |  19145005/250792 | 8.737172
127 | 13 | -6 | 6 | 7 |  307/4 |  |  | 8.760708
128 | 12 | -4 | 5 | 9 |  |  369/5; 2/5 |  795830295/10506727 | 8.703151
129 | 13 | -5 | 5 | 8 |  313/4 |  |  | 8.845903
130 | 13 | 0 | 1 | 10 |  $\frac{3182571+312741\sqrt{97}}{84872}$; $\frac{27\sqrt{97}-143}{412}$ |  $\frac{18088+16\sqrt{1192555}}{469}$ | 8.707601
131 | 13 | -1 | 1 | 10 |  |  3094091/40804; 81/202 |  192224305737/2452589252 | 8.853026
132 | 12 | 0 | 0 | 11 |  $\frac{165700+10000\sqrt{258}}{4107}$; $\frac{5\sqrt{258}-30}{111}$ |  $f(132)$ | 9.033127

Continuation of Table 1.

$k$ | $i_{1}$ | $j_{1}$ | $i_{2}$ | $j_{2}$ |  $d^{2}(k)$ [; $r$] | $d(k)$
---|---|---|---|---|---|---
|  |  |  |  |  regular |  semi-regular |  rectilinear | approx.
133 | 12 | -1 | 1 | 11 |  331/4 |  |  | 9.096703
134 | 12 | -1 | 2 | 11 |  |  323/4; 9/20 |  6558879037/79690552 | 9.072174
135 | 12 | -1 | 3 | 11 |  |  787/10; 2/5 |  $\frac{10272242758901}{127501682004}$ | 8.975831
136 | 13 | -2 | 3 | 10 |  |  181289/2209; 27/47 |  20689231/248261 | 9.128889
137 | 13 | -4 | 5 | 9 |  |  2598661/31684; 73/178 |  $\frac{74222844365281}{892932274564}$ | 9.117159
138 | 12 | -2 | 3 | 11 |  |  27389/338; 19/52 |  1122833/13505 | 9.118225
139 | 13 | -3 | 3 | 10 |  343/4 |  |  | 9.260130
140 | 14 | -6 | 7 | 7 |  |  134518/1587; 157/276 |  80581/949 | 9.214743
141 | 14 | -3 | 5 | 9 |  |  929347/11236; 67/106 |  196014091823/2324718468 | 9.182448
142 | 13 | -3 | 4 | 10 |  |  117688/1369; 21/37 |  15150070069/174437335 | 9.319392
143 | 13 | 0 | 1 | 11 |  $\frac{399015+4455\sqrt{7553}}{8978}$; $\frac{5\sqrt{7553}-99}{536}$ |  | 9.357805
144 | 12 | 0 | 0 | 12 |  363/4 |  |  | 9.526279
145 | 13 | -1 | 2 | 11 |  |  59675/676; 15/26 |  3279780625/36423236 | 9.489277
146 | 14 | -4 | 5 | 9 |  |  119853/1352; 59/104 |  16400637/183488 | 9.454238
147 | 13 | -2 | 2 | 11 |  367/4 |  |  | 9.578622
148 | 14 | -6 | 6 | 8 |  91 |  |  | 9.539392
149 | 13 | -2 | 3 | 11 |  |  995/11; 5/11 |  1164694451/12649484 | 9.595544
150 | 13 | -4 | 5 | 10 |  |  112529/1280; 103/160 |  286692851539/3160973281 | 9.523531
151 | 14 | -5 | 5 | 9 |  373/4 |  |  | 9.656604
152 | 14 | -3 | 4 | 10 |  |  405821/4374; 46/81 |  7532633043/80257313 | 9.687932
153 | 15 | -6 | 8 | 7 |  |  1030243/11236; 199/318 |  $\frac{894213901557691}{9724451350300}$ | 9.589328
154 | 14 | -6 | 7 | 8 |  |  771407/8192; 109/256 |  1311592/13897 | 9.714912
155 | 14 | -1 | 1 | 11 |  |  133856/1445; 53/85 |  $\frac{15557710498703}{164015315380}$ | 9.739365
156 | 13 | 0 | 0 | 12 |  97 |  |  $f(156)$ | 9.904296
157 | 13 | -1 | 1 | 12 |  397/4 |  |  | 9.962429
158 | 13 | -1 | 2 | 12 |  |  97; 5/11 |  75425702171/763451750 | 9.939599
159 | 13 | -1 | 3 | 12 | 12176361/128164; 223/358 |  184097879615/1897621332 | 9.849623
160 | 14 | -2 | 3 | 11 |  |  677416/6889; 189/332 |  9726226252/97433683 | 9.991199
161 | 15 | -7 | 8 | 7 |  |  9284209/93636; 173/306 |  | 9.957516
162 | 15 | -6 | 7 | 8 |  |  301273/3025; 31/55 |  195650535281/1958502269 | 9.994901
163 | 14 | -3 | 3 | 11 |  409/4 |  |  | 10.11187
164 | 14 | -4 | 6 | 10 |  |  292/3; 5/12 |  267614/2691 | 9.972351
165 | 15 | -5 | 6 | 9 |  |  2487641/24500; 197/350 |  $\frac{33252082489597}{325566557812}$ | 10.10624
166 | 14 | -3 | 4 | 11 |  |  813/8; 11/24 |  544579139/5269928 | 10.16549
167 | 15 | -1 | 2 | 11 |  |  6025357/60516; 77/246 |  7169589821/70929460 | 10.05388
168 | 14 | 0 | 1 | 12 |  $\frac{662311+5070\sqrt{16226}}{12482}$; $\frac{\sqrt{16226}-30}{158}$ |  | 10.23727
169 | 13 | 0 | 0 | 13 |  108 |  |  | 10.39230
170 | 14 | -1 | 2 | 12 |  |  10317/98; 4/7 |  1618378002/15084233 | 10.35806
171 | 15 | -6 | 6 | 9 |  427/4 |  |  | 10.33199
172 | 14 | -2 | 2 | 12 |  109 |  |  | 10.44031
173 | 15 | -2 | 4 | 11 | 15153631/145924; 239/382 |  $\frac{37802632989571}{354544269100}$ | 10.32585
174 | 14 | -2 | 3 | 12 |  |  2579/24; 11/24 |  1893749/17324 | 10.45531
175 | 15 | -5 | 5 | 10 |  439/4 |  |  | 10.47616

3 |  | 7 |
---|---|---|---
4 |  | 8 |
5 |  | 9 |
6 |  | 10 |
*6 |  | 11 |
Figure 4: Optimal hexagonal (* and pentagonal) tilings for $3\leq k\leq 11$.  12 |  | 17 |
---|---|---|---
13 |  | 18 |
14 |  | 20 |
15 |  | 22 |
16 |  | 23 |
Figure 5: Optimal hexagonal tilings for $12\leq k\leq 23$.

Tilings with the maximum $d(k)$ that we obtained for various small $k$ are
shown in Figures 4 and 5. All tilings are given in the same scale. Note that
the hexagons with $k=17,22$ and $23$ are not precisely semi-regular.

### Classification.

Other than the Löschian class $k=3n^{2}$ ($n\geq 1$), for which
$d(k)=\frac{3n}{2}-1$ in the regular setting (except for $k=147$, where
another colouring of regular hexagons wins), $d(k)$ is not often rational; the
only examples we found are regular $k=84$ and semi-regular $k=5$, 8, 22, 33,
40, 55, 56, 85, 94, 95, 96, 106 and 119 (some of which are beaten by
rectilinear examples).

However, $d^{2}(k)$ is usually rational. The obtained values can be classified
according to the degree of the polynomial of which $d^{2}(k)$ is a root:

i) When $d^{2}(k)$ is rational, in the semi-regular setting the denominator is
always a square or a small multiple of a square, and is seldom much greater
than $k^{2}$. (Also, the denominator of $r$ does not exceed $3k$.) In the
rectilinear setting, however, the denominator is often far larger, reaching
the tens of millions even for $k$ as small as $29$.

ii) In all cases where the semi-regular $d^{2}(k)$ is irrational, it has the
form $\mathbb{Q}[\sqrt{a}]$ for some integer $a$. All such solutions appear at
$k=n(n+m)$, $0<m<<n$. In the rectilinear setting, there are a few examples
where the best $d^{2}(k)$ has the quadratic form. In the range of $k$ that we
explored, the only such examples are $k=18$ and $130$ in which there is a
rectilinear tiling that beats the semi-regular one, and $k=35,99$ and $143$ in
which the quadratic-form semi-regular $d(k)$ is unbeaten.

iii) The largest irrational class is $k=n(n+1)$ for $n\geq 2$, with the
exception of the Löschian case $k=12$. Here, the rectilinear $d^{2}(k)$,
denoted in the table by $f(k)$, is the largest real root333The other two roots
of $f(k)$ tend to $-1/3$ and $+1$ as $k$ increases. of the cubic polynomial
$a_{3}x^{3}+a_{2}x^{2}+a_{1}x+a_{0}=0$, where $p=n(n-1)=k+1-\sqrt{4k+1}$,
$a_{3}=4p(p^{2}+3p+1)$, $a_{2}=-3p^{4}-8p^{3}+2p^{2}+4p+1$,
$a_{1}=2p^{2}(p^{2}-2p-1)$, $a_{0}=p^{4}$. We do not know why.

iv) Only four values that we checked, namely $k=11,23,45$ and $187$, have an
optimal $x=d^{2}(k)$ in the rectilinear setting that is a root of a polynomial
of degree greater than 3, in all cases quartic. For completeness we list those
polynomials here (take the smaller of the two real roots):

$25600x^{4}-1459616x^{3}+8840377x^{2}-18735876x+13623552$

$97344x^{4}-14493200x^{3}+287680857x^{2}-2041299600x+4968218624$

$548800x^{4}-70030800x^{3}+2721532527x^{2}-44612184348x+279110188800$

$\;\;391936098304x^{4}-214522652834752x^{3}+43720132398169569x^{2}-\\\
-4062662189783485600x+145700736445997574400$

### Observations.

As expected, with Löschian values of $k$ we found that regular hexagons are
hard to beat. The smallest Löschian $k$ for which regular hexagons are beaten
by another class (in this case a rectilinear one but not a semi-regular one)
is $k=112$.

By contrast, we obtained substantial improvements on earlier results with many
non-Löschian $k$. For example, we achieved $d(8)=7/5$, a considerable
improvement on the value of 1.37542 reported in [1].

In the case of $k=156$ not only do rectilinear hexagons beat regular ones, but
both are record breaking compared to the previous $d(155)$. We checked $f(k)$
for other cases where $k=n(n+1)$ is Löschian: $\\{756,1332,$ $2352,\dots\\}$.
For all subsequent $k<10^{10}$ rectilinear hexagons beat regular ones. It is
natural to conjecture that this will be true for all $k>12$.

Rectilinear hexagons can usually beat semi-regular ones, especially in the
case where semi-regular hexagons beat regular ones. The only $k\leq 175$ for
which a semi-regular hexagon beats regular ones, is not beaten by a
rectilinear one, and beats all smaller $k$, are $k=8,15,33,96,99,143,168$.

Except for the few cases $k=11,15,18,23,45,77$ and $187$, the optimal tiling
has each hexagon equidistant from six same-coloured ones, i.e.
$d_{01}=d_{02}=d_{12}$. Note that the exceptions completely include the class
associated with the fourth degree polynomial.

For $k=80$ and $120$, a quadratic-form semi-regular $d^{2}(k)$ is beaten by a
rectilinear example in which $d^{2}(k)$ is rational.

For some values of $k$ there are strikingly large values of $n$ for which
$d(k)<d(k-n)$. For example, $d(29)<d(25)$, $d(77)<d(69)$, $d(98)<d(91)$,
$d(187)<d(169)$.

##  4 Avenues for future work

### More general hexagons.

Clearly this work can be extended to the analysis of more general tilings. The
stipulation that the hexagon must be rectilinear can be relaxed, for example.
More ambitious would be to allow hexagons of multiple shapes, or to assign
colours to sublattices that are not translations of each other. In general, we
do not expect such tilings to beat the ones we describe here – though the fact
that there remain so many cases where $d(k)<d(k-1)$ may be viewed as evidence
to the contrary.

### Beyond hexagons.

Looking beyond hexagons, however, the outlook is much rosier. We found a
tiling that can be 8-coloured with $d\approx 1.444157$, a remarkable
improvement over the value of $7/5$ found in this study. We also found 14- and
15-tilings that beat hexagons by similarly impressive margins: $d(14)\approx
2.260808$, $d(15)\approx 2.346969$. Further details will be provided in a
forthcoming article (see also the discussion within the Polymath16 project
[5]).

### Towards d(6) = 1.

Finally we come to the only report of which we are aware that has discussed
$d(k)$ in the context of tiling with identical but non-regular hexagons.
Ivanov reported [3] that $d(6)$ can be as much as $\approx 0.991445$ using
rectilinear hexagons. As shown in Table 1, we obtained a fractionally better
value with the same design, $\approx 0.992076$ (see Fig. 4); we do not know
how Ivanov derived the parameters leading to his value. For comparison, the
best pentagonal tiling we know has $d(6)=\sqrt{55/56}\approx 0.991031$ (marked
with an asterisk in Fig. 4).

Since a block-structured $6$-tiling of the plane excluding distance 1 can only
use tiles with at most five edges, simplistic application of the methods used
here cannot reach $d(6)=1$. However, we cannot exclude such possibility with
more general polygons. The remarkable improvements to $d(8),d(14)$ and $d(15)$
noted above may inspire optimism in this regard.

### Redundant colours?

As we already noted, the nonmonotonicity of the function $d(k)$ is associated
with the restrictions of the study area. However, if we remove all additional
restrictions, then for some $k$ we still get $d(k)=d(k-1)$ and even
$d(k)=d(k-2)$ (minimal examples with $d>0$ here are $k=5$ and $k=11$
respectively). In other words, adding one or two colours does not change the
maximum distance between tiles of the same colour. And this state of affairs
is extremely surprising.

The question arises whether the function $d(k)$ is strictly monotonic or are
some colours actually redundant? The available data leave this question open.

## References

  * [1] J. Chybowska-Sokół, K. Junosza-Szaniawski, and K. Węsek. Coloring distance graphs on the plane. arXiv:2201.04499v1, 12 Jan 2022,  https://arxiv.org/abs/2201.04499
  * [2] M. Gardner. Tilings with Convex Polygons. Time Travel and Other Mathematical Bewilderments. New York: W. H. Freeman, 1988, ch. 13, pp. 162–176.
  * [3] L.L. Ivanov. On the chromatic numbers of $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ with an interval of forbidden distances. Discrete and Continuous Models and Applied Computational Science, no. 1, 2011, pp. 14–29. (In Russian.)
  * [4] J. Parts. On upper bounds for the multi-fold chromatic numbers of the plane. Geombinatorics, vol. 30, no. 4, 2021, pp. 177–189.
  * [5] D.H.J. Polymath. Comments in Polymath16, thread 17, Dec 2021 – Feb 2022,  https://asone.ai/polymath/index.php?title= Hadwiger-Nelson_problem
  * [6] M. Rao. Exhaustive search of convex pentagons which tile the plane. arXiv: 1708.00274, 1 Aug 2017, https://arxiv.org/abs/1708.00274
  * [7] A. Soifer. The mathematical coloring book. Mathematics of coloring and the colorful life of its creators. Springer, New York, 2009.