# Superconductivity in doped Weyl semimetal Mo0.9Ir0.1Te2 with broken
inversion symmetry

Manasi Mandal  Department of Physics, Indian Institute of Science Education
and Research Bhopal, Bhopal, 462066, India  Chandan Patra  Department of
Physics, Indian Institute of Science Education and Research Bhopal, Bhopal,
462066, India  Anshu Kataria  Department of Physics, Indian Institute of
Science Education and Research Bhopal, Bhopal, 462066, India  Suvodeep Paul
Department of Physics, Indian Institute of Science Education and Research
Bhopal, Bhopal, 462066, India  Surajit Saha  Department of Physics, Indian
Institute of Science Education and Research Bhopal, Bhopal, 462066, India  R.
P. Singh<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Science Education and Research Bhopal, Bhopal, 462066, India

###### Abstract

This work presents the emergence of superconductivity in Ir - doped Weyl
semimetal Td \- MoTe2 with broken inversion symmetry. Chiral anomaly induced
planar Hall effect and anisotropic magneto-resistance confirm the topological
semimetallic nature of Mo1-xIrxTe2. Observation of weak anisotropic,
moderately coupled type-II superconductivity in Td -Mo1-xIrxTe2 makes it a
promising candidate for topological superconductor.

##  I INTRODUCTION

Topological superconductor (TSC) is an exotic state of matter where the
nontrivial topology of bulk promotes the emergence of Majorana bound states
within the bulk superconducting gap TI2 ; TI3 . The realization of topological
superconductors (TSCs) is a focus of research interest in condensed-matter
physics due to their potential applications in fault-tolerant topological
quantum computation. Theoretical studies have proposed that topological
superconductivity can be induced at the interface of a topological material
and a conventional superconductor TI4 ; TI5 . However, the possible complexity
at the heterostructure fabrication extends the idea to search for
superconductors with topological surface states. Recently, superconductivity
in Weyl semimetal offers a new platform towards understanding topological
superconductivity since the bulk nodal points have definite chirality and the
Fermi surfaces are intrinsically topological TM1 . An interesting theoretical
proposal suggests that the doped Weyl semimetals with four Weyl points are
natural candidates to realize higher-order topological superconductors that
exhibit a fully gapped bulk state, and the surface hosts robust gapless chiral
hinge states. A finite-range attractive interaction in these materials favours
p+ip pairing symmetry, and such a chiral pairing state is identified as a
higher-order topological superconductor HOTI1 .

MoTe2 is a promising member of this family, having an abundant collection of
structural variation with unusual attractive electronic properties, such as
charge density wave (CDW) CDW , an edge supercurrent edge_supercurrent  and
surface superconductivity surface_superconductivity . Noncentrosymmetric
orthorhombic Td-MoTe2 was predicted to be type-II Weyl semimetal (WSM), which
was further confirmed by experimental calculations like angle-resolved
photoemission spectroscopy (ARPES) and scanning tunneling microscopy (STM)
MoTe2_topology1 ; manasi_raman ; MoTe2_topology7 ; MoTe2_STM . Td-phase MoTe2
shows superconductivity at transition temperature TC = 0.10 K and
superconducting transition temperature is sensitive to the chemical pressure
MTS ; MTSe1 ; MTSe2 . It is worthy to mention that S-doped MoTe2 shows two
band superconductivity with s$\pm$ pairing, and the quasi-particle
interference patterns with band-structure calculations reveal the existence of
Fermi arcs SMoTe2 ; SMoTe2_2 . Scanning tunnelling spectroscopy (STS) studies
on the surface of doped MoTe2 reveal a much larger superconducting gap that
could be interpreted as topological nontrivial superconductivity based on the
pairing of Fermi arc surface states SMoTe2 . Re doping in Mo-site in MoTe2
enhances TC up to 4.1 K with new unconventional electronic properties like
CDW, and strain-induced pseudo magnetic field and anomalous phonon behaviour
manasi1 ; manasi2 ; jpcm ; manasi_raman . Also, the broken inversion symmetry
in Td phase MoTe2 can generate anti-symmetric spin-orbit coupling and lift the
spin degeneracy. This allows an admixture of spin-singlet and spin-triplet
superconducting ground states, which may lead to host exotic time-reversal
symmetry breaking ground state ReZr ; ReHf ; Ti ; ReTi2 ; La7Ir ; La7Rh ;
La7Ni . Therefore, it is essential to study superconductivity in MoTe2 through
chemical doping to identify the exotic topological states and determine the
role of non-trivial topological surface states in the superconducting pairing
mechanism.

In this paper, we report the emergence of superconductivity in Td-Mo1-xIrxTe2
single crystal. It is characterised by Powder X-ray diffraction (XRD), energy-
dispersive X-ray (EDAX) and temperature dependent Raman measurements to
confirm the phase purity and structural phase transition. Magnetization, AC
transport, specific heat measurements were performed to probe
superconductivity in Ir-doped MoTe2 system. It is observed that Ir-
substitution in Mo-site enhances the superconducting transition temperature up
to 2.65(3) K. Planner Hall effect and anisotropic magneto-resistance confirm
the Weyl semimetallic nature after the doping. Realisation of
superconductivity with intact topological surface state suggests that
Mo1-xIrxTe2 as a new candidate of topological superconductor.

Figure 1: (a) Room-temperature crystal structure of the Mo0.9Ir0.1Te2 material
(space group P21/m). (b) Rietveld refinement of powder XRD pattern for
Mo1-xIrxTe2 (x = 0.1) recorded at room temperature. (c) Room temperature
single crystal XRD shows the orientation of crystal growth along (00n)
direction whereas inset shows the crystal image.

##  II EXPERIMENTAL DETAILS

Pure phase polycrystalline samples of composition Mo1-xIrxTe2 (x = 0 and 0.1)
were prepared from stoichiometric mixtures of Mo (99.9 % pure), Ir (99.99 %
pure) and Te (99.99% pure) by the standard solid-state reaction process
manasi1 . We tried to dope Ir with a higher percentage, but 10% was the
highest possible concentration successfully achieved. Single crystals of
Mo1-xIrxTe2 (x = 0.1) were grown by horizontal flux method with excess Te as
flux flux_method . Stoichiometric mixtures of Mo, Ir, and Te powders were
sealed together in an evacuated quartz tube with charge and flux ratio 1:18.
Crystallization was carried out at 1100$\degree$C for seven days in a box
furnace, followed by the furnace switched off to avoid the hexagonal phase
formation. Single crystals were trapped in the Te flux, which was taken out
mechanically.

Powder x-ray diffraction (XRD) was carried out at room temperature on a
PANalytical diffractometer equipped with Cu-Kα radiation ($\lambda$ = 1.54056
Å)for the characterization of crystal structure and phase purity. Rietveld
refinement was performed using Full Prof Suite Software. Sample compositions
were checked by a scanning electron microscope (SEM) equipped with an energy-
dispersive X-ray (EDAX) spectrometer, which confirms Mo, Ir, and Te in the
samples. The Raman spectra were measured in a Horiba JY LabRam HR Evolution
Raman spectrometer and temperature-dependent Raman measurements were performed
using a Linkam stage. To study the superconducting state, we measured DC and
AC susceptibility by Superconducting Quantum Interference Device (SQUID MPMS,
Quantum Design) at different fields. Specific heat measurements were performed
by the two tau time-relaxation method in zero fields. Using the conventional
four-probe technique with horizontal Rotator in the Physical Property
Measurement System (PPMS, Quantum Design, Inc.), AC transport measurements
were performed.

##  III RESULTS AND DISCUSSION

###  a Sample characterization

Polycrystalline sample was crushed in fine powder to perform XRD. The Rietveld
refinement of XRD patterns confirms that the polycrystalline samples
crystallized into single phase centrosymmetric monoclinic CdI2 type structure
having space group P21/m as shown in Fig. 1 (b). XRD pattern of Mo0.9Ir0.1Te2
single-crystal indicates the crystal orientation in (00n) direction (Fig. 1
(c)) and inset shows the crystal image. Raman spectra further confirms the
phase purity (see Appendix B). EDAX analysis shows the presence of Ir in the
single crystal though the quantity of Ir is (4$\%$) lesser than the nominal
one (see Appendix A).

Figure 2:  (a) and (b) Temperature dependence of the dc magnetic moment on
polycrystal and single crystal sample at H = 1 mT.

###  b Superconductivity

Superconductivity in Mo0.9Ir0.1Te2 polycrystal and the single crystal was
further confirmed by magnetization measurements. Due to the small size of the
crystals and less magnetic moments, a number of crystals were stacked together
for magnetization measurements. The magnetic moment of the single crystal was
checked in zero fields cooled (ZFC) and field cooled (FC) with 1 mT applied
field. Fig. 2 exhibits a clear signature of type-II superconducting state with
a diamagnetic signal. The polycrystalline sample shows superconducting
transition temperature (TC) at 2.32(2) K. In contrast, single-crystal yields
TC at 2.65(3) K.

Figure 3: (a) The low-field magnetization curves for Mo0.9Ir0.1Te2 single
crystal at different temperature. (b) Temperature dependence of the lower
critical field HC1 was fitted using Ginzburg-Landau relation.

A sharp drop of electrical resistivity, $\rho$(T) further confirms the
superconductivity at temperature, TC = 2.70(2) K (see Appendix C). An anomaly
with hysteresis in the $\rho$(T), which is associated with the first-order
structural phase transition (TS) from 1T′ phase to Td phase is observed for
Mo0.9Ir0.1Te2 single crystal at 150 K, which is consistent with Raman study
(see Appendix B).

To find the lower and upper critical fields, magnetic measurements were done
on single crystal. We estimated the lower critical field HC1(0) from the low
field magnetization curves M(H) in the range of 0-5 mT took at different
temperatures, as shown in (Fig. 3). The lower critical field, HC1(T), is
determined from the first deviation from the linearity of the initial slope as
the field is increased. HC1(0) was estimated to be 3.45(4) mT by fitting the
data in accordance with Ginzburg-Landau (GL) equation
$H_{C1}(T)=H_{C1}(0)\left(1-\left(\frac{T}{T_{C}}\right)^{2}\right)$.

The temperature dependence of the magnetic moment under various magnetic
fields was also measured to calculate the second-order transition field
HC2(0), as shown in Fig. 4 (a). Since, in the case of a type- II material, the
magnetic field can penetrate the sample and reduce the gap function, TC shifts
to lower temperature with the increment of the applied field. We estimated the
value of the upper critical field by using Ginzburg-Landau (GL) formula
$H_{C2}(T)=H_{C2}(0)\frac{(1-t^{2})}{(1+t^{2})}$, where t = T/TC. The
estimated value of HC2(0) is around 1.01(3) T as shown in Fig. 4(b).

Figure 4: (a) The $\frac{M}{H}$ curves for various applied magnetic fields.
(c) HC2(0) was calculated by Ginzburg-Landau fits.

The Ginzburg Landau coherence length $\xi_{GL}(0)$ is related to H${}_{C2}(0)$
by the relation $H_{C2}(0)$ = $\frac{\Phi_{0}}{2\pi\xi_{GL}^{2}}$ where
$\Phi_{0}$ is the magnetic flux quantum having value 2.07 $\times$10-15 T m2
Tinkham . The another superconducting characteristics parameter Penetration
depth $\lambda_{GL}$(0) associated with HC1(0) by the relation (1):

$H_{C1}(0)=\frac{\Phi_{0}}{4\pi\lambda_{GL}^{2}(0)}\left(ln\frac{\lambda_{GL}(0)}{\xi_{GL}(0)}+0.497\right)$
(1)

The estimated values of $\xi_{GL}(0)$ and $\lambda_{GL}$(0) are 180.6 Å and
2819 Å respectively. GL parameter $k_{GL}$ =
$\frac{\lambda_{GL}(0)}{\xi_{GL}(0)}$ = 21.2 ($\gg\frac{1}{\sqrt{2}}$),
indicating that Mo0.9Ir0.1Te2 is a type-II superconductor. Thermodynamic
critical field HC is evaluated around 0.034(1) T using the relation
$H_{C1}(0)H_{C2}(0)$ = $H_{C}^{2}lnk_{GL}$. Breaking of Cooper pair depends on
the combine consequence of Pauli paramagnetic and orbital diamagnetic effect.
The Pauli limiting field within the BCS theory is given by H${}_{C2}^{p}$(0) =
1.83 TC, which gives H${}_{C2}^{p}$(0) = 4.85(1) T. In accordance to
Werthamer-Helfand-Hohenberg (WHH) EH ; NRW  expression the upper critical
field in orbital limit, H${}_{C2}^{orbital}$(0) can be calculated by the
relation $H_{C2}^{orbital}(0)=-\alpha
T_{C}\left.\frac{dH_{C2}(T)}{dT}\right|_{T=T_{C}}$. The obtained vale of
H${}_{C2}^{orbital}$(0) is 0.67(2) T using $\alpha$ = 0.69 for dirty limit
superconductors and the initial slope
$\left(\frac{-dH_{C2}(T)}{dT}\right)|_{T=T_{C}}$= 0.37 calculated from
$H_{C2}$ -T phase diagram. The relative strengths of the orbital and Pauli
limiting values of HC2 is given by the Maki parameter
$\alpha_{M}=\sqrt{2}H_{C2}^{orb}(0)/H_{C2}^{p}(0)$ maki . We obtain
$\alpha_{M}$ = 0.20 indicating that orbital effect is negligible in breaking
of Cooper pairs.

Figure 5:  Temperature dependence of electrical resistivity $\rho$(T) at
different field angle, $\theta$ = (a) 0$\degree$ and (b) 90$\degree$. (c)
shows field dependence of electrical resistivity $\rho$(T) at different
$\theta$ with fixed temperature, T = 1.8 K. (d) Angle dependent $\rho$(T) at T
= 1.8 K whereas solid red line is 3D GL-fit (Eq. (2)).

We have further characterized the superconducting property by angle-dependent
resistivity measurements. Fig. 5(a) and (b) show the temperature dependence of
resistivity for the magnetic fields H $\parallel$ c and H $\perp$ c, which
clearly indicates anisotropy in the upper critical field, HC2. Fig. 5(c) shows
the field dependence of the resistivity at different angles of H with c-axis
at lowest possible temperature T = 1.8 K ($\simeq$ 0.6 TC). The angular
dependence of HC2($\theta$) determined from the onset value of the curve. In
order to clarify the origin of the characteristic angular dependence of
HC2($\theta$), the experimental data were fitted with GL-model. According to
the three dimensional (3D) anisotropic mass GL model HC2($\theta$) can be
described by the following relation CeIr3 ; NbSe2 ; NbSe2 ; NbSe2_2 ; NbSe2_3
:

$\left(\frac{H_{C2}(T,\theta)cos\theta}{H_{C2}^{\perp}(T)}\right)^{2}+\left(\frac{H_{C2}(T,\theta)sin\theta}{H_{C2}^{\parallel}(T)}\right)^{2}=1$
(2)

The experimental data ware well fitted (red solid in Fig. 5(d)) with the
relation Eq. (2) with low anisotropy factor $\gamma$ =
$\frac{H_{C2}(90\degree)}{H_{C2}(0\degree)}$ = 2.1. The nature of
HC2($\theta$) is in sharp contrast to pressure-induced two-dimensional
superconductivity though the $\gamma$ value is comparable MoTe2 . Similar type
of feature in HC2($\theta$) has observed in 3D-anisotropic superconductor
2H-NbSe2 NbSe2 ; NbSe2 ; NbSe2_2 ; NbSe2_3 , WTe2 WTe2_hc2 , and gated MoS2
MoS2 . Furthermore, the in-plane, $\xi_{GL}^{\parallel}$ and out-of-plane
coherence lengths $\xi_{GL}^{\perp}$ at T = 1.8 K can be extracted from the
HC2 data, giving 247 and 166 Å. The coherence length along the c direction
($\xi_{GL}^{\perp}$) is much larger than the lattice parameter along the same
direction (13.7 Å).
Heat capacity measurements were performed at zero fields, where the
superconducting transition is manifested by a jump in the heat capacity data
at 2.4 K (Fig. 6). The low temperature normal-state specific heat data was
analyzed in the temperature range 1.9 K $\leq$ T $\leq$ 8 K by using the
relation $\frac{C}{T}$ = $\gamma_{n}+\beta_{3}T^{2}+\beta_{5}T^{4}$ where,
$\gamma_{n}$ is Sommerfeld coefficient, $\beta_{3}$ is Debye constant and
$\beta_{5}$ is the anharmonic contribution to the specific heat. The
extrapolation of normal state behavior below TC, to the T $\to$ 0 limits,
allows the determination of normal state coefficients. The linear fits to the
C/T vs T2 gave $\gamma_{n}$ = 5.05 $\pm$ 0.06 mJ mol-1 K-2, $\beta_{3}$ = 1.22
$\pm$ 0.01 mJ mol-1 K-4 and $\beta_{5}$ = (5 $\pm$ 0.1)*10-7 $\mu$J mol-1 K-6
for Mo0.9Ir0.1Te2. We calculated the Debye temperatures $\theta_{D}$ of the
compounds using the formula $\theta_{D}$ =
$\left(\frac{12\pi^{4}RN}{5\beta_{3}}\right)^{\frac{1}{3}}$ where N (= 3) is
the number of atoms per formula unit, R is the molar gas constant (= 8.314 J
mol-1 K-1), which was based on the simple Debye model for the phonon
contribution to the specific heat. The estimated $\theta_{D}$ value was 168.5
K. For non-interacting fermions, the Sommerfeld coefficient, $\gamma_{n}$ is
proportional to the density of states $D_{C}(E_{F})$ at the Fermi level which
was calculated to be 2.14 $\frac{states}{eVf.u}$ from the relation
$\gamma_{n}$ = $\left(\frac{\pi^{2}k_{B}^{2}}{3}\right)D_{C}(E_{f})$, where kB
$\simeq$ 1.38 $\times$ 10-23 J K-1. The strength of the attractive interaction
between electron and phonon can be estimated by the electron-phonon coupling
constant, $\lambda_{e-ph}$ by McMillan equation Mc1
$\lambda_{e-ph}=\frac{1.04+\mu^{*}ln(\theta_{D}/1.45T_{C})}{(1-0.62\mu^{*})ln(\theta_{D}/1.45T_{C})-1.04}$,
where $\mu^{*}$ (= 0.13 for many superconductors) is the repulsive screened
coulomb parameter. The calculated values of $\lambda_{e-ph}$ = 0.61 suggesting
that Mo0.9Ir0.1Te2 is a moderately coupled superconductor such as Mo1-xRexTe2
manasi1 , Mg10Ir19B16 MgIrB , and LaPtGe LPG . The estimated parameters for
Mo0.9Ir0.1Te2 single crystal are comparable to Mo1-xRexTe2. We have summarized
all the superconducting and normal state parameters in Table 1. By
substituting Ir for Mo in MoTe2, facilitate the enhancement of the electron-
phonon coupling and DOS at the Fermi energy similar to Re doped MoTe2 compound
manasi1 . So, we may expect unusual physical property in this system like
Mo1-xRexTe2 compound. More investigation on a single crystal upto very low
temperature is needed to explore the effect of Ir doping and the actual nature
of the superconducting gap in this interesting system. To explore the
topological nature in the Ir-doped MoTe2 compound, we have done AC transport
measurements.

Figure 6:  (a) C/T vs T2 at 0 Oe field. (b) shows the data fitted with $\frac{C}{T}$ = $\gamma_{n}+\beta_{3}T^{2}+\beta_{5}T^{4}$. Table 1: Normal and superconducting parameters of Mo0.9Ir0.1Te2 Parameters | Unit | Mo0.9Ir0.1Te2 | Mo0.7Re0.3Te2 manasi1
---|---|---|---
T${}_{C}^{res}$ | K | 2.70(2) | 4.1
H${}_{C1}(0)$ | mT | 3.45(4) | 5
HC2(0) | T | 1.01(3) | 7.2
H${}_{C}(0)$ | T | 0.034(1) |
$\xi_{GL}$ | Å | 180.6 | 72.2
$\lambda_{GL}$ | Å | 3819 |
$k_{GL}$ |  | 21.2 |
$\gamma_{n}$ | mJ/mol K2 | 5.05(6) | 9
$\beta_{3}$ | mJ/mol K4 | 1.22(1) |
$\theta_{D}$ | K | 165.5 |
$\lambda_{e-ph}$ |  | 0.61 | 0.64
DC(Ef) | states/eV f.u. | 2.14 | 3.83

###  c Topological Semimetallic Property

In transport experiments, WSMs usually distinctively attribute the chiral
anomaly-induced negative longitudinal magnetoresistance (LMR) PHE1 ; PHE2 ;
PHE3 ; PHE4 . This feature indicates the nonconservation of chiral charge
around the Weyl nodes due to the coplanar electric and magnetic fields
($\boldsymbol{E}.\boldsymbol{B}$ = $\mathrm{E}\mathrm{B}\mathrm{cos\phi}$)
PHE1 ; PHE2 . The experimental observation of negative LMR is very critical in
type-II Weyl semimetal due to several extrinsic effects, such as the current
jetting effect. Negative LMR has not been reported in the MoTe2 system so far
within our knowledge LMR ; LMR_2 . Planar Hall effect (PHE) is also a
practical tool to probe the chiral charge pumping effect that directly
correlates the physics of topological nature and the transport properties.
Experimentally PHE is observed in Td-MoTe2 MoTe2_topology2 ; MoTe2_topology4
as well as other Dirac semimetals such as WTe2 WTe2 ; WTe2_2 , ZrTe5 ZrTe5
and Cd3As2 CdAs ; CdAs_2 .

Figure 7:  The extracted intrinsic planar Hall resistivity $\rho_{xy}^{PHE}$
(a) at T = 5 K under different fields (b) at H = 9 T under different
temperatures. The variation of anisotropic resistivity,
$\Delta\rho_{xy}^{PHE}$ with field and temperature is shown in (c) and (d).

We have measured planar Hall resistivity on Mo0.9Ir0.1Te2 crystal where the
applied field is rotated with angle $\phi$ to the current direction. A
schematic of planar Hall measurements is shown in the inset of Fig. 7. The
obtained Hall resistivity has the total contribution of regular Hall and
planar one. Normal Hall resistivity is antisymmetric under the opposite
applied field, whereas the chiral-anomaly induced Hall resistivity is
symmetric under the antisymmetric applied field PHE3 ; PHE4 . So planar hall
resistivity can be extracted by average data of opposite applied field:

$\rho_{xy}^{PHE}=[\rho_{xy}(B)+\rho_{xy}(-B)]/2$ (3)

Fig. 7 (a) shows the angular dependence of planar Hall resistivity at T = 5 K
under different applied magnetic fields. Fig. 7 (b) shows the same at H = 9 T
under different temperatures. In WSMs, chiral anomaly induced planar Hall
effect (PHE) can be expressed as PHE3 ; PHE4

$\rho_{xy}^{PHE}=\Delta\rho_{xy}^{chiral}sin\phi cos\phi$ (4)

Where $\Delta\rho_{xy}^{chiral}$ is the chiral anomaly induced anisotropic
resistivity that can be expressed as
$\Delta\rho_{xy}^{chiral}=\rho_{\perp}-\rho_{\parallel}$. Where $\rho_{\perp}$
and $\rho_{\parallel}$ are the resistivity corresponding to the magnetic field
perpendicular to and along the direction of the current flow (I).

Figure 8:  (a) Angle dependent $\rho_{xx}$ at T = 5 K under different fields
whereas (b) shows the angle dependent ($\rho_{xx}-\rho_{\perp}$) at H = 9 T
under different temperatures. The solid curve is a fit to equation. The
variation of anisotropic resistivity, $\Delta\rho_{xx}$ with field and
temperature extracted from fitting is shown in (c) and (d).

Fig. 7 shows that angular dependence of $\rho_{xy}^{PHE}$ has a period of
$\pi$ with maximums appear at 45$\degree$ and 135$\degree$ similar to MoTe2
and several other topological semimetals. We have fitted the data at T = 5 K
with Eq. (4) and plotted the field dependence of $\Delta\rho_{xy}^{chiral}$ in
Fig. 7(c), which follows quadratic nature as expected for chiral anomaly
induced PHE. Temperature dependence of $\Delta\rho_{xy}^{chiral}$ at H = 9 T
is also shown in Fig. 7(d) which decreases with increases temperature similar
to the nature of DSMs and WSMs.

Anisotropic magnetoresistivity (AMR), $\Delta\rho_{xx}$ can also be observed
by the angular dependence ($\phi$) of longitudinal resistivity $\rho_{xx}$. In
WSMs, chiral anomaly related AMR can be expressed as PHE3 ; PHE4 :

$\rho_{xx}=\rho_{\perp}-\Delta\rho_{xx}^{chiral}cos^{2}\phi$ (5)

Fig. 8(a) shows the the angular dependence ($\phi$) of $\rho_{xx}$ at T = 5 K
under different applied fields. The data were fitted by the Eq. (5), and
$\Delta\rho_{xx}^{chiral}$ was extracted for a different field. We have shown
the field variation of $\Delta\rho_{xx}^{chiral}$ in Fig. 8(c) that fitted
well with quadratic equation ($\Delta\rho_{xx}^{chiral}$ $\propto$ H2). The
nature of the curve is consistent to Fig. 7(c). The slight variation in value
of $\Delta\rho_{xx}^{chiral}$ and $\Delta\rho_{xy}^{chiral}$ be due to the
sample dimensions used for the purpose of the measurements.

Fig. 8(b) shows the temperature dependence of $\rho_{xx}-\rho_{\perp}$ taken
at various temperature at constant applied magnetic filed H = 9 T. The data
were fitted with Eq. (5) which yields temperature dependence of
$\Delta\rho_{xx}^{chiral}$ as shown in Fig. 8(d). Similar nature was observed
in Fig. 8(d) validating our measurements and analysis.

It is observed that nontrivial Berry-curvature-induced PHE is prounced in Ir-
doped Td-MoTe2 and still visible upto 150 K. The quadratic field dependence of
$\Delta\rho_{xy}^{chiral}$ indicates the weak chiral anomaly induced coupling
strength between electric and chiral charge whereas parent MoTe2 shows
intermediate coupling strength with high field at T $\leq$ 50 K
MoTe2_topology4 . The change in coupling strength may be due to the small
change in lattice constant or carrier concentration due to Ir-doping similar
to the theoretically prediction theory_MoTe2 ; theory_MoTe2_2 .

##  IV Conclusion

In summary, we have successfully prepared polycrystal and single crystal of
Mo1-xIrxTe2 (x = 0 and 0.1) compound crystallized into a centrosymmetric
monoclinic structure having space group P21/m. Magnetization measurements
reveal that Mo0.9Ir0.1Te2 is a type-II superconductor with a transition
temperature of 2.65(3) K, which is 26 times enhanced with respect to parent
MoTe2. Resistivity measurements indicate that structural transition is
decreased with Ir doping similar to the Re-doping effect in MoTe2 compound.
Raman study confirms the phase purity and 1T′ to Td structural phase
transition at TS $\simeq$ 150 K. Angle-dependent higher critical field,
H${}_{C2}(\theta)$ follows 3D anisotropic GL-model with factor 2.1. The
enhanced density of states, electron-phonon coupling, and other estimated
values of superconducting and normal state parameters make Mo0.9Ir0.1Te2
compound interesting to study further to explore unusual attractive electronic
state like Re-doped MoTe2 system. Above all, exotic chiral anomaly induced PHE
is observed up to 150 K in Mo1-xIrxTe2 (x = 0.1). We can conclude this system
as a possible topological superconductor with an exotic angle-dependent higher
critical field. Further microscopic investigations like muon spectroscopy and
scanning tunneling microscopy required to probe superconducting gap symmetry
and ground state.

##  V Acknowledgments

R. P. S. acknowledges the Science and Engineering Research Board, Government
of India for the Core Research Grant CRG/2019/001028. S. S. acknowledges the
Science and Engineering Research Board, Government of India for the Research
Grant ECR/2016/001376. Authors thank Mr. Devesh Negi [CSIR fellow, File No.:
09/1020(0139)/2018-EMR-I] and Ms. Bommareddy Poojitha for their help with the
Raman measurements.

##  Appendix A Characterization

Figure 9:  (a) Room temperature powder XRD pattern for x = 0 and 0.1 (b) EDX
pattern for x = 0.1.

EDAX analysis is summarized in Table 2.

Table 2: Average elemental concentration (average of 10 points) obtained from the EDS measurements for Mo1-xIrxTe2 (x = 0.1) samples. | Nominal composition | From EDS
---|---|---
Polycrystal | Mo0.9Ir0.1Te2 | Mo0.9Ir0.08Te2
Single crystal | Mo0.9Ir0.1Te2 | Mo0.9Ir0.04Te2

##  Appendix B Raman Study

Micro-Raman measurements are very effective in confirming the phase purity of
investigated samples down to micron length scales. The micro-Raman
measurements reported here were performed at a laser spot size of 0.67 $\mu$m,
thereby allowing the investigation of any impurity phase with a spatial
resolution of $<1$ $\mu$m. Fig. 10(a) compares the room-temperature Raman
spectra of MoTe2 and Mo0.9Ir0.1Te2. The Raman modes observed in pure MoTe2 are
labelled as P1-P9, which correspond to out-of-plane vibrations with symmetries
Ag or Bg raman1 ; raman2 ; raman3 ; raman4 . We also observe that the room-
temperature Raman spectrum obtained for Mo0.9Ir0.1Te2 shows similar modes as
in for the pure phase, thus, confirming the absence of any impure phases in
the investigated flakes. In addition to the confirmation of the phase purity,
our temperature-dependent Raman measurements provide evidence of the
topological phase transition in Mo0.9Ir0.1Te2 at $\simeq$150 K from a
crystallographic point of view. Fig. 10(b) shows the Raman spectra obtained
for Mo0.9Ir0.1Te2 at various temperatures. We clearly observe a splitting of
the P6 mode into two bands (P6A and P6B) below $\simeq$150 K. A similar
observation of mode-splitting was previously reported for pure manasi_raman ;
raman1 ; raman2 ; raman3 ; raman4  and Re-doped MoTe2 manasi_raman ,
attributed to the topological phase transition from the high-temperature 1T′
to low-temperature Td phase. Importantly, the topological phase transition is
associated with a structural transformation characterized by a modification of
the layer stacking order. The high-temperature phase (1T′) has a monoclinic
structure possessing a center of inversion symmetry, while the low-temperature
phase (Td) is orthorhombic without the inversion symmetry. Therefore, the
topological phase transition is associated with a breaking of inversion
symmetry, which is also a fundamental requirement for the appearance of the
Weyl semimetallic phase, that can be probed using Raman spectroscopy. The
splitting of the P6 mode is a result of the renormalization of phonons induced
by the breaking of inversion symmetry at $\simeq$ 150 K. The observation also
matches well with the AC transport measurements (Appendix C), which show a
hysteretic change in $\rho$(T) around $\simeq$150 K (Fig. 11).

Figure 10:  (a) A comparison of the room temperature Raman spectra of MoTe2
and Mo0.9Ir0.1Te2. (b) Raman spectra of Mo0.9Ir0.1Te2 taken at different
temperature revealing the splitting of the P6 mode to P6A and P6B at low
temperature.

##  Appendix C AC Transport Measurements

Figure 11:  Normal state resistivity data $\rho$(T) over the temperature range
1.9 K $\leq$ T $\leq$ 300 K. A hysteresis in $\rho$(T) observed under the
cooling and warming modes indicates polar structural transition temperature as
shown by arrow symbol. Inset of (b) highlights the superconducting phase
transition for x = 0.1.

Electrical resistivity, $\rho$(T) were measured on Mo1-x IrxTe2 (x = 0 and
0.1) single crystal samples over the temperature range 1.9 K $\leq$ T $\leq$
300 K in zero magnetic field in the cooling and warming modes. An anomaly with
hysteresis in the $\rho$(T), which is associated with the first-order
structural phase transition (TS) from 1T′ phase to Td phase is observed for
both the samples (x = 0 and 0.1). The polar structural transition temperature,
TS, is found at around 240 K in parent MoTe2 whereas the same has decreased to
150 K for Mo0.9Ir0.1Te2, which is consistent with Raman study. The structural
transition temperature is indicated by a arrow symbol in Fig. 11. Similar type
behaviour of resisitivity was observed in Re-doped MoTe2 system manasi1 . A
sharp drop of $\rho$(T) is observed in Mo0.9Ir0.1Te2 sample at temperaute, TC
= 2.70(2) K. This type of nature indicates the superconducting phase
transition.

Figure 12:  Regular Hall measurements performed on Mo1-x IrxTe2 (x = 0 and
0.1) at 10 K and 2 K upto 9 T.

In order to find the nature of carrier we have performed regular Hall effect
on the ab plane with the applied filed along the c-direction (Hall
configuration is shown in Fig. 12). The linear field dependence of hall
resistivity shows negative slope indicating the dominating electron carrier in
Mo0.9Ir0.1Te2 system as shown in Fig. 12. Ir doping may enhance the electron
concentration in the system similar to Re doping in Mo-site manasi1 .

Figure 13:  MR data taken at 5 K and 2 K on Mo0.9Ir0.1Te2 crystal at different
angles, $\theta$ between the direction of c-axis and applied magnetic field.

We have also measured magneto-resistance (MR) on Ir-doped crystal at 5 K and 2
K at different angles, $\theta$ between the c-axis and magnetic field
direction, H. The the magnetic field dependence of the MR is defined as
$\frac{\rho_{xx}(H)-\rho_{xx}(0)}{\rho_{xx}(0)}$. As shown in Fig. 13(a) and
(b), MR decreases with increasing angle $\theta$ from 0$\degree$ to
90$\degree$ at 5 K and 2 K, suggesting an anisotropy of transport behavior. MR
data at T = 2 K indicates the superconducting nature at low field region. MR
data suggests that magneto-resistance has been suppressed with Ir-doping in
MoTe2 system which may be due to the emergence of superconductivity
MoTe2_topology4 .

## References

  * (1) X. L. Qi and S. C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (2) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * (3) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * (4) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * (5) J. Wang, Natl. Sci. Rev. 6, 2 (2019).
  * (6) A. Jahin, A. Tiwari, and Y. Wang, arXiv:2103.05010 (2021).
  * (7) L. Wang, Y. Wu, Y. Yu, A. Chen, H. Li, W. Ren, S. Lu, S. Ding, H. Yang, Q.-K. Xue, F.-S. Li, and G. Wang, ACS Nano 14, 8299 (2020).
  * (8) W. Wang, S. Kim, M. Liu, F. Cevallos, R. Cava, and N. Ong, Science 368, 534 (2020).
  * (9) Y. Naidyuk, O. Kvitnitskaya, D. Bashlakov, S. Aswartham, I. Morozov, I. Chernyavskii, G. Fuchs, D. stefan Ludwig, R. Hühne, K. Nielsch, B. Buchner, and D. Efremov, 2D Mater. 5 (2018).
  * (10) K. Deng, G. Wan, P. Deng, K. Zhang, S. Ding, E. Wang, M. Yan, H. Huang, H. Zhang, Z. Xu, J. Denlinger, A. Fedorov, H. Yang, W. Duan, H. Yao, Y. Wu, S. Fan, H. Zhang, X. Chen, and S. Zhou, Nat. Phys. 12, 1105 (2016).
  * (11) S. Paul, S. Karak, M. Mandal, A. Ram, S. Marik, R. P. Singh, and S. Saha, Phys. Rev. B 102, 054103 (2020).
  * (12) J. Jiang, Z. K. Liu, Y. Sun, H. F. Yang, C. R. Rajamathi, Y. P. Qi, L. X. Yang, C. Chen, H. Peng, C. C. Hwang, S. Z. Sun, S. K. Mo, I. Vobornik, J. Fujii, S. S. P. Parkin, C. Felser, B. H. Yan, and Y. L. Chen, Nat. Commun. 8, 13973 (2017).
  * (13) P. Deng, Z. Xu, K. Deng, K. Zhang, Y. Wu, H. Zhang, S. Zhou, and X. Chen, Phys. Rev. B 95, 245110 (2017).
  * (14) F. Chen, X. Luo, R. Xiao, W. Lu, B. Zhang, H. Yang, J. Li, Q. Pei, D.-F. Shao, R. Zhang, L. Ling, C. Xi, W. Song, and Y. Sun, Appl. Phys. Lett. 108, 162601 (2016).
  * (15) H. Takahashi, T. Akiba, K. Imura, T. Shiino, K. Deguchi, N. K. Sato, H. Sakai, M. S. Bahramy, and S. Ishiwata, Phys. Rev. B 95, 100501 (2017).
  * (16) P. Li, J. Cui, J. Zhou, D. Guo, Z. Zhao, J. Yi, J. Fan, Z. Ji, X. Jing, F. Qu, C. Yang, L. Lu, J. Lin, Z. Liu, and G. Liu, Adv. Mater. 31, e1904641 (2019).
  * (17) Y. Lia, Q. Gua, C. Chend, J. Zhangd, Q. Liud, X. Hua, J. Liug, Y. Liua, L. Lingg, M. Tiang, Y. Wangh, N. Samarthc, S. Lid, T. Zhangd, J. Feng, and J. Wanga, PNAS 115, 9503 (2018).
  * (18) P. Hosur, X. Dai, Z. Fang, and X.-L. Qi, Phys. Rev. B 90, 045130 (2014).
  * (19) M. Mandal, S. Marik, K. P. Sajilesh, Arushi, D. Singh, J. Chakraborty, N. Ganguli, and R. P. Singh, Phys. Rev. Mater. 2, 094201 (2018).
  * (20) S. Kamboj, P. S. Rana, A. Sirohi, A. Vasdev, M. Mandal, S. Marik, R. P. Singh, T. Das, and G. Sheet, Phys. Rev. B 100, 115105 (2019).
  * (21) A. Vasdev, S. Kamboj1, A. Sirohi, M. Mandal, S. Marik, R. P. Singh, and G. Sheet, J. Phys.: Condens. Matter 33, 255401 (2021).
  * (22) R. P. Singh, A. D. Hillier, B. Mazidian, J. Quintanilla, J. F. Annett, D. M. Paul, G. Balakrishnan, and M. R. Lees, Phys. Rev. Lett. 112, 107002 (2014).
  * (23) D. Singh, J. A. T. Barker, A. Thamizhavel, D. M. Paul, A. D. Hillier, and R. P. Singh, Phys. Rev. B 96, 180501 (2017).
  * (24) D. Singh, K. P. Sajilesh, J. A. T. Barker, D. M. Paul, A. D. Hillier, and R. P. Singh, Phys. Rev. B 97, 100505 (2018).
  * (25) D. Singh, Sajilesh K. P., J. A. T. Barker, D. McK. Paul, A. D. Hillier, and R. P. Singh, Phys. Rev. B 97, 100505(R) (2018).
  * (26) P. K. Biswas, H. Luetkens, T. Neupert, T. Sturzer, C. Baines, G. Pascua, A. P. Schnyder, M. H. Fischer, J. Goryo, M. R. Lees, H. Maeter, F. Bruckner, H. H. Klauss, M. Nicklas, P. J. Baker, A. D. Hillier, M. Sigrist, A. Amato, and D. Johrendt, Phys. Rev. B 87, 180503(R) (2013).
  * (27) D. Singh, M. S. Scheurer, A. D. Hillier, D. T. Adroja, and R. P. Singh, Phys. Rev. B 102, 134511 (2020).
  * (28) Arushi, D. Singh, A. D. Hillier, M. S. Scheurer, and R. P. Singh, Phys. Rev. B 103, 174502 (2021).
  * (29) Y. Tanaka, Y. Mizuno, T. Yokoyama, K. Yada, and M. Sato, Phys. Rev. Lett. 105, 097002 (2010).
  * (30) X. L. Qi and S. C. Zhang, Rev. Mod. Phys. 83, 1057 (2011)
  * (31) J.-Q. Yan, B. C. Sales, M. A. Susner, and M. A. McGuire, Phys. Rev. Mater. 1, 023402 (2017).
  * (32) M. Tinkham, ed., Introduction to Superconductivity, 2nd ed. (McGraw - Hill, New York, 1996).
  * (33) E. Helfand and N. R. Werthamer, Phys. Rev. 147, 288 (1966).
  * (34) N. R. Werthamer, E. Helfand, and P. C. Hohenberg, Phys. Rev. 147, 295 (1966).
  * (35) K. Maki, Phys. Rev. 148, 362 (1966).
  * (36) Y. J. Sato, F. Honda, Y. Shimizu, A. Nakamura, Y. Homma, A. Maurya, D. Li, T. Koizumi, and D. Aoki, Phys. Rev. B 102, 174503 (2020).
  * (37) R. C. Morris, R. V. Coleman, and R. Bhandari, Phys. Rev. B 5, 895 (1972).
  * (38) S. Yasuzuka, S. Uji, S. Sugiura, T. Terashima, Y. Nogami, K. Ichimura, S. Tanda, J. Supercond. Nov. Magn. 33, 953 (2020).
  * (39) H. Matsuoka, M. Nakano, T. Shitaokoshi, T. Ouchi, Y. Wang, Y. Kashiwabara, S. Yoshida, K. Ishizaka, M. Kawasaki, Y. Kohama, T. Nojima, and Y. Iwasa, Phys. rev. res. 2, 012064(R) (2020).
  * (40) Y. J. Hu, Y. Tai Chan, Kwing To Lai, K. O. Ho, X. Guo, H.-P. Sun, K. Y. Yip, Dickon H. L. Ng, H.-Z. Lu, and Swee K. Goh, Phys. Rev. Mater. 3, 034201 (2019).
  * (41) Y. Tai Chan, P. L. Alireza, K. Y. Yip, Q. Niu, K. T. Lai, and Swee K. Goh, Phys. Rev. B 96, 180504(R) (2017).
  * (42) J. M. Lu, O. Zheliuk, I. Leemakers, N. F. Q. Yuan, U. Zeltler, K. T. Law, J. T. Ye, Science, 350, 1353 (2015).
  * (43) W. L. McMillan, Phys. Rev. 167, 331 (1968).
  * (44) T. Klimczuk, F. Ronning, V. Sidorov, R. J. Cava, and J. D. Thompson, Phys. Rev. Lett. 99, 257004 (2007).
  * (45) S. K. P., D. Singh, P. K. Biswas, A. D. Hillier, and R. P. Singh, Phys. Rev. B 98, 214505 (2018).
  * (46) C. L. Zhang, S. Y. Xu, I. Belopolski, Z. Yuan, Z. Lin, B. Tong, G. Bian, N. Alidoust, C. C. Lee, S. M. Huang, T. R. Chang, G. Chang, C. H. Hsu, H. T. Jeng, M. Neupane, D. S. Sanchez, H. Zheng, J. Wang, H. Lin, C. Zhang, H. Z. Lu, S. Q. Shen, T. Neupert, M. Z. Hasan, and S. Jia, Nat. Commun. 7, 10735 (2016).
  * (47) Y. Y. Lv, X. Li, B. B. Zhang, W. Y. Deng, S. H. Yao, Y. B. Chen, J. Zhou, S. T. Zhang, M. H. Lu, L. Zhang, M. Tian, L. Sheng, and Y. F. Chen, Phys. Rev. Lett. 118, 096603 (2017).
  * (48) A. A. Burkov, Phys. Rev. B 96, 041110 (2017).
  * (49) S. Nandy, G. Sharma, A. Taraphder, and S. Tewari, Phys. Rev. Lett. 119, 176804 (2017).
  * (50) Y. Sun, S. C. Wu, M. N. Ali, C. Felser, and B. Yan, Phys. Rev. B 92, 161107 (2015).
  * (51) Z. Wang, D. Gresch, A. A. Soluyanov, W. Xie, S. Kushwaha, X. Dai, M. Troyer, R. J. Cava, and B. A. Bernevig, Phys. Rev. Lett. 117, 056805 (2016).
  * (52) X. Luo, F. Chen, J. Zhang, Q. Pei, G. Lin, W. Lu, Y. Han, C. Xi, W. Song, and Y. Sun, Appl. Phys. Lett. 109, 102601 (2016).
  * (53) F. C. Chen, X. Luo, J. Yan, Y. Sun, H. Y. Lv, W. J. Lu, C. Y. Xi, P. Tong, Z. G. Sheng, X. B. Zhu, W. H. Song, and Y. P. Sun, Phys. Rev. B 98, 041114(R) (2018).
  * (54) A. Soluyanov, D. Gresch, Z. Wang, Q.-S. Wu, M. Troyer, X. Dai, and B. Bernevig, Nature 527, 495 (2015).
  * (55) P. Li, C. Zhang, Y. Wen, L. Cheng, G. Nichols, David G. Cory, G.-X. Miao, and Xi-Xiang Zhang, Phys. Rev. B 100, 205128 (2019).
  * (56) P. Li, C. H. Zhang, J. W. Zhang, Y. Wen, and X. X. Zhang, Phys. Rev. B 98, 121108 (2018).
  * (57) H. Li, H. Wang, H. He, J. Wang, and S. Q. Shen, Phys. Rev. B 97, 201110 (2018).
  * (58) M. Wu, G. Zheng, W. Chu, W. Gao, H. Zhang, J. Lu, Y. Han, J. Yang, H. Du, W. Ning, Y. Zhang, and M. Tian, Phys. Rev. B 98, 161110 (2018).
  * (59) Y. Sun, S. C. Wu, M. N. Ali, C. Felser, and B. Yan, Phys. Rev. B 92, 161107 (2015).
  * (60) Z. Wang, D. Gresch, A. A. Soluyanov, W. Xie, S. Kushwaha, X. Dai, M. Troyer, R. J. Cava, and B. A. Bernevig, Phys. Rev. Lett. 117, 056805 (2016).
  * (61) K. Zhang, C. Bao, Q. Gu, X. Ren, H. Zhang, K. Deng, Y. Wu, Y. Li, J. Feng, and S. Zhou, Nat. Commun. 7, 13552 (2016).
  * (62) Y. Y. Lv, L. Cao, X. Li, B. B. Zhang, K. Wang, B. Pang, L. Ma, D. Lin, S. H. Yao, J. Zhou, and Y. B. Chen, Sci. Rep. 7, 44587 (2017).
  * (63) S. M. Oliver, R. Beams, S. Krylyuk, I. Kalish, A. K. Singh, A. Bruma, F. Tavazza, J. Joshi, I. R. Stone, S. J. Stranick, and A. V. Davydov, 2D Mater. 4, 045008 (2017).
  * (64) X. J. Yan, Y. Y. Lv, L. Li, X. Li, S. H. Yao, Y. B. Chen, X. P. Liu, H. Lu, M. H. Lu, and Y. F. Chen, npj Quantum Mater. 2, 31 (2017).


# Probe the regolith characteristics of asteroids from 9-years infrared
observations of WISE/NEOWISE: A case study of the Main-Belt Object (656)
Beagle

Liang-Liang Yu  Institute of Science and Technology for Deep Space
Exploration, Nanjing University-Suzhou Campus, Suzhou, 215163, China

(Received MM DD, YY; Revised MM DD, YY; Accepted MM DD, YY)

###### Abstract

This work presents data processing, fitting procedure, modelling and analyzing
of 9-years infrared light curves provided by the WISE/NEOWISE telescope, by
which the regolith characteristics of Main-Belt Object (656) Beagle is
studied. We determine Beagle’s effective diameter $D_{\rm
eff}=57.3^{+4.5}_{-2.2}$ km, geometric albedo $p_{\rm
v}=0.05^{+0.004}_{-0.007}$, mean roughness $\theta_{\rm RMS}=44\pm 4^{\circ}$,
mean grain size $b=100^{+350}_{-90}~{}\mu$m, mean specific heat capacity
$c_{\rm p}=173\sim 516\rm~{}JKg^{-1}K^{-1}$, mean thermal conductivity
$\kappa=0.7\sim 1.3\times 10^{-3}\rm~{}Wm^{-1}K^{-1}$ and mean thermal inertia
$\Gamma=14\sim 32\rm~{}Jm^{-2}s^{-0.5}K^{-1}$. The albedo of Beagle is a
little anomalous that the albedos of Beagle’s neighbouring asteroids are more
close to Themis, rather than Beagle itself. The W1-band near-infrared light
curves don’t reveal significant heterogeneous NIR features on the surface of
Beagle, being inconsistent with the expectation of a family parent that has
members with diverse NIR spectral types. These results add new clues of Beagle
probably being an interloper or a sister, rather than the parent of its
neighbouring asteroids including the first main-belt comet (MBC) 133P, hence
may lead to new scenarios about the origin of famous MBC 133P. Besides, we
found that asteroidal shape models from inversion of optical light curves are
imperfect for modeling infrared lightcurves, thus could mislead evaluations of
both the heterogeneity of regolith reflectivity at near infrared and
thermophysical characteristics at thermal infrared.

Light curves — Infrared photometry — Small Solar System bodies — Main belt
asteroids — Asteroid surfaces

revtex4-1.clsrevtex4-2.cls

##  1 Introduction

The main-belt asteroid (656) Beagle (hereafter Beagle for short) came into the
sight of astronomers since it was proposed to be the parent of the famous
main-belt comet (MBC) 133P/Elst-Pizarro (Nesvorný et al., 2008). And the
Beagle family is proposed to be a young family with an age of $<10$ Myr
(Nesvorný et al., 2008) or $\sim 14$ Myr (Carruba, 2019). Then the story about
the origin of 133P is widely accepted as the following scenario: 133P is a
young icy fragment of Beagle, which is a daughter of the icy main-belt object
(24) Themis, on which absorption feature of water ice was observed by the IRTF
spectra (Campins et al., 2010; Rivkin & Emery, 2010). The formation scenario
was somewhat supported by the work of (Fornasier et al., 2016), which reported
Beagle to have a geometric albedo of $p_{\rm v}=0.0782\pm 0.0222$, being
similar to the geometric albedos of its neighbouring asteroids (mean $p_{\rm
v}=0.0941\pm 0.0055$) and potential parent (24) Themis ($p_{\rm v}=0.0641\pm
0.0157$).

Besides, Fornasier et al. (2016) also found that the neighbouring asteroids of
Beagle exhibit diverse spectral types based on the reflectance spectrum at
near infrared (NIR), implying that the composition of these asteroids differs
from each other. If these neighbouring asteroids of Beagle were formed from
the same impact event, then their diverse compositions may originate from the
parent, or from the impactor, or produced by the high temperature and pressure
during impact process. According to Michel et al. (2015), if a parent body is
heterogeneous in composition, then the resulting family is expected to show a
variety of spectral properties within its members. On the other hand, although
the original parent may not have spectral variability on its surface, the
impact that produced heterogenous members will make the survived parent to
show surface heterogeneity. For example, Vesta is observed to show
heterogeneous surface features due to impact exposing heterogeneous subsurface
and upper crust materials (Rousseau, De Sanctis & Raponi et al, 2021).
Therefore, if Beagle is the survived parent of these neighbouring asteroids
with diverse compositions, then Beagle is expected to contain diverse
compositions on its surface and show heterogeneous NIR features across its
surface.

However, Mainzer et al. (2016) reported Beagle to have a very low geometric
albedo $p_{\rm v}=0.045\pm 0.005$, being even lower than that of (24) Themis.
On the other hand, the recent work of Yu et al. (2020) found that 133P is more
likely to have an old age $>100$ Myr, being inconsistent with the young age
$<10$ Myr of the Beagle family. These new observations and theoretical
predictions thus raise the question of whether Beagle having genetic
connections with 133P and (24) Themis.

So in this work, we attempt to use more data and new method to investigate
whether asteroid Beagle has a geometric albedo as low as $\sim 0.05$ and
significant heterogeneous reflectivity at near infrared. This objective can be
achieved by analyzing the multi-year infrared light curves of Beagle from
WISE/NEOWISE with the well-tested thermophysical model — RSTPM (Yu & Ip,
2021). As the WISE/NEOWISE observations of main-belt objects at band-W1 are
dominated by sunlight reflection at near infrared. In addition, heterogeneity
of regolith thermophysical characteristics can also be evaluated from the
other three bands of W2, W3, and W4.

##  2 The Radiometric model

###  2.1 Observations and data processing

The Wide-field Infrared Survey Explorer (WISE) mission has mapped entire sky
in four bands around 3.4 (W1), 4.6 (W2), 12 (W3), and 22 (W4) $\mu$m with
resolutions from $6.1^{\prime\prime}$ to $12^{\prime\prime}$ (Wright et al.,
2010). All four bands were imaged simultaneously, and the exposure times were
7.7 s in W1 and W2, and 8.8 s in W3 and W4 (Wright et al., 2010). The four-
bands survey started from 2010 January 7, and ended on 2010 August 6 after the
outer cryogen tank was exhausted, making the W4 channel be no longer able to
be used to obtain survey data. The W3 channel continued operation until 2010
September 29 when the inner cryogen reserve was exhausted, while the W1 and W2
channel kept working until the telescope was set into hibernation on 2011
February 1 (Mainzer et al., 2011; Cutri et al., 2012). The two-band survey was
then resumed on 2013 December 13 (known as NEOWISE) (Mainzer et al., 2014;
Cutri et al., 2015);, and is still in service, which has obtained nearly
10-year observations.

We found 9-years observations of Beagle from the WISE archive (see the website
of the NASA/IPAC Infrared Science Archive http://irsa.ipac.caltech.edu/). The
datasets are summerized in Table 1.

Table 1: The collected WISE/NEOWISE observations of (656) Beagle from 2010 to 2019. UTC | Bands | Obs | $r_{\rm helio}$ | $\Delta_{\rm obs}$ | $\alpha$
---|---|---|---|---|---
|  | Num | (AU) | (AU) | (∘)
Jan. 31-01, 2010 | W1-4 | 14 | 2.82 | 2.66 | 20.42
Jul. 26-28, 2010 | W1-4 | 16 | 3.03 | 2.78 | -19.51
Feb. 15-16, 2014 | W1-2 | 11 | 3.01 | 2.81 | -19.17
Dem. 05-06, 2014 | W1-2 | 14 | 2.73 | 2.52 | 21.11
May. 19-20, 2015 | W1-2 | 16 | 2.75 | 2.40 | -21.28
Mar. 09-11, 2016 | W1-2 | 11 | 3.06 | 2.88 | 18.95
Aug. 20-21, 2016 | W1-2 | 11 | 3.27 | 2.88 | -17.52
May. 14-15, 2017 | W1-2 | 14 | 3.52 | 3.39 | 16.70
Oct. 19-20, 2017 | W1-2 | 12 | 3.57 | 3.12 | -15.32
Jul. 17-18, 2018 | W1-2 | 11 | 3.49 | 3.33 | 16.95
Dem. 12-13, 2018 | W1-2 | 13 | 3.35 | 2.84 | -15.66
Oct. 03-04, 2019 | W1-2 | 12 | 2.99 | 2.81 | 19.58
$r_{\rm helio}$: heliocentric distance;
$\Delta_{\rm obs}$: observation distance;
$\alpha$: solar phase angle.

####  2.1.1 Flux colour corrections

All the four-band data need colour corrections, especially for the W1 and
W2-band data, the infrared fluxes of which contain not only self-thermal
emission but also sunlight diffused from the object’s surface, thus the color
corrections should be done separately for the thermal component and the
reflection component, as the color correction factors are different. Besides,
the color correction factor for the thermal component is temperature
dependent, thus will have different value if the helio-centric distance was
different for each observation epoch. So we implement the color-correction
procedure as follows:

First, convert the data-base magnitude to flux without any color correction,
giving the total integrated flux $F(\lambda\pm\Delta\lambda)_{\rm tot,obs}$
for each band. The derived band-integrated fluxes together with the
observation geometry are listed in Table 2 and 3, and each flux should at
least have an associated uncertainty of $\pm$10 percent according to Wright et
al. (2010);

Table 2: WISE observations of (656) Beagle. The derived fluxes are band-integrated fluxes of each band without color correction. Each flux should have an associated uncertainty of ±10 percent according to Wright et al. (2010). MJD | W1 | W2 | W3 | W4 | $r_{\rm helio}$ | $\Delta_{\rm obs}$ | $\alpha$
---|---|---|---|---|---|---|---
| (mJy) | (mJy) | (Jy) | (Jy) | (AU) | (AU) |  (∘)
55227.504 | 1.49 | 6.64 | 1.26 | 3.18 | 2.82 | 2.67 | 20.42
55227.637 | 1.44 | 6.71 | 1.16 | 2.99 | 2.82 | 2.67 | 20.42
55227.637 | 1.51 | 6.57 | 1.22 | 3.22 | 2.82 | 2.67 | 20.42
55227.769 | 1.42 | 5.44 | 0.94 | 2.58 | 2.82 | 2.67 | 20.42
55227.901 | 0.96 | - | 0.74 | 2.16 | 2.82 | 2.66 | 20.42
55228.034 | 0.77 | 2.95 | 0.54 | 1.64 | 2.82 | 2.66 | 20.42
55228.100 | 1.53 | 6.93 | 1.16 | 3.37 | 2.82 | 2.66 | 20.42
55228.166 | 0.73 | 2.84 | 0.57 | 1.69 | 2.82 | 2.66 | 20.42
55228.232 | 1.58 | 7.36 | 1.23 | 3.34 | 2.82 | 2.66 | 20.42
55228.298 | 0.87 | 4.14 | 0.75 | 2.25 | 2.82 | 2.66 | 20.42
55228.364 | 1.29 | 5.95 | 1.06 | 2.93 | 2.82 | 2.66 | 20.42
55228.497 | 1.17 | 4.68 | 0.87 | 2.44 | 2.82 | 2.66 | 20.42
55228.629 | 0.95 | 3.90 | 0.73 | 1.95 | 2.82 | 2.65 | 20.42
55228.762 | - | 3.23 | 0.55 | 1.59 | 2.82 | 2.65 | 20.42
55403.908 | 1.17 | 3.00 | 0.68 | 2.18 | 3.03 | 2.77 | -19.50
55404.040 | 0.80 | - | 0.50 | 1.70 | 3.03 | 2.77 | -19.50
55404.172 | 0.58 | 1.73 | 0.40 | 1.34 | 3.03 | 2.77 | -19.51
55404.173 | 0.58 | 1.68 | 0.41 | 1.34 | 3.03 | 2.77 | -19.51
55404.305 | 0.62 | 1.59 | 0.36 | 1.19 | 3.03 | 2.77 | -19.51
55404.371 | 1.21 | 3.68 | 0.88 | 2.54 | 3.03 | 2.78 | -19.51
55404.437 | 0.76 | 2.05 | 0.41 | 1.30 | 3.03 | 2.78 | -19.51
55404.503 | 1.01 | 2.70 | 0.70 | 2.32 | 3.03 | 2.78 | -19.51
55404.569 | 0.90 | 2.67 | 0.58 | 1.76 | 3.03 | 2.78 | -19.51
55404.636 | 0.98 | 2.33 | 0.63 | 2.01 | 3.03 | 2.78 | -19.51
55404.768 | 0.81 | 2.02 | 0.48 | 1.63 | 3.03 | 2.78 | -19.51
55404.900 | 0.38 | - | 0.35 | 1.20 | 3.03 | 2.78 | -19.51
55405.032 | 0.50 | 1.43 | 0.33 | 1.02 | 3.03 | 2.79 | -19.52
55405.032 | 0.53 | 1.26 | 0.33 | 1.04 | 3.03 | 2.79 | -19.52
55405.165 | 0.68 | 2.09 | 0.45 | 1.41 | 3.03 | 2.79 | -19.52
55405.165 | 0.65 | 2.31 | 0.50 | 1.53 | 3.03 | 2.79 | -19.52
Table 3: NEOWISE observations of (656) Beagle. The derived fluxes are band-integrated fluxes of each band without color correction. Each flux should have an associated uncertainty of ±10 percent according to Wright et al. (2010). MJD | W1 | W2 | $r_{\rm helio}$ | $\Delta_{\rm obs}$ | $\alpha$ |  | MJD | W1 | W2 | $r_{\rm helio}$ | $\Delta_{\rm obs}$ | $\alpha$
---|---|---|---|---|---|---|---|---|---|---|---|---
| (mJy) | (mJy) | (AU) | (AU) |  (∘) |  |  | (mJy) | (mJy) | (AU) | (AU) |  (∘)
56703.138 | 0.99 | 2.23 | 3.01 | 2.80 | -19.16 |  | 57887.678 | 0.32 | 0.61 | 3.52 | 3.40 | 16.70
56703.270 | 0.75 | 1.69 | 3.01 | 2.80 | -19.16 |  | 57887.809 | 0.25 | - | 3.52 | 3.40 | 16.70
56703.401 | 0.54 | 1.21 | 3.01 | 2.80 | -19.16 |  | 57887.940 | 0.36 | 0.65 | 3.52 | 3.39 | 16.70
56703.533 | 0.66 | 1.68 | 3.01 | 2.81 | -19.17 |  | 57887.940 | 0.41 | 0.52 | 3.52 | 3.39 | 16.70
56703.599 | 1.18 | 3.22 | 3.01 | 2.81 | -19.17 |  | 57888.071 | 0.54 | 0.80 | 3.52 | 3.39 | 16.70
56703.665 | 1.01 | 2.80 | 3.01 | 2.81 | -19.17 |  | 57888.136 | 0.47 | 0.67 | 3.52 | 3.39 | 16.70
56703.731 | 1.14 | 2.62 | 3.01 | 2.81 | -19.17 |  | 57888.202 | 0.64 | 0.86 | 3.52 | 3.39 | 16.70
56703.797 | 1.05 | 3.13 | 3.01 | 2.81 | -19.17 |  | 57888.267 | 0.34 | 0.47 | 3.52 | 3.39 | 16.70
56703.994 | 0.60 | 1.36 | 3.01 | 2.81 | -19.18 |  | 57888.333 | 0.68 | 1.23 | 3.52 | 3.39 | 16.71
56704.126 | 0.70 | 1.54 | 3.01 | 2.81 | -19.18 |  | 57888.333 | 0.71 | 1.38 | 3.52 | 3.39 | 16.71
56704.258 | 0.72 | 2.64 | 3.01 | 2.82 | -19.18 |  | 57888.398 | 0.28 | 0.51 | 3.52 | 3.39 | 16.71
56996.277 | 1.71 | 8.14 | 2.73 | 2.53 | 21.12 |  | 57888.529 | 0.36 | 0.45 | 3.52 | 3.39 | 16.71
56996.409 | 1.51 | 7.25 | 2.73 | 2.53 | 21.12 |  | 57888.660 | 0.56 | 0.82 | 3.52 | 3.38 | 16.71
56996.540 | 1.17 | 5.58 | 2.73 | 2.53 | 21.11 |  | 57888.791 | 0.60 | 0.98 | 3.52 | 3.38 | 16.71
56996.672 | 0.84 | 3.77 | 2.73 | 2.52 | 21.11 |  | 58045.491 | 0.40 | 0.53 | 3.57 | 3.12 | -15.28
56996.737 | 1.91 | 9.87 | 2.73 | 2.52 | 21.11 |  | 58045.622 | 0.38 | 0.46 | 3.57 | 3.12 | -15.29
56996.803 | 0.66 | 3.16 | 2.73 | 2.52 | 21.11 |  | 58045.622 | 0.31 | 0.53 | 3.57 | 3.12 | -15.29
56996.869 | 1.84 | 9.27 | 2.73 | 2.52 | 21.11 |  | 58045.753 | 0.44 | 0.59 | 3.57 | 3.12 | -15.30
56996.934 | 0.94 | 4.62 | 2.73 | 2.52 | 21.11 |  | 58045.884 | 0.56 | 0.97 | 3.57 | 3.12 | -15.31
56997.000 | 1.64 | 8.15 | 2.73 | 2.52 | 21.11 |  | 58046.015 | 0.57 | 1.06 | 3.57 | 3.12 | -15.32
56997.066 | 1.53 | 7.41 | 2.73 | 2.52 | 21.11 |  | 58046.080 | 0.50 | 0.69 | 3.57 | 3.13 | -15.32
56997.132 | 1.37 | 6.12 | 2.73 | 2.52 | 21.11 |  | 58046.146 | 0.93 | 1.08 | 3.57 | 3.13 | -15.33
56997.263 | 0.93 | 4.10 | 2.73 | 2.52 | 21.11 |  | 58046.211 | 0.39 | 0.42 | 3.57 | 3.13 | -15.34
56997.394 | 0.76 | 3.28 | 2.73 | 2.51 | 21.11 |  | 58046.342 | 0.39 | 0.54 | 3.57 | 3.13 | -15.35
56997.526 | 0.96 | 4.15 | 2.73 | 2.51 | 21.10 |  | 58046.473 | 0.52 | 0.57 | 3.57 | 3.13 | -15.36
57161.421 | 1.82 | 7.32 | 2.75 | 2.39 | -21.25 |  | 58046.604 | 0.67 | 0.94 | 3.57 | 3.13 | -15.37
57161.552 | 1.41 | 5.31 | 2.75 | 2.40 | -21.26 |  | 58316.418 | 0.33 | 0.60 | 3.49 | 3.34 | 16.95
57161.552 | 1.65 | 5.52 | 2.75 | 2.40 | -21.26 |  | 58316.549 | 0.44 | 0.84 | 3.49 | 3.34 | 16.95
57161.683 | 1.20 | 3.57 | 2.75 | 2.40 | -21.26 |  | 58316.680 | 0.57 | 0.95 | 3.49 | 3.34 | 16.95
57161.814 | 0.87 | 2.76 | 2.75 | 2.40 | -21.27 |  | 58316.811 | 0.67 | 1.26 | 3.49 | 3.33 | 16.95
57161.880 | 1.87 | 7.80 | 2.75 | 2.40 | -21.27 |  | 58316.942 | 0.62 | 1.20 | 3.49 | 3.33 | 16.95
57161.946 | 0.89 | 3.42 | 2.75 | 2.40 | -21.28 |  | 58317.007 | 0.36 | 0.53 | 3.49 | 3.33 | 16.95
57161.946 | 0.87 | 3.51 | 2.75 | 2.40 | -21.28 |  | 58317.073 | 0.63 | 1.03 | 3.49 | 3.33 | 16.95
57162.011 | 1.98 | 8.11 | 2.75 | 2.40 | -21.28 |  | 58317.138 | 0.37 | 0.71 | 3.49 | 3.33 | 16.95
57162.077 | 1.22 | 5.61 | 2.75 | 2.40 | -21.28 |  | 58317.269 | 0.59 | 0.85 | 3.49 | 3.33 | 16.95
57162.143 | 1.63 | 5.65 | 2.75 | 2.40 | -21.29 |  | 58317.400 | 0.72 | 1.18 | 3.49 | 3.32 | 16.95
57162.274 | 1.19 | 4.16 | 2.75 | 2.41 | -21.29 |  | 58317.531 | 0.57 | 0.97 | 3.49 | 3.32 | 16.95
57162.405 | 0.95 | 3.04 | 2.75 | 2.41 | -21.30 |  | 58464.862 | 0.47 | 0.90 | 3.35 | 2.83 | -15.60
57162.536 | 1.02 | 3.38 | 2.75 | 2.41 | -21.31 |  | 58464.992 | 0.62 | 1.29 | 3.35 | 2.84 | -15.62
57162.536 | 0.89 | 3.22 | 2.75 | 2.41 | -21.31 |  | 58464.992 | 0.71 | 1.27 | 3.35 | 2.84 | -15.62
57162.667 | 1.04 | 5.02 | 2.75 | 2.41 | -21.31 |  | 58465.123 | 0.86 | 1.63 | 3.35 | 2.84 | -15.63
57456.914 | 1.39 | 3.76 | 3.06 | 2.89 | 18.95 |  | 58465.254 | 1.03 | 1.87 | 3.35 | 2.84 | -15.65
57457.045 | 1.02 | 3.40 | 3.06 | 2.89 | 18.95 |  | 58465.320 | 0.55 | 0.77 | 3.35 | 2.84 | -15.66
57457.176 | 1.19 | 3.26 | 3.06 | 2.89 | 18.95 |  | 58465.385 | 1.13 | 2.01 | 3.35 | 2.84 | -15.66
57457.307 | 0.98 | 2.63 | 3.06 | 2.89 | 18.95 |  | 58465.450 | 0.40 | - | 3.35 | 2.84 | -15.67
57457.373 | 1.24 | 2.93 | 3.06 | 2.89 | 18.95 |  | 58465.451 | 0.42 | 0.78 | 3.35 | 2.84 | -15.67
57457.438 | 0.93 | 2.10 | 3.06 | 2.88 | 18.95 |  | 58465.516 | 0.97 | 1.61 | 3.35 | 2.84 | -15.68
57457.635 | - | 3.98 | 3.06 | 2.88 | 18.94 |  | 58465.581 | 0.67 | 1.19 | 3.35 | 2.84 | -15.69
57457.766 | 1.18 | 3.48 | 3.06 | 2.88 | 18.94 |  | 58465.712 | 0.84 | 1.45 | 3.35 | 2.84 | -15.70
57457.897 | 1.08 | 2.94 | 3.06 | 2.88 | 18.94 |  | 58465.843 | 0.98 | 1.83 | 3.35 | 2.85 | -15.72
57458.028 | 0.68 | 2.09 | 3.06 | 2.88 | 18.94 |  | 58759.742 | 0.50 | 1.67 | 2.99 | 2.82 | 19.58
57458.159 | 0.68 | 1.62 | 3.06 | 2.88 | 18.94 |  | 58759.742 | 0.55 | 1.67 | 2.99 | 2.82 | 19.58
57620.359 | 0.44 | - | 3.27 | 2.88 | -17.49 |  | 58759.872 | 0.59 | 2.11 | 2.99 | 2.82 | 19.58
57620.490 | 0.48 | 0.99 | 3.27 | 2.88 | -17.50 |  | 58760.134 | 0.92 | 3.20 | 2.99 | 2.81 | 19.58
57620.621 | 0.81 | 1.47 | 3.27 | 2.88 | -17.51 |  | 58760.200 | 0.70 | 2.54 | 2.99 | 2.81 | 19.58
57620.752 | 0.92 | 2.06 | 3.27 | 2.88 | -17.52 |  | 58760.265 | 1.06 | 4.34 | 2.99 | 2.81 | 19.58
57620.817 | 0.56 | 1.06 | 3.27 | 2.88 | -17.52 |  | 58760.330 | 0.76 | 2.04 | 2.99 | 2.81 | 19.58
57620.883 | - | 2.43 | 3.27 | 2.88 | -17.52 |  | 58760.396 | 1.21 | 4.10 | 2.99 | 2.81 | 19.58
57620.948 | 0.58 | 0.92 | 3.27 | 2.88 | -17.53 |  | 58760.461 | 0.66 | 1.62 | 2.99 | 2.81 | 19.58
57621.014 | 0.84 | 2.00 | 3.27 | 2.89 | -17.53 |  | 58760.592 | 0.68 | 2.35 | 2.99 | 2.81 | 19.58
57621.079 | 0.49 | 1.03 | 3.27 | 2.89 | -17.53 |  | 58760.723 | 0.90 | 3.02 | 2.99 | 2.80 | 19.58
57621.210 | 0.71 | 1.56 | 3.27 | 2.89 | -17.54 |  | 58760.854 | 1.23 | 4.36 | 2.98 | 2.80 | 19.58
57621.341 | 0.76 | 1.93 | 3.27 | 2.89 | -17.55 |  |  |  |  |  |  |

Second, $F(\lambda\pm\Delta\lambda)_{\rm tot,obs}$ contains both integrated
thermal emission and sunlight reflection, and can be modelled as

$F(\lambda\pm\Delta\lambda)_{\rm tot,model}=F(\lambda)_{\rm rl}*f_{\rm
c,rl}+F(\lambda)_{\rm th}*f_{\rm c,th},$ (1)

where $F(\lambda)_{\rm rl}$ and $F(\lambda)_{\rm th}$ represent the
monochromatic sunlight reflection and thermal emission that can be calculated
from reflection and thermal models; $f_{\rm c,rl}$ and $f_{\rm c,th}$ stand
for the color correction factor of the reflection component and the thermal
component respectively; here $f_{\rm c,rl}$ can be chosen to be the color
correction factor of G2V star, thus fixing

$f_{\rm c,rl}=1.0049,1.0193,1.0024,1.0012$

for the four band — W1, W2, W3 and W4 respectively (Wright et al., 2010);

Third, since $f_{\rm c,th}$ is temperature dependent, so its value for each
band will be obtained from an interpolation procedure according to the
effective temperature

$T_{\rm eff}=\left[\frac{(1-A_{\rm B})L_{\odot}}{\varepsilon\sigma\pi
d_{\odot}^{2}}\right]^{1/4}\approx\frac{300\rm~{}K}{\sqrt{d_{\odot}}},$ (2)

of the asteroid at the time of observation, on the basis of the Table 1 of
Wright et al. (2010). In Equation (2), $d_{\odot}$ means the heliocentric
distance in AU; $L_{\odot}$ is the Solar constant, about
$1361.5\rm~{}Wm^{-2}$, $A_{\rm B}$ is the surface bond albedo; $\varepsilon$
is the averaged thermal emissivity over the entire emission spectrum of the
surface and $\sigma$ is the Stefan Boltzmann constant;

Finally, calculate the theoretical monochromatic sunlight reflection
$F(\lambda)_{\rm rl}$ and thermal emission $F(\lambda)_{\rm th}$ via our RSTPM
code, and then obtain the theoretical total flux
$F(\lambda\pm\Delta\lambda)_{\rm tot,model}$ via Equation (1), so as to do
comparison with the observed total flux $F(\lambda\pm\Delta\lambda)_{\rm
tot,obs}$ at the first step.

####  2.1.2 Anomalous data removal

Due to some noise pollution, e.g. cosmic rays and zodiacal light, some of the
observed fluxes may deviate too much from the actual emission of the target.
Such anomalous data need to be removed from the parameter inversion procedure,
so as to make the results about the physical characteristics of the target get
to be more close to reality.

The basic way to remove anomalous data is on the basis of ’signal to noise
ratio ($snr$)’. For the WISE/NEOWISE observations of asteroids, generally we
remove the data with $snr<3$.

For asteroids with shape model, there is also an additional way to remove
anomalous data. As the multi-year data of WISE and NEOWISE make it possible to
generate infrared light curves, so by doing comparisons between the
observational and theoretical infrared light curves, the data sets deviate too
much from the theoretical curves can be treated as anomalous data and remove
them can make the reversion procedure get better results.

###  2.2 Themophysical model

As has been mentioned above, this work uses the thermophysical model for
realistic surface layers on airless small bodies — RSTPM (Yu & Ip, 2021). The
model considers not only real shape and rough surface, but also real orbital
cycle, rotational cycle, and even temperature dependent thermal parameters in
the simulation process, as well as contribution of sunlight reflection in the
infrared radiometric procedure.

In comparison to previous models, the advantages of the RSTPM are: (1) A
different mathematical technique is used to solve the influence of surface
roughness on the energy balance equation of surface boundary; (2) For the aim
to remove the degeneracy of thermal inertia and roughness by interpreting
multi-epoch thermal light-curves, variation of thermal parameters due to
temperature variation caused by orbital cycle and rotation cycle is taken into
consideration; (3) A combination model of simultaneously computing thermal
emission and sunlight-reflection under the same surface topography is proposed
to fit infrared data in case of the data containing significant sunlight
reflection.

So RSTPM has advantages for the small bodies which has regolith on the surface
and has large orbital eccentricity or an obliquity close to 90 degrees. And
for the target asteroid Beagle, we assume it have fine regolith on its
surface, as it has a much larger size than Eros, Ryugu and Bennu.

##  3 Fitting procedure and Results

###  3.1 The Input and Free Parameters

In order to interpret the multi-year observations, RSTPM needs several input
parameters, including observation geometry, shape model, spin orientation,
rotation phase $ph$, scattering weight-factor $w_{\rm f}$, geometric albedo
$p_{\rm v}$, effective diameter $D_{\rm eff}$, mean grain radius $b$, and
roughness fraction $f_{\rm r}$, which can be related to the root-mean-square
slope as $\theta_{\rm RMS}=\sqrt{f_{\rm r}}\times 50^{\circ}$ (Yu & Ip, 2021).

The observation geometry at the time of each observation can be easily
obtained according to the orbit of Beagle and WISE. The spin period and spin
orientation together with shape model of Beagle are also known by light curve
inversion method from the DAMIT database. The utilized shape model of Beagle
is shown in Figure 1, where $ph$ represents rotational phase defined as
$ph=1-\varphi/(2\pi)$, in which $\varphi$ is the local longitude of the
observer in the body-fixed coordinate system.

Figure 1: The light-curve inversion shape model of (656) Beagle from the DAMIT
database. $ph$ represents rotational phase defined as $ph=1-\varphi/(2\pi)$,
in which $\varphi$ is the local longitude of the observer in the body-fixed
coordinate system.

The scattering weight-factor $w_{\rm f}$ is introduced in the combined
Lambert-Lommel-Seeliger law via

$C_{\rm L}(\psi_{i},\psi_{\rm o,i},\alpha,w_{\rm f})=f(\alpha)\left(w_{\rm
f}+\frac{1}{\psi_{i}+\psi_{o,\rm i}}\right),$ (3)

where $C_{\rm L}$ represents a correction coefficient to the Lambertian
reflection, $\psi_{i}$ and $\psi_{\rm o,i}$ are the cosines of the incident
angle and emergence angle on facet $i$ respectively, $\alpha$ is the solar
phase angle; $f(\alpha)$ is the phase correction function, according to
(Kaasalainen, Torppa & Muinonen, 2001),

$f(\alpha)\sim 0.5\exp(-\alpha/0.1)-0.5\alpha+1.$

Parameter $w_{\rm f}$ represents the weight of Lambertian term in the
scattering law, and is an artificial factor to interpret sunlight refection.
Its physical significance isn’t that clear, thus we only need a scattering
weight-factor $w_{\rm f}$ that could achieve best-fitting degree to the
observations.

The effective diameter $D_{\rm eff}$, defined by the diameter of a sphere with
the same area to that of the shape model, can be related to the geometric
albedo $p_{v}$ and absolute visual magnitude $H_{v}$ via:

$D_{\rm eff}=\frac{1329\times 10^{-H_{v}/5}}{\sqrt{p_{v}}}~{}(\rm km)~{}.$ (4)

In addition, the geometric albedo $p_{v}$ is related to the effective Bond
albedo $A_{\rm eff,B}$ by

$A_{\rm eff,B}=p_{v}(0.290+0.684G)~{},$ (5)

in which $G$ is the slope parameter in the $H,G$ magnitude system of Bowell et
al. (1989). As for Beagle, this work uses $H_{v}=10.8$ and $G=0.15$ (IAU MPC).

So parameters including geometric albedo $p_{\rm v}$, roughness $\theta_{\rm
RMS}$, and mean grain radius $b$, are the free parameters that would be
evaluated from the fitting procedures.

Table 4: Fixed input parameters of (656) Beagle. Properties | Value
---|---
Facets of shape model | 1144
Nodes of shape model | 574
Spin ecliptic orientation | ($237^{\circ}$,$86^{\circ}$)
Spin period | 7.0331 hr
Absolute visual magnitude $H_{v}$ | 10.8
Slope parameter $G$ | 0.15
Mean thermal emissivity $\varepsilon$ | 0.9
Regolith grain density $\rho_{\rm g}$ | $3110~{}{\rm kgm^{-3}}$
Regolith porosity $\phi$ | 0.5

###  3.2 Wavelength-dependent Thermal Emissivity

In general cases, an constant thermal emissivity $\epsilon\approx 0.9$ is
assumed for each band in the model to calculate the thermal emission
component. However, in the case of Beagle, we find that, if using wavelength-
dependent thermal emissivity, the model results match the observations better,
as shown in Figure 2.

Figure 2: The observation/model ratios as a function of wavelength (upper
panel) and solar phase angle (under panel) for the case of best-fit
parameters. Panel A (left): best-fit model results are obtained by using
constant input thermal emissivity $\epsilon\approx 0.9$; Panel B (right):
using wavelength-dependent thermal emissivity $\epsilon\approx
0.72,0.72,0.72,0.95$ for the W1, W2, W3, W4 band respectively.

After a process of adjusting the thermal emissivity to fit the observations,
we find that, the fitting degree can be the best by using $\epsilon\approx
0.72,0.72,0.72,0.95$ for the W1, W2, W3, W4 band respectively. The thermal
emissivity of band W4 is almost close to 1, indicating that band W4 is near
the Wien peak of the thermal emission spectra of Beagle.

On the other hand, with RSTPM, the fraction of sunlight-reflection in the
observed flux of Beagle for each band of WISE/NEOWISE can be evaluated, the
results are plotted in Figure 3. The W1-band is reflection dominated with
reflection ratio $\sim 95\%$, the W2-band becomes thermal emission dominated
but have significant reflection ratio between $10\sim 45\%$, finally the W3
and W4 are almost all thermal emission with reflection ratio $<0.1\%$.

Figure 3: The fraction of sunlight-reflection in the observed flux of Beagle
for each band of WISE/NEOWISE.

###  3.3 Results of Roughness, Grain Size and Albedo

By adopting different thermal emissivity for the four band, we then fit the
observations by scanning roughness $\theta_{\rm RMS}$ in the range of $0\sim
50^{\circ}$ and mean grain radius $b$ in the range of $1\sim 1000~{}\mu$m.
With each pair of ($\theta_{\rm RMS}$,$b$), a best-fit geometric albedo
$p_{\rm v}$ together with effective diameter $D_{\rm eff}$ is found to compute
the reduced $\chi^{2}_{\rm r}$. The results are presented in Figure 4 as a
contour of $\chi^{2}_{\rm r}$($\theta_{\rm RMS}$,$b$).

Figure 4: The contour of $\chi_{\rm r}^{2}(\theta_{\rm RMS},b)$, which is
obtained by fitting the observations with two free parameters: the roughness
RMS slope $\theta_{\rm RMS}$ and the mean grain radius $b$. The blue region
stands for the $1\sigma$-level constraint, and the red region represents the
$3\sigma$-level constraint.

According to Figure 4, a well constrained $1\sigma$-level limit is derived for
roughness $\theta_{\rm RMS}$ and mean grain radius $b$, giving $\theta_{\rm
RMS}=44\pm 2^{\circ}$, $b=100^{+130}_{-70}~{}\mu$m respectively.
$3\sigma$-level constraint for roughness is derived as $\theta_{\rm RMS}=44\pm
4^{\circ}$, whereas for mean grain radius, a relatively wide $3\sigma$-level
limit is obtained as $b=100^{+350}_{-90}(10\sim 450)~{}\mu$m.

According to the above derived $1\sigma$ and $3\sigma$ ranges of roughness and
mean grain radius, the corresponding geometric albedo $p_{\rm v}$ and
$\chi^{2}_{\rm r}$ are picked out, leading to the $p_{\rm v}\sim\chi^{2}_{\rm
r}$ relation as shown in Figure 5. In this way, we obtain the $1\sigma$ and
$3\sigma$-level limits of geometric albedo as $p_{\rm
v}=0.05^{+0.002}_{-0.003}$ and $p_{\rm v}=0.05^{+0.004}_{-0.007}$
respectively, and simultaneously the effective diameter of Beagle is obtained
as $D_{\rm eff}=57.3^{+1.8}_{-1.1}$ km ($1\sigma$) and $D_{\rm
eff}=57.3^{+4.5}_{-2.2}$ km ($3\sigma$).

Figure 5: $p_{\rm v}\sim\chi^{2}_{\rm r}$ profiles fit to the observations in
consideration of the derived $1\sigma$ and $3\sigma$ ranges of roughness and
mean grain radius.

###  3.4 Regolith thermophysical characteristics

Thermal inertia is a strong function of temperature. Now with the above
derived profile of mean grain radius, we can evaluate the change of surface
thermal inertia of Beagle due to the influence of seasonal temperature
variation according to the relationships between thermal inertia, thermal
conductivity, specific heat capacity and temperature given in Yu & Ip (2021).

In the left panel of Figure 6, a map of surface temperature of Beagle is
plotted as a function of local latitude and orbital mean anomaly. Each
temperature has been averaged over one rotational period. We can clearly see
that temperature on each local latitude can reach maximum (summer) or minimum
(winter) at different orbital positions as a result of seasonal effect.
Temperature on the poles can vary from $\sim$26 K to $\sim$183 K.

Figure 6: Left panel: Seasonal variation of the diurnal-averaged surface
temperature as a function of local latitude. Right panel: Seasonal average
temperature of each local latitude. The so-called local latitude is defined as
the complementary angle of the angle between the local normal vector and the
rotation axis.  Figure 7: The regolith specific heat capacity $c_{\rm p}$,
thermal conductivity $\kappa$, and thermal inertia $\Gamma$ of Beagle in
consideration of seasonal temperature variation and $3\sigma$-level range of
mean grain radius $b$. $\Gamma=\sqrt{(1-\phi)\rho_{\rm g}c_{\rm p}\kappa}$,
and porosity $\phi\sim 0.5$, regolith grain density $\rho_{\rm
g}=3110\rm~{}kgm^{-3}$.

Figure 7 shows that, if considering seasonal temperature variation and
$3\sigma$-level range of mean grain radius $b$, the regolith thermal inertia,
thermal conductivity, specific heat capacity of Beagle can be estimated to
vary between $2\sim 45\rm~{}Jm^{-2}s^{-0.5}K^{-1}$, $0.3\sim 2.3\times
10^{-3}\rm~{}Wm^{-1}K^{-1}$, $7\sim 562\rm~{}JKg^{-1}K^{-1}$.

Despite the fact that thermal parameters is temperature dependent, the
majority of relevant works ignore such temperature dependence and only
estimate the average values. For comparison with such existing results, we
estimate the seasonal average thermal parameters of Beagle with inputs of the
derived mean grain radius and seasonal averaged temperature. The seasonal
average temperature would be a function of local latitude,
$\tilde{T}(\theta)$, and can be estimated as

$(1-A_{\rm eff,B})\tilde{L}_{\rm
s}(\theta)=\varepsilon\sigma\tilde{T}(\theta)^{4},$ (6)

where $A_{\rm eff,B}$ is the bond albedo, $\varepsilon\sim 0.9$ is the average
thermal emissivity, $\tilde{L}_{\rm s}(\theta)$ is the annual average incoming
solar flux on each latitude. The results are presented in The right panel of
Figure 6, giving the seasonal average temperature of Beagle to be $90\sim
170\rm~{}K$, and accordingly the average specific heat capacity

$c_{\rm p}=173\sim 516\rm~{}JKg^{-1}K^{-1},$

the average thermal conductivity

$\kappa=0.7\sim 1.3\times 10^{-3}\rm~{}Wm^{-1}K^{-1},$

and the average thermal inertia

$\Gamma=\sqrt{(1-\phi)\rho_{\rm g}c_{\rm p}\kappa}=14\sim
32\rm~{}Jm^{-2}s^{-0.5}K^{-1},$

where the porosity $\phi$ uses 0.5, and the regolith grain density $\rho_{\rm
g}$ uses $3110\rm~{}kgm^{-3}$ (Opeil et al., 2010).

###  3.5 Infrared light curves

With the above derived physical and thermophysical parameters of Beagle, now
we are able to determine the rotation phase of Beagle at each time of
WISE/NEOWISE observation by doing comparisons between the observational and
theoretical light curves. To do so, the 3D shape model is used to define the
local body-fixed coordinate system, where the z-axis is chosen to be the
rotation axis. Moreover, if we define the view angle of one observation with
respect to the body-fixed coordinate system to be $(\varphi,\theta)$, where
$\varphi$ stands for local longitude, and $\theta$ means local latitude, then
the rotational phase $ph$ of this observation can be related to the local
longitude $\varphi$ via

$ph=1-\varphi/(2\pi),$ (7)

and ”zero” rotational phase is chosen to be the ”Equatorial view
($0^{\circ}$)” as shown in Figure 1.

If selecting a reference epoch, and assuming the rotational phase at this
epoch to be $zph$, then all the rotational phases of other data could be
derived in consideration of the observation time and geometry. Furthermore,
for some particular epoch, light curves can be derived for each band by
correcting the observed flux at various epochs into one rotation period at
this epoch, where the correction is implemented via

$F_{i,\rm corr}=F_{i}\left(\frac{r_{i,\rm helio}}{r_{0,\rm
helio}}\right)^{2}\left(\frac{\Delta_{i,\rm obs}}{\Delta_{0,\rm
obs}}\right)^{2},$ (8)

in which $F_{i,\rm corr}$ is the flux after correction, $F_{i}$ is the
original observed flux, $r_{i,\rm helio}$ and $r_{0,\rm helio}$ are the
heliocentric distance of epoch $i$ and the reference epoch, while
$\Delta_{i,\rm obs}$ and $\Delta_{0,\rm obs}$ are the observation distance.

Following the above method, firstly we select ’2010-01-31 12:06’ as the
reference epoch for deriving the reference rotational phase $zph$. But for
correction of flux, in order to reduce flux errors caused by correction
Equation (8), we select 28 separate reference epochs, so as to use data close
to each reference epoch (data within three days) to generate infrared light
curves. By this way, we obtain 28 infrared light curves. Then for each of the
reference epoch, theoretical infrared light curves are simulated by RSTPM to
fit the above generated observational light curves. The best-fit results are
plotted in Figure 8. Roughly speaking, our modeled infrared light curves match
well with all the four-band observational light curves of WISE/NEOWISE in
consideration of the observation errors, indicating that our model results are
reliable.

Figure 8: Best-fit results to the infrared light curves of WISE/NEOWISE for
each observation epoch; Panel A,B,C: thermal light curves at band W4, W3, W2
respectively; Panel D: reflection light curves at band W1.  Figure 9: Ratios
of observation/model as a function of rotation phase for each band
respectively.

To investigate whether Beagle has surface heterogeneity, the ratios of
observation/model are plotted as a function of rotation phase for each band
respectively in Figure 9, where a clear rotation phase dependent feature is
observed in the data of band-W4, W3, and W2, showing that the
observation/model ratios appear to be significantly higher around $ph=0.2$ and
0.7, whereas the band-W1 data don’t show such feature. Since both irregular
shape and surface heterogeneity can contribute to the rotational variation of
light-curves, and bands of W4, W3, and W2 are thermal emission dominated, so
the heterogeneous features of the thermal light curves of WISE/NEOWISE may
imply that: 1, the light-curve inversion shape model is not perfect for
modeling thermal emission; or 2, the regolith of Beagle may have heterogeneous
thermophysical characteristics along longitude.

##  4 Discussion

The multi-epoch infrared data of WISE/NEOWISE make it possible to probe the
regolith characteristics of asteroids in detail. As an example, this work uses
9 years of infrared data from WISE/NEOWISE to study the regolith of Main-Belt
Object (656) Beagle by the well-tested thermophysical model — RSTPM (Yu & Ip,
2021). Details about the data processing, fitting procedure, modelling and
analyzing of infrared light curves are presented. The results show that Main-
Belt Object Beagle has a low mean thermal inertia of $\Gamma=14\sim
32\rm~{}Jm^{-2}s^{-0.5}K^{-1}$, being consistent with that of most large main
belt objects (MacLennan & Emery, 2021, 2022; Huang et al., 2022); and Beagle
has a geometric visual albedo as low as $0.043\sim 0.054$, being far lower
than the geometric albedos of its neighbouring asteroids (mean $p_{\rm
v}=0.0941\pm 0.0055$) (Fornasier et al., 2016), even lower than the geometric
albedo of (24) Themis ($p_{\rm v}=0.064^{+0.008}_{-0.011}$) (Yu & Ip, 2021).
Although the low albedo $\sim 0.05$ of Beagle still lies on the edge of
$1\sigma$ range of the Beagle family and the Themis family, it is a little
anomalous that the albedos of Beagle’s neighbouring asteroids are more close
to Themis, rather than Beagle itself.

In addition, Figure 9 shows that Beagle exhibits significant rotation-phase
dependent features at band-W4, W3 and W2, but not at W1-band. According to
Figure 3, band-W4, W3 and W2 are thermal emission dominated, whereas W1-band
is dominated by NIR reflection. So the result implies that the surface of
Beagle doesn’t have significant heterogeneous NIR reflectivity. As mentioned
before, Fornasier et al. (2016) shows that the neighbouring asteroids of
Beagle have diverse NIR spectral types, and Beagle is expected to show
heterogeneous NIR features across its surface if Beagle is the parent of those
neighbouring asteroids. However, here our result shows that Beagle has no
significant heterogeneous NIR reflectivity, indicating that Beagle may be not
the parent of its neighbouring asteroids, or surface heterogeneity of Beagle
at NIR are eliminated for some reasons, for example, the utilization of
imperfect light-curve inversion shape model, which doesn’t remove the
degeneracy between regional NIR reflectivity and slope.

On the other hand, for band-W4, W3 and W2, if doing comparisons between Figure
9 and 1, we can find that the higher observation/model ratios appear around
$ph=0.2$ and 0.7, the phases of which Beagle happens to have the minimum
cross-section area, implying that the rotation-phase dependent feature of
band-W4, W3 and W2 has a strong relationship with the utilized shape model.
Also we have done several tests by adjusting the rotation period within
$7.033\pm 0.001$ to do the fitting procedure, but the above rotation-phase
dependent feature remains unchange, indicating that the utilized shape rather
than rotation period is more likely to be cause of the rotation-phase
dependent feature. Hence the light-curve inversion shape model from optical
light curves is not perfect for infrared light curves, similar phenomenons
have been found by Hanus et al. (2015) and Ďurech, et al. (2017).
Nevertheless, the possibility of heterogeneous regolith thermophysical
characteristics across the surface of Beagle remains positive, because it
would be reasonable to expect different physical characteristics (e.g.
roughness) between the body part and the smaller head part, as such
phenomenons have been observed on other asteroids by in-situ space missions,
for example (25143) Itokawa (Abe et al., 2006).

##  5 Conclusion

By analyzing 9-years infrared light curves of Beagle from WISE/NEOWISE, we get
the following results:

1). Beagle has an effective diameter $D_{\rm eff}=57.3^{+4.5}_{-2.2}$ km,
geometric albedo $p_{\rm v}=0.05^{+0.004}_{-0.007}$, mean roughness
$\theta_{\rm RMS}=44\pm 4^{\circ}$, mean grain size $b=100^{+350}_{-90}(10\sim
450)~{}\mu$m, mean specific heat capacity $c_{\rm p}=173\sim
516\rm~{}JKg^{-1}K^{-1}$, thermal conductivity $\kappa=0.7\sim 1.3\times
10^{-3}\rm~{}Wm^{-1}K^{-1}$ and thermal inertia $\Gamma=14\sim
32\rm~{}Jm^{-2}s^{-0.5}K^{-1}$.

2). We confirm that the Beagle has an anomalous low albedo that the albedos of
Beagle’s neighbouring asteroids are more close to Themis, rather than Beagle
itself. In addition, the W1-band NIR light curves of Beagle don’t reveal
significant heterogeneous NIR reflectivity across the surface of Beagle. This
result does not support Beagle to be the parent of its neighbouring asteroids
which have diverse NIR spectral types (Fornasier et al., 2016). A possible
explanation to these results is that Beagle may be an interloper or a sister,
rather than the parent of its neighbouring asteroids including the first MBC
133P. Considering Beagle’s albedo anomaly, we prefer to surmise Beagle to be
an interloper to the Themis family, although it still remains possible for
Beagle being an anomalous member of the Themis family. So our results may add
new clues of Beagle probably having no genetic connection with its
neighbouring asteroids and even (24) Themis, and may lead to new scenarios
about the origin of the famous MBC 133P.

3). Asteroidal shape models from inversion of optical light curves could have
obvious imperfections. As a result, the heterogeneity of surface reflectivity
at near infrared is difficult to be discovered if using such shapes to model
the theoretical infrared light curves. Besides, the imperfect shapes can also
lead to a rotation dependent features in thermal light curves (e.g. band W2,
W3, W4 of WISE and NEOWISE), which would mislead the evaluation of the
heterogeneity of regolith thermophysical characteristics. Therefore, to obtain
truthful information about the surface heterogeneity of asteroids, it is
necessary to input shapes that are more close to reality.

## Acknowledgments

We would like to thank the WISE teams for providing public data. This work was
supported by the grants from The Science and Technology Development Fund,
Macau SAR (File No. 0051/2021/A1). and Faculty Research Grants of The Macau
University of Science and Technology (File).

## References

  * Abe et al. [2006] Abe, M., Takagi, Y., Kitazato, K., Abe, S., et al., 2006. Near-Infrared Spectral Results of Asteroid Itokawa from the Hayabusa Spacecraft, Science 312: 1334-1338
  * Bowell et al. [1989] Bowell, E., Hapke, B., Domingue, D., et al., 1989. Application of photometric models to asteroids. In _Asteroids II_ , pp. 524-556
  * Campins et al. [2010] Campins, H., Hargrove, K., Pinilla-Alonso, N., et al., 2010, Nature, 464(7293): 1320-1321
  * Carruba [2019] Carruba, V. 2019, P&SS, 166, 90
  * Cutri et al. [2015] Cutri, R., Mainzer, A., Conrow, T., et al. 2015, Explanatory Supplement to the NEOWISE Data Release Products, http://wise2.ipac.caltech.edu/docs/release/neowise/expsup/
  * Cutri et al. [2012] Cutri, R., Wright, E., Conrow, T., et al. 2012, Explanatory Supplement to the WISE All-Sky Data Release Products, http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/index.html
  * Ďurech, et al. [2017] Ďurech, J., Delbo’, M., Carry, B., Hanuš, J., & Alí-Lagoa, V., 2017, A&A, 604, A27
  * MacLennan & Emery [2021] MacLennan, E. M., & Emery, J. P., 2021, Planet. Sci. J., 2(161): 1-12
  * MacLennan & Emery [2022] MacLennan, E. M., & Emery, J. P., 2022, Planet. Sci. J., 3(47): 1-23
  * Fornasier et al. [2016] Fornasier, S., Lantz, C., Perna, D., & Campins, H., et a., 2016, Icarus, 269: 1-14
  * Hanus et al. [2015] Hanuš, J., Delbo’, M., Ďurech, J., & Alí-Lagoa, V., 2015, Icarus, 256, 101
  * Huang et al. [2022] Huang, D., Hanuš, J., Masiero, J. R., & Tholen, D. J., 2022, Planet. Sci. J., 3(56): 1-25
  * Kaasalainen, Torppa & Muinonen [2001] Kaasalainen, M., Torppa, J., & Muinonen, K., 2001, 153, 37
  * Mainzer et al. [2011] Mainzer, A., Bauer, J., Grav, T., et al. 2011, ApJ, 731, 53
  * Mainzer et al. [2014] Mainzer, A., Bauer, J., Cutri, R., et al., 2014, ApJ, 792, 30
  * Mainzer et al. [2016] Mainzer, A. K., Bauer, J. M., Cutri, R. M., et al. 2016, PDSS, EAR-A-COMPIL-5-NEOWISEDIAM-V1.0
  * Michel et al. [2015] Michel P., Richardson D. C., Durda D. D., Jutzi M., and Asphaug E., 2015, Collisional formation and modeling of asteroid families. In Asteroids IV (P. Michel et al., eds.), pp. 341–354. Univ. of Arizona, Tucson
  * Nesvorný et al. [2008] Nesvorný, D., Bottke, W. F., Vokrouhlický, D., et al., 2008, ApJL, 679: 143-146
  * Opeil et al. [2010] Opeil, C. P., Consolmagno, G. J., Britt, D. T., 2010, Icarus, 208, 449-454
  * Rivkin & Emery [2010] Rivkin, A.S., & Emery, J. P., 2010, Nature, 464(7293): 1322-1323
  * Rousseau, De Sanctis & Raponi et al [2021] Rousseau, B., De Sanctis, M. C., Raponi, A., & Ciarniello, M., et al., 2021, A&A, 653, 118
  * Wright et al. [2010] Wright, E. L., Eisenhardt, P. R. M., & Mainzer, A. K., et al., 2010, AJ, 140, 1868
  * Yu et al. [2020] Yu, L. L., Hsia, C.H., & Ip, W. H., 2020, AJ, 159(66): 1-10
  * Yu & Ip [2021] Yu L. L., & Ip W. H., 2021, ApJ, 913(96): 1-22