# Lagrangian particle simulation of hydrogen pellets and SPI into runaway
electron beam in ITER

Shaohua Yuan1, Nizar Naitlho1, Roman Samulyak1,2, Bernard Pégourié3, Eric
Nardon3, Eric Hollmann4, Paul Parks5, Michael Lehnen6 1 Department of Applied
Mathematics and Statistics, Stony Brook University, Stony Brook, NY, USA  2
Computational Science Initiative, Brookhaven National Laboratory, Uptown, NY,
USA  3 CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France  4 University of
California - San Diego, La Jolla, CA, USA  5 General Atomics, San Diego, CA,
USA  6 ITER Organization, Route de Vinon sur Verdon, 13115 St Paul Lez
Durance, France<EMAIL_ADDRESS>
###### Abstract

Numerical studies of the ablation of pellets and shattered pellet injection
(SPI) fragments into a runaway electron beam in ITER have been performed using
a time-dependent pellet ablation code [R. Samulyak at el., Nucl Fusion, 61
(4), 046007 (2021)]. The code resolves detailed ablation physics near pellet
fragments and large-scale expansion of ablated clouds. The study of a single
fragment ablation quantifies the influence of various factors, in particular
the impact ionization by runaway electrons and cross-field transport models,
on the dynamics of ablated plasma and its penetration into the runaway beam.
Simulations of SPI performed using different numbers of pellet fragments study
the formation and evolution of ablation clouds and their large-scale dynamics
in ITER. The penetration depth of ablation clouds is found to be of the order
of 50 cm.

##  1 Introduction

Recent experimental results from DIII-D and JET [15, 16] suggest that pure D2
(or H2) shattered pellet injection (SPI) into a Runaway Electron (RE) beam may
be a promising RE mitigation strategy for ITER. In DIII-D and JET,
sufficiently large D2 injection is observed to cause rapid ($\sim$20 ms or
less) thermal background plasma volume recombination and a simultaneous radial
transport of background high-Z (argon) impurities out of the RE beam. Upon
subsequent compression against the plasma wall, these “purged” RE beams are
then observed to be lost to the wall in a single large MHD instability with
large RE wetted area, low local heat loads, and nearly full conversion of RE
current into thermal plasma current [15, 16]. Presently, the physics of this
unusually large, single RE loss event is not understood. Experimentally,
however, the large MHD event appears to depend on achieving background plasma
volume recombination. In present devices, background plasma volume
recombination is understood as resulting from rapid neutral diffusion of D and
D2 into the vessel volume and subsequent molecular recombination of plasma
ions [9]. The initial D2 deposition profile from the pellet ablation does not
appear to be important in present devices [20]. In ITER, however, due to the
larger plasma minor radius, it is expected that RE plateau diffusive
timescales will be longer than in present devices [8], so the initial pellet
deposition profile may become important for determining whether or not plasma
volume recombination can be achieved on timescales less than the RE plateau
vertical drift time ( 100 ms). Preliminary work based on a simple 1D model
suggests that gas cloud expansion in the vicinity of SPI shards ablated by REs
may improve injected D penetration into the RE beam [11]. The main objective
of this paper is to improve on this work and provide a numerical estimate of
penetration of ablated mass into a RE beam resulting from the expanding SPI D2
shard ablation plumes.

Simulations of pellets and SPI using tokamak MHD codes are restricted to
relatively coarse meshes that cannot resolve SPI fragments and fine physics
details of the ablation process. In this work, we use a highly adaptive 3D
Lagrangian particle pellet code [19] based on the Lagrangian particle method
[18] for hydrodynamic equations. The code is capable of resolving a large
spectrum of spatial scales ranging from sub-millimeter scales associated with
phase transition processes on pellet surfaces and steep density gradients in
cold ablation clouds to tokamak-scale expansion of ablation clouds. By
resolving relevant physics processes near pellet / fragment surfaces and in
the far field, the LP code self-consistently computes pellet ablation rates
and shows their dependence on 3D effects. The code was validated using
experimental data on the injection of hydrogen fueling pellets and DIII-D
experiment with small neon pellets [10].

The paper is organized as follows. In Section 2, we describe the problem setup
and main physics parameters. Main models, governing equations and their
numerical implementation are presented in Section 3. Results of numerical
simulations of a single fragment injection and shattered pellet injection
(SPI) into a runaway electron beam in ITER are presented in Section 4. We
conclude the paper with a summary of our results and plans for the future
work.

##  2 Problem setup

Consider a large, cylindrical, pure hydrogen pellet with 28.5 mm diameter and
57 mm length that is presently the foreseen size of large pellets in the ITER
disruption mitigation system. The pellet is injected with the velocity of
$V=500$ m/s into a runaway electron beam in ITER. Before entering the plasma,
the pellet is shattered into a large number of fragments creating a cloud with
a 20 degree half angle and the injection time of about 2.2 ms.

Schematic of the simulation setup is shown in Figure 1. The left image in
Figure 1(a) shows an ITER cross-section containing a runaway electron beam
region outlined by the blue line. The red ellipse denotes the core RE beam
with constant RE density and the layer between the red and blue ellipses is an
RE density ramp with 60 cm thickness. The right image in Figure 1(a) shows an
enlarged distribution of SPI fragments. The real-size fragment distribution is
injected into the RE region in the direction of the red arrow. Figure 1(b)
depicts the profile of relative (with respect to the core) RE density in the
ramp between the RE core and the edge. The x-coordinate in this plot is
consistent with the x-coordinate in the ITER cross-section plot.

(a) Schematic of SPI into RE beam in ITER

(b) Profile of relative RE density

Figure 1: (a) Schematic of SPI into RE beam in ITER. The left image shows ITER
cross-section with a runaway electron beam region. The right image
schematically shows enlarged distribution of SPI fragments. The injection
direction is shown by the red arrow. (b) Profile of the relative RE density in
the ramp between the RE core and the edge.

The RE beam and its companion plasma parameters are approximated as follows:
the RE beam is monoenergetic with $E_{RE}=20MeV$ and the core of the RE beam
has uniform density of $n_{RE}=5\times 10^{16}m^{-3}$ which corresponds to
$10$ MA of runaway current. The companion plasma consists of a mixture of
deuterium and neon with densities $n_{D}=1\times 10^{20}m^{-3}$ and
$n_{Ne}=5\times 10^{19}m^{-3}$, where the neon density corresponds to pure Ne
injection resulting in $50$ ms current quench time. Electron densities are
based on coronal equilibrium ionization stages exclusive of ionization through
the RE beam. The background plasma temperature is $T=20eV$ and the electron
density is $n_{e}=3\times 10^{20}m^{-3}$.

The goal of this simulation study is to characterize the ablation of
fragments, heating of the ablated material by the RE beam, and evolution of
the ablated cloud as it ionizes, channels along magnetic fields lines and
experiences cross-field transport by various drift forces. We are interested
in computing the penetration depth of the ablated cloud into the RE beam and
clarifying the dependence of ablation cloud evolution on the number and size
of pellet fragments.

##  3 Main models, governing equations and their numerical implementation

###  3.1 MHD in low Magnetic Reynolds number approximation

As a cryogenic pellet or an SPI fragment is injected into a cold tokamak
plasma in the presence of a runaway beam, the pellet experiences an intense
volumetric heating by REs and a rapid ablation (sublimation). A cold, neutral,
high-pressure ablation cloud is formed around the pellet that initially
expands spherically by pressure gradients. The expanding ablation material is
heated up by REs, gradually ionizes via the thermal Saha mechanism as well as
by direct impact of fast runaway electrons, and starts flowing predominantly
along magnetic field lines by the Lorentz force action. To describe the
evolution of such a cold, partially ionized ablation cloud, we use the system
of MHD equations in the low Magnetic Reynolds number approximation, ${\delta
B}/B\sim R_{m}<<1$, where $\delta B$ is the eddy current induced magnetic
field. In Lagrangian coordinates, the equations are

$\displaystyle\frac{d\rho}{dt}=-\rho\nabla\mathbf{u},$ (1)
$\displaystyle\rho\frac{d\mathbf{u}}{dt}=-\nabla P+\mathbf{J\times
B}+\mathbf{f}_{D},$ (2)
$\displaystyle\rho\frac{de}{dt}=-P\nabla\cdot\mathbf{u}+\frac{1}{\sigma}\mathbf{J}^{2}+Q_{RE},$
(3) $\displaystyle P=P(\rho,e),$ (4)

where $d/dt$ is the Lagrangian time derivative, $\mathbf{u}$, $\rho$ and $e$
are the velocity, density and specific internal energy, respectively, $P$ is
the pressure, $\mathbf{B}$ is the magnetic field induction, $\mathbf{J}$ is
the current density, and $\sigma$ is the fluid conductivity. The runaway
electron heat flux is represented by an external heat source $Q_{RE}$. Viscous
forces in the ablation cloud are neglected. As we show in Section 4, the
ablated material heating by cold background plasma via heat conduction is
negligibly small in the energy balance in the ablation cloud, hence excluded
in the equations. Finally, the ablation cloud experiences various drift forces
in the transverse direction to the magnetic field represented by the term
$\mathbf{f}_{D}$. The corresponding forces are explained in detail in Section
3.3, dedicated to the cross-field transport.

The current density for (2 \- 3) is obtained from Ohm’s law

$\mathbf{J}=\sigma(-\nabla\phi+\mathbf{u\times B}),$ (5)

where $\phi$ is the electric potential in the cloud. For a general 3D problem,
the electric potential must be found from the following Poisson equation

$\nabla\sigma\nabla\phi=\nabla\cdot(\mathbf{u\times B})$

that follows from the charge conservation equation $\nabla\cdot\mathbf{J}=0$,
subject to an appropriate boundary condition. In this work, we assume that the
ablation cloud is uniformly charged resulting in a constant value of $\phi$.
The current density becomes

$\mathbf{J}=\sigma_{\perp}\mathbf{u\times B},$

where ${\mathbf{B}}=(0,0,B_{\parallel})$ and the transverse conductivity
$\sigma_{\perp}$ is approximated by a model derived for hydrogenic species in
[17]. Similar approximations for ablation cloud were used in simulations of
pellet ablation by hot thermal plasma electrons [19, 3]

The EOS model used in this study implements two sources of ionization. First,
it employs the local thermodynamic equilibrium (LTE) Saha model [21]. The
second component of the model is based on the direct impact ionization of the
ablated hydrogen by energetic REs. The electron impact ionization rates are
computed as

$R=N_{e}(r)\int v\sigma(E)f_{e}(E)dE,$ (6)

where $v$ is the electron impact velocity, $\sigma(E)$ is the collision cross-
section, $f_{e}$ is the electron energy distribution function, and $N_{e}(r)$
is the RE density that depends on the radial location in the tokamak as shown
in Figure 1(b). Electron impact ionization rates for the parameter space of
our problems of interest have been obtained courtesy of N. Garland [7].

We would like to comment on omitting hydrogen radiation losses in equation
(3). Using simulation data presented in later sections, we extracted time-
dependence of hydrodynamic states on Lagrangian particles at various phases of
the ablation cloud evolution and used the collisional-radiative code
PrismSPECT [1] to compute hydrogen radiation power density. The results show
that hydrogen radiation power density is always much smaller compared to the
RE heating power density. Only at the end of the process when the cloud
density is very low, hydrogen radiation becomes of the same order but still
smaller compared to the RE heating. The radiation of REs is always much higher
compared to the radiation of hydrogen. Since a small fraction of radiation is
absorbed by the hydrogen cloud, the RE radiation energy source partially
compensates the hydrogen radiation energy sink. Since adding using full
radiation hydrodynamics algorithms is not feasible for the present problems,
neglecting radiation seems to be a reasonable 1st approximation.

###  3.2 Runaway electron heat flux

The energy of RE in our model is assumed to be in the range of 20 MeV.
Consider an incident RE beam passing through a pellet or ablated material with
atomic density $n_{0}$ until it is stopped. The following expression [4] gives
the stopping power of the ultrarelativistic electrons:

$\frac{d}{dt}\frac{\epsilon_{RE}}{m_{e}c^{2}}=-4\pi
r^{2}_{e}c\ln\Lambda_{free}(n_{free}+\frac{1}{2}n_{bound}),$ (7)

where $\epsilon_{RE}$ is the RE energy, $r_{e}\equiv e^{2}/mc^{2}$ is the
classic electron radius, $\ln\Lambda_{free}\approx 21$ is the Coulomb
logarithm, and $n_{free}$ and $n_{bound}$ are the free and the bound electron
number densities, respectively. In the solid state, nearly all the electrons
are bound so we have $n_{free}=0$ and $n_{bound}=Zn_{0}$, where $Z$ is the
atomic number of the pellet material. Equation (7) gives the stopping length
of RE [11] as

$L\approx\frac{\epsilon_{RE}}{m_{e}c^{2}}\frac{1}{2\pi r_{e}^{2}\ln\Lambda
Zn_{0}}.$ (8)

For the solid hydrogen density $n_{0}=4.5\times 10^{28}\,1/m^{-3}$, the
stopping length of RE is $L_{H}\approx 78$ cm. Therefore, we conclude that REs
pass through the pellet / ablated clouds multiple times before being stopped
and provide volumetric heating. The heat source $Q_{RE}$ from REs streaming
into the pellet and the ablation cloud can be computed from (7) as

$Q_{RE}=4\pi
m_{e}r^{2}_{e}c^{3}n_{RE}\ln\Lambda_{free}(n_{free}+\frac{1}{2}n_{bound})$ (9)

###  3.3 Transverse flow dynamics

A cold and weakly ionized ablation cloud polarizes as it travels across
magnetic field lines and the resulting polarization $\mathbf{E\times B}$ drift
attempts to maintain the cloud motion with the initial injection velocity. In
addition, the cloud polarization evolves, due to (a) the cloud-plasma pressure
gradient which tends to maintain it and (b) the Alfvén wave emission and
internal and external cloud self-connections along field lines which tend to
reduce it [12]. The resulting equation governing the $\mathbf{E\times B}$
drift velocity is:

$\displaystyle\frac{dv_{D}}{dt}=\frac{1}{1+(1-P_{Alf}-P_{con})(L^{\prime}_{con}\rho_{\infty}/Z_{0}\rho)}\times$
(10)
$\displaystyle\left[\frac{qR}{Z_{0}}\sin\left(\frac{Z_{0}}{qR}\right)\frac{2<P(1+\frac{M_{\parallel}^{2}}{2})-P_{\infty}>}{R<\rho>}-v_{D}\frac{B_{||}^{2}}{<\rho>}H(\tau_{Bdiff}-\tau_{adv})\times\right.$
$\displaystyle\left.\left\\{P_{Alf}\frac{2}{\mu_{0}v_{A}}+P_{con}\frac{\left[1-\exp(-t/(\tau_{coll}^{e}+\tau_{L}))\right]\pi
r_{0}^{2}\\-\sigma_{\infty}}{L_{con}}\right\\}\right].$

In this expression,

$<A>\equiv\int_{-Z_{0}}^{Z_{0}}Adz$, the integral is along a magnetic field
line, $A$ is any quantity,

$Z_{0}$: half-cloud length,

$v_{D}$: drift velocity,

$P_{con}$: proportion of cloud section self-connected (+ vs. -) at time $t$,

$P_{Alf}$: proportion of cloud section not self-connected at time $t$,

$L^{\prime}_{con}$: arithmetic average of the length of the field lines self-
connecting the cloud,

$L_{con}$: harmonic average of the length of the field lines self-connecting
the cloud,

$B_{||}$: background magnetic field,

$\rho$: cloud density,

$P$: cloud pressure,

$\rho_{\infty}$: background plasma density,

$P_{\infty}$: background plasma pressure,

$r_{0}$: cloud radius,

$q=2$: tokamak safety factor,

$R$: tokamak major radius,

$M_{\parallel}$: Mach number of the plasmoid parallel expansion,

$\sigma_{\infty}$: plasma parallel conductivity,

$\tau_{coll}^{e}$: plasma electron collision time,

$\tau_{L}=\sigma_{\infty}\mu_{0}P_{con}r_{0}^{2}$: damping current delay time,

$H$: the Heaviside function that depends on

$\tau_{Bdiff}=\mu_{0}\sigma_{\perp}L_{char}^{2}$: characteristic time for
magnetic field diffusion over characteristic length $L_{char}=5$ cm,

$\tau_{adv}=L_{char}/v_{D}$: characteristic time for cloud motion across
magnetic field lines.

After the ablation cloud ionizes and expands along magnetic field lines, it
starts drifting across the magnetic field in the direction of the outward
major radius $R$. This motion, described by the left term in the second row of
equation (10), has been attributed to $\mathbf{E\times B}$ drift resulting
from a vertical curvature and grad-B drift polarization induced inside the
ionized ablated material by the $1/R$ toroidal magnetic field variation [2,
14, 13]. As the cloud continues to heat-up and ionize, the characteristic
magnetic field diffusion time $\tau_{Bdiff}$ becomes longer compared to the
characteristic time for cloud motion across magnetic field lines,
$\tau_{adv}$, i.e. magnetic flux becomes frozen into the cloud. At this time,
the Alfven wave drag and the external connection (Pégourié) current - induced
drag, described by the two terms in the third line of equation (10), become
essential [12, 6]. Formally, these terms are added as the Heaviside function
in the second line of equation (10) changes from zero to one. As we mentioned
above, the plasmoid polarizes in the transverse direction as it moves across
the magnetic field lines. This creates a rapidly expanding Alfven-wave induced
polarization in the ambient plasma. After time
$\tau_{con}=L^{\prime}_{con}/v_{A}$, the plasmoid becomes self-connected by
this polarization perturbation. Proportion of the cloud section not self-
connected at this time, $P_{Alf}$, contributes to the Alfven wave drag force.
Proportion of the cloud section connected to the region of opposite polarity,
$P_{con}$, contributes to the external connection current drag term.
Proportion of the cloud section connected to the region of the same polarity,
$P_{same}=1-P_{con}-P_{Alf}$, increases the effective moving mass as indicated
in the first term on the right-hand-side in equation (10). Further details,
including the topology of self-connection paths in a tokamak, can be found in
[12, 6].

###  3.4 Numerical implementation

The system of equations (1) - (4) is solved by the Lagrangian Particle (LP)
pellet ablation code [19] based on the Lagrangian particle method [18] for
hydrodynamic equations. The code is applicable to pellets and SPI fragments
ablated by hot plasma and REs. The LP method represents fluid cells with
Lagrangian particles and is suitable for the simulation of complex free
surface or multiphase flows. It is highly adaptive to density changes, a
critically important property for 3D simulations of the ablation of pellets
and, especially, SPI fragments. The LP-based pellet ablation code implements
MHD equations in the low magnetic Reynolds number approximation and the
physics models described above. The LP code is optimized for massively
parallel supercomputers by using p4est (”parallel forest of K-trees” ) [5], a
parallel library that implements a dynamic management of a collection of
adaptive K-trees on distributed memory supercomputers.

Lagranigan particle algorithms for pellet / SPI fragment ablation by hot
plasma electrons were developed in [19] and validated using experimental data
on the injection of hydrogen fueling pellets and DIII-D experiment with small
neon pellets [10]. Physics models for pellet ablation clouds in REs beams were
implemented in LP within the present work. In addition, the cross-field
transport model for ablated clouds has been significantly improved. Since [19]
dealt with relatively small geometric domains around ablating pellets,
sufficient for computing pellet ablation rates, it implemented simplified
grad-B drift and Alfven wave drag models responsible for establishing pellet
shielding lengths. The improved model described in Section 3.3 is necessary
for large-scale simulation of ablated clouds reported in this work.

##  4 Simulation results

###  4.1 The simulation of a single fragment into RE beam

For better understanding of main processes in the ablation cloud, we simulate
first a single $5mm$ diameter hydrogen pellet injected with initial velocity
of $500m/s$ into the runaway beam from the low-field-side, as shown in the
schematic in Section 2. At the start of the simulation, the pellet is located
at the RE beam edge since the pellet ablation in low temperature plasma is
negligibly small: the loss of pellet mass due to the ablation by thermal
electrons was estimated to be 2 %.

(a) Density (1/cc) at $t=40\mu s$

(b) Density (1/cc) at $t=80\mu s$

(c) Pressure (bar) at $t=40\mu s$

(d) Pressure (bar) at $t=80\mu s$

(e) Temperature (eV) at $t=40\mu s$

(f) Temperature (eV) at $t=80\mu s$

Figure 2: Distribution of ablation cloud states on 2D slices at $t=40\mu s$
(left images) and $t=80\mu s$ (right images). Horizontal axis is along
magnetic filed lines. There are two vertical scales on both sides of images:
penetration depth (right axis), global coordinates (left axis). RE beam region
is above 70 cm in vertical direction. The pellet is injected from the bottom
side of the plot.

As the pellet enters the RE beam ramp, it rapidly sublimates and forms a cold,
dense ablation cloud. The short-time-scale sublimation process is not
explicitly resolved - the pellet is converted instantaneously into gas state
at the RE beam edge in our numerical model. Figures 2(a)-(f) show states of
the ablation cloud in early stages of expansion, at 40$\mu s$ (left images)
and 80$\mu s$ (right images). The expansion at 40$\mu s$ remains almost
spherical, with some deformation at the top due to extra heating in the non-
uniform RE ramp. Ionization and MHD effects become visible at 80$\mu s$ as the
top part of the cloud expands along magnetic field lines and experiences the
grad-B drift.

(a) Connection coefficients

(b) Ratio of drift terms

Figure 3: (a) Evolution of connection coefficients $P_{con}$, $P_{Alf}$, and
$P_{same}$. (b) Ratio of magnitudes of the external connection damping term
and the Alfven wave drag term.

As the cloud heats-up and ionizes, the magnetic field diffusion time exceeds
the cloud advection time and the Alfven wave and external connection current
drag terms turn on by the Heaviside function in equation (10). Figure 3
explains the dynamics and relative magnitude of these terms. During the first
15 microseconds of the action of drag terms, the magnitudes of both connection
coefficients are close, $P_{con}\sim 0.3$ and $P_{Alf}\sim 0.37$, and the
Alfven wave drag dominates. The Alfven wave drag connection coefficient
reduces to small values within the next 100$\mu s$ and the external connection
current drag becomes the dominant factor. These terms rapidly change the
penetration dynamics of the ablation cloud as shown in Figure 4. At 0.75 ms
after injection, the whole cloud still moves with the transverse velocity
close to the injection velocity of 500 m/s (Figure 4(a)). This is explained by
a very low ionization level of the ablation cloud at early time resulting in a
weak grad-B drift force. Ionization of front layers in the ablation cloud
increases after 0.75 ms and the grad-B drift significantly slows them down.
Figure 4(b) shows that top layers of the cloud slow down by 5 times while the
bottom layers are still largely unaffected at 1 ms. At this time, the
Heaviside function turns on and the cloud starts experiencing the Alfven wave
and the external connection current drag forces. The penetration rapidly stops
within the next 50 $\mu s$ as we see in Figure 4(c): the top layers slowly
drift backward by the grad-B drift while the bottom layers very slowly
continue the forward motion. The penetration depth evolution is shown in
Figure 5(a). We would like to note that we turned on the Heaviside function at
exactly 1 ms for simplicity of analysis; this is within 4% from the exact
timing when $\tau_{Bdiff}=\tau_{adv}$. Due to a small diameter of the ablation
cloud, the Heaviside function turns on the drift forces uniformly across the
cloud. We employ a more realistic approximation for SPI simulations in the
next section.

(a) Drift velocity (m/s) at 0.75 ms

(b) Drift velocity (m/s) at 1 ms

(c) Drift velocity (m/s) at 1.05 ms

Figure 4: Distribution of the ablation cloud drift velocity at 0.75ms, 1ms and
1.05ms. The longitudinal axis (m) is along magnetic field lines. There are two
vertical axes: the left axis (m) is the global coordinate and the right axis
(cm) is the penetration depth. The pellet is injected from the bottom side of
the plot. The positive direction of drift velocity is the pellet injection
direction.

(a) Penetration depth

(b) Cloud length

Figure 5: (a) Evolution of the penetration depth of the ablation cloud into RE
beam. (b) Evolution of the ablation cloud length.

The longitudinal cloud expansion remains almost linear in time throughout
simulation (Figure 5(b)). Since the drag terms also prevent strong grad-B
drift, the cloud remains at approximately the same depth in the RE beam. At
about 3 ms, the cloud reconnects longitudinally at the safety factor $q=2$
(Figure 6).

Figure 6: Long-term dynamics of the ablation cloud. Velocity distribution
(m/s) is shown.

Since our simulations used the direct electron impact ionization in addition
to the thermal Saha ionization, it is important to clarify the role of the two
ionization processes. We do this in simplified geometry and physics
approximations to avoid the influence of other factors. In particular, we
artificially turn off the cross-field terms (10) and place the pellet in a
uniform RE beam. In 1MA RE current, the cloud expands spherically and no MHD
effects are observed if only thermal Saha ionization is used. At $0.2ms$, the
ionization level in the bulk of the cloud is around $10^{-7}$ (Figure 7(a)).
After adding the electron ionization by impact, we observed significant MHD
effects in the ablation cloud and the channeling of the ablated material along
magnetic field lines. The ionization level in the bulk of the cloud was around
$0.004$ (Figure 7(b)) and the cloud length reached $190cm$ after $0.2ms$. If
the RE current is increased to 10 MA, the Saha ionization alone is sufficient
to cause the longitudinal motion along magnetic field lines, but the addition
of impact ionization changes quantitative characteristics of the ablation
cloud dynamics. At 200$\mu s$, the cloud was 3.2 m long and 40 cm in diameter
with only Saha ionization and it expanded to 6.6 m length and reduced to 20 cm
diameter when the impact ionization was added. Therefore, the electron impact
ionization changes the cloud dynamics at all conditions, and it is especially
important at low values of the RE density in the ramp.

(a) Saha ionization

(b) Saha + impact ionization

Figure 7: Temperature distribution (eV) of the ablation cloud at 0.2ms with
(a) thermal Saha ionization and (b) both Saha ionization and ionization by
impact.

We would like to add a note on the heat diffusion in the ablation cloud. Our
simulations show that the longitudinal heat diffusion provides negligibly
small factor in the energy balance. Even for the most favorable case for the
longitudinal heat diffusion - 20 eV background plasma, only small regions at
the front and back of the ablation cloud heat up to the background temperature
at times relevant to our simulations; the bulk ablation plasma temperatures
are mostly defined by the RE energy deposition and the expansion dynamics
remain practically the same. Since the transverse heat conductivity is several
orders of magnitude lower than the parallel conductivity, the transverse heat
diffusion does not have any noticeable effect on the cloud dynamics despite
the difference in scales: the cloud radius is much smaller than the cloud
length. We verified our results by comparing 3D simulations to simplified 1D
simulations. Considering the small effect of the thermal diffusion, 3D
simulations reported in this work were performed with the heat diffusion
turned off to reduce computing time. This explains why the temperature of the
ablated plasma edge is not at the background plasma temperature.

###  4.2 Simulation of SPI into RE beam

In this section, we report numerical simulation results of SPI into the RE
beam in ITER. We assume uniform size of pellet fragments and compare
simulations using 555 fragments of 5 mm diameter and 5000 fragments of 2.4 mm
diameter. Each cloud of SPI fragments forms a cone with a 20 degree semi-
angle. The distribution of fragments is uniform with random perturbations. The
injection time for clouds of both fragment sizes is 2.2 ms.

As first fragments enter the RE beam ramp, they rapidly ablate and form dense,
high-pressure individual clouds (called here cloudlets). The initial stage of
the injection and ablation process (at 100 $\mu s$) is illustrated in Figure 8
which shows the ablated material velocity distribution in the direction of
injection. At this time, we observe small spherical cloudlets created from
most recent fragments as well as larger cloudlets from earlier fragments that
became sufficiently ionized and expanded along magnetic field lines by MHD
forces. Cloudlets from larger size fragments (Figure 8(a)) evolve individually
at this stage while cloudlets from smaller fragments (Figure 8(b)) merge into
a global cloud. In addition, higher initial pressure in cloudlets from 5 mm
diameter fragments push some fraction of cold and neutral ablated material
away from the RE region. Similar process is also present for smaller SPI
fragments but to a less extent, resulting in more material being deposited in
the RE region at earlier stages of the ablation process. At later time,
cloudlets created by large fragments also merge into a single albeit less
uniform ablation cloud.

(a) 555 fragments

(b) 5000 fragments

Figure 8: Initial phase of SPI: state of fragments after 100 $\mu s$. Velocity
distribution (m/s) in ablated clouds is shown in the injection direction (view
from the LFS) for (a) 555 fragments and (b) 5000 fragments.

The ablation clouds continue to heat-up by absorbing the energy of REs, expand
along magnetic fields lines, penetrate into the RE beam ramp by the $E\times
B$ drift, and partially drift in the direction of large tokamak radius by the
grad-B drift force. At some point in time, the front layers of the clouds in
the injection direction become sufficiently ionized that the magnetic field
diffusion time exceeds the advection time, causing the Alfven wave drag and
external connection currents forces. Because of a large transverse size of the
ablation clouds, we do not turn-on the drag forces uniformly in the entire
cloud, as we did in simulations of a single fragment, but instead check the
Heaviside function argument in equation (10) individually for small tubes of
the ablated material along magnetic field lines and turn-on the drag forces
accordingly. These forces decelerate the clouds and prevent them from
penetrating deeply into the RE region.

Figures 9 shows temperature distributions in both clouds at the time of their
maximum penetration into the RE beam. Plots show that the longitudinal
expansion (approximately 36 m) and temperature values are very similar for
both clouds obtained from large and small size fragments. Close to the pellet
fragment injection region, we observe cold clouds of neutral material,
approximately spherical in shape, expelled from the RE beam. The fraction of
expelled cold particles is about 15% and 10% for the 555 and 5000 fragment
SPI, respectively. At this time, a significant fraction of hot ablated
material has also drifted outside of the RE region: the corresponding
fractions for the 555 and 5000 fragment SPI are 50% and 65%. Depositing more
ablated material into the RE beam at the initial stages of ablation in the
5000 fragment SPI leads to higher pressure values in the ablation cloud and
stronger grad-B drift (Figure 10(b)). As a result, less material remains in
the RE beam over longer time.

(a) 555 fragments

(b) 5000 fragments

Figure 9: State of SPI ablation cloud near maximum penetration depth, at time
1 ms. Temperature distribution (eV) in ablated clouds is shown in the
injection direction (view from the HFS) for (a) 555 fragments and (b) 5000
fragments.

(a) 555 fragments

(b) 5000 fragments

Figure 10: Transverse velocity distribution (m/s) in ablated clouds for (a)
555 fragments and (b) 5000 fragments near maximum penetration depth at time 1
ms. The line of view is along magnetic field lines and the vertical axis is
along the injection direction, giving the penetration depth.

The global dynamics of ablation clouds is shown in Figure 11. At approximately
1.7 ms, the clouds re-connect in the longitudinal direction assuming the
tokamak safety factor of 2. Both clouds contain wide layers of ablated
hydrogen drifting beyond the RE region. As new fragments are still injected
into the RE beam, the fraction of material drifted outside of the RE region
slowly increases: it is a few percent higher at the re-connection time
compared to the reported values at 1.1 ms.

The penetration depth of the leading edge of both clouds in the RE beam is
similar; each case is plotted with different color dots in Figure 12 and the
line corresponds to a polynomial fit using both sets of data. The penetration
depth is slightly smaller compared to that of the single cloud fragment
reported in the previous section.

(a) 555 fragments

(b) 5000 fragments

Figure 11: State of SPI ablation cloud after longitudinal reconnection at
$q=2$ at time 1.7 ms. Velocity distribution (m/s) in ablated clouds is shown
for (a) 555 fragments and (b) 5000 fragments. Figure 12: Evolution of
penetration depth of SPI ablated cloud into the runaway beam.

##  5 Conclusions and future work

Lagrangian particle simulation of the evolution of ablation clouds of hydrogen
pellets and SPI fragments heated by a runaway electron beam in ITER have been
performed. To accomplish this task, capabilities of the Lagrangian particle
pellet / SPI code [19] have been extended by including physics models for the
injection of hydrogen pellets and SPI into a runaway electron beam. In
particular, ionization of ablated materials by direct electron impact and the
volumetric heating by REs have been implemented. The transverse flow model has
been significantly improved by including the internal (poloidal) and external
Pégourié connection currents and the corresponding numerical algorithms have
been developed and implemented.

Simulations of the ablation of a single pellet fragment in the RE region in
ITER demonstrated that the thermal Saha ionization alone is not sufficient for
ionizing the cloud at the edge of the RE beam. Within the RE density ramp, the
pellet cloud evolution strongly depends on the impact ionization by REs which
is crucial for causing strong MHD effects. Both 3D and 1D simulations showed
that the longitudinal heat diffusion provides negligibly small factor to the
total energy balance. Single fragment simulations quantified various terms in
the transverse dynamics model. The penetration depth of the ablated cloud into
the RE region was about 55 cm, reached at 1.0 ms. At about 3 ms , the cloud
reconnected longitudinally at q=2.

SPI simulations show details of the evolution of ablated clouds, in particular
the expansion along magnetic field lines, cross-field transport due to grad-B
drift, Alfven wave drag and connection current damping terms, global dynamics
of the ablated cloud at tokamak scales as well as the maximum penetration of
the leading edge of ablated clouds into the RE region. We observe that a
fraction of the ablated hydrogen gas (10 - 15%) is expelled from the RE beam
by pressure gradients soon after each fragment injection. The expelled cold
gas slowly expands in the background plasma. In the RE beam, the ablated
material becomes ionized by both the thermal Saha process and by direct RE
impact, absorbs energy of RE’s, rapidly propagates along magnetic field lines
and experiences various cross-field drift effects. The SPI cloud penetrates
approximately 40 cm into the RE beam for the given RE ramp profile. The
deepest penetration was achieved at 1.1 - 1.2ms. For $q=2$, the cloud
reconnects longitudinally at 1.7ms. At these stages of the cloud evolution,
more than 50% of the ablated gas is moved outside of the RE beam by the grad-B
drift.

We also studied the dependence of ablation cloud evolution processes on the
number of fragments in the SPI plume. More material was initially deposited
into the RE beam by the SPI of 5000 fragments compared to 555 fragments: for
the case of smaller number of large fragments, more neutral gas was expelled
from the RE region at the initial stage of ablation. But depositing more
material into the RE beam led to higher ablation cloud pressure values and
stronger grad-B drift. As a result, a larger fraction of hot ablated material
drifted outside the RE beam in the case of 5000 fragment SPI. Despite these
differences, the global dynamics and the cloud reconnection time were similar.
In particular, both clouds achieved very similar penetration depth at any
time.

In the future, we will improve physics models for REs. In particular, we will
add an evolution model for the RE beam that accounts for REs decay in the
ablated material.

Acknowledgments This research has been partially supported by the U.S.
Department of Energy grants SciDAC Center for Tokamak Transient Simulations
and DE-FG02-07ER54917. The authors are grateful to N. Garland for providing
electron impact ionization rates. Permission to use PrismSPECT from I.
Golovkin is gratefully acknowledged.

Data availability The data that support the findings of this study are
available from the corresponding author upon reasonable request.

## References

  * [1] J. E. Bailey, G. A. Rochau, R. C. Mancini, C. A. Iglesias, J. J. MacFarlane, I. E. Golovkin, C. Blancard, Ph. Cosse, and G. Faussurier.  Experimental investigation of opacity models for stellar interior, inertial fusion, and high energy density plasmas.  Physics of Plasmas, 16(5):058101, 2009.
  * [2] L.R. Baylor, T.C. Jernigan, P. Gohil, G.L. Schmidt, K. H. Burrell, S.K. Combs, D.R. Ernst, C.M. Greenfield, R.J. Groebnerand, W.A. Houlberg, et al.  Improved fueling and transport barrier formation with pellet injection from different locations on DIII-D.  in Proc. 18th Fusion Energy Conference, IAEA-CN-77, 01 2000.
  * [3] Nicolas Bosviel, Paul Parks, and Roman Samulyak.  Near-field models and simulations of pellet ablation in tokamaks.  Physics of Plasmas, 28(1):012506, 2021.
  * [4] B.N. Breizman, P. Aleynikov, E. M. Hollmann, and M. Lehnen.  Physics of runaway electrons in tokamaks.  Nuclear Fusion, 59(8):083001, 2019.
  * [5] Carsten Burstedde, Lucas C. Wilcox, and Omar Ghattas.  p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees.  SIAM Journal on Scientific Computing, 33(3):1103–1133, 2011.
  * [6] N. Commaux, B. Pégourié, L.R. Baylor, F. Kochl, P.B. Parks, T.C. Jernigan, A. Géraud, and H. Nehme.  Influence of the low order rational $q$ surfaces on the pellet deposition profile.  Nucl. Fusion, 50:025011, 2010.
  * [7] Nathan Garland.  Private communication.
  * [8] E. M. Hollmann, A. Bortolon, F. Effenberg, N. Eidietis, D. Shiraki, I. Bykov, B. E. Chapman, J. Chen, S. Haskey, J. Herfindal, A. Lvovskiy, C. Marini, A. McLean, T. O’Gorman, M. D. Pandya, et al.  Dynamic measurement of impurity ion transport in runaway electron plateaus in DIII-D.  Physics of Plasmas, 29(2):022503, 2022.
  * [9] E.M. Hollmann, I. Bykov, N.W. Eidietis, J. L. Herfindal, A. Lvovskiy, R.A. Moyer, P.B. Parks, C. Paz-Soldan, , A. Yu. Pigarov, D.L. Rudakov, et al.  Study of argon expulsion from the post-disruption runaway electron plateau following low-Z massive gas injection in DIII-D.  Phys. Plasmas, 27:042515, 2020.
  * [10] E.M. Hollmann, N. Naitlho, S. Yuan, R. Samulyak, P.B. Parks, D. Shiraki, J. Herfindal, and C. Marini.  Measurement and simulation of small cryogenic neon pellet Ne-I 640 nm photon efficiency during ablation in DIII-D plasma.  Phys. Plasmas, 2022.  Submitted.
  * [11] D.I. Kiramov and B. N. Breizman.  Pellet sublimation and expansion under runaway electron flux.  Nuclear Fusion, 60(8):084004, 2020.
  * [12] F. Koechl, B. Pégourié, A Matsuyama, H. Nehme, D. Frigione, L. Garzotti, G. Kamelander, and V. Parail.  Modelling of pellet particle ablation and deposition: The hydrogen pellet injection code HPI2.  EFDA–JET–PR(12)57:1–86, 2012.
  * [13] A. Matsuyama, F. Koechl, B. Pégourié, R. Sakamoto, G. Motojima, and H. Yamada.  Modeling of the pellet depostion profile and $\nabla$B-induced drift displacament in non-axisymmetric configurations.  Nucl. Fusion, 52:123017, 2012.
  * [14] P.B. Parks and L.R. Baylor.  Effect of parallel flows and toroidicity on cross-field transport of pellet ablation matter in tokamak plasmas.  Physical Review Letters, 94:125002, 2005.
  * [15] C. Paz-Soldan, C. Reux, K. Aleynikova, P. Aleynikov, V. Bandaru, M. Beidler, N. Eidietis, Y.Q. Liu, C. Liu, A. Lvovskiy, et al.  A novel path to runaway electron mitigation via deuterium injection and current-driven MHD instability.  Nuclear Fusion, 61(11):116058, oct 2021.
  * [16] Cédric Reux, Carlos Paz-Soldan, Pavel Aleynikov, Vinodh Bandaru, Ondrej Ficker, Scott Silburn, Matthias Hoelzl, Stefan Jachmich, Nicholas Eidietis, Michael Lehnen, et al.  Demonstration of safe termination of megaampere relativistic electron beams in tokamaks.  Phys. Rev. Lett., 126:175001, Apr 2021.
  * [17] R. Samulyak, T. Lu, and P.B. Parks.  A magnetohydrodynamic simulation of pellet ablation in the electrostatic approximation.  Nuclear Fusion, 47(2):103–118, Jan 2007.
  * [18] R. Samulyak, W. Xingyu, and H.-C. Chen.  Lagrangian particle method for compressible fluid dynamics.  Journal of Computational Physics, 362:1–19, June 2018.
  * [19] R. Samulyak, S. Yuan, N. Naitlho, and P. B. Parks.  Lagrangian particle model for 3d simulation of pellets and spi fragments in tokamaks.  Nuclear Fusion, 61(4):046007, 2021.
  * [20] D. Shiraki, N. Commaux, L.R. Baylor, C.M. Cooper, N.W. Eidietis, E.M. Hollmann, C. Paz-Soldan, S.K. Combs, and S.J. Meitner.  Dissipation of post-disruption runaway electron plateaus by shattered pellet injection in DIII-d.  Nuclear Fusion, 58(5):056006, mar 2018.
  * [21] Y.B. Zel’dovich and Y.P. Raiser.  Physics of shock waves and high temperature hydrodynamic phenomena.  Dover, 2002.


# Albedo heterogeneity on the surface of (1943) Anteros

Joseph Masiero11affiliation: Institute for Astronomy, University of Hawaii,
2680 Woodlawn Dr, Honolulu, HI 96822 22affiliation: Jet Propulsion Lab,
California Institute of Technology, 4800 Oak Grove Dr., MS 264-767, Pasadena,
CA 91106<EMAIL_ADDRESS>
###### Abstract

We have investigated the effect of rotation on the polarization of scattered
light for the near-Earth asteroid (1943) Anteros using the Dual Beam Imaging
Polarimeter on the University of Hawaii’s $2.2~{}$m telescope. Anteros is an
L-type asteroid that has not been previously observed polarimetrically. We
find weak but significant variations in the polarization of Anteros as a
function of rotation, indicating albedo changes across the surface.
Specifically, we find that Anteros has a background albedo of $p_{v}=0.18\pm
0.02$ with a dark spot of $p_{v}<0.09$ covering $<2\%$ of the surface.

##  1 Introduction

As the last remnants from an epoch of accretive formation, asteroids provide
us windows into the composition and history of the inner Solar System. Except
for the few largest bodies, asteroids did not heat up enough via decay of
short-lived radionuclides or dissipation of gravitational potential energy to
undergo complete differentiation. As such the minerals observed on their
surfaces capture the elemental and temperature history of the local region of
the protoplanetary disk at the time of their formation. By understanding
asteroid surfaces we can directly probe those early disk conditions.

As the illuminated cross section of an asteroid changes the observed
brightness fluctuates. Given a large enough sample of data a full shape model
of a rotating body can be constructed (Kaasalainen, et al., 2001) even though
it is unresolved. Photometric surveys for asteroid light curves have set
limits on the composition and density of asteroids as a population (Pravec, et
al., 2002) and have estimated the average shape distribution of small Main
Belt asteroids (Masiero, et al., 2009a). All of these results however assume
that the light curve is dominated by the object’s shape and that the entire
surface has a uniform composition and albedo.

It is possible to test for albedo variations using optical imaging
polarimetry, color variations, or even simply photometric variations under the
assumption of a regular shape (Akimov, et al., 1983). In the case of
polarimetry there are strong empirical correlations between the albedo of an
asteroid and both the slope of the polarization-phase curve and the location
of the minimum (negative) polarization (most recently: Cellino, et al., 1999).
The polarization of light scattered off of an atmosphereless body as a
function of phase angle depends on the distance between scattering elements
and their index of refraction (Muinonen, 1989; Muinonen, et al., 2002c). Index
of refraction is an inherent mineralogical property and recent work has shown
that the inter-element scattering distance is likewise determined by the
surface chemistry (Masiero, et al., 2009b). It is not unexpected then that
asteroids of different spectral classification show different polarization-
phase curves (Muinonen, et al., 2002b), or that a differentiated-then-broken
object like (4) Vesta would show polarization variations with rotation.

In almost every way investigated so far the asteroid Vesta stands out as an
interesting and unique object, and this is similarly the case for albedo
variation studies. Although a handful of other asteroids have weak detections
of rotational modulation of their polarization and thus albedos (e.g. Broglia
& Manara, 1992, etc.), Vesta represents one clear case of an object with
polarization changes across its surface caused by changes in composition
(Degewij, et al., 1979; Lupishko, et al., 1988; Broglia & Manara, 1989). This,
along with photometry, spectroscopy and adaptive optic imaging, has lead to
the current interpretation of Vesta as a differentiated body that has
undergone a nearly-catastrophic impact event leaving a giant crater in its
southern hemisphere. The crater reveals a now-solid mantle distinctly
different in color and composition from the original crust material (Cellino,
et al., 1987; Thomas, et al., 1997).

Nakayama, et al. (2000) found rotational modulation of the polarization for
the asteroid (9) Metis with amplitude similar to what is observed for Vesta.
Metis is a $D\sim 180~{}$km asteroid that may have two large spots of
significantly higher albedo ($p_{v}\sim 0.24-0.28$) than the background
material ($p_{v}=0.14$). The authors find that both bright areas are on the
leading (for prograde rotation) or trailing (for retrograde rotation) faces of
the model ellipsoid (Mitchell, et al., 1995). The cause of albedo
heterogeneity across the surface of objects smaller than Vesta is still
undetermined. It is possible that non-disruptive collisions with impactors of
different composition can leave localized deviations from the average
mineralogy. Alternatively, a late formation with a history free of melting may
preserve the varied composition of the protoplanetary disk. However this
theory is complicated by recent work showing that asteroids likely were born
big, and most objects smaller than a few hundred kilometers in diameter should
be collisionally-created fragments (Morbidelli, et al., 2009).

Identifying albedo variations for small asteroids allows us to evaluate the
accuracy of the assumption that flux changes are solely dependent on shape.
This has important implications for results based on this assumption,
especially shape models. Additionally we can also quantify the effect of
collisions between small bodies in determining an asteroid’s local regolith
properties.

All asteroid polarization-phase curves follow the same general trend with
increasing phase angle: zero polarization at zero phase, becoming negative to
some minimum value and then increasing in an approximately linear fashion.
Note that as is standard for Solar system polarimetry the reference direction
for the angle of polarization is aligned with the vector perpendicular to the
plane of scattering such that “positive” and “negative” polarization are
defined as perpendicular and parallel to the scattering plane, respectively.
The results presented here follow this convention. Each polarization-phase
curve displays three distinguishing values used to classify its properties:
the minimum negative polarization ($P_{min}$), the phase angle at which the
polarization returns to zero (the inversion angle, $\alpha_{0}$) and the
linear slope of the polarization-phase relation beyond the inversion angle
($h$). Making use of the albedo-polarization relation from Cellino, et al.
(1999),

$\displaystyle\log p_{v}$ $\displaystyle=$ $\displaystyle(-1.12\pm 0.07)\log
h-(1.78\pm 0.06)$ (1)

(where $p_{v}$ is the geometric V-band albedo) we can use imaging polarimetry
to test for changes in polarization that directly indicate albedo
heterogeneity across an asteroid’s surface.

##  2 Observations and Discussion

Changes in the polarization of the scattered light across the surface of an
asteroid will be small even in the best-case scenarios. To obtain a
significant measurement of the largest of these variations we require an
instrument that can attain better than $0.1\%$ polarization accuracy. Our
study made use of the Dual Beam Imaging Polarimeter (DBIP) located on the
University of Hawaii’s $2.2~{}$m telescope on Mauna Kea, Hawaii (Masiero, et
al., 2007). DBIP uses a double-calcite Savart plate in series with a quarter-
wave and a half-wave retarder to simultaneously measure linear and circular
polarization with accuracy better than $0.1\%$ (Masiero, et al., 2008). DBIP
uses a $g^{\prime}+r^{\prime}$ filter with a bandpass of $400-700~{}$nm. While
asteroid polarization does depend on color (Cellino et al., 2005) changes are
usually small in this wavelength range and typically within measurement
errors.

Observation of our target asteroid were supplemented with polarized and
unpolarized standards to confirm consistency of setup, stability of the
instrument, and accuracy of the measurements. Standards were taken from
Fossati, et al. (2007) as well as the standard list for
Keck/LRISp111http://www2.keck.hawaii.edu/inst/lris/polarimeter/polarimeter.html
which includes the Hubble Space Telescope polarimetric standards (Schmidt, et
al., 1992). These measurements all verified that the errors were within the
range expected from previous calibrations.

As albedo is related to the polarization-phase slope $h$, for a given albedo
variation the respective polarization change will be larger when observed at
higher phase angles (for $\alpha>\alpha_{0}$). Geometric restrictions prevent
Main Belt asteroids (MBAs) from ever reaching phases angles larger than
$\alpha\sim 30^{\circ}$, but near-Earth asteroids (NEAs) pass closer to Earth
and so can reach much larger phase angles. For this reason, NEAs are preferred
targets when looking for albedo variations. At high phases the polarization of
scattered light takes on a linear trend that increases up to the level of
$\sim 5-10\%$ polarized depending on surface mineralogy. These large
polarizations mean that any variation with rotation at high phases can be
easily interpreted as changes in the integrated surface albedo using the
slope-albedo relation (Eq 1) and that the absolute value for the range of the
albedo can be determined.

From July to September of 2009 the NEA (1943) Anteros passed through phase
angles of $16-40^{\circ}$ all while brighter than $V=17~{}$mag presenting a
prime opportunity to measure the polarization, slope and albedo with high
accuracy. The optical/NIR spectrum of Anteros displays a spectral slope
comparable to S-types but with a muted $1~{}\mu$m absorption band resulting in
a classification of L-type (Binzel, et al., 2004). With a measured period of
$P=2.8695\pm 0.0002~{}$hr and a single-peaked photometric light curve with
amplitude $A=0.09~{}$mag (Pravec, et al., 1998), Anteros is an excellent
target to test for rotational variation in polarization and albedo in a few
nights of observing. In particular, a single-peaked low-amplitude light curve
indicates a shape very close to spherical. (Pravec, et al. (1998) found an
amplitude of $0.09~{}$mag across phase angles ranging from
$19^{\circ}<\alpha<32^{\circ}$, resulting in a shape approximation of
$a/b<1.1$).

In Table 1 we present our polarimetric observations of Anteros. Included for
each night are the V magnitude, exposure time, number of 6-exposure
polarimetry measurements acquired ($n_{meas}$), phase angle, ecliptic
longitude, summed linear polarization of all measurements and linear
polarization angle with respect to the vector perpendicular to the scattering
plane ($\theta_{p}$). No significant circular polarization was detected on any
of the nights. The average nightly polarizations are shown in Fig 1 along with
generic model polarization-phase curves for typical S-type asteroids (dotted)
and C-type asteroids (dashed). The model curves were made using the linear-
exponential modeling technique presented by Muinonen, et al. (2002a) and
fitted by-eye to the data shown in Fig 1 of Muinonen, et al. (2002b), to act
as useful approximations. The constants used in this case were, for the C-type
model: $P_{a}=5.5$, $P_{d}=6$, $P_{k}=0.3$ and for the S-type model:
$P_{a}=4.3$, $P_{d}=12$, $P_{k}=0.17$.

The polarization of Anteros is clearly most closely related to an S-type
polarization curve as expected from its spectral features. We measure for
Anteros an inversion angle of $\alpha_{0}=20.3\pm 0.3^{\circ}$ and a slope
beyond the inversion angle of $h=0.122\pm 0.001$. Both $\alpha_{0}$ and $h$
(and their respective errors) were found by conducting a least-squares
minimization fit of a line to the four nights of data. As the data span a
large range of phases and have small individual errors, the resultant error on
$h$ is small. Following Eq 1 we derive a bulk albedo of $p_{v}=0.175\pm
0.002\pm 0.02$ (errors relative and absolute, respectively). Note that the
limiting error on albedo ($\sigma_{abs}=0.02$) derives from the uncertainty on
the constants in Eq 1 and will affect absolute albedo measurements. This does
not apply to relative comparisons between measured albedos, which have an
error of $\sigma_{rel}=0.002$ in the above case. For all observations reported
here, the calibration error on Eq 1 is greater than the noise error by nearly
an order of magnitude, and thus dominates the final error on the measured
albedos. (All errors reported in this paper are $1~{}\sigma$.)

Following (Muinonen, et al., 2002c) polarization can be approximated as

$P\sim\frac{\alpha^{2}}{2n}-\left(\frac{n-1}{n+1}\right)^{2}\frac{(k~{}d~{}\alpha)^{2}}{2~{}[1+(k~{}d~{}\alpha)^{2}]}$

and solving this for the case of zero polarization at the inversion angle we
find

$\displaystyle\alpha_{0}\sim\sqrt{n\left(\frac{n-1}{n+1}\right)^{2}-\frac{1}{(kd)^{2}}}$
(2)

where $k=2\pi/\lambda$, $d$ is the inter-element scattering distance, and $n$
the index of refraction. For a central wavelength of $0.55~{}\mu$m for DBIP
and a typical scattering distance for NEAs of $d\sim 4~{}\mu$m (Masiero, et
al., 2009b) the second term on the right in Eq 2 is negligible, and the index
of refraction can be determined for a given inversion angle. For Anteros we
find an index of refraction of $n\sim 1.74$.

In Fig 2 we show the nightly polarization light curves for Anteros.
Observations have been wrapped to the measured $2.8695~{}$hr photometric
period and all measurements within a $0.1~{}$phase bin have been co-added to
reduce measurement error. The zero point for rotation phase was chosen
arbitrarily on the first night. When all nights of data are wrapped to this
zero-point the error on the period translates to a phase error of $\pm
0.01~{}$rotations between each observing night. Thus, features at specific
phases can be compared across nights.

We find weak variation in the polarization at a $4\sigma$ significance level
for the night of 2009-07-22 with an amplitude of $0.3\%$ when comparing the
data in the rotation phase range of $0.6-0.8$ to the data between rotation
phases $0.8-0.4$. These polarization changes most likely indicate a variation
in the albedo across the surface of Anteros. The amplitude of the polarization
variation scales with the absolute polarization so it is not surprising that
the other observing nights at lower phases show no clear variation, e.g. a
variation with amplitude of $P=0.3\%$ at a phase of $\alpha=40.5^{\circ}$
would be expected to have an amplitude of $P=0.1\%$ at a phase of
$\alpha=28.0^{\circ}$ which is below our threshold for significant detection.
Additionally, changes in the observing geometry between observations could
cause an area that was observable on the night of 2009-07-22 to have a reduced
visibility on subsequent observations, or even be beyond the horizon. However
even in the most extreme case, where the rotation axis is parallel the to
plane of Earth’s orbit, the line-of-sight vector only moves a total of $\sim
18^{\circ}$ with respect to the rotation axis over the dates of the
observations (this is equivalent to the change in ecliptic longitude). Though
this could account for the changes in polarization under specific
circumstances, it is not the most likely scenario.

It has been suggested that the change in polarization alternatively could be
due to an extreme topographical feature that deviates significantly from the
surrounding area. At the high phases at which we observed Anteros first-order
scattering is dominant and so we may simply apply basic scattering properties
to the surface (e.g. the angle of incidence and angle of scattering are equal,
etc). Thus the light we observe necessarily must have been scattered by planes
on the surface normal to the vector that bisects the phase angle. Even on
unusual surfaces, microroughness will provide the appropriate scattering
facets. In the extreme argument, should the plane be a perfectly flat surface
(e.g. a mirror) it will scatter no light to the observer when away from a
perfect alignment. This will decrease the overall flux but not change the
percent of the received flux polarized by the surrounding area. If aligned
perfectly, the plane will still behave polarimetrically as the underlying
material it is made of, polarizing the same fraction of the light as
determined by its albedo. Thus even in case of extreme topography, percent
polarization measurements will only be sensitive to the underlying material
composition.

Changes in the absolute level of polarization as the phase angle changes
nightly precludes wrapping all four nights of observing onto a single rotation
phase. However we can correlate features at similar rotation phases across
nights and from this fit different polarization slopes to different locations
on the asteroid. We use the phase range of $0.6-0.8$ to represent the peak
polarization and the phase range of $0.8-0.4$ (wrapped) to represent the
baseline background polarization (as determined from the 2009-07-22 observing
night). We interpret the phase range of $0.4-0.6$ as a transition region (see
below) and do not include it in either measurement. Using the summed maximum
and baseline values for the first night as location benchmarks we determine
the absolute change in albedo across the surface of Anteros.

We find a background surface albedo of $p_{v}=0.181\pm 0.002\pm 0.02$ with a
single spot of much lower albedo. Note that this value does not vary
significantly from the albedo of $p_{v}=0.13\pm 0.03$ found from radiometric
modeling (Veeder, et al., 1981) or the one published in the compilation by
Chapman, et al. (1994) of $0.17$ (no error given). We measure an upper limit
to the albedo for the dark area of $p_{v}=0.160\pm 0.004\pm 0.02$ however this
feature is unresolved and thus the albedo measurement assumes coverage of the
maximal possible area allowable by the observing geometry ($44\%$ of the total
surface area for 2009-07-22, corresponding to a projected area of
$3.7~{}$km2). It is likely that the dark spot covers a much smaller area with
a much lower reflectance.

Reflected polarization from a mottled surface will mix when unresolved to give
a value between the two extremes. For the case of Anteros on the night of
2009-07-22, taking a background polarization for $\alpha=40.5^{\circ}$ of
$P_{bkg}=2.46\%$ (the mean of the baseline polarization value for the night)
and a peak measured polarization of $2.70\%$ we find that the true percent
polarization of the dark spot ($P_{dark}$) on that night can be described as

$P_{dark}*C+P_{bkg}(1-C)=2.7$

where $C$ is the fraction of the projected illuminated surface that the dark
region covers. This is simply because the polarized light from the dark area
is diluted by the signal from any background material also visible. This
relation can be simplified to:

$\displaystyle P_{dark}=\frac{0.24}{C}+2.46$ (3)

The 2009-07-22 data plotted in Fig 2 show a gradual build up in the
polarization value with a rapid falloff from the peak level to the base level
over slightly more than one tenth of a rotation. This is likely due to the
dark region rising over the horizon as seen from Earth and then passing from
the lit side quickly across the terminator and thus out of illumination. This
scenario would require Anteros to have a prograde rotation state. From the
rapid falloff we can calculate that the maximum size of the dark feature along
the direction of rotation is $<0.7~{}$km assuming it is located on the equator
of the asteroid. If the spot is not on the equator this argument would derive
a smaller value for the size.

Using this size as the diameter of a circular crater this corresponds to a
total projected surface coverage of the dark spot of $C<11\%$ (equivalent to
$<2\%$ of the total unprojected surface area, assuming the rotation pole is
parallel to the plane of the sky), giving a polarization value for this area
of $P_{dark}=4.7\%$. If instead the rotation axis is inclined to the plane of
the sky, the projected area would be smaller. This would increase the required
value of $P_{dark}$ which would cause the calculated value for the albedo (see
below) to decrease. Thus the assumptions made here represent the ’brightest-
case’ scenario for the spot.

As the polarization for the other three nights is consistent at all phases we
can only set an upper limit for the value of $P_{dark}$ on those nights based
on a maximum unresolved change in percent polarization of $0.1\%$. Using the
analouges of Eq 3 for the other nights we find limits of $P_{dark}<1.85\%$ for
the night of 2009-08-11, $P_{dark}<0.88\%$ for 2009-08-26, and
$P_{dark}<-0.14\%$ for 2009-09-10. We show the polarization for the dark
region separated out from the background material in Fig 3 as well as the same
S- and C-type generic models from Fig 1. Using these values we calculate a
lower limit on the slope and thus an upper limit on albedo. The albedo of the
dark spot is limited to $p_{v}<0.09$.

##  3 Conclusions

Using DBIP on the University of Hawaii’s $2.2$ m telescope we have
investigated the NEA (1943) Anteros for surface heterogeneity. We find that
Anteros shows significant polarimetric variation as a function of rotation at
high phase angles, implying an albedo gradient and corresponding surface
composition variations. We determine that Anteros has a base albedo of
$p_{v}=0.181\pm 0.002\pm 0.02$ (errors are relative and absolute,
respectively) consistent with literature values as well as a dark spot of
albedo $p_{v}<0.09$ covering $<2\%$ of its surface.

A single small asteroid showing albedo variations does not invalidate the
assumption that shape dominates the light curves of these bodies. Indeed it is
clear from spacecraft visits that shape does play an important role in
determining the reflected flux. Additionally, albedo variations for most
asteroids appear to be small and localized. However the potential for albedo
variations for small asteroids can not be discounted outright.

It is currently unclear what processes could cause localized albedo changes
across the surfaces of small asteroids. Further polarimetric studies of small
NEAs are thus necessary to constrain the frequency of albedo and composition
changes across the surfaces of these bodies. Once an account the population is
established evolutionary pathways to the creation of these features can be
explored.

## Acknowledgments

J.M. was supported under NASA PAST grant NNG06GI46G. The author would like to
thank Rob Jedicke and Alan Tokunaga for providing comments on the manuscript,
as well as V. Rosenbush and an anonymous referee for helpful reviews that
improved the paper. The author wishes to recognize and acknowledge the very
significant cultural role and reverence that the summit on Mauna Kea has
always had within the indigenous Hawaiian community. I am most fortunate to
have the opportunity to conduct observations from this sacred mountain.

## References

  * Akimov, et al. (1983) Akimov, L.A., Lupishko, D.F. & Belskaya, I.N., 1983, “Photometric heterogeneity of asteroid surfaces”, Soviet Astronomy, 27, 577.
  * Binzel, et al. (2004) Binzel, R.P., Birlan, M., Bus, S.J., Harris, A.W., Rivkin, A.S. & Fornasier, S., 2004, “Spectral observations for near-Earth objects including potential target 4660 Nereus : Results from Meudon remote observations at the NASA Infrared Telescope Facility (IRTF)”, P&SS, 52, 291.
  * Broglia & Manara (1992) Broglia, P. & Manara, A., 1992, “A study of the polarimetric lightcurve of the asteroid 16 Psyche”, A&A, 257, 770.
  * Broglia & Manara (1989) Broglia, P. & Manara, A., 1989, “Rotational variations in the optical polarization of 4 Vesta”, A&A, 214, 389.
  * Cellino et al. (2005) Cellino, A., Gil Hutton, R., di Martino, M., Bendjoya, Ph., Belskaya, I.N. & Tedesco, E.F., 2005, “Asteroid polarimetric observations using the Torino UBVRI photopolarimeter”, Icarus, 179, 304.
  * Cellino, et al. (1999) Cellino, A., Gil Hutton, R., Tedesco, E.F., di Martino, M. & Brunini, A., 1999, “Polarimetric Observations of Small Asteroids: Preliminary Results”, Icarus, 138, 129.
  * Cellino, et al. (1987) Cellino, A., Zappala, V., di Martino, M., Farinella, P. & Paolicchi, P., 1987, “Flattening, pole, and albedo features of 4 Vesta from photometric data”, Icarus, 70, 546.
  * Chapman, et al. (1994) Chapman, C.R., Harris, A.W. & Binzel, R., 1994, “Physical properties of near-Earth asteroids: implications for the hazard issue”, Hazards Due to Comets and Asteroids, ed. Gehrels (Univ of Arizona Press), 537.
  * Degewij, et al. (1979) Degewij, J., Tedesco, E.F. & Zellner, B., 1979, “Albedo and color contrasts on asteroid surfaces”, Icarus, 40, 364.
  * Fossati, et al. (2007) Fossati, L., Bagnulo, S., Mason, E. & Landi Del’Innocenti, E., 2007, “Standard Stars for Linear Polarization Observed with FORS1”, ASP Conf., 364, 503.
  * Kaasalainen, et al. (2001) Kaasalainen, M., Torppa, J. & Muinonen, K., 2001, “Optimization Methods for Asteroid Lightcurve Inversion. II. The Complete Inverse Problem”, Icarus, 153, 37.
  * Lupishko, et al. (1988) Lupishko, D.F., Bel’Skaia, I.N., Kvaratskheliia, O.I., Kiselev, N.N. & Morozhenko, A.V, 1988, “Polarimetry of Vesta in 1986 opposition”, Astron.Vestnik, 22, 142.
  * Masiero, et al. (2007) Masiero, J., Hodapp, K.-W., Harrington, D. & Lin, H., 2007, “Commissioning of the Dual-Beam Imaging Polarimeter for the UH 88-inch telescope”, PASP, 119, 1126.
  * Masiero, et al. (2008) Masiero, J., Hodapp, K.-W., Harrington, D. & Lin, H., 2008, “Extended Commissioning and Calibration of the Dual-Beam Imaging Polarimeter”, to appear in ASP Conf. Series for Astronomical Polarimetry 2008; arXiv:0809.4313.
  * Masiero, et al. (2009a) Masiero, J., Jedicke, R., Ďurech, J., Gwyn, S., Denneau, L. & Larsen, J., 2009a, “The Thousand Asteroid Light Curve Survey”, Icarus, 204, 145.
  * Masiero, et al. (2009b) Masiero, J., Hartzell, C. & Scheeres, D.J., 2009b, “The effect of the dust size distribution on asteroid polarization”, AJ, 138, 1557.
  * Morbidelli, et al. (2009) Morbidelli, A., Bottke, W., Nesvorný, D. & Levison, H.F., 2009, “Asteroids were born big”, Icarus, in press.
  * Mitchell, et al. (1995) Mitchell, D.L., Ostro, S.J., Rosema, K.D., Hudson, R.S., Campbell, D.B., Chandler, J.F. & Shapiro, I.I., 1995, “Radar observations of asteroids 7 Iris, 9 Metis, 12 Victoria, 216 Kleopatra, and 654 Zelinda”, Icarus, 118, 105.
  * Muinonen (1989) Muinonen, K., 1989, “Electromagnetic Scattering by Two Interacting Dipoles”, Proc. URSI International Symp. on Electromagnetic Theory, 428.
  * Muinonen, et al. (2002a) Muinonen, K., Piironen, J., Kaasalainen, S. & Cellino, A., 2002a, “Asteroid photometric and polarimetric phase curves: empirical modeling”, Mem. S.A.It., 2002, 73, 716.
  * Muinonen, et al. (2002b) Muinonen, K., Piironen, J., Shkuratov, Y., Ovcharenko, A. & Clark, B., 2002b, “Asteroid Photometric and Polarimetric Phase Effects”, Asteroids III, ed. Bottke, Cellino, Paolicchi & Binzel (Univ of Arizona Press), 123.
  * Muinonen, et al. (2002c) Muinonen, K., Videen, G., Zubko, E. & Shkuratov, Yu., 2002b, “Numerical Techniques for Backscattering by Random Media”, Optics of Cosmic Dust, ed. Videen & Kocifaj (NATO Science Series), 79, 261.
  * Nakayama, et al. (2000) Nakayama, H., Fujii, Y., Ishiguro, M., et al., 2000, “Observations of Polarization and Brightness Variations with the Rotation for Asteroids 9 Metis, 52 Europa, and 1036 Ganymed”, Icarus, 146, 220.
  * Pravec, et al. (2002) Pravec, P., Harris, A. W., Michaĺowski, T., 2002, “Asteroid Rotations”, Asteroids III (W. F. Bottke, Jr. et al., eds.), 113.
  * Pravec, et al. (1998) Pravec, P., Wolf, M. & Sarounova, L., 1998, “Lightcurves of 26 Near-Earth Asteroids”, Icarus, 136, 124.
  * Schmidt, et al. (1992) Schmidt, G.D., Elston, R., & Lupie, O.L., 1992, “The Hubble Space Telescope Northern-Hemisphere Grid of Stellar Polarimetric Standards”, AJ, 104, 1563.
  * Thomas, et al. (1997) Thomas, P.C., Binzel, R.P., Gaffey, M.J., Storrs, A.D., Wells, E.N. & Zellner, B.H., 1997, “Impact excavation on asteroid 4 Vesta: Hubble Space Telescope results”, Science, 277, 1492.
  * Veeder, et al. (1981) Veeder, G.J., Tedesco, E.F., Tholen, D.J., Tokunaga, A., Kowal, C.T., Matthews, K., Neugebauer, G. & Soifer, B. T., 1981, “The diameter and albedo of 1943 Anteros”, Icarus, 46, 281.
  * Zellner, Gehrels & Gradie (1974) Zellner, B., Gehrels, T., & Gradie, J., 1974, “Minor planets and related objects. XVI. Polarimetric diameters”, AJ, 79, 1100.

Table 1: Polarimetry of (1943) Anteros

UT Obs Date | V mag |  $t_{exp}$ (sec) | $n_{meas}$ | $\alpha$ | Ecliptic Long | Lin $\%~{}$Pol | $\theta_{p}$
---|---|---|---|---|---|---|---
2009-07-22 | 16.6 | 130 | 21 | $40.5^{\circ}$ | $356.9^{\circ}$ | $2.53\pm 0.02$ | $178.8\pm 0.3^{\circ}$
2009-08-11 | 16.3 | 130 | 26 | $28.0^{\circ}$ | $352.9^{\circ}$ | $0.83\pm 0.02$ | $0.8\pm 0.6^{\circ}$
2009-08-26 | 16.1 | 90 | 27 | $19.1^{\circ}$ | $346.3^{\circ}$ | $-0.14\pm 0.02$ | $86\pm 4^{\circ}$
2009-09-10 | 16.3 | 100 | 30 | $16.4^{\circ}$ | $339.3^{\circ}$ | $-0.41\pm 0.02$ | $94\pm 1^{\circ}$

Figure 1:  New observations of the polarization of (1943) Anteros along with
model curves for S-type (dotted) and C-type (dashed), for reference. Note the
errors on the percent polarization are comparable to the size of the points.
Figure 2:  Polarization of Anteros as a function of rotation phase.
Observations within $0.1~{}$phase-wide bins have been co-added to reduce
errors. UT date of each observation set is listed.  Figure 3:  Polarization of
the background and the dark region of (1943) Anteros. The values of the
polarization for the dark region for phase angles below $40.5^{\circ}$ are
upper limits. Note that the error bars on the points are comparable to their
size.