Current address: ]CrystMat Company, CH-8037 Zurich, Switzerland

# Field-angle dependent vortex lattice phase diagram in MgB2

A. W. D. Leishman  Department of Physics, University of Notre Dame, Notre
Dame, Indiana 46556, USA  A. Sokolova  Australian Centre for Neutron
Scattering, Australian Nuclear Science and Technology Organization, NSW 2234,
Australia  M. Bleuel  NIST Center for Neutron Research, National Institute of
Standards and Technology, 100 Bureau Drive, Gaithersburg, Maryland 20899-8562,
USA  Department of Materials Science and Engineering, University of Maryland,
College Park, Maryland 20742-2115, USA  N. D. Zhigadlo  [  Laboratory for
Solid State Physics, ETH, CH-8093 Zurich, Switzerland  M. R. Eskildsen
<EMAIL_ADDRESS>Department of Physics, University of Notre Dame, Notre Dame,
Indiana 46556, USA

###### Abstract

Using small-angle neutron scattering we have studied the superconducting
vortex lattice (VL) phase diagram in MgB2 as the applied magnetic field is
rotated away from the $c$ axis and towards the basal plane. The field rotation
gradually suppresses the intermediate VL phase which exists between end states
aligned with two high symmetry directions in the hexagonal basal plane for
$\bm{H}\parallel\bm{c}$. Above a critical angle, the intermediate state
disappears, and the previously continuous transition becomes discontinuous.
The evolution towards the discontinuous transition can be parametrized by a
vanishing twelvefold anisotropy term in the VL free energy.

##  I Introduction

Vortex matter in type-II superconductors is highly sensitive to the
environment provided by the host material. An example is the vortex lattice
(VL) symmetry and orientation, which is governed by anisotropies in the
screening current plane perpendicular to the applied field and the associated
nonlocal vortex interactions.Kogan (1981); Kogan  _et al._ (1997) Such
anisotropies may arise from the Fermi surface or, in non-$s$ wave
superconductors, from nodes in the superconducting gap. A rich VL phase
diagram often arises when this anisotropy is incommensurate with an
equilateral triangular VL, as seen in Nb with the applied field along the
[100] crystalline direction.Laver  _et al._ (2006, 2009); Mühlbauer  _et al._
(2009) However, structural transitions between different VL configurations can
also arise when the field is applied perpendicular to a sixfold symmetric
crystal plane. In such cases higher-order contributions to the screening
current plane anisotropy become relevant, affecting the orientation of the
triangular VL relative to the crystalline axes as seen in UPt3Gannon  _et al._
(2015); Avers  _et al._ (2020) and MgB2.Cubitt  _et al._ (2003); Das  _et al._
(2012)

In the hexagonal superconductor MgB2, three different triangular VL phases
(labeled F, L and I) were observed for $\bm{H}\parallel\bm{c}$, distinguished
by their orientation relative to the crystalline axes.Das  _et al._ (2012) The
multiple phases arise from a competition between the superconducting $\pi$\-
and $\sigma$-bands which have opposite in-plane anisotropies,Hirano  _et al._
(2013) and are shown in the field-temperature ($H\\!-\\!T$) plane of Fig.
1(a).

Figure 1:  (a) Qualitative MgB2 VL phase diagram in the $H\\!-\\!T$
($\Omega=0$) and $H\\!-\\!\Omega$ ($T=0$) planes. The former was determined
previously,Das  _et al._ (2012) and the latter is discussed in this work. The
VL L phase vanishes at the critical angle, $\Omega_{0}$. (b) Experimental
geometry, indicating the direction and rotation axis of the magnetic field
relative to the crystalline axes.

For low fields, the $\pi$ band anisotropy causes the VL to align with the
$\bm{a}$ crystal axis, producing the F phase. At high fields, where
superconductivity on the $\pi$-band is suppressed,Eskildsen  _et al._ (2002);
Cubitt  _et al._ (2003) the $\sigma$ band anisotropy aligns the VL with the
$\bm{a^{*}}$ crystal axis, producing the I phase. At intermediate fields, the
VL undergoes a continuous $30^{\circ}$ rotation transition between the two
extremes in the L phase. Since clockwise and counterclockwise rotations are
degenerate, the VL fractures into oppositely rotated domains. This domain
formation leads to robust metastable states which inhibit transitions between
the three phases and demonstrate activated behavior.Das  _et al._ (2012);
Rastovski  _et al._ (2013); Louden  _et al._ (2019a, b)

Due to the hexagonal crystal structure of MgB2 and the $s$-wave pairing, the
VL free energy for $\bm{H}\parallel\bm{c}$ can be expanded in term with
anisotropies that are multiples of six.Zhitomirsky and Dao (2004) Moreover,
the continuous rotation in the L phase implies that at least the six- and
twelvefold terms are sufficiently strong to influence the VL orientation, as
the transition would otherwise be discontinuous. Here, we have sought to
explore the evolution of the VL phase diagram as the twelvefold anisotropy is
suppressed by rotating the applied field away from the $c$ axis. We find that
the twelvefold anisotropy decreases linearly as the rotation angle is
increased, reducing the size of the L phase until it disappears entirely from
the phase diagram at a critical value. Above the critical angle, the VL
undergoes a first order phase transition directly from the F to I phase.

##  II Experimental Details

The VL was studied using small-angle neutron scattering (SANS) Mühlbauer  _et
al._ (2019). The SANS measurements were performed at the Bilby
instrumentSokolova  _et al._ (2019) at the Australian Nuclear Science and
Technology Organization. Preliminary SANS measurements were carried out at the
NG7 beam lineGlinka  _et al._ (1998) at the National Institute of Standards
and Technology Center for Neutron Research.

The experimental geometry used for the SANS measurements is shown in Fig.
1(b). Here the magnetic field is applied along the horizontal neutron beam
direction, and at an angle $\Omega$ relative to the crystalline $c$ axis
achieved by rotating the sample about the vertical axis in situ. The SANS
measurements used a neutron wavelength $\lambda_{n}=0.6$ nm and bandwidth
$\Delta\lambda_{n}/\lambda_{n}=10\%$. All measurements were performed at 2 K.

Measurements were performed on a 200 $\mu$g single crystal of MgB2 grown using
a high pressure cubic anvil technique.Karpinski  _et al._ (2003) The crystal
has a platelet geometry, roughly 1 mm $\times$ 1 mm wide and 50 $\mu$m thick
along the $c$ axis, and is isotopically enriched with 11B to reduce neutron
absorption. The superconducting critical temperature of the sample is
$T_{\text{c}}\approx 38$ K, and the upper critical field increases from
$H_{\text{c2}}=3.1$ T to $\sim\\!18$ T as the field is rotated from the $c$
axis ($\Omega=0$) to the basal plane ($\Omega=90^{\circ}$).Karpinski  _et al._
(2003)

##  III Results

Figure 2(a) shows a VL diffraction pattern obtained for
$\bm{H}\parallel\bm{c}$, with all six Bragg peaks lying on a circle of radius
$q=2\pi(2\mu_{0}H/\sqrt{3}\Phi_{0})^{1/2}$ where $\Phi_{0}=h/2e=2069$ T nm2 is
the flux quantum.

Figure 2:  Vortex lattice diffraction patterns at $\mu_{0}H=0.5$ T and (a)
$\bm{H}\parallel\bm{c}$ and (b) $\Omega=45^{\circ}$. The orientation of the
crystalline axes is indicated in (a). Both measurements were performed
following a field cooling from above $T_{\text{c}}$. Background scattering
near the detector center is masked off.

Here, the field and crystal have been rotated together through the Bragg
condition of all six VL peaks. The crystal was deliberately mounted with the
$\bm{a}$ axis roughly 10∘ from the vertical to investigate whether the
degeneracy of the two VL domain orientations in the L phase can be lifted by
rotating the applied field away from the $c$ axis. The field rotation
introduces a distortion of the VL due to the different penetration depth
within the basal plane vs perpendicular to it.Campbell  _et al._ (1988) This
is seen in the SANS diffraction patterns in Fig. 2(b) as a relocation of the
six Bragg peaks such that they lie on an ellipse with the minor axis parallel
to the axis of rotation. We note that the diffraction patterns in Fig. 2 were
obtained following a field cooling, which for $\bm{H}\parallel\bm{c}$ left the
VL in a metastable F phase.Das  _et al._ (2012)

To test the effect of the field rotation on the VL phase diagram, measurements
were made with $\Omega=0^{\circ},30^{\circ},45^{\circ},60^{\circ}$, and
$70^{\circ}$ and fields between 0.3 T and 1.2 T. Prior to each measurement,
the magnitude of the applied field was cycled through a damped oscillation
about the desired value with large enough amplitude (+5 mT, -4 mT, +3 mT, -2
mT, +1 mT) to overcome the activated behavior in MgB2 and ensure an
equilibrium VL configuration. For $\Omega=0$, this produces an equilibrium L
phase configuration at 0.5 T shown in Fig. 3(a) as opposed to the metastable F
phase configuration shown in Fig. 2(a).Levett  _et al._ (2002); Das  _et al._
(2012); Louden  _et al._ (2019a) As seen in Figs. 3(a)-3(d), the elliptical
distortion increases with $\Omega$.

Figure 3:  Vortex lattice diffraction patterns at (a) $H=0.5$ T and
$\Omega=0^{\circ}$, (b) 0.5 T and $45^{\circ}$, (c) 0.85 T and $45^{\circ}$,
and (d) 1.2 T and $70^{\circ}$. In all cases only a single (in (a) and (c)
split) peak fully satisfies the Bragg condition. The other peaks are indicated
by open circles in their predicted locations. For each diffraction pattern the
fitted geometric anisotropy ($\varepsilon$) is indicated. All measurements
were performed following a damped field oscillation. Background scattering
near the detector center is masked off. (e) Rotation of the VL Bragg peaks
around $q=0$ as a function of field for different $\Omega$. Angles are
measured with respect to the VL orientation at 0.3 T (F phase). Error bars
represent one standard deviation.

It is useful to quantify this distortion by the geometric anisotropy of the
ellipse $\varepsilon$, defined as the ratio of its major and minor axes. At
$\Omega=90^{\circ}$, $\varepsilon$ is expected to reach the penetration depth
anisotropy Campbell  _et al._ (1988). To conserve beam time not all VL peaks
were rocked through the Bragg condition for all measurements; however, their
location can be determined from symmetry and the analysis discussed below, and
are indicated by open circles in Figs. 3(a)-3(d). Within the L phase, split
Bragg peaks corresponding to the two degenerate VL orientations are observed
as seen in Figs. 3(a) and 3(c).

The angular rotation $\Delta\varphi$ of the VL Bragg peak as a function of
field and $\Omega$ is summarized in Fig. 3(e). At each $\Omega$ the rotation
is measured relative to the peak position in the F phase at 0.3 T,
corresponding to the projection of the $\bm{a}$ axis onto the scattering
plane, as shown in Figs. 3(c) and 3(d) and given by
$\tan\varphi^{\prime}=\tan\varphi\cdot\cos\Omega$. For $\Omega=0$, we observe
the same VL rotation reported earlier,Das  _et al._ (2012) associated with the
progression through the L phase. In the L phase, where the VL Bragg peaks are
split, we show the orientation of the clockwise rotating domain. As $\Omega$
is increased, the VL distortion allows the rotation of peaks which start near
the minor axis to exceed the $30^{\circ}$ range for the F to I transition when
$\bm{H}\parallel\bm{c}$. This is seen most clearly at $\Omega=70^{\circ}$,
where the VL rotation approaches $50^{\circ}$ at 1.2 T.

The location of the VL peak positions in the SANS data are governed by two
separate effects: the VL rotation transition within the L phase, and the
geometric distortion due to the penetration depth anisotropy discussed above.
To analyze the progression of the rotation transition it is useful to first
remove the effect of the geometric distortion, which corresponds to mapping
the VL Bragg peaks positions from lying on an ellipse back onto a circle. Due
to flux quantization, the area of the circle in reciprocal space must be the
same as the original ellipse, and the transformation can treated as a squeeze
mapping of all points $(q_{x}^{\prime},q_{y}^{\prime})$ in the circle to all
points $(q_{x},q_{y})$ in the ellipse:

$\left(\begin{array}[]{c}q_{x}^{\prime}\\\
q_{y}^{\prime}\end{array}\right)=\left(\begin{array}[]{cc}\varepsilon^{-1/2}&0\\\
0&\varepsilon^{1/2}\end{array}\right)\left(\begin{array}[]{c}q_{x}\\\
q_{y}\end{array}\right).$ (1)

Converting to polar coordinates $(q^{\prime},\varphi^{\prime})$ and
$(q,\varphi)$, and solving for $\varepsilon$ yields a transcendental equation
that can be solved numerically given any point on the ellipse:

$\varepsilon^{1/2}=\frac{q^{\prime}}{q\sin{\varphi}}\sin\left(\arccos\left[\frac{q\,\cos{\varphi}}{q^{\prime}\,\varepsilon^{1/2}}\right]\right).$
(2)

Values of $q$ and $\varphi$ are determined from 2D Gaussian fits to VL Bragg
peaks on the SANS detector, and $q^{\prime}$ is determined from the
measurement performed at the same applied field with $\Omega$ = 0. Once
$\varepsilon$ is determined, the re-scaled VL orientation $\varphi^{\prime}$
is given by

$\varphi^{\prime}=\arccos\left[\frac{q\,\cos{\varphi}}{q^{\prime}\,\varepsilon^{1/2}}\right].$
(3)

Finally, the rescaled VL rotations can be calculated by subtracting the VL
orientation at 0.3 T,

$\Delta\varphi^{\prime}(H,\Omega)=\varphi^{\prime}(H,\Omega)-\varphi^{\prime}(0.3~{}T,\Omega),$
(4)

and are shown in Fig. 4(a).

Figure 4:  Evolution of VL rotation transition versus applied field for
increasing field angle: (a) rescaled VL rotation corresponding to Fig. 3(e)
using Eqs. (2) – (4), (b) square root of the geometric distortion,
$\varepsilon$, including linear fits, and (c) anisotropy ratio, $\kappa$, (d)
inverse slope of the anisotropy ratio near $\kappa$ = 0 as a function of
$\Omega$. In all panels, error bars represent one standard deviation.

The transformation rescales the VL rotation so that the phase transition
always occurs in the range $\Delta\varphi^{\prime}=0^{\circ}$ to $30^{\circ}$
for all $\Omega$, consistent with the L phase behavior for
$\bm{H}\parallel\bm{c}$.

From Fig. 4(a) it is clear that the VL phase diagram changes dramatically with
increasing $\Omega$. Firstly, the onset of the rotation transition moves to
higher fields, indicating that the F phase
($\Delta\varphi^{\prime}=0^{\circ}$) is expanding. Secondly, once the rotation
does begin, it occurs more rapidly at high $\Omega$. Thirdly, the L phase
($0^{\circ}<\Delta\varphi^{\prime}<30^{\circ}$) either becomes very narrow or
vanishes entirely for $\Omega>70^{\circ}$. This would correspond to a
discontinuous phase transition directly from the F to the I phase
($\Delta\varphi^{\prime}=30^{\circ}$). The increasingly abrupt rotation
transition is also evident in Figs. 4(b), which shows an abrupt jump in the
slope of $\sqrt{\varepsilon}$ versus field for $\Omega\geq 60^{\circ}$. Here,
we note that the field dependence of $\varepsilon$ at constant $\Omega$ is due
to the gradual suppression of superconductivity on the $\pi$-band.Cubitt  _et
al._ (2003) This rise in $\varepsilon$ with respect to field is similar to
phenomena observed in V3Si,Christen  _et al._ (1985) 2$H$-NbSe2,Gammel  _et
al._ (1994) and Sr2RuO4.Kuhn  _et al._ (2017)

##  IV Discussion

The VL phase diagram in MgB2 can be modeled by a free energy

$\delta
F(\Delta\varphi^{\prime})=K_{6}\cos{(6\,\Delta\varphi^{\prime})}+K_{12}\cos{(12\,\Delta\varphi^{\prime})},$
(5)

containing six- and twelve-fold anisotropy terms. This form of the free energy
was originally proposed for $\bm{H}\parallel\bm{c}$ where
$\varphi^{\prime}=\varphi$,Zhitomirsky and Dao (2004); Das  _et al._ (2012);
Olszewski  _et al._ (2020) but here it is generalized to non-zero $\Omega$. In
this model $K_{12}$ is positive, and the continuous rotation transition in the
L phase occurs when the anisotropy ratio $\kappa=K_{6}/K_{12}$ is varied
between $-4$ and 4. Note that if $K_{12}$ was negative, the VL would undergo a
discontinuous reorientation transition between the F and I phases when $K_{6}$
changes sign. For each measurement it is possible to calculate the anisotropy
ratio

$\kappa=-4\cos{(6\,\Delta\varphi^{\prime})},$ (6)

obtained from a minimization of the free energy in Eq. (5). This is shown in
Fig. 4(c), and is seen to mimic the behavior of the re-scaled rotation angle
$\Delta\varphi^{\prime}$. We note that actual values of $\kappa$ most likely
extend outside the $\pm 4$ range allowed by Eq. (6). However, this does not
affect the analysis below as it focuses on the region where $\kappa$ is close
to zero.

The increasingly abrupt transition seen in Figs. 4(a) and 4(c) suggests a
reduction and possible sign change of $K_{12}$ as $\Omega$ increases. It is
possible to model the increasingly rapid reorientation transition by Taylor
expanding $K_{6}$ to linear order around the field $H_{0}$ where $K_{6}$, and
thus $\kappa$, vanish:

$K_{6}(H)\approx\alpha_{6}(H-H_{0}).$ (7)

Here, $\alpha_{6}$ is an undetermined constant. The slope of $\kappa$ with
respect to field near the transition at $H=H_{0}$ is then an indirect measure
of $K_{12}$:

$\left.\frac{\partial\kappa}{\partial
H}\right|_{H=H_{0}}=\frac{\alpha_{6}}{K_{12}},$ (8)

where we have taken $K_{12}$ to be field independent. Figure 4(d) shows the
inverse slope $(\partial\kappa/\partial H)^{-1}$ measured near $\kappa=0$ for
each $\Omega$. The inverse slope depends linearly on $\Omega$, suggesting a
functional form for $K_{12}$:

$K_{12}(\Omega)=\alpha_{12}(\Omega-\Omega_{0}).$ (9)

Extrapolating to $K_{12}/{\alpha_{6}}=0$, we find an estimate for
$\Omega_{0}=\left(83^{+3}_{-6}\right)^{\circ}$, where the error bars both
represent one standard deviation. For values of $\Omega>\Omega_{0}$, $K_{12}$
will be negative and the VL will undergo a first-order transition between the
F and I phases. As a result, the L phase vanishes and $H_{0}$ becomes a
critical field corresponding to an equal F and I phase energy (rather than
simply the field where $\Delta\varphi^{\prime}=15^{\circ}$ which is the case
for $\Omega<\Omega_{0}$). We note that $\Omega=70^{\circ}$ may already be
above $\Omega_{0}$, as the uncertainly on $(\partial\kappa/\partial H)^{-1}$
at this angle extends to negative values. Behavior qualitatively similar to
MgB2 has also been observed in YBCO,Yethiraj  _et al._ (1993) where a
triangular to square VL transition switches from second to first order at a
critical angle of $10^{\circ}$. The much smaller $\Omega_{0}$ for YBCO is
likely due to the highly two-dimensional nature of this material, making it
relatively more susceptible to field rotations away from the $c$ axis.

Finally, we return to the question of whether rotating the applied field away
from the $c$ axis will break the VL domain degeneracy within the L phase. Such
an effect was previously observed in TmNi2B2C where the VL undergoes a
triangular to square transition with a degenerate intermediate rhombic
phase,Eskildsen  _et al._ (1998) and where rotating the field away from the
$c$ axis by $\approx 10^{\circ}$ is sufficient to suppress one of the two
rhombic domains.Dewhurst  In an ideal uniaxial superconductor with an
isotropic basal plane, London theory predicts that the two-fold anisotropy
introduced by field rotation will favor a VL orientation with Bragg peaks on
the minor axis of the ellipse in reciprocal space.Campbell  _et al._ (1988)
However, any real material will exhibit some basal plane anisotropy which may
compete with the uniaxial effect, and the relative strength of the two will
determine the VL orientation. As an example, one can consider NbSe2, where the
triangular VL is oriented in a manner corresponding to the maximal energy
according to the London theory.Gammel  _et al._ (1994)

In our SANS measurements the MgB2 crystal was mounted such that within the L
phase Bragg peaks for one of the domains lie near the vertical axis, as seen
in Fig. 3(a), and are thus favored by the uniaxial effect. However, peaks
corresponding to both L phase domains were clearly observed at
$\Omega=45^{\circ}$and 0.6 T and at $60^{\circ}$ and 0.9 T, indicating that
the degeneracy is not readily lifted in MgB2. This suggests that at these
values of $\Omega$ the uniaxial anisotropy is weaker than $K_{12}$, and that a
suppression of the degeneracy will only occur just below $\Omega_{0}$. It is
also consistent with the inability of the uniaxial effect to cause a VL
reorientation in the F phases to have the Bragg peaks on the minor axis of the
ellipse, since $K_{6}$ and $K_{12}$ are of the same order of magnitude except
in the immediate vicinity of $H_{0}$ and $\Omega_{0}$. That said, the uniaxial
anisotropy does appear to lower the free energy of the F phase relative to the
I phase, causing $H_{0}$ to shift to higher fields in Figs. 4(a) and 4(c) as
$\Omega\rightarrow\Omega_{0}$. This conclusion is supported by previous SANS
measurements where the magnetic field was rotated about the $\bm{a}^{*}$.Pal
_et al._ (2006) Here the L phase (rather than the F phase) is favored by the
uniaxial anisotropy, and at 0.4 T the critical angle is shifted to a lower
value $\Omega_{0}\sim 70^{\circ}$.

##  V CONCLUSION

In summary, we have demonstrated how the continuous (second order) VL rotation
transition observed in MgB2 for $\bm{H}\parallel\bm{c}$ gradually evolves
towards and finally becomes discontinuous (first order) as the applied
magnetic field is rotated away from the $c$ axis by $\approx 83^{\circ}$. For
rotation angles below this critical value, domain formation in the
intermediate L phase persists. We speculate that the long lived metastable VL
phases, attributed to domain formation in the L phase, will thus no longer
occur above the critical angle, and should be the subject of further SANS
studies.

## Acknowledgements

This work was supported in part by the U.S. Department of Energy, Office of
Basic Energy Sciences, under Award No. DE-SC0005051 (A.W.D.L, M.R.E). We
acknowledge the support of the National Institute of Standards and Technology,
U.S. Department of Commerce, in providing the neutron research facilities used
in this work. We acknowledge useful conversations with G. Longbons and A.
Francisco in preparation for these experiments, as well as the assistance of
J. Krzywon, A. Qiang, and C. Baldwin in completing them.

## References

  * Kogan (1981) V. G. Kogan, Phys. Rev. B 24, 1572 (1981).
  * Kogan  _et al._ (1997) V. G. Kogan, M. Bullock, B. Harmon, P. Miranovic, L. Dobrosavlevic-Grujic, P. L. Gammel, and D. J. Bishop, Phys. Rev. B 55, R8693 (1997).
  * Laver  _et al._ (2006) M. Laver, E. M. Forgan, S. P. Brown, D. Charalambous, D. Fort, C. Bowell, S. Ramos, R. J. Lycett, D. K. Christen, J. Kohlbrecher, C. D. Dewhurst, and R. Cubitt, Phys. Rev. Lett. 96, 167002 (2006).
  * Laver  _et al._ (2009) M. Laver, C. J. Bowell, E. M. Forgan, A. B. Abrahamsen, D. Fort, C. D. Dewhurst, S. Mühlbauer, D. K. Christen, J. Kohlbrecher, R. Cubitt, and S. Ramos, Phys. Rev. B 79, 014518 (2009).
  * Mühlbauer  _et al._ (2009) S. Mühlbauer, C. Pfleiderer, P. Böni, M. Laver, E. M. Forgan, D. Fort, U. Keiderling, and G. Behr, Phys. Rev. Lett. 102, 136408 (2009).
  * Gannon  _et al._ (2015) W. J. Gannon, W. P. Halperin, C. Rastovski, K. J. Schlesinger, J. Hlevyack, M. R. Eskildsen, A. B. Vorontsov, J. Gavilano, U. Gasser, and G. Nagy, New J. Phys. 17, 023041 (2015).
  * Avers  _et al._ (2020) K. E. Avers, W. J. Gannon, S. J. Kuhn, W. P. Halperin, J. A. Sauls, L. DeBeer-Schmitt, C. D. Dewhurst, J. Gavilano, G. Nagy, U. Gasser, and M. R. Eskildsen, Nat. Phys. 16, 531 (2020).
  * Cubitt  _et al._ (2003) R. Cubitt, M. R. Eskildsen, C. D. Dewhurst, J. Jun, S. M. Kazakov, and J. Karpinski, Phys. Rev. Lett. 91, 047002 (2003).
  * Das  _et al._ (2012) P. Das, C. Rastovski, T. R. O’Brien, K. J. Schlesinger, C. D. Dewhurst, L. DeBeer-Schmitt, N. D. Zhigadlo, J. Karpinski, and M. R. Eskildsen, Phys. Rev. Lett. 108, 167001 (2012).
  * Hirano  _et al._ (2013) T. Hirano, K. Takamori, M. Ichioka, and K. Machida, J. Phys. Soc. Japan 82, 063708 (2013).
  * Eskildsen  _et al._ (2002) M. R. Eskildsen, M. Kugler, S. Tanaka, J. Jun, S. M. Kazakov, J. Karpinski, and Ø. Fischer, Phys. Rev. Lett. 89, 187003 (2002).
  * Rastovski  _et al._ (2013) C. Rastovski, K. J. Schlesinger, W. J. Gannon, C. D. Dewhurst, L. DeBeer-Schmitt, N. D. Zhigadlo, J. Karpinski, and M. R. Eskildsen, Phys. Rev. Lett. 111, 107002 (2013).
  * Louden  _et al._ (2019a) E. R. Louden, C. Rastovski, S. J. Kuhn, A. W. D. Leishman, L. DeBeer-Schmitt, C. D. Dewhurst, N. D. Zhigadlo, and M. R. Eskildsen, Phys. Rev. B 99, 060502(R) (2019a).
  * Louden  _et al._ (2019b) E. R. Louden, C. Rastovski, L. DeBeer-Schmitt, C. D. Dewhurst, N. D. Zhigadlo, and M. R. Eskildsen, Phys. Rev. B 99, 144515 (2019b).
  * Zhitomirsky and Dao (2004) M. E. Zhitomirsky and V. H. Dao, Phys. Rev. B 69, 054508 (2004).
  * Mühlbauer  _et al._ (2019) S. Mühlbauer, D. Honecker, È. A. Pèrigo, F. Bergner, S. Disch, A. Heinemann, S. Erokhin, D. Berkov, C. Leighton, M. R. Eskildsen, and A. Michels, Rev. Mod. Phys. 91, 015004 (2019).
  * Sokolova  _et al._ (2019) A. Sokolova, A. E. Whitten, L. de Campo, J. Christoforidis, A. Eltobaji, J. Barnes, F. Darmann, and A. Berry, J. Appl. Crystallogr. 52, 1 (2019).
  * Glinka  _et al._ (1998) C. J. Glinka, J. G. Barker, B. Hammouda, S. Krueger, J. J. Moyer, and W. J. Orts, J. Appl. Cryst. 31, 430 (1998).
  * Karpinski  _et al._ (2003) J. Karpinski, M. Angst, J. Jun, S. M. Kazakov, R. Puzniak, A. Wisniewski, J. Roos, H. Keller, A. Perucchi, L. Degiorgi, M. R. Eskildsen, P. Bordet, L. Vinnikov, and A. Mironov, Supercond. Sci. Technol. 16, 221 (2003).
  * Campbell  _et al._ (1988) L. J. Campbell, M. M. Doria, and V. G. Kogan, Phys. Rev. B 38, 2439 (1988).
  * Levett  _et al._ (2002) S. J. Levett, C. D. Dewhurst, and D. M. Paul, Phys. Rev. B 66, 014515 (2002).
  * Christen  _et al._ (1985) D. K. Christen, H. R. Kerchner, S. T. Sekula, and Y. K. Chang, Physica B 135, 369 (1985).
  * Gammel  _et al._ (1994) P. L. Gammel, D. A. Huse, R. N. Kleiman, B. Batlogg, C. S. Oglesby, E. Bucher, D. J. Bishop, T. E. Mason, and K. Mortensen, Phys. Rev. Lett. 72, 278 (1994).
  * Kuhn  _et al._ (2017) S. J. Kuhn, W. Morgenlander, E. R. Louden, C. Rastovski, W. J. Gannon, H. Takatsu, D. C. Peets, Y. Maeno, C. D. Dewhurst, J. Gavilano, and M. R. Eskildsen, Phys. Rev. B 96, 174507 (2017).
  * Olszewski  _et al._ (2020) M. W. Olszewski, M. R. Eskildsen, C. Reichhardt, and C. J. O. Reichhardt, Phys. Rev. B 101, 224504 (2020).
  * Yethiraj  _et al._ (1993) M. Yethiraj, H. A. Mook, G. D. Wignall, R. Cubitt, E. M. Forgan, S. L. Lee, D. M. Paul, and T. Armstrong, Phys. Rev. Lett. 71, 3019 (1993).
  * Eskildsen  _et al._ (1998) M. R. Eskildsen, K. Harada, P. L. Gammel, A. B. Abrahamsen, N. H. Andersen, G. Ernst, A. P. Ramirez, D. J. Bishop, K. Mortensen, D. G. Naugle, K. D. D. Rathnayaka, and P. C. Canfield, Nature 393, 242 (1998).
  * (28) C. D. Dewhurst, private communication.
  * Pal  _et al._ (2006) D. Pal, L. DeBeer-Schmitt, T. Bera, R. Cubitt, C. D. Dewhurst, J. Jun, N. D. Zhigadlo, J. Karpinski, V. G. Kogan, and M. R. Eskildsen, Phys. Rev. B 73, 012513 (2006).


Influence of asteroid shape elongation during a distant planetary encounter
Y Kim et al.

mode = title]The surface sensitivity of rubble-pile asteroids during a distant planetary encounter: Influence of asteroid shape elongation

1]Yaeji Kim[type=editor,
[1]

[1]Department of Aerospace Engineering, Auburn University, 211 Davis Hall, Auburn, AL 36849-5338, USA

1]Masatoshi Hirabayashi

2]Richard P. Binzel

3]Marina Brozović

4]Daniel J. Scheeres

5]Derek C. Richardson

[2]Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139 USA

[3]Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Mail Stop 301-121, Pasadena, CA 91109-8099, USA

[4]Aerospace Engineering Department, University of Colorado Boulder, 429 UCB, Boulder, CO 80309, USA

[5]Department of Astronomy, University of Maryland, College Park, MD 20742, USA

* Corresponding author

We numerically investigate how an asteroid's elongation controls the sensitivity of its surface to tidal effects during a distant planetary encounter beyond the Roche limit. We analyze the surface slope and its variation by considering the shape elongation, as well as the spin period and orbital conditions. A more elongated asteroid tends to have a higher slope variation, while there may not be a monotonic increase of the total area having such a variation.

* We numerically investigate the influence of an asteroid’s elongation on the sensitivity of its surface to tidal effect during a distant planetary encounter beyond the Roche limit.

* Objects spinning rapidly, but below the spin barrier, tend to experience a dramatic change of the surface slope during the distant encounters.

* A more elongated asteroid tends to have a higher slope variation, while the asteroid elongation causes a non-linear increase in the total area having such a variation.

Asteroids, dynamicsAsteroids, surfaces Geological processes Tides, solid body

§ INTRODUCTION

S- and Q-type asteroids exhibit the compositional features of ordinary chondrites, while their spectral slope and absorption band are slightly different. S-type asteroids have steeper spectral slopes, weaker absorption bands near wavelengths between 1 $\mu$m and 2 $\mu$m, and lower albedos than Q-type asteroids [3, 7, 33, 48]. Space weathering driven by ion implantation and micrometeorite bombardments cause optical variations in material compositions [7, 38]. For S- and Q-type asteroids, space weathering is considered to redden their surface materials. The timescale of space weathering may range between 10 ka and 1 Ma in the inner Solar System [16, 45]. For near-Earth asteroids (NEAs), the dynamical lifetime is 2 Ma [5] and thus longer than the space weathering timescale. Without additional processes, Q-type asteroids should be altered to S-type asteroids and depleted continuously. However, this contradicts the high fraction of Q-type asteroids among NEAs ($\sim25\%$). One possible hypothesis is that redder S-type asteroids are resurfaced to become bluer Q-type asteroids <cit.>.

Possible resurfacing mechanisms include 1) planetary encounters with terrestrial planets [3, 6, 11, 23, 30, 29], 2) thermal fatigue by cyclic diurnal temperature variations [9, 25, 26], and 3) rotational instability driven by the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect [15]. Among these mechanisms, the tidal effect may be a critical driver that resurfaces rubble-pile asteroids. The distribution of Q-type asteroids correlates with the perihelion distance and the minimum orbit intersection distance, suggesting that encounters with massive planets may resurface asteroids enough to change their surface spectral properties [3, 23, 29]. [29] proposed that Q-type asteroids may be resurfaced if the close encounter is within $\sim5$ planetary radii, thus outside the Roche limit ($\sim3.4$ planetary radii for an asteroid with a bulk density of 2.0 g cm$^{-3}$). These studies statistically showed that the resurfacing process of S-type asteroids is possibly related to their distant planetary encounters. However, the resurfacing mechanism is still not well understood.

Numerical studies have been reported for investigations of the tidal effects on rubble-pile asteroids during planetary encounters. They are in general divided in two categories; catastrophic, where an object is destroyed by the tides, and non-catastrophic, where either shape or surface or both are altered. For catastrophic disruptions, research has shown how small bodies that closely encounter planetary bodies are broken up due to strong tidal forces. The breakup conditions depend on small bodies' mechanical strength, material compositions, and orbital parameters <cit.>. In 1994, comet Shoemaker-Levy 9 encountered Jupiter and was broken into approximately 20 similar-sized fragments [42]. [41] and [1] modeled this event to better quantify catastrophic disruption processes during close tidal encounters.

Tidal processes during distant encounters outside the Roche limit have also been studied in the literature. [19, 20] applied the theory of hill slope stability to evaluate the stability of asteroid regolith during the distant planetary flyby. They had two conclusions: rapidly rotating asteroids are highly likely to experience resurfacing, and the resurfacing process on asteroids may occur at larger flyby distance ($\sim$10 Earth Radii) than the previous estimation ($\sim$5 Earth Radii). In this paper, we extend their works to gain a more comprehensive understanding of resurfacing by adding an element that had not been discussed in the earlier works: an asteroid's elongation. We finally note that studies have shown limited changes in the shape and surface conditions of Apophis, which is going to closely flyby the Earth within 6 Earth radii [10, 40, 52]. Our study will give further insight into the mechanism of mild tidal effects on a rubble-pile asteroid during a distant planetary encounter, which may be a critical source of Q-type asteroids.

In this study, we investigate the influence of a rubble-pile asteroid's elongation on surface sensitivity to tidal effects when it approaches the Earth outside the Roche limit by parameterizing the elongation with different planetary encounter conditions. We first demonstrate the resurfacing mechanism by considering the surface slope and its variation in Section <ref>. Our investigation and simulation settings are described in Sections <ref> and <ref>. In Section <ref>, the results of all conducted simulations are presented. We discuss our findings, potential issues of the current numerical model, and future work in Section <ref>.

§ RESURFACING MECHANISM

Resurfacing processes remove a space-weathered, redder surface layer on S-type asteroids by exposing fresh materials beneath it. In situ observations by spacecraft revealed that the space weathered layer is likely very thin and correlated with the surface topography. Despite limited knowledge of space weathering on S-type asteroids, the sampled particles from NEA (25143) Itokawa suggested that the weathered thickness is only $\sim80$ nm [32, 31], indicating that only a very-thin top surface layer had been affected by space weathering. We thus speculate that even mass movements at a tiny scale can induce color variations in top surfaces. Furthermore, [46, 28] reported on bluer, likely unweathered regions around the pole regions and the equatorial ridges of NEA (162171) Ryugu (C-type asteroid), although the space weathering mechanism may be different from that on S- and Q-type asteroids [22].

We use two geophysical parameters to describe the surface sensitivity to possible granular flows: the surface slope and its variation. The surface slope describes how the surface element is tilted from the direction of its total acceleration combined with local gravity, the tidal effect, and the rotational effect on each surface element. The slope variation shows how the surface is affected by time-varying acceleration during the encounter. Using these parameters, we estimate the locations of where granular flows can occur. When the surface slope reaches its critical slope, i.e., the angle of repose, which is $35^\circ$ for a typical geological material without cohesion [21], granular flows may start to occur [8, 37]. Furthermore, if the slope variation is high enough to cause surface mobility, granular flows may occur for sub-critical slopes. [2] showed this mechanism operating on the surface of the Martian moon Phobos, leading to the resurfacing process; however, this mechanism may not be proper for the planetary encounter resurfacing because it is a non-periodic event. In our analysis, we use the slope variation to see how sensitive surfaces become during the distant planetary encounter.

Mass movements in and on asteroids may be influenced by many different elements. Electrostatic forces may induce mass mobility [17]. Particle ejection may also occur by thermal fatigue [27]. Impact cratering can excavate fresh materials. Granular convection such as the Brazil-nut effect can transport fresh materials from the interior [47, 50]. Seismic wave propagation may enhance resurfacing [36, 51]. The consideration of these effects on mass movements is beyond our scope. Instead, we focus on surface flows driven by the tidal force, as well as the gravitational force and the centrifugal force to better understand how an asteroid's elongation affects the tidally induced resurfacing mechanism.

§ MODELING

Our model computes the surface slope ($\theta$) by using the following equation,

\begin{equation}
\theta = \arccos \left ( \frac{- \vec{n} \cdot \vec{a_g}}{\| \vec{n} \| \|\vec{a_g} \|} \right ).
\end{equation}

where $\vec{a_g}$ is the net acceleration vector, and $\vec{n}$ is a normal vector to a surface element. The surface slope variation ($\delta\theta$) is the angle difference of a given element between the maximum surface slope during the flyby and the initial slope before the flyby. This quantity was used to analyze resurfacing on the martian moon, Phobos [2].

Next, we discuss the net acceleration vector, $\vec{a_g}$, during a planetary encounter. Since we focus on the onset of the grain motion on an asteroid surface, particles are assumed to rest initially. The assumption yields $\dot {\vec{r}} = 0$ and thus excludes the Coriolis effect. This setting may provide a conservative condition for the occurrence of resurfacing. $\vec{a_g}$ is given as

\begin{equation}
\vec{a}_g = -G\rho_A\int_{V_A}\frac{\Vec{r}}{r^3} dV_A - \frac{GM_E}{{R_c}^3}(\Vec{r}-\frac{3({\Vec{R}_c}+\Vec{r})\cdot\Vec{r}}{R_c^2}{\Vec{R}}_c) - \dot{\vec{{\omega}}}\times\vec{r} - \Vec{\omega}\times(\Vec{\omega}\times\Vec{r}).
\label{Eq:acc}\end{equation}

In this equation, $G$ is the gravitational constant, $\rho_A$ is the asteroid's bulk density, $V_A$ is the asteroid's volume, $\vec{r}$ is the position of a surface element relative to the center of mass (COM) of the asteroid, $M_E$ is the mass of the planet, $\vec{R_c}$ is the position of the COM of the asteroid relative to that of the planet, $\Vec{\omega}$ is the asteroid's angular velocity vector, and $\dot{\Vec{\omega}}$ is the asteroid's angular acceleration vector. The first term is the acceleration by self-gravity, the second term describes the tidal acceleration, and the third and fourth terms represent the rotational effect. For self-gravity, we use an elliptical integral for computing the gravity field around a biaxial ellipsoid [18].

For the rotational motion, we solve the following equation:

\begin{equation}
[I]\Dot{\Vec{\omega}} + \Vec{\omega}\times[I]\Vec{\omega} = \frac{3GM_E}{R^5_c}\Vec{R_c}\times[I]\Vec{R_c}. \label{Eq:rot1}
\end{equation}

where $[I]$ is the moment of inertia matrix and computed by [13]. We also note the changes in the asteroid's orientation during the encounter. The transformation matrix, $[A]$, is given as

\begin{equation}
[\Dot{A}] = -[\Tilde{\omega}][A]. \label{Eq:rot2}
\end{equation}

The initial condition of [A] is fixed as an identity matrix. The evolution of [A] is eventually used to convert $\Vec{R_c}$, which is a vector in the planet-centered inertia reference frame, to the one in the body-fixed frame in Equation (<ref>). Here, $[\tilde \omega]$ is the skew matrix of $\vec{\omega}$, which is given as

\begin{equation}
\Tilde{\omega} = \begin{bmatrix}
0 & -\omega_3 & \omega_2 \\
\omega_3 & 0 & -\omega_1 \\
-\omega_2 & \omega_1 & 0
\end{bmatrix}
\label{Eq:rot3}
\end{equation}

We use a Runge-Kutta 4th order integrator with a constant time-step to propagate $\vec{\omega}$ and $[A]$ over a considered time frame. $\vec{R_c}$ is described in the body-fixed frame centered at the COM of the planetary body, where the axes are aligned with the asteroid's moment of inertia axes. To determine $\vec{R_c}$, we design a hyperbolic trajectory with a given periapsis and a constant encounter velocity and obtain a list of $\Vec{R_c}$ within a considered time span with a proper time step (see Section <ref>). This list is incorporated into the integrator that solves Equations (<ref>) and (<ref>) to describe the rotational motion. At last, all computed parameters are applied to Equation (<ref>) to obtain $\vec{a}_g$ of each element on the asteroid. We developed a MATLAB program package for this simulation process.

§ SIMULATION SETTINGS

Physical parameters and shape model information used in our simulations.

Parameter     Values

Earth-like planet mass ($M_E$)     5.97 $\times$ $10^{24}$ $kg$

Asteroid's volume ($V_A$)     1.98 $\times$ $10^{7}$ $m^3$

Asteroid's bulk density ($\rho_A$)     2.0 $g$ $cm^{-3}$

Simulation length     < 20 $Earth Radii$

Simulation time step     5 $sec$

Encounter speed     15 $km$ $s^{-1}$

Angle of repose of top surface layers     35$^\circ{}$

Number of facets      10440

Mean facet size     35 $m^{3}$ ($\sim$0.01$\%$ of the entire surface)

We focus on the slope variation on uniformly rotating biaxial ellipsoids with different aspect ratios ($AR$). Here, $AR$ is defined as a ratio of the major axis to the minor axis. The volume of the shape models is set identical to that of a sphere having an equivalent diameter of 340 m, which is a typical size of rubble pile asteroids. A bulk density is assumed to be a constant value of 2.0 g cm$^{-3}$, which is consistent with that of Itokawa (S-type)  [14]. In each hyperbolic trajectory, the initial and final points of a test asteroid are 20 Earth-like planet radii ($\sim$6371 km for one Earth-like planet radius) away from the planet center, and thus the tidal effect is negligible at those points. The encounter speed is set to 15 km s$^{-1}$ for all the cases [4], while we consider three hyperbolic trajectories with different periapses in the range of plausible resurfacing distances (3.5, 5, and 10 planet radii) [19]. The total length of simulation is approximately 4.5 hours with a time step of 5 seconds. In all test cases, the body is set to rotate along their maximum moment of inertia axes, and its spin axis is perpendicular to the orbital plane. The angle of the repose is set to be 35$^\circ$.

§ RESULTS

We introduce the nomenclature used in this session. To characterize the area affected by high slope variation ($\delta \theta$), we introduce two characteristic regions: Regions I and II. Region I is an area whose surface slope is higher than the angle of repose of $35^\circ$ before the encounter. Region II, on the other hand, is an area that has the slope exceeding the angle of repose during the encounter. We also use the area difference between Regions I and II to indicate the surface regions that reach above the angle of repose due to the tidal effect. A total of 33 simulations are performed to analyze how $\delta \theta$ is affected by the asteroid elongation during distant encounters.

We consider three cases of the spin periods, 1.5, 2.8, and 3.1 hr. These spin periods are selected to demonstrate high $\delta \theta$ variations around the spin barrier, where gravitational aggregates may fall apart if they do not have tensile strength [34]. Given the spin periods defined, the test body may have high slopes even before a tidal encounter (i.e., Region I) in some cases. The 1.5 hr spin period case causes all the $AR$ cases to have Region I across the entire surface. Highly elongated shapes with lower $AR$s have Region I at even longer spin periods. At a spin period of 2.8 hr (see Figure <ref>), shapes with $AR=0.8$ and $AR=0.57$ have Region I at middle latitudes around the longest axis. At a spin period of 3.1 hr, the $AR=0.57$ shape still has Region I, which would be sensitive to granular flow. We emphasize that the reason for this selection is to illustrate the mechanisms of the slope variation due to the elongation and show its transition around the spin barrier. Thus, these conditions allow for visualizing how shape elongation controls surface slope variations, enhanced by rotation, during a distant encounter.

To show how the elongation affects $\delta \theta$, Figure <ref> illustrates sample cases in which asteroids with $AR$s of 1.0, 0.8, and 0.57 reach 5 Earth radii with a spin period of 2.8 hr. The results show that $\delta \theta$ becomes higher as $AR$ is lower. When the asteroid has $AR=1$, high $\delta\theta$ occurs at a latitude of $\sim$18 deg. The high $\delta\theta$ regions always face the planet because of its axisymmetric shape. When $AR=0.8$, high $\delta \theta$ regions tend to be distributed widely around the edges along the longest axis. When $AR=0.57$, higher $\delta \theta$ concentrates on the edges.

Maximum surface slope and slope variations for each shape. The top panels (a, c, and e) show the surface slopes, which are affected by the tidal effect, during the planetary encounter with the periapsis of 5 Earth radii and the spin period of 2.8 hr. The bottom panels (b, d, and f) illustrate the maximum surface slope variation during this encounter. The aspect ratio of each column is 1, 0.8, and 0.57, respectively. All plots are projected from a longitude of $130^\circ$ and a latitude of $15^\circ$. The black arrows denote the spin axes.

Surface slope data. $\delta\theta$ shows the minimum and maximum value at Region II, excluding Region I. If this value is unchanged (i.e., Region I = Region II), $\delta\theta$ shows the slope variation across the entire surface with a parenthesis. The increase rate defines the area difference between Region II and Region I. For the 2.8 hr case, five $AR$ cases are introduced to illustrate the variations in the increase rate.

2CSimulation Setting     2CPre-encounter Time
3CPlanetary Encounter Time
1-2 3-4 5-7
1C$AR$   1CPeriapsis ($R_{E}$)   1CMean slope ($^\circ$)   1CRegion I ($\%$)   1CRegion II
($\%$)   1CIncrease Rate
($\%$)   1C$\delta\theta$ ($^\circ$)

2LSpin period - 1.5 hr

1.0     3.5     83.2     97     97.4     0.4     0.47 - 1.94

0.8     3.5     82.7     96.2     96.2     0     (0.20 - 4.68)

0.57     3.5     82.3     96.6     96.6     0    2.50 - 3.17

1.0     5     83.2     97     97.4     0.4     0.31 - 0.69

0.8     5     82.7     96.2     96.2     0     (0.06 - 1.62)

0.57     5     82.3     96.6     96.6     0     0.63 - 0.80

1.0     10     83.2     97     97.4     0.4     0.06 - 0.09

0.8     10     82.7     96.2     96.2     0     (0.01 - 0.23)

0.57     10     82.3     96.6     96.6     0     (0.01 - 0.30)

2LSpin period - 2.8 hr

1.0     3.5     22.3     0     22.7     22.7     4.12 - 16.23

0.9     3.5     22.7     0     31.5     31.5     3.13 - 19.89

0.8     3.5     24     10.2     46.9     36.7     5.12 - 34.48

0.7     3.5     26.7     22.1     55     32.9     6.89 - 148.83

0.57     3.5     33.1     33.7     63.4     29.7     9.14 - 38.41

1.0     5     22.3     0     0.7     0.7     3.96 - 4.22

0.9     5     22.7     0     6.9     6.9     1.62 - 4.67

0.8     5     24     10.2     20.9     10.7     2.21 - 6.92

0.7     5     26.7     22.1     33.2     11.1     2.88 - 34.24

0.57     5     33.1     33.7     43.6     9.9     3.80 - 11.22

1.0     10     22.3     0     0     0     (0.01 - 0.53)

0.9     10     22.7     0     0     0     (0.01 - 0.74)

0.8     10     24     10.2     11.8     1.6     0.31 - 0.75

0.7     10     26.7     22.1     23.5     1.4     0.34 - 2.66

0.57     10     33.1     33.7     34.8     1.1     0.39 - 0.98

2LSpin period - 3.1 hr

1.0     3.5     16.6     0     0     0      (0.04 - 13.08)

0.8     3.5     16.9     0     0     0      (0.10 - 11.08)

0.57     3.5     22.5     15.4     29.7     14.3     7.13 - 21.88

1.0     5     16.6     0     0     0      (0.02 - 3.72)

0.8     5     16.9     0     0     0      (0.02 - 4.73)

0.57     5     22.5     15.4     17.7     2.3     2.41 - 8.36

1.0     10     16.6     0     0     0      (0.01 - 0.39)

0.8     10     16.9     0     0     0      (0.01 - 0.60)

0.57     10     22.5     15.4     16.8     1.4     0.78 - 1.72

Table <ref> shows the results of the cases simulations. We first discuss the cases of 1.5 hr and 3.1 hr. When the spin period is 1.5 hr, since the centrifugal force is already dominant, Region I spreads over the almost entire regions for any shapes. This extreme condition does not cause the slope to change significantly. When the spin period is 3.1 hr, self-gravity starts to be dominant, and tidal acceleration does not influence the slope variation remarkably. The $AR=0.57$ case still shows non-zero area difference between Region II and Region I at the asteroid's edge along the longest axis, where self-gravitational acceleration and rotational acceleration are comparable.

When the spin period is 2.8 hr, rotational acceleration becomes comparable to self-gravitational acceleration while some locations on elongated shapes may still be dominantly influenced by centrifugal forces. We show more simulation cases for this spin period case in Table <ref>, compared to other spin periods, to examine the slope variation depending on the elongation. Region II becomes wider than Region I, strongly depending on the elongation. With the decrease of $AR$, the area difference between Region II and Region I (see the Increase Rate in Table <ref>) increases in the range of $AR$ from 1 to $\sim$0.8 but decreases from $\sim$0.8 to 0.57. Consider when the periapsis is 3.5 Earth radii. For $AR = 1$, Region II is $\sim$23% broader than Region I. The area difference becomes $\sim$37% for $AR=0.8$ and then gradually decreases to $\sim$30% for $AR=0.57$, which is still higher than the spherical case. This similar trend is observed when the periapsis is 5 and 10 planet radii. The decrease of the area difference from $AR=0.8$ to $AR=0.57$ results from the fact that a more elongated shape has a limited region in which a high slope variation occurs because of its narrower edge. As the periapsis becomes distant, tidal acceleration becomes smaller, leading to the decrease of the slope variation.

Our results are consistent with [19, 20, 52] in terms of the influence of the spin period and orbit on $\delta \theta$. Importantly, we newly address how the elongation impacts the slope variation during distant encounters. If the shape is moderately elongated, broader areas may have reasonable slope variations, which may cause the surface slope to reach the angle of repose. If the shape is highly elongated, the slope variation becomes higher although the affected area may be limited. While our parameter analysis is still coarse, the 2.8 hr case implies that $AR=0.8$ has the most noticeable change between Regions II and I because the total acceleration of rotation and self-gravity becomes small across the broader regions, the tidal acceleration can easily change their slope conditions.

§ DISCUSSION

We numerically showed how the shape elongation affects the slope variations and the location of sensitive regions during distant encounters. We observed two critical features. First, the shape elongation enhances the slope variation. High slope variations are observed when the shape is elongated. When an asteroid spins at a spin period at which centrifugal and gravity acceleration are comparable, tidal acceleration can dramatically change the surface slope. Thus, an elongated body tends to be more exposed to such a dramatic change. Second, the increase of the area whose surface slope exceeds 35$^\circ$ during a distant encounter (Region II - Region I) depends on shape elongation. As the shape is more elongated from the spherical condition, the area that reaches the angle of repose of 35$^\circ$ increases. However, after the shape reaches a certain elongation (in our simulation, it was $\sim$0.8), the area difference, Region II - Region I, starts decreasing due to the limited area of high slope variations. This trend is reasonable because the slope condition on a highly elongated body is clearly separated into gravitational and centrifugal acceleration-dominant regions along the longest axis, causing "fewer" areas to have high slope variations. From these two trends, we conclude that elongation is a strong contributor to the surface slope variation during a distant encounter.

The shape elongation of NEAs has been reported by observational investigations. The statistical characterizations of light-curve amplitudes by [24] showed that the shape elongation of NEAs may be between 0.6 and 0.8 if the objects are assumed to be a population of prolate spheroids. Furthermore, there is an apparent excess of fast rotators in the NEAs with D > 200 m near the spin barrier. Our study implies that if the fast spinning asteroids experience distant encounters, high slope variations may occur, and their elongation may control the affected area. This interpretation is consistent with the arguments of earlier works <cit.>, and we emphasize that our study newly analyzes the influence of the elongation on slope variations.

We point out potential issues with our current numerical model. First, this model does not take into account cohesion, which may prevent resurfacing at fast spin but enhance its magnitude if resurfacing occurs. However, the surface slope without cohesion is still a meaningful parameter based on recent observations <cit.> that thin top-surface layers of an asteroid may be covered with weak, cohesionless materials. Second, the shape of an asteroid is assumed to be rigid in our model; in other words, its deformation is ignored. If deformation occurs, the location of resurfaced regions and the magnitude of resurfacing are likely to change due to the variation of $[I]$ and $\vec{r}$. Third, in the performed simulations, we did not parameterize the asteroid's spin orientation. Depending the direction of the spin axis, the tidal effect may affect the rotational motion differently [39], causing variations in surface acceleration. Thus, the slope distribution may evolve differently due to the rotational conditions during distant encounters. Fourth, we only considered biaxial ellipsoids for our asteroid shapes, which only represent an approximate subset of true shapes of NEAs. Lastly, this study only estimated incipient conditions when surface regolith could be set in motion but did not give the magnitude of the surface mobility. We leave these issues as critical elements to be solved in future work.

§ ACKNOWLEDGEMENT

Y.K. and M.H. acknowledge support from Auburn University's Intramural Grant Program.

[1]
Asphaug, E. and Benz, W. (1994).
Density of comet Shoemaker–Levy 9 deduced by modelling breakup of
the parent 'rubble pile'.
Nature, 370(6485):120–124.

[2]
Ballouz, R.-L., Baresi, N., Crites, S. T., Kawakatsu, Y., and Fujimoto, M.
Surface refreshing of martian moon phobos by orbital
eccentricity-driven grain motion.
Nature Geoscience, 12(4):229.

[3]
Binzel, R. P., Morbidelli, A., Merouane, S., DeMeo, F. E., Birlan, M.,
Vernazza, P., Thomas, C. A., Rivkin, A. S., Bus, S. J., and Tokunaga, A. T.
Earth encounters as the origin of fresh surfaces on near-earth
Nature, 463(7279):331.

[4]
Bottke, W., Nolan, M. C., Greenberg, R., and Kolvoord, R. A. (1994).
Collisional lifetimes and impact statistics of near-earth asteroids.
Hazards due to comets and asteroids, 337.

[5]
Bottke, W. F. (2002).
Asteroids III.
University of Arizona Press.

[6]
Carry, B., Solano, E., Eggl, S., and DeMeo, F. E. (2016).
Spectral properties of near-earth and mars-crossing asteroids using
sloan photometry.
Icarus, 268:340–354.

[7]
Chapman, C. R. (2004).
Space weathering of asteroid surfaces.
Annu. Rev. Earth Planet. Sci., 32:539–567.

[8]
Culling, W. (1960).
Analytical theory of erosion.
The Journal of Geology, 68(3):336–344.

[9]
Delbo, M., Libourel, G., Wilkerson, J., Murdoch, N., Michel, P., Ramesh, K.,
Ganino, C., Verati, C., and Marchi, S. (2014).
Thermal fatigue as the origin of regolith on small asteroids.
Nature, 508(7495):233.

[10]
DeMartini, J. V., Richardson, D. C., Barnouin, O. S., Schmerr, N. C., Plescia,
J. B., Scheirich, P., and Pravec, P. (2019).
Using a discrete element method to investigate seismic response and
spin change of 99942 apophis during its 2029 tidal encounter with earth.
Icarus, 328:93–103.

[11]
DeMeo, F. E., Binzel, R. P., and Lockhart, M. (2014).
Mars encounters cause fresh surfaces on some near-earth asteroids.
Icarus, 227:112–122.

[12]
Dobrovolskis, A. R. (1990).
Tidal disruption of solid bodies.
Icarus, 88(1):24–38.

[13]
Dobrovolskis, A. R. (1996).
Inertia of any polyhedron.
Icarus, 124(2):698–704.

[14]
Fujiwara, A., Kawaguchi, J., Yeomans, D., Abe, M., Mukai, T., Okada, T., Saito,
J., Yano, H., Yoshikawa, M., Scheeres, D., et al. (2006).
The rubble-pile asteroid itokawa as observed by hayabusa.
Science, 312(5778):1330–1334.

[15]
Graves, K. J., Minton, D. A., Hirabayashi, M., DeMeo, F. E., and Carry, B.
Resurfacing asteroids from yorp spin-up and failure.
Icarus, 304:162–171.

[16]
Hapke, B. (2001).
Space weathering from mercury to the asteroid belt.
Journal of Geophysical Research: Planets, 106(E5):10039–10073.

[17]
Hartzell, C. M. (2019).
Dynamics of 2d electrostatic dust levitation at asteroids.
Icarus, 333:234–242.

[18]
Holsapple, K. (2001).
Equilibrium configurations of solid cohesionless bodies.
Icarus, 154(2):432–448.

[19]
Keane, J. T. and Matsuyama, I. (2014).
Rejuvenating neos: the efficiency of asteroid resurfacing via
planetary flybys.
In AAS/Division for Planetary Sciences Meeting Abstracts# 46,
volume 46.

[20]
Keane, J. T. and Matsuyama, I. (2015).
Rejuvenating asteroids during planetary flybys: Applications to
(99942) apophis and other near-earth asteroids.
In 46th Lunar and Planetary Science Conference, page Abstract
#2996, Houston. Lunar and Planetary Institute.

[21]
Lambe, T. W. and Whitman, R. V. (1969).
Soil mechanics.
Wiley, New York (553 pp.).

[22]
Lantz, C., Binzel, R., and DeMeo, F. (2018).
Space weathering trends on carbonaceous asteroids: A possible
explanation for bennu’s blue slope?
Icarus, 302:10–17.

[23]
Marchi, S., Magrin, S., Nesvornỳ, D., Paolicchi, P., and Lazzarin, M.
A spectral slope versus perihelion distance correlation for
planet-crossing asteroids.
Monthly Notices of the Royal Astronomical Society: Letters,

[24]
McNeill, A., Hora, J. L., Gustafsson, A., Trilling, D. E., and Mommert, M.
Constraining the shape distribution of near-earth objects from
partial light curves.
The Astronomical Journal, 157(4):164.

[25]
Molaro, J. L., Byrne, S., and Langer, S. A. (2015).
Grain-scale thermoelastic stresses and spatiotemporal temperature
gradients on airless bodies, implications for rock breakdown.
Journal of Geophysical Research: Planets, 120(2):255–277.

[26]
Molaro, J. L., Byrne, S., and Le, J.-L. (2017).
Thermally induced stresses in boulders on airless body surfaces, and
implications for rock breakdown.
Icarus, 294:247–261.

[27]
Molaro, J. L., Hergenrother, C. W., Chesley, S., Hanna, R. D., Haberle, C. W.,
Ballouz, R.-L., Schwartz, S. R., Bottke, W., Walsh, K. J., Campins, H., and
Lauretta, D. (2019).
Thermal fatigue as a driving mechanism for activity on asteroid
Earth and Space Science Open Archive.

[28]
Morota, T., Sugita, S., Cho, Y., Kanamaru, M., Tatsumi, E., Sakatani, N.,
Honda, R., Hirata, N., Kikuchi, H., Yamada, M., et al. (2020).
Sample collection from asteroid (162173) ryugu by hayabusa2:
Implications for surface evolution.
Science, 368(6491):654–659.

[29]
Nesvornỳ, D., Bottke, W. F., Vokrouhlickỳ, D., Chapman, C. R., and
Rafkin, S. (2010).
Do planetary encounters reset surfaces of near earth asteroids?
Icarus, 209(2):510–519.

[30]
Nesvornỳ, D., Jedicke, R., Whiteley, R. J., and Ivezić, Ž.
Evidence for asteroid space weathering from the sloan digital sky
Icarus, 173(1):132–152.

[31]
Noguchi, T., Kimura, M., Hashimoto, T., Konno, M., Nakamura, T., Zolensky,
M. E., Okazaki, R., Tanaka, M., Tsuchiyama, A., Nakato, A., et al. (2014).
Space weathered rims found on the surfaces of the itokawa dust
Meteoritics & Planetary Science, 49(2):188–214.

[32]
Noguchi, T., Nakamura, T., Kimura, M., Zolensky, M., Tanaka, M., Hashimoto, T.,
Konno, M., Nakato, A., Ogami, T., Fujimura, A., et al. (2011).
Incipient space weathering observed on the surface of itokawa dust
Science, 333(6046):1121–1125.

[33]
Pieters, C. M., Taylor, L. A., Noble, S. K., Keller, L. P., Hapke, B., Morris,
R. V., Allen, C. C., McKAY, D. S., and Wentworth, S. (2000).
Space weathering on airless bodies: Resolving a mystery with lunar
Meteoritics & Planetary Science, 35(5):1101–1107.

[34]
Pravec, P., Harris, A. W., and Warner, B. (2006).
Nea rotations and binaries.
Proceedings of the International Astronomical Union,

[35]
Richardson, D. C., Bottke Jr, W. F., and Love, S. G. (1998).
Tidal distortion and disruption of earth-crossing asteroids.
Icarus, 134(1):47–76.

[36]
Richardson Jr, J. E., Melosh, H. J., Greenberg, R. J., and O'Brien, D. P.
The global effects of impact-induced seismic activity on fractured
asteroid surface morphology.
Icarus, 179(2):325–349.

[37]
Roering, J. J., Kirchner, J. W., and Dietrich, W. E. (1999).
Evidence for nonlinear, diffusive sediment transport on hillslopes
and implications for landscape morphology.
Water Resources Research, 35(3):853–870.

[38]
Sasaki, S., Nakamura, K., Hamabe, Y., Kurahashi, E., and Hiroi, T. (2001).
Production of iron nanoparticles by laser irradiation in a simulation
of lunar-like space weathering.
Nature, 410(6828):555.

[39]
Scheeres, D. (2001).
Changes in rotational angular momentum due to gravitational
interactions between two finite bodies.
Celestial Mechanics and Dynamical Astronomy, 81(1-2):39–44.

[40]
Scheeres, D., Benner, L., Ostro, S., Rossi, A., Marzari, F., and Washabaugh, P.
Abrupt alteration of asteroid 2004 mn4's spin state during its 2029
earth flyby.
Icarus, 178(1):281–283.

[41]
Scotti, J. and Melosh, H. (1993).
Estimate of the size of comet shoemaker-levy 9 from a tidal breakup
Nature, 365(6448):733.

[42]
Sekanina, Z., Chodas, P. W., and Yeomans, D. K. (1994).
Tidal disruption and the appearance of periodic comet shoemaker-levy
Astronomy and Astrophysics, 289:607–636.

[43]
Sharma, I., Jenkins, J. T., and Burns, J. A. (2006).
Tidal encounters of ellipsoidal granular asteroids with planets.
Icarus, 183(2):312–330.

[44]
Sridhar, S. and Tremaine, S. (1992).
Tidal disruption of viscous bodies.
Icarus, 95(1):86–99.

[45]
Strazzulla, G., Dotto, E., Binzel, R., Brunetto, R., Barucci, M., Blanco, A.,
and Orofino, V. (2005).
Spectral alteration of the meteorite epinal (h5) induced by heavy ion
irradiation: A simulation of space weathering effects on near-earth
Icarus, 174(1):31–35.

[46]
Sugita, S., Honda, R., Morota, T., Kameda, S., Sawada, H., Tatsumi, E., Yamada,
M., Honda, C., Yokota, Y., Kouyama, T., et al. (2019).
The geomorphology, color, and thermal properties of ryugu:
Implications for parent-body processes.
Science, 364(6437):eaaw0422.

[47]
Tancredi, G., Maciel, A., Heredia, L., Richeri, P., and Nesmachnow, S. (2012).
Granular physics in low-gravity environments using discrete element
Monthly Notices of the Royal Astronomical Society,

[48]
Vernazza, P., Binzel, R., Thomas, C., DeMeo, F., Bus, S., Rivkin, A., and
Tokunaga, A. (2008).
Compositional differences between meteorites and near-earth
Nature, 454(7206):858.

[49]
Willner, K., Shi, X., and Oberst, J. (2014).
Phobos' shape and topography models.
Planetary and Space Science, 102:51–59.

[50]
Yamada, T. M. and Katsuragi, H. (2014).
Scaling of convective velocity in a vertically vibrated granular bed.
Planetary and Space Science, 100:79–86.

[51]
Yasui, M., Matsumoto, E., and Arakawa, M. (2015).
Experimental study on impact-induced seismic wave propagation through
granular materials.
Icarus, 260:320–331.

[52]
Yu, Y., Richardson, D. C., Michel, P., Schwartz, S. R., and Ballouz, R.-L.
Numerical predictions of surface effects during the 2029 close
approach of asteroid 99942 apophis.
Icarus, 242:82–96.