# Stable Object Placing using Curl and Diff Features
of Vision-based Tactile Sensors

Kuniyuki Takahashi1, Shimpei Masuda1, Tadahiro Taniguchi2 1K.Takahashi and S.
Masuda are with Preferred Networks, Inc. {takahashi,
<EMAIL_ADDRESS>is with Ritsumeikan University, College of
Information Science and Engineering<EMAIL_ADDRESS>
###### Abstract

Ensuring stable object placement is crucial to prevent objects from toppling
over, breaking, or causing spills. When an object makes initial contact to a
surface, and some force is exerted, the moment of rotation caused by the
instability of the object’s placing can cause the object to rotate in a
certain direction (henceforth referred to as direction of corrective
rotation). Existing methods often employ a Force/Torque (F/T) sensor to
estimate the direction of corrective rotation by detecting the moment of
rotation as a torque. However, its effectiveness may be hampered by sensor
noise and the tension of the external wiring of robot cables. To address these
issues, we propose a method for stable object placing using GelSights, vision-
based tactile sensors, as an alternative to F/T sensors. Our method estimates
the direction of corrective rotation of objects using the displacement of the
black dot pattern on the elastomeric surface of GelSight. We calculate the
Curl from vector analysis, indicative of the rotational field magnitude and
direction of the displacement of the black dots pattern. Simultaneously, we
calculate the difference (Diff) of displacement between the left and right
fingers’ GelSight’s black dots. Then, the robot can manipulate the objects’
pose using Curl and Diff features, facilitating stable placing. Across
experiments, handling 18 differently characterized objects, our method
achieves precise placing accuracy (less than 1-degree error) in nearly 100% of
cases. 333An accompanying video is available at the following link:
https://youtu.be/fQbmCksVHlU

Figure 1:  Concept of the proposed method: The displacement of the black dots
on GelSight is analyzed in two patterns corresponding to an object’s
directions of corrective rotation: roll and pitch. By calculating the Curl and
Diff of the displacement of the black dot to estimate the object’s direction
of corrective rotation and manipulate the object accordingly, the robot
ensures its stable placement.

##  I Introduction

Stable placing of grasped objects is important for many reasons, such as 1)
preventing objects from toppling over, 2) preventing damage to the objects
themselves or their surroundings, and 3) keeping their contents, such as
liquids or small particles, from spilling. Consider settings like chemistry
and biology laboratories, for instance, where beakers and flasks are handled
routinely. With the growing demand for laboratory automation, the risk of
these objects toppling over with resulting breakage or spillage of their
contents has become a significant concern.

We define _stable placing_ as a state where a flat-bottomed object’s bottom
surface is parallel to a table or a non-flat object contacts the table at
three or more points. Conversely, _unstable placing_ is defined as the state
where a flat-bottom object’s bottom surface isn’t parallel to the table or a
non-flat object contacts the table at less than three points. Unstable
placement can result from many factors, such as inaccurate object pose
estimation.

The moment of rotation caused by the instability of the object’s placing can
cause the object to rotate in a certain direction (henceforth referred to as
direction of corrective rotation). The existing method for stable object
placing is to use a Force/Torque (F/T) sensor mounted on the robot’s wrist by
estimating the object’s direction of corrective rotation from the torque
applied to the grasped object [1]. However, several issues when using F/T
sensors make this estimation uncertain. 1) Tension of external wiring of robot
cables: The F/T sensor needs to be calibrated to compensate for the weight of
the robot’s wrist or gripper. However, changes in the robot’s pose can cause
variable tension in the cables, thereby complicating accurate calibration and
preventing precise force and torque detection. 2) Sensor noise: The F/T sensor
contains noise, which makes the measured sensor values uncertain. The F/T
sensors attached to the wrist often require a larger rated capacity. In
general, increasing the sensor’s rated capacity sacrifices the ability to
measure changes in a small range.

To address these issues, typical approaches include mounting a F/T or tactile
sensor on the fingertip. As precision manipulation needs detailed information
like fingertip contact area, interest in the use of tactile sensors,
particularly vision-based ones like GelSight [2], are growing [3, 4, 5]. In
this context, this study uses GelSight tactile sensors. Unlike F/T sensors
attached to the wrist, tactile sensors are attached to the end-effector’s
fingertip and are not affected by cable tension from other sensors (such as a
camera’s USB cable). Moreover, since detecting small magnitudes of force and
torque is required, they are typically designed to minimize the impact of
noise.

GelSight uses a camera to measure the deformation of an attached elastomer
during contact with a surface. Some GelSight devices feature black dots on
their surface to track their displacement when they make contact with a
surface. The displacement of GelSight’s black dots provides significant
information, and some previous studies used it for a learning-based approach
in peg-insertion task [6, 5].

Our work proposes utilizing these displacements as an alternative to torque
measurements traditionally given by F/T sensors to achieve stable object
placing. Our observations of how tactile sensors respond to various objects
during unstable placement suggest that the displacement patterns of GelSight’s
black dots can be categorized into two types with two directions of corrective
rotation: roll and pitch (Fig. 1). After analyzing each classification, the
displacement of GelSight’s black dots can be expressed as a one-dimensional
feature for each direction of corrective rotation. The direction of corrective
rotation of the object is estimated by calculating the Curl from vector
analysis for the roll axis, indicative of the rotational field magnitude and
direction of the displacement of the black dots pattern, and the difference
(Diff) of displacement for the pitch axis between the left and right fingers’
GelSight’s black dots. These features allow accurate stable object placing on
the robot. The novelty and contribution of this method are as follows:

  * •

Estimation of the direction of corrective rotation of objects for stable
placing using a training-free method employing Curl and Diff features based on
the GelSight’s black dots displacement. To the best of our knowledge,
achieving stable placement with proposed Curl and Diff is a novel and non-
trivial achievement.

  * •

The use of these features yielded a nearly 100% success rate of stable
placing, achieving a consistently high rate across various objects with less
than a 1-degree error. Despite the apparent simplicity of our proposed method,
it is versatile enough to be applied to many objects.

##  II Related Works

###  II-A Object Pose Estimation

When considering a placing task, a typical first step is the use of object
pose estimation methods. Object pose estimation is a well-studied problem in
computer vision and is important for robotic tasks. While many approaches are
image-based, they often face issues associated with occlusion and coarse
accuracy [7, 8]. It is worth noting that there exist methods that reduce
uncertainty and increase object pose estimation accuracy by having the grasped
object contact other objects [9, 10, 11]. However, these require prior
knowledge of the object’s shape, limiting their use for arbitrary objects.

Research utilizing tactile sensors has demonstrated their ability to estimate
pose with greater accuracy [3, 12]. In particular, vision-based tactile
sensors, such as GelSight, have been used in recent research, and some have
achieved success in tasks with milli-order to sub-milli-order accuracy [3, 5,
4, 13, 14]. However, these methods only estimate an object’s relative or
absolute pose and do not consider whether the object can be placed stably on a
desk or other surface. If there is an error in the estimation result, the
object may be unstable and may topple over. Therefore, the focus of this paper
is the stable placing of objects.

###  II-B Object Placing

Object placing can be approached from both hardware and software perspectives.
A well-known method for the hardware approach is to incorporate passive
compliance in the robot’s wrist [15]. This mechanical compliance allows the
robot to adjust to its environment and achieve stable placing without
intricate control manipulations. However, hardware approaches are often
specialized for specific objects, making them less flexible in accommodating a
variety of objects.

As for the software approach, conventional methods involve control of applied
force using robot joint torque or F/T sensors [1]. Yet, as discussed in the
introduction, these methods often require substantial engineering effort to
address problems like cable tension and sensor noise.

Several studies have proposed the use of tactile sensors, either solely or
combined with F/T sensors, to estimate and control the force and torque on an
object [14, 16, 5, 6]. These studies can be broadly divided into learning-
based and analytical approaches. The learning-based approach, such as [14,
16], predicts where and how much external forces are applied to an object and
performs object manipulation based on these predictions. In [5], a policy for
controlling the object pose directly from raw images of GelSight is proposed.
These studies indicate that tactile sensors can be used as F/T sensors, either
directly or indirectly, by learning from raw images of vision-based tactile
sensors and flow images showing tactile sensor displacement. While these
approaches can generally adapt to various objects, the prerequisite of
extensive training can be time-consuming and may result in out-of-distribution
problems if data is insufficient.

In an analytical approach, [6] successfully measured the force and torque
applied to a single GelSight sensor by using the Helmholtz-Hodge Decomposition
to calculate the relationship between multiple forces and torques applied to
different objects and GelSight’s black dots displacement. This method does not
require the collection of large-scale data on various objects. However, to
fully capture the potential of these approaches, a conversion to force and
torque is attempted, which may necessitate calibration with dedicated
equipment. In our work, we propose an analytical approach that avoids the need
for calibration and shows promise in achieving stable object placing across
various objects.

Figure 2: Robot setup and coordinate system used for the experiment

##  III Preliminaries

As Preliminaries, this section describes the estimation of the direction of
corrective rotation of an object when the object is placed unstably by
pressing the object against a desk. Then, the robot control method for stable
placing based on sensor information is described.

Three coordinate systems used in this study are described below (Fig. 2). The
first one is the robot coordinate system, which is set at the base of the
robot arm. The second is the end-effector (EE) coordinate system, which is set
at the center of the two robot gripper’s fingertips. The robot manipulates
objects by applying force and torque to them for stable placing. The point
where the robot arm generates force is the origin of the EE coordinate system.
In this paper, ‘the point where the robot arm generates force’ is referred to
as the force application point. The third is the tactile sensor coordinate
system, in which the robot coordinate system is rotated around the z-axis so
that the x-axis direction of the tactile sensor coordinate system is parallel
to the x-axis direction of the EE coordinate system.

Figure 3: The estimation of the direction of corrective rotation of an object
when the object is pressed against a desk.

###  III-A Estimation of Direction of Corrective Rotation

The idea behind the estimation of the direction of corrective rotation can be
explained by considering the moment around the point of contact between the
object and the environment [17]. The direction of corrective rotation is
determined based on whether the force application point is to the right or
left of the vertical line from the object’s contact point (Fig 3 (1)). The
direction of corrective rotation is clockwise when the force application point
is on the right side of the vertical line and counterclockwise when the force
application point is on the left. If the force application point is on the
vertical line, the object does not rotate.

The sign of the torque at the force application point indicates the direction
of corrective rotation. The torque is positive, indicating clockwise rotation,
if the force application point is to the right side of the vertical line;
negative, for counterclockwise rotation, if on the left. We note that
depending on the convention used, positive torque can indicate either
clockwise or counterclockwise rotation but should stay consistent across
future experiments. If the torque is zero, the force application point is on
the line, indicating that the object does not rotate. Based on the sign of
this torque detected by sensors, the robot is controlled to perform stable
placing of the object. This control method of the robot is described in
Section III-B.

Sensors measure these torques to estimate the direction of corrective
rotation. While the method using the F/T sensor (referred to in this paper as
the F/T sensor-based method, which is the baseline for this paper) can
estimate the direction of corrective rotation, it is prone to inaccuracies
primarily due to the cable tension and sensor noise. Especially near the
vertical line, misrecognition of the direction of corrective rotation may
occur as the measured sign of torque is likely to invert. Our proposed method
uses tactile sensors to overcome these issues, providing a more accurate
estimation of the direction of corrective rotation. Although tactile sensors
measure different physical quantities than torque sensors, their measurements
can still be interpreted as positive or negative with respect to the vertical
line, providing a consistent approach to the estimation of the direction of
corrective rotation.

###  III-B Robot Control for Stable Placing

The previous section described that the direction of corrective rotation of
the object can be estimated based on the sign of the torque. This section
describes how to control the robot for stable placing based on the torque.

If the torque measured at the force application point is not zero, the object,
which is in an unstable state, can be rotated in the direction of corrective
rotation. Stable placing can be performed by controlling the object by
applying torque to the force application point so that the torque becomes 0.
We make assumptions for our model, which are detailed below: As shown in Fig.
3 (2), the force application point is in the area where the object can be
placed stably. We note that the area where this stable placement is feasible
is not constant but varies depending on the position of the vertical line from
the object’s contact point at the time the object is unstably placed. Under
this assumption, manipulating the object so that the detected torque is zero
results in a stable placing.

To realize stable placing, the robot arm is controlled by admittance control
so that force
$F_{target}=(F_{target-x},F_{target-y},F_{target-z})=(0,0,F_{tz})$ and torque
$\tau_{target}=(\tau_{target-x},\tau_{target-y},\tau_{target-z})=(0,0,0)$ for
the force application point exerts the force and torque in the robot
coordinate system. The combination of force and torque is called a wrench. Let
$w_{target}=(F_{target},\tau_{target})$ be the target wrench and
$w_{current}=(F,\tau)$ be the current wrench. For this target wrench, when the
object is pushed against the desk, the reaction force from the desk is
$F_{tz}$, and the other forces and torques are 0. In this experiment,
$F_{tz}=5\mathrm{\,N}$.

The end-effector velocity of the robot arm, $v_{ee}$ from the EE coordinate
system, is calculated as follows to achieve the target wrench $w_{target}$:

$v_{ee}=k(w_{current}-w_{target})$ (1)

where $k$ is a gain. In this paper, the same value for this gain is used for
all objects.

In the case of the F/T sensor-based approach, the current wrench $w_{current}$
uses an F/T sensor mounted on the wrist. Let’s denote the force measured by
the F/T sensor on the wrist as $F^{\prime}$ and the torque as $\tau^{\prime}$.
The force $F$ and torque $\tau=(\tau_{x},\tau_{y},\tau_{z})$ at the force
application point can be calculated as follows:

$\displaystyle F$ $\displaystyle=R\cdot F^{\prime}$ $\displaystyle\tau$
$\displaystyle=R\cdot\tau^{\prime}+d\times(R\cdot F^{\prime})$ (2)

Here, $R$ is a 3x3 rotation matrix, and $d$ is a 3x1 displacement vector from
the F/T sensor coordinate to the EE coordinate in the F/T sensor coordinate.

In our method, we use the Curl and Diff features values calculated by
GelSight, instead of the $\tau_{x}$ and $\tau_{y}$ values for roll and pitch
that are calculated from the F/T sensor.

##  IV Method: Tactile-based Placing

The previous section described that the direction of corrective rotation can
be estimated based on the sign of the torque. This section describes how to
use tactile sensors as an alternative to F/T sensors to estimate the direction
of corrective rotation.

Our observations of how tactile sensors respond to various objects during
unstable placement suggest that the displacement patterns of GelSight’s black
dots can be categorized into two types with two directions of corrective
rotation: roll of the tactile sensor coordinate system and pitch of the EE
coordinate system (Fig. 1). The following section describes how to calculate
the rotations for pitch and roll based on the displacement of the black dots
on the GelSight.

Figure 4: Schematic of torque applied to fingertip

###  IV-A Calculation of Rotation in Pitch Direction by Curl

We describe a method where the Curl features calculated from the displacement
of the black point of the tactile sensor, used as an alternative to $\tau_{y}$
(pitch rotation at the force application point) calculated from the F/T sensor
by eq. 2.

When the object is pressed against the desk, a torque is generated that
attempts to rotate the object (Fig. 4 (1)). This torque is transmitted to the
GelSight and results in a shift in the position of the black dots on the
GelSight, which represents the applied torque. This torque at the contact
point on GelSight can be observed as the displacement of the black dots on
GelSight (Fig. 1 (1) Pitch rotation). From the red arrows in Fig. 1 (1), when
the left side of the object contacts the desk (with the gripper considered as
the center), we observe a displacement of the arrows representing the black
dots on the GelSight moving in a specific rotational direction. If the right
side of the object is in contact instead, the arrows exhibit the opposite
rotational movement. Specifically, in the former case, the arrows move in a
clockwise direction, whereas in the latter, they move in a counterclockwise
direction.

When the displacement of these arrows, which represent the movement of the
black dots on the GelSight, becomes zero, the object is in a stable placing
state with respect to the pitch rotation. To achieve a stable placement of the
object by minimizing the torque, the robot should rotate around the force
application point in the direction indicated by these arrows, which in turn
reduces the displacement of the black dots. We note that the contact point is
the fingertip, and the force application point is the middle of the contact
point between the two fingers. Thus, this rotation center is considered the
origin of the EE coordinate system.

From the arrows, the direction of rotation and its degree of rotation are
calculated by Curl of vector analysis. The Curl is indicative of the
rotational field magnitude and direction of the displacement of the black dots
pattern. The displacements of the black dots on the GelSight give a vector
field $A_{tactile}=(A_{x},A_{z})$. In the EE coordinate system, the Curl of
this field can be computed as follows:

$\text{Curl}\>A_{tactile}=\partial A_{z}/\partial x-\partial A_{x}/\partial z$
(3)

We use the mean of the calculated Curl of all black dots as an alternative to
$\tau_{y}$ calculated from the F/T sensor by eq. 2, the torque at the force
application point. The goal is to control the end-effector speed by eq. 1 of
the robot so that this value becomes $0$.

###  IV-B Calculation of Rotation in Roll Direction by Diff

We describe a method where the Diff features calculated from the displacement
for the pitch axis between the left and right fingers’ GelSight’s black dots,
used as an alternative to $\tau_{x}$ (roll rotation at the force application
point) calculated from the F/T sensor by eq. 2.

When the object is pressed against the desk, a torque is generated that
attempts to rotate the object (Fig. 4 (2)). This torque is transmitted to the
GelSight, and is exerted as shear forces in opposite directions on the right
and left GelSight. This makes a difference in the displacement of the black
dots in GelSights between the left and right fingers. The shear force arising
at the contact point on GelSights can be observed by the displacement of the
black dots on the GelSights (Fig. 1 (2) Roll rotation). Additionally, there
are also forces generated by the action of pushing the object against the
table that is exerted on the GelSights, whose magnitudes may change the
overall direction of displacement. Fig. 1 (2) demonstrates the displacement of
the black dot of GelSight when the tilted object makes contact with the desk.
On the side where the object is in contact, the displacement of the black dot
is represented as an arrow pointing vertically upward (in the z-axis direction
of the tactile sensor coordinate system). This displacement is noticeably
larger compared to the side where the object is not making contact.

When there is no longer any imbalance between the left and right of the
displacement of the black dot of GelSight, the object is in a state of stable
placing with respect to the direction of rotation of the roll. In order to
reduce the imbalance of the left and right of the displacement of the black
dot of GelSight, the robot should rotate around the force application point
with respect to the direction of the object’s rotation. We note that the
contact point is the fingertip, and the force application point is the middle
of the contact point between the two fingers. Thus, this rotation center is
considered the origin of the EE coordinate system.

In this study, we use two GelSight sensors, GelSight1 and GelSight2.
Considering the black dots’ displacement in the z-axis direction of the
tactile sensor coordinate system (normal facing upwards to the desk) helps
calculate the GelSights’ displacement affecting the roll axis. This makes it
possible to be irrespective of the pose of the robot and GelSights. We compute
the average displacement values of the black dots on each of the two GelSights
along the z-axis of the tactile sensor coordinate system. We denote these
values as $\bar{A}_{z1}$ and $\bar{A}_{z2}$ respectively. We define a metric
Diff as the difference between $\bar{A}_{z1}$ and $\bar{A}_{z2}$:

$Diff=\bar{A}_{z1}-\bar{A}_{z2}$ (4)

Diff serves as the torque at the force application point in the tactile sensor
coordinate system. The goal is to control the robot’s end-effector speed to
reduce Diff to zero. We note that the end-effector speed calculated from Diff
corresponds to the speed measured from the tactile sensor coordinate system.
However, as $\tau_{y}$ in eq. 2 is calculated in the EE coordinate system, we
must convert the value of Diff from the tactile sensor to the EE coordinate
system. This can be achieved using standard coordinate transformation methods.

##  V Experimental Setup

###  V-A Robot Setup

Our robotic system, as shown in Fig. 2, features a 7-DoF Franka Emika Panda
Arm. The end-effector command to the robot is updated at $200\mathrm{\,Hz}$.
This robotic arm is equipped with a custom-designed parallel gripper. The
gripper features vision-based tactile sensors, known as GelSight Mini [18].
The sampling rate of GelSights is $10\mathrm{\,Hz}$. A Leptrino force/torque
sensor (FFS055YA501U6) is mounted between the robot arm and the gripper. The
rated capacity of force of F/T sensor is $F_{xyz}=\pm 500\mathrm{\,N}$, torque
is $T{xyz}=\pm 4\mathrm{\,Nm}$, resolution is $\pm 1/2000$, and sampling rate
is $200\mathrm{\,Hz}$. An Intel RealSense Depth Camera D435i is used to
capture the pose of the arUco marker [19]. The system runs on a PC equipped
with 32 GB RAM and an Intel Core i7-8700 CPU, running Ubuntu 20.04.6 LTS with
ROS Noetic.

Figure 5: Experiment process for object placing Figure 6: Target objects used
in the experiment and GelSight images while grasping the target objects

###  V-B Task Setup and Process of Object Placing

We conduct experiments where a robot grasps objects and places them on a desk
to evaluate our proposed system. To aid evaluation, an arUco marker is placed
on the object, and the object’s pose is measured when the object is placed on
the desk. The displacement is then calculated by comparing the object’s pose
initially in its manually assigned stable state on the table with the final
pose achieved through the robot’s stable placing process. The detailed steps
of the experimental procedure are graphically illustrated in Fig. 5.

This experiment begins by manually placing an object with an attached arUco
marker of size $25~{}\mathrm{\,mm}$ onto a desk. We note that the experiment
is repeated 10 times for each of the proposed and baseline methods but for the
same object placed in the same location. The robot moves to a predetermined
pre-grasping position (Fig. 5 (1)). The pre-grasping positions are specific to
each object.

The robot grasps the object, and then, the camera (RealSense) in front of the
robot (Fig. 2), with a sampling rate of $30\mathrm{\,Hz}$, is used to capture
the pose of the arUco marker (Fig. 5 (2)). The initially placed object is
defined as being in a stably placed state in this experiment. To ensure
accuracy and reduce potential noise during marker recognition, we use the
average pose from the marker across 90 static images captured by the camera
over three seconds.

After capturing the arUco marker, the robot lifts it off the ground by moving
along the z-axis, and tilts it (Fig. 5 (3)-(4)). To evaluate the stable
placing, the grasped object is tilted (either fixed or randomly sampled) in
this pre-placing process. The way the tilt is specified is different for each
experiment and will be described in the section VI. In the tilted state, the
value read from the F/T sensor is recorded as an offset. During control
procedures, we will use the value derived from the F/T sensor after removing
this offset. The displacement of the black dots of the tactile sensor is
calculated from the position of the black dots of the tactile sensor at this
stage.

The robot then moves under velocity control at $0.005~{}\mathrm{\,m/s}$ in the
-z-axis direction in the robot coordinate system (Fig. 5 (5)). When the force
value of the F/T sensor along the z-axis of the robot coordinate system
exceeds $1\mathrm{\,N}$, the process of stable placing begins using both the
proposed and baseline methods (Fig. 5 (5)).

After sufficient time has passed to ensure the completion of stable placing,
the placing process is stopped, and the pose of the arUco marker is captured
again (Fig. 5 (6)). As with the acquisition of the arUco marker taken in the
beginning, the mean value of the 6d pose of the arUco marker was calculated
from 90 static images captured over three seconds. Under the control of stable
placing, the displacements in the roll and yaw axes of the arUco marker are
evaluated. This evaluation is done by comparing the marker’s pose before and
after the stable placing procedure is completed (Fig. 2).

Finally, the robot releases the object and returns to the home position. We
performed this experimental process ten times for each object and calculated
the mean and variance of the displacements of the arUco marker at that time.

###  V-C Target Objects

For the evaluation, we chose 18 objects considering various factors such as
size (lengths of the short and long axes and thickness), softness, surface
friction, texture, left-right asymmetrical shape and texture, and whether the
object had a flat contact surface or not (Fig. 6).

##  VI Experiment Results

This evaluation focused on two key aspects: 1) Evaluation of stable placing,
considering changes in the object’s grasp positions and the size of its
support polygon \- the area defined by the contact points between the object
and the surface upon which it rests. 2) Comparison of the proposed method,
Tactile-based method, with baseline, F/T sensor-based method across a variety
of objects.

TABLE I:  The accuracy of stable placing when the grasping position and the size of the support polygon of the object are different.  Object |  Grasp height | Tactile-based (Ours) | F/T-based (Baseline)
---|---|---|---
Roll [deg] | Pitch [deg] | All [deg] | $<1$ deg | $<2$ deg | Roll [deg] | Pitch [deg] | All [deg] | $<1$ deg | $<2$ deg
Large rectangular | $1.0~{}\mathrm{\,cm}$ | $0.4\pm 0.2$ | $0.2\pm 0.2$ | $0.3\pm 0.2$ | $\bf{10/10}$ | $\bf{10/10}$ | $0.4\pm 0.2$ | $0.2\pm 0.2$ | $0.3\pm 0.2$ | $\bf{10/10}$ | $\bf{10/10}$
$2.5~{}\mathrm{\,cm}$ | $0.4\pm 0.3$ | $0.2\pm 0.1$ | $0.3\pm 0.2$ | $\bf{10/10}$ | $\bf{10/10}$ | $0.3\pm 0.2$ | $0.3\pm 0.2$ | $0.3\pm 0.2$ | $\bf{10/10}$ | $\bf{10/10}$
Small rectangular | $1.0~{}\mathrm{\,cm}$ | $0.3\pm 0.1$ | $0.1\pm 0.1$ | $0.2\pm 0.2$ | $\bf{10/10}$ | $\bf{10/10}$ | - | - | - | $6/10$ | $6/10$
$2.5~{}\mathrm{\,cm}$ | $0.3\pm 0.2$ | $0.1\pm 0.1$ | $0.2\pm 0.2$ | $\bf{10/10}$ | $\bf{10/10}$ | - | - | - | $0/10$ | $0/10$
TABLE II:  The accuracy of stable placing.  Object type | Object | Tactile-based (Ours) | F/T-based (Baseline)
---|---|---|---
Roll [deg] | Pitch [deg] | All [deg] |  $<1$ deg |  $<2$ deg | Roll [deg] | Pitch [deg] | All [deg] |  $<1$ deg |  $<2$ deg
Elongated shape | Small rectangular | $0.4\pm 0.2$ | $0.1\pm 0.1$ | $0.2\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $2/10$ |  $2/10$
Large rectangular | $0.4\pm 0.3$ | $0.3\pm 0.2$ | $0.4\pm 0.3$ |  $\bf{10/10}$ |  $\bf{10/10}$ | $0.4\pm 0.2$ | $0.4\pm 0.3$ | $0.4\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$
Wooden rectangular | $0.4\pm 0.2$ | $0.2\pm 0.1$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $7/10$ |  $7/10$
Wooden cylinder | $0.2\pm 0.2$ | $0.1\pm 0.1$ | $0.2\pm 0.1$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $6/10$ |  $7/10$
Beaker | $0.2\pm 0.2$ | $0.3\pm 0.1$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $4/10$ |  $4/10$
Metal rectangular | $0.2\pm 0.1$ | $0.2\pm 0.1$ | $0.2\pm 0.1$ |  $\bf{10/10}$ |  $\bf{10/10}$ | $0.3\pm 0.2$ | $0.1\pm 0.1$ | $0.2\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$
V shape | V block | $0.4\pm 0.3$ | $0.2\pm 0.1$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $1/10$ |  $1/10$
Joint | $0.7\pm 0.5$ | $0.2\pm 0.1$ | $0.4\pm 0.5$ |  $\bf{8/10}$ |  $\bf{10/10}$ | - | - | - |  $0/10$ |  $0/10$
Asymmetry | LEGO blocks | $0.3\pm 0.2$ | $0.9\pm 0.1$ | $0.6\pm 0.3$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $6/10$ |  $6/10$
Metal object | $0.2\pm 0.2$ | $0.2\pm 0.1$ | $0.2\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $0/10$ |  $0/10$
Non-surface contact | Metal pedestal | $0.2\pm 0.2$ | $0.2\pm 0.1$ | $0.2\pm 0.1$ |  $\bf{10/10}$ |  $\bf{10/10}$ | $0.2\pm 0.2$ | $0.1\pm 0.1$ | $0.2\pm 0.1$ |  $\bf{10/10}$ |  $\bf{10/10}$
Conversion connector | $0.4\pm 0.2$ | $0.2\pm 0.1$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $2/10$ |  $2/10$
Soft material | Wood bond tube | $0.3\pm 0.2$ | $0.2\pm 0.1$ | $0.2\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $1/10$ |  $1/10$
Hand cream tube | $0.4\pm 0.3$ | $0.4\pm 0.2$ | $0.4\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $4/10$ |  $4/10$
Fewer contact | Test tube stand | $0.1\pm 0.1$ | $0.6\pm 0.2$ | $0.3\pm 0.3$ |  $\bf{10/10}$ |  $\bf{10/10}$ | $0.2\pm 0.2$ | $0.5\pm 0.2$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$
Metal part | $0.4\pm 0.2$ | $0.1\pm 0.1$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | $0.3\pm 0.2$ | $0.1\pm 0.1$ | $0.2\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$
Round shape | Round bottom flask | $0.4\pm 0.2$ | $0.2\pm 0.2$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $0/10$ |  $0/10$
Liquid | Bottle with liquid | $0.4\pm 0.3$ | $0.2\pm 0.2$ | $0.3\pm 0.2$ |  $\bf{10/10}$ |  $\bf{10/10}$ | - | - | - |  $4/10$ |  $5/10$

###  VI-A Different Grasp Positions and Support Polygon Size

We evaluate the influence of the object’s grasping position and the size of
the object’s support polygon on stable placing. To investigate these effects,
we utilize two rectangular objects, small rectangular and large rectangular,
which have the same height and different sizes of support polygon (Fig. 6).
The dimensions of small and large rectangular are, respectively,
$30~{}\mathrm{\,mm}\times 30~{}\mathrm{\,mm}\times 70~{}\mathrm{\,mm}$,
$50~{}\mathrm{\,mm}\times 50~{}\mathrm{\,mm}\times 70~{}\mathrm{\,mm}$ in
terms of width, depth, and height. These two objects were made by 3D printing,
and they both weigh $142\mathrm{\,g}$. The two grasping positions are at
heights of $1.0~{}\mathrm{\,cm}$ and $2.5~{}\mathrm{\,cm}$ from the desk. To
compare under the same noise conditions, the stable placing experiment,
following the process outlined in section V-B, was conducted 10 times with the
object tilted 10 degrees to the roll and pitch in the tactile-sensor
coordinate in the pre-placing process (Fig. 5 (4)). Table I presents the mean
and standard deviation of the roll and pitch misalignment and their combined
values from ten iterations of stable placing experiments. However, when the
error magnitude exceeds 2 degrees, we do not compute the mean and standard
deviation as the means appear exceedingly large when the object topples over.
We also report the number of times when roll and pitch errors were below 1 and
2 degrees. We chose the thresholds of 1 degree and 2 degrees based on our
experimental considerations. Specifically, the 1-degree threshold was selected
due to the detection accuracy of the arUco marker. The 2-degree threshold was
chosen because our experimental experience suggested that exceeding this value
often leads to the object toppling over.

Using our proposed tactile-based method, stable placing was achieved under all
conditions regardless of the grasping position and the size of the support
polygon. In the case of the baseline F/T-based method, for the large
rectangular object, stable placing was achieved at both grasp heights. For the
small rectangular object, while there were six successful outcomes with a low
grasping height ($1.0~{}\mathrm{\,cm}$), there were no successes when the
grasping height was high ($2.5~{}\mathrm{\,cm}$).

The reasons for this are explained below. The direction of corrective rotation
can be estimated through the sign of the torque. For a rectangular object, as
illustrated in Fig. 3 (2), bias in estimating the direction of corrective
rotation depends on the grasping position and the size of the supporting
polygon. A low grasping position or large support polygon may still enable
accurate estimation despite sensor inaccuracies since the area in which stable
placement can be achieved is larger. Conversely, a high grasping position or
small support polygon decreases the area in which stable placing can be
feasible and increases the risk of sensor inaccuracies negatively affecting
the direction of corrective rotation estimation. F/T sensors’ measured values
have low accuracy due to their noisy nature. We note that a low pass filter is
used to remove noise from the F/T sensor, but noise that cannot be removed
still has an effect. Therefore, even under the same conditions of conducting
10 stable placings with a 10-degree tilt in roll and pitch, failures may
occasionally occur due to the inherent uncertainty of the F/T sensor.
Additionally, cable tension can cause a shift in the torque value, which could
potentially lead to misinterpretations of the object’s direction of corrective
rotation.

Figure 7: Examples of stable placing. For the sake of visibility, the length
of the arrows is eight times the actual displacement.

###  VI-B Evaluation across Diverse Objects

We evaluate the proposed method and a baseline F/T-based method for the stable
placing of 18 different objects. This stable placing experiment, which follows
the process outlined in section V-B, was conducted 10 times with the object
tilted. The tilt in the pre-placing process (Fig. 5 (4)) is calculated from a
normal distribution with a mean of $\pm 10$ degrees (and the direction of the
tilt is uniformly sampled to be positive or negative) and a variance of 1 in
both roll and pitch directions in the tactile-sensor coordinate. For objects
prone to toppling over, we specify their tilts as follows: The conversion
connector, which is tilted even in the stable placing state, has a normally
distributed noise with a mean of $\pm 10$ degrees in the roll direction and a
mean of $\pm 3$ degrees in the pitch direction, while the Joint object has a
mean of -5 degrees for noise in the pitch direction and a mean of 10 degrees
in the roll direction. For hand cream, which has a small support polygon, a
mean of $\pm 5$ degrees. The variance of tilt for all items tested is 1. The
Table II shows the results.

In the Tactile-based method, stable placing with an error of less than 1
degree was achieved for all objects except for the Joint object. The method
was successful for objects with asymmetrical shapes and textures, objects with
a point of contact with the desk instead of a surface, soft objects, objects
with a small size of support polygon, and objects with a changing center of
gravity, such as liquids. Fig. 7 shows examples of the stable placing process
using the proposed method. The Curl and Diff values are nearly zero for a
variety of objects, and stable placement is achieved.

Additionally, the method was also robust against cracks in the sensor of the
surface and manufacturing differences in the size and shape of the black dots.
The image of the tactile sensor in Fig. 6 shows that the surface of the
tactile sensor is cracked. The tactile sensor was completely damaged and was
replaced by another GelSight, which was also handled with a large black dot
due to individual differences in the GelSight (Metal object and Test tube
stand in Fig. 7).

For the Joint object, while it did not topple over, an error of less than 2
degrees remained. We analyze this occurrence as follows: The displacement of
the black dots on GelSight is calculated based on their positions before
stable placing when noise is added to the object’s pose. In the Joint object
experiment, the object was grasped at the tip of the V-shape, which was far
from the center of gravity. As gravity causes the object to rotate, the black
dots on GelSight deviate from their positions in the stable state. Even when
the Joint object reaches a stable placement, some displacement of the black
point remains present. Thus, the robot arm is controlled to apply force to the
force application point to reduce the displacement of the black point. Unlike
objects with a larger support polygon that retain their pose, the Joint
object, with its smaller support polygon, appears to tilt due to the force
intended for its rotation. However, the tilt is small, less than 2 degrees,
and the risk of toppling over is small.

The Tactile-based method achieved high accuracy for a variety of objects,
while the F/T-based method often failed. A common feature of the objects with
a 100% success rate in the F/T-based method is that the support polygon is
large, such as large rectangular, metal rectangular, metal pedestal, test tube
stand, and metal part. Objects with small support polygon failed with
probability or did not succeed at all. As explained in Section VI-A, the F/T
sensor is affected by sensor noise and the cable tension, causing uncertainty
in sensor values and misrecognizing the direction of rotation of the object.
As a result, there were many failures that resulted in the object toppling
over in the F/T sensor-based approach. There were also failures that the
object’s tilt was kept after contact between the object and the desk.
Moreover, some were seen where the robot arm applied force in the direction
that the object was to be upright, but the object passed through the upright
state and toppled over.

In summary, the F/T sensor method fails in stable placing due to sensor noise
and the cable tension, causing deviations in the measured torque values, which
result in misrecognition of the object’s direction of corrective rotation. On
the other hand, the Tactile-based method succeeds in stable placing with high
accuracy as a result of its precise estimation of the object’s direction of
corrective rotation.

##  VII Conclusion

We propose a method for stable object placing using tactile sensors. Under
unstable object placing conditions, we observed that the displacement pattern
of the black dots in GelSight aligns with the two possible directions of
corrective rotation. By estimating these two object directions of corrective
rotation from Curl and Diff, the robot can manipulate the object’s pose for
stable placing. This method of estimating the direction of corrective rotation
is robust against sensor failures such as cracks in the surface, and
manufacturing differences in the size and shape of the black dots. In
experiments with 18 objects with different characteristics, the F/T sensor-
based method fails because it misrecognizes the direction of rotation for
objects other than those with a large support polygon due to sensor noise and
cable tension. In contrast, the proposed Tactile sensor-based approach can
achieve stable placing with high accuracy (less than 1-degree error) in nearly
100% of cases across all objects, highlighting its potential as an effective
F/T sensor alternative.

## ACKNOWLEDGMENT

The authors thank Avinash Ummadisingu and Dr. Naoki Fukaya for the many
discussions about this research. This work was supported by JST [Moonshot
R&D][Grant Number JPMJMS2033].

## References

  * [1] M. H. Raibert and J. J. Craig, “Hybrid position/force control of manipulators,” 1981.
  * [2] W. Yuan, _et al._ , “GelSight: High-Resolution Robot Tactile Sensors for Estimating Geometry and Force,” _Sensors_ , vol. 17, no. 12, p. 2762, 2017.
  * [3] T. Anzai and K. Takahashi, “Deep gated multi-modal learning: In-hand object pose changes estimation using tactile and image data,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 9361–9368.
  * [4] M. Bauza, _et al._ , “Tac2pose: Tactile object pose estimation from the first touch,” _The International Journal of Robotics Research_ , vol. 42, no. 13, pp. 1185–1209, 2023.
  * [5] S. Dong, _et al._ , “Tactile-rl for insertion: Generalization to objects of unknown geometry,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 6437–6443.
  * [6] Y. Zhang, _et al._ , “Effective estimation of contact force and torque for vision-based tactile sensors with helmholtz–hodge decomposition,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 4094–4101, 2019.
  * [7] J. Sun, _et al._ , “Onepose: One-shot object pose estimation without cad models,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 6825–6834.
  * [8] B. Li, _et al._ , “VoxDet: Voxel Learning for Novel Instance Detection,” in _Proceedings of the Advances in Neural Information Processing Systems (NeurIPS)_ , 2023.
  * [9] N. Chavan-Dafle and A. Rodriguez, “Stable prehensile pushing: In-hand manipulation with alternating sticking contacts,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 254–261.
  * [10] F. Von Drigalski, _et al._ , “Contact-based in-hand pose estimation using bayesian state estimation and particle filtering,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 7294–7299.
  * [11] J. Pankert and M. Hutter, “Learning Contact-Based State Estimation for Assembly Tasks,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2023, pp. 5087–5094.
  * [12] L. Lach, _et al._ , “Placing by Touching: An empirical study on the importance of tactile sensing for precise object placing,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2023, pp. 8964–8971.
  * [13] S. Pai, _et al._ , “Laboratory Automation: Precision Insertion with Adaptive Fingers utilizing Contact through Sliding with Tactile-based Pose Estimation,” _arXiv preprint arXiv:2309.16170_ , 2023.
  * [14] K. Ota, _et al._ , “Tactile Estimation of Extrinsic Contact Patch for Stable Placement,” _arXiv preprint arXiv:2309.14552_ , 2023.
  * [15] D. E. Whitney _et al._ , “Quasi-static assembly of compliantly supported rigid parts,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 104, no. 1, pp. 65–77, 1982.
  * [16] C. Higuera, _et al._ , “Neural contact fields: Tracking extrinsic contact with tactile sensing,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2023, pp. 12 576–12 582.
  * [17] M. Shimizu and K. Kosuge, “Spatial parts mating with fiction using structured compliance with compliance center,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , vol. 2. IEEE, 2002, pp. 1585–1590.
  * [18] “Gelsight mini,” https://www.gelsight.com/gelsightmini/, 2024.
  * [19] S. Garrido-Jurado, _et al._ , “Automatic generation and detection of highly reliable fiducial markers under occlusion,” _Pattern Recognition_ , vol. 47, no. 6, pp. 2280–2292, 2014.


# Statistical study of the optimal local sources for cosmic ray nuclei and
electrons

Qing Luo  Hebei Normal University, Shijiazhuang, Hebei 050024, China  Key
Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese
Academy of Sciences, Beijing 100049, China  Bing-qiang Qiao  Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Wei Liu  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  Shu-wang Cui  Hebei Normal University, Shijiazhuang,
Hebei 050024, China  Yi-qing Guo  Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China  University of Chinese Academy of Sciences, Beijing 100049, China

###### Abstract

The local sources, such as the Geminga supernova remnant (SNR), may play an
important role in the anomaly of protons, electrons and anisotropy in the past
works. In fact, there exist twelve SNRs around the solar system within $1$
kiloparsec (kpc). One question is that can other SNRs also possibly contribute
to the spectra of nuclei and electrons and explain the special structure of
the anisotropy? In this work, under the spatial-dependent propagation, we
systematically study the contribution of all local SNRs within 1 kpc around
the solar system to the spectra of nuclei and electrons, as well as the energy
dependence of the anisotropy. As a result, only the Geminga, the Monogem, and
the Vela SNRs have quantitative contribution to the nuclei and electron
spectra and the anisotropy. Here, the Geminga SNR is the sole optimal
candidate and the Monogem SNR is controversial due to the tension of the
anisotropy between the model calculation and the observations. The Vela SNR
contributes to a new spectral structure beyond TeV energy, hinted by the HESS,
the VERITAS, the DAMPE, and the CALET measurements. More interestingly, the
electron anisotropy satisfies the Fermi-LAT limit below TeV energy, but rises
greatly and reaches $10\%$ at several TeV. This novel structure will shed new
light on verifying our model. We hope that the new structure of the electron
spectrum and anisotropy can be observed by the space-borne DAMPE and HERD, and
the ground-based HAWC and LHAASO experiments in the near future.

##  1 Introduction

It is well known that supernova remnants (SNRs) are the dominant sources of
galactic cosmic rays (GCRs) (Boulares, 1989; Blasi, 2013; Ackermann et al.,
2013). In this scenario, the expanding diffusive shocks accelerate cosmic rays
(CRs) to very high energy (VHE). The electrons and nuclei are concomitant and
can both be accelerated to VHE simultaneously. This means that the nuclei and
electrons should have some common origin (Yuan & Bi, 2013). The combined study
for such multi-messenger topic is important to unveil the enigma of the
origins of the CRs. And the typical properties are required to be observed to
support this view of point.

The measurements of CRs are stepping into a high precision era with the new
generation of space-borne and ground-based experiments. A series of new
phenomena are revealed by those precise measurements. Firstly, a fine
structure of spectral hardening at 200 GeV for nuclei was discovered by the
ATIC-2, the CREAM and the PAMELA experiments (Panov et al., 2007, 2009; Ahn et
al., 2010; Yoon et al., 2017; Adriani et al., 2011). Lately, the AMS-02
experiment confirmed it with unprecedented precision (Aguilar et al., 2015b,
a). More interestingly, the spectral break-off around $\sim$14 TeV was
observed by the CREAM, the NUCLEON and the DAMPE experiments (Yoon et al.,
2017; Atkin et al., 2017, 2018; An et al., 2019; Alemanno et al., 2021). Three
categories of models such as the local sources model (Sveshnikova et al.,
2013; Liu et al., 2017, 2019; Vladimirov et al., 2012), the combined effects
from different group sources model (Malkov & Moskalenko, 2020) and the
spatial-dependent propagation (SDP) model (Guo et al., 2016; Jin et al., 2016;
Guo & Yuan, 2018; Liu et al., 2018) are proposed to explain the these new
structures. Similar to that of the nuclei, the electrons should also exist
such a component of spectral break-off above $\sim$100 GeV. The recent precise
spectrum measurement shows sharp drop-off at 284 GeV for positrons, but no
obvious change until TeV for the total spectrum of electrons and positrons
(Aguilar et al., 2019). The deficit of positrons above 284 GeV requires the
compensation of primary electrons. This means that the excess of electrons
exists in a way similar to the fine structure of nuclei, which is possibly
accelerated by the SNRs scenario (Mertsch, 2011; Bernard et al., 2012;
Serpico, 2012; Di Mauro et al., 2014; Fang et al., 2018b; Tang & Piran, 2019).

What more important is the anisotropy evolution with energy. Though the
arrival directions of GCRs are highly isotropic due to their diffusive
propagation in the galactic magnetic field, a weak dipole-like anisotropy is
consistently observed, with difference in intensity of up to $\sim
10^{-4}-10^{-3}$. So far, a large amount of observations with anisotropies
ranging from TeV to PeV have been carried out by the ground-based experiments,
for example, the Super-Kamiokande (Guillian et al., 2007), the Tibet
(Amenomori et al., 2005, 2006, 2010, 2017), the Milagro (Abdo et al., 2008,
2009), the IceCube/Ice-Top (Abbasi et al., 2010, 2011, 2012; Aartsen et al.,
2013, 2016), the ARGO-YBJ (Bartoli et al., 2013, 2015) and the HAWC
(Abeysekara et al., 2014). It is clear that the phase of the anisotropy below
$100$ TeV energy roughly directs to the galactic anticenter, which is totally
paradoxical with the conventional propagation model (CPM). However, above
$100$ TeV, the phase gradually turns to the direction of the galactic center,
until about 100 PeV energies. This is consistent with the expectation of the
CPM. Correspondingly, the amplitude has experienced similar transition at the
critical energy of $100$ TeV. In addition, what most important is that there
exists a common transition energy scale between the structures of the energy
spectra and the anisotropy. The local sources possibly play a very important
role in resolving the conjunct problems of the spectra and the anisotropy (Liu
et al., 2019; Qiao et al., 2019).

Based on the above discussion, the local sources are required to reproduce the
multi-messenger anomaly for the spectra of protons and electrons and the
nuclear anisotropy. In our recent work, we propose a local source under the
SDP model to reproduce the co-evolution of the spectra and the anisotropy. We
find that the Geminga SNR at its birth place could be a preferred candidate
(Liu et al., 2019; Qiao et al., 2019; Zhang et al., 2021). One natural
question follows is that there exist dozens of SNRs near the solar system
within 1 kiloparsec (kpc), as shown in Figure.1, then, how is about the
contribution from the other SNRs? In this work, we systematically study the
contribution of all the local SNRs within 1 kpc around the solar system to the
spectra of nuclei and electrons, the detailed parameters of these local SNRs
are shown in Table 1. The paper is organized as follows, Section 2 describes
the model description, Section 3 presents the calculated results, and Section
4 gives the conclusion.

Fig. 1: Face-on view of the galaxy showing positions of the local SNRs within
$1$ kpc from the solar system. The blue solid squares mark the SNRs. The
solide circle at the center represents the location of the solar system. The
solid red, green and pink pentagrams represent the Geminga, the Monogem, and
the Vela SNRs, respectively, and the hollow ones represent the homologous
pulsars. The gray arrow points to the galactic center and the brown arrow
indicates the direction of CRs anisotropy below $\sim 100$ TeV.

$Number$ | SNR | Other name |  Distance |  Radio index $\alpha_{r}$ |  Age | Pulsar |  References
---|---|---|---|---|---|---|---
| G+long+lat |  | [kpc] |  | [kyr] |  |
1 | G194.3-13.1 | Geminga | 0.33 |  | 345 | J0633+1746 | (1)
2 | G65.3+5.7 |  | $0.9\pm 0.1$ | $0.58\pm 0.07$ | $26\pm 1$ |  | (2-5)
3 | G65.7+1.2 | DA 495 | $1.0\pm 0.4$ | $0.45\pm 0.1$ | $16.75\pm 3.25$ | unkown | (6,7)
4 | G74.0-8.5 | Cygnus Loop | $0.54_{-0.08}^{+0.10}$ | $0.4\pm 0.06$ | $10\pm 1$ | J1952+3252 | (5,8,9)
5 | G114.3+0.3 |  | 0.7 | $0.49\pm 0.25$ | $7.7\pm 0.1$ |  | (5,6,10,11)
6 | G127.1+0.5 | R5 | $1\pm 0.1$ | $0.43\pm 0.1$ | $25\pm 5$ |  | (5,6,11-13)
7 | G156.2+5.7 |  | $0.8\pm 0.5$ | $2.0_{-0.7}^{+1.1}$ | $10\pm 1$ | B0450+55 | (5,6,14-17)
8 | G160.9+2.6 | HB 9 | $0.8\pm 0.4$ | $0.48\pm 0.03$ | $5.5\pm 1.5$ | B0458+46 | (5,6,11,18,19)
9 | G203.0+12.0 | Monogem Ring | $0.288_{-0.027}^{+0.033}$ |  | $86\pm 1$ | B0656+14 | (20-21)
10 | G263.9-3.3 | Vela(XYZ) | $0.295\pm 0.075$ | variable | $11.2\pm 0.1$ | B0833-45 | (5,22-25)
11 | G266.2-1.2 | Vela Jr | $0.75\pm 0.01$ |  | $3.5\pm 0.8$ | J0855-4644 | (5,26-29)
12 | G347.3-0.5 | SN393 | $1\pm 0.3$ |  | 4.9 |  | (5,30-31)

1(Smith et al., 1994);2(Gorham et al., 1996);3(Mavromatakis et al.,
2002);4(Xiao et al., 2009);5(Green, 2014);6(Kothes et al., 2006);7(Kothes et
al., 2008);8(Blair et al., 2005);9(Sun et al., 2006);10(Yar-Uyaniker et al.,
2004);11(Reich et al., 2014);12(Joncas et al., 1989);13(Leahy & Tian,
2006);14(Reich et al., 1992);15(Yamauchi et al., 2000);16(Xu et al.,
2007);17(Katsuda et al., 2009);18(Reich et al., 2003);19(Leahy & Tian,
2007);20(Plucinsky et al., 1996);21(Brisken et al., 2003);22(Cha et al.,
1999);23(Alvarez et al., 2001);24(Caraveo et al., 2001);25(Miceli et al.,
2008);26(Aschenbach, 1998);27(Iyudin et al., 1998);28(Redman & Meaburn,
2005);29(Katsuda et al., 2008);30(Lazendic et al., 2004);31(Morlino et al.,
2009).

Tab. 1:  Characteristic parameters of the nearby SNRs within 1 kpc

##  2 Model Description

After entering the interstellar space, CRs undergo the random walks within the
galactic magnetic field by bouncing off the magnetic waves and
magnetohydrodynamic turbulence. The diffusive region, which is called the
magnetic halo, is approximated as a flat cylinder with radius of $R=20$ kpc,
equivalent to the galactic radius. The half-thickness $z_{h}$ is unknown, and
is typically constrained by fitting the B/C ratio. The galactic disk, where
both the CR sources and the interstellar gas are mainly spread across, is
located in the middle of the magnetic halo. The width of the galactic disk is
approximated to be invariant spatially and equals to $200$ pc. In addition to
diffusion, CRs may also go through convection, diffusive-reacceleration,
fragmentation, radioactive decay and other energy losses before arriving at
the solar system. In fact, the process of convection is ignored in this work.
The CR nuclei lose their energy principally via ionization, Coulomb scattering
and adiabatic expansion. For electrons and positrons, their major energy loss
mechanisms are the bremsstrahlung, synchrotron radiation and the inverse
Compton scattering. This comprehensive process can be described by the
propagation equation as,

$\begin{array}[]{lcll}\frac{\partial\psi(\boldsymbol{r},p,t)}{\partial
t}&=&q(\boldsymbol{r},p,t)+\boldsymbol{\nabla}\cdot\left(D_{xx}\boldsymbol{\nabla}\psi-\boldsymbol{V_{c}}\psi\right)\\\
&+&\frac{\partial}{\partial p}p^{2}D_{pp}\frac{\partial}{\partial
p}\frac{1}{p^{2}}\psi-\frac{\partial}{\partial
p}\left[\dot{p}\psi-\frac{p}{3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{V_{c}}\psi\right)\right]\\\
&-&\frac{\psi}{\tau_{f}}-\frac{\psi}{\tau_{r}}\end{array}~{},$ (1)

where $\rm q(\boldsymbol{r},p,t)$ is the acceleration sources,
$\rm\psi(\boldsymbol{r},p,t)$ the density of CR particles per unit momentum
$\rm p$ at position $\rm\boldsymbol{r}$, $\rm\boldsymbol{V_{c}}$ the
convection velocity, $\rm\dot{p}\equiv dp/dt$ the momentum loss rate,
$\rm\tau_{f}$ and $\rm\tau_{r}$ the characteristic time scales for
fragmentation and radioactive decay, respectively, and $\rm D_{xx}$ and $\rm
D_{pp}$ the diffusion coefficients in the coordinate and the momentum space,
respectively.

###  2.1 Spatially-dependent propagation

The SDP of CRs has received a lot of attention in recent years. It is first
introduced as a Two Halo model (THM) (Tomassetti, 2012) to explain the
spectral hardening of both protons and helium above $200$ GeV (Adriani et al.,
2011). Afterwards, it is further applied to the secondary and heavier
components (Tomassetti, 2015; Feng et al., 2016; Guo et al., 2016; Liu et al.,
2018; Tian et al., 2020; Yuan et al., 2020a), diffuse gamma-ray distribution
(Guo & Yuan, 2018) and large-scale anisotropy (Liu et al., 2019; Qiao et al.,
2019). For a comprehensive introduction, one can refer to (Guo et al., 2016)
and (Liu et al., 2018).

In the SDP model, the whole diffusive halo is divided into two parts. The
galactic disk and its surrounding area are called the inner halo (IH) region,
in which the diffusion coefficient is spatially-dependent and relevant to the
radial distribution of background CR sources. The extensive diffusive region
outside the IH is named as the outer halo (OH) region, where the diffusion is
regarded as only rigidity dependent. The spatially dependent diffusion
coefficient $D_{xx}$ is thus parameterized as,

$D_{xx}(r,z,{\cal R})=D_{0}F(r,z)\beta^{\eta}\left(\dfrac{\cal R}{{\cal
R}_{0}^{{}^{\prime}}}\right)^{\delta(r,z)}~{},$ (2)

where $r$ and $z$ are cylindrical coordinate, ${\cal R}$ the particle’s
rigidity, $\beta$ the particle’s velocity in unit of light speed, and $D_{0}$
and $\eta$ are constants. For the parameterization of $F(r,z)$ and
$\delta(r,z)$, one can refer to (Tian et al., 2020). The total half-thickness
of the propagation halo is $z_{h}$, and the half-thickness of the inner halo
is $\xi z_{h}$.

In this work, we adopt the common diffusion-reacceleration (DR) model, with
the diffusive-reacceleration coefficient $D_{pp}$ coupled to $D_{xx}$ by
$D_{pp}D_{xx}=\dfrac{4p^{2}v_{A}^{2}}{3\delta(4-\delta^{2})(4-\delta)}$, where
$v_{A}$ is the Alfvén velocity, $p$ is the momentum, and $\delta$ the rigidity
dependence slope of the diffusion coefficient (Seo & Ptuskin, 1994). The
numerical package DRAGON is used to solve the SDP equation to obtain the
distribution of the CR positrons, electrons and protons. Lower than tens of
GeV, the CR fluxes are impacted by the solar modulation. The well-known force-
field approximation (Gleeson & Axford, 1968; Perko, 1987) is applied to
describe such an effect, with a modulation potential $\phi$ adjusted to fit
the low energy data.

###  2.2 Background sources

The SNRs are regarded as the most likely sites for the acceleration of GCRs by
default, in which the charge particles are accelerated to a power-law
distribution through the diffusive shock acceleration. The distribution of the
SNRs are approximated as axisymmetric, which is usually parameterized as,

$f(r,z)=\left(\dfrac{r}{r_{\odot}}\right)^{\alpha}\exp\left[-\dfrac{\beta(r-r_{\odot})}{r_{\odot}}\right]\exp\left(-\dfrac{|z|}{z_{s}}\right)~{},$
(3)

where $r_{\odot}\equiv 8.5$ kpc represents the distance from the solar system
to the galactic center. The parameters $\alpha$ and $\beta$ are taken as
$1.69$ and $3.33$, respectively, in this work (Green, 2015). The density
distribution of the SNRs decreases exponentially along the vertical height
from the galactic plane, with $z_{s}=200$ pc.

Past studies, for example, (Trotta et al., 2011; Boschini et al., 2017; Yuan
et al., 2020b), have pointed out that the assumption of a single power-law
injection spectrum could not reproduce the observations very well, especially
in the diffusion-reacceleration model. Furthermore, this low-energy break has
also been indirectly supported by the observations of nearby molecular clouds,
see A. Neronov et al. (Neronov et al., 2012). Therefore, we introduce a break
at low energy of the injection spectrum,

$q({\cal R})=q_{0}\left\\{\begin{array}[]{lll}\left(\dfrac{{\cal R}}{{\cal
R}_{\rm br1}}\right)^{\nu_{1}},~{}{{\cal R}\leq{\cal R}_{\rm br1}}\\\ \\\
\left(\dfrac{{\cal R}}{{\cal R}_{\rm
br1}}\right)^{\nu_{2}}\exp\left[-\dfrac{\cal R}{{\cal R}_{\rm
c}}\right],~{}{{\cal R}>{\cal R}_{\rm br1}}\end{array}~{},\right.$ (4)

where $q_{0}$ is the normalization factor, $\nu_{1,2}$ the spectral incides,
${\cal R}_{\rm br1}$ break rigidities, and ${\cal R}_{\rm c}$ the cutoff
rigidity.

###  2.3 Local SNRs

At $\sim$ TeV energies, CR electrons are found to originate from local sources
within $\sim 1$ kpc around the solar system (Yuan & Feng, 2018). In this small
region, the hypothesis of continuous distribution may not be valid any more.
Studies show that the discrete effect of the nearby CR sources could induce
large fluctuations, especially at high energies (Mertsch, 2011; Bernard et
al., 2012; Fang et al., 2017). The contribution of the nearby sources to the
CR nuclei and electrons have been studied in the past works (see e.g.,
(Serpico, 2012; Di Mauro et al., 2014; Liu et al., 2017; Fang et al., 2018a)).
In this work, we assume that a nearby SNR accounts for the excess of electrons
(protons) above $\sim 100$ (or 200) GeV, respectively. The propagation of CRs
injected instantaneously from a point source is described by a time-dependent
propagation equation (Atoyan et al., 1995). The injection rate as a function
of time and rigidity is assumed to be,

$Q^{\rm snr}({\cal R},t)=Q^{\rm snr}_{0}(t)\left(\dfrac{\cal R}{{\cal
R}_{0}}\right)^{-\gamma}\exp\left[-\dfrac{\cal R}{{\cal R}^{\rm e^{-},p}_{\rm
c}}\right]~{},$ (5)

where ${\cal R}^{\rm e^{-},p}_{\rm c}$ is the cutoff rigidity of its
accelerated CRs.

The local SNRs may accelerate primary nuclei and electrons during their early
evolution stage. This contribution of primary electrons from local sources may
be necessary given the different spectral behaviors of the positrons and the
electrons (Zhang et al., 2021). The injection process of the SNRs is
approximated as burst-like. The source injection rate is assumed as the
following,

$Q_{0}^{\rm snr}(t)=q_{0}^{\rm snr}\delta(t-t_{0})~{},\\\ $ (6)

where $t_{0}$ is the time of the supernova explosion. The propagated spectrum
from the local SNRs is thus a convolution of the Green’s function and the
time-dependent injection rate $Q_{0}(t)$ (Atoyan et al., 1995),

$\varphi(\boldsymbol{r},{\cal
R},t)=\int_{t_{i}}^{t}G(\boldsymbol{r}-\boldsymbol{r}^{\prime},t-t^{\prime},{\cal
R})Q_{0}(t^{\prime})dt^{\prime}.$ (7)

##  3 Results

In this work, both the propagation and the injection parameters are fitted
manually. By fitting the B/C ratio, the diffusion coefficient parameters are
$D_{0}=4.87\times 10^{28}$ cm2 s-1, $\delta_{0}=0.55$, $N_{m}=0.6$, $n=4$,
$\xi=0.1$, and $\eta=0.05$. The half thickness of the propagation halo is
$z_{h}=5$ kpc, and the Alfvénic velocity is $v_{A}=6$ km s-1. The comparison
of the B/C ratio between the model prediction and the observation data is
given in Figure 2, which indicates that the relevant parameters are reasonably
matched. In order to determine which nearby SNRs are the appropriate
candidates on the basis of their contributions to the proton and electron
spectra, respectively, we carry out more detailed analyses in the following.

Fig. 2: Model prediction of the B/C ratio compared with the AMS-02 measurement
(Aguilar et al., 2015b, 2017).

###  3.1 Spectra of all the local SNRs

Figure 3 presents the calculated proton and electron spectra for all the SNRs
as listed in Table 1. Here, it can be seen that only the Geminga, the Monogem,
and the Vela SNRs have significant contributions to the proton and electron
spectra, respectively. Whereas most others have very limited contribution to
the spectra. To appear in the same figure, we multiply a scale factor in the
fluxes for the faint SNRs to show properly, as shown by the dashed line in
Figure 3. The detailed parameters of the local sources are listed in Table 2.
In the following subsection, we give the detailed studies on the contribution
to the spectra and anisotropy by the Geminga, the Monogem, and the Vela SNRs.
By doing so, we can expect to pinpoint the optimal local SNR.

Fig. 3: The calculated energy spectra of protons (left) and electrons (right)
from the different local SNRs. The red, blue and purple solid lines represent
the energy spectra of the Geminga, the Monogem, and the Vela SNRs,
respectively. The dashed lines with different colors represent the energy
spectra of other local SNRs, respectively. For comparing the results
conveniently, the fluxes of some local SNRs are multiplied by a scale factor.
The red-square, green-triangle and magenta-square data points are from the
AMS-02, the DAMPE, and the CREAM experiments, respectively (Aguilar et al.,
2015b, a; Yoon et al., 2017; An et al., 2019).

###  3.2 Spectra and anisotropy of the Geminga SNR

Figure 4 illustrates the calculated energy spectra and anisotropy of protons
and primary electrons, respectively. The blue solid lines are the fluxes from
the background sources. Taking into account the contributions of the Geminga
SNR as shown in Figure 3, the model predictions of the proton and primary
electron spectra are in good agreement with the experimental data, as shown in
the upper panel of Figure 4. Table 2 shows the relevant background and local
sources parameters. In fact, the anisotropy is a more effective factor to
determine the dominant local SNRs from their direction, age and distance. The
bottom panel of Figure 4 shows the model calculated anisotropy for CRs and
electrons, respectively. It is obvious that the model calculated anisotropy is
consistent with that of the observations of CRs. And the electron anisotropy
is roughly at the same level and far below the Fermi-LAT limit. Therefore, the
results support our previous conclusion of the co-evolved spectrum and
anisotropy for the Geminga SNR (Liu et al., 2019; Qiao et al., 2019; Zhang et
al., 2021).

Fig. 4: The calculated energy spectra (top) and anisotropy (bottom) for
protons (left) and primary electrons (right) from the Geminga SNR,
respectively. For the energy spectra, the blue lines represent the fluxes from
the background sources, the red lines the fluxes from the local Geminga SNR,
and the black solid lines their corresponding sum. For the anisotropies, the
black solid lines are the model expected results.

Bckgroud | Normalization† |  $\nu 1$ |  $\mathcal{R}_{br1}$ |  $\nu 2$ |  $\mathcal{R}_{c}$ |
---|---|---|---|---|---|---
| $[({\rm m}^{2}\cdot{\rm sr}\cdot{\rm s}\cdot{\rm GeV})^{-1}]$ |  | [GV] |  | [GV] |
P | $4.62\times 10^{-2}$ | 2 | $6$ | 2.41 | $7\times 10^{6}$ |
Electron | $2.80\times 10^{-1}$ | 1.50 | 5.10 | 2.80 | 300 |
SNR | $r_{snr}$ | $t_{snr}$ | $\gamma_{snr}$ | $\mathcal{R}_{c}$ | $q_{0}^{prot}$ | $q_{0}^{elec}$
| [kpc] | [yrs] |  | [TV] | $[\rm GeV^{-1}]$ | $[\rm GeV^{-1}]$
Geminga | 0.33 | $3.45\times 10^{5}$ | 2.13 | 55 | $2.50\times 10^{52}$ | $1.20\times 10^{50}$
Monogem | 0.29 | $8.60\times 10^{4}$ | 2.13 | 27 | $1.20\times 10^{52}$ | $2.10\times 10^{49}$
Others | … | … | $2\alpha_{r}+1$ | 55 | $2.50\times 10^{52}$ | $1.20\times 10^{50}$

†The proton and electron normalization are set at kinetic energy per nucleon
$E_{k}=100$ GeV/n and $10$ GeV/n, respectively.

Tab. 2: The background and local SNRs injected parameters for protons and
electrons. For the age and distance parameters of the local sources, refer to
Table 1

###  3.3 Spectra and anisotropy of the Monogem SNR

Similar to the results of the Geminga SNR, the model calculations of the
Monogem SNR are presented in Figure 5. The proton spectrum is also consistent
with the observations. However, the calculated anisotropy greatly overthrows
the observations. Obviously, there exists tension between the spectrum and the
anisotropy. In addition, the cut-off energy of the electron spectrum is also
beyond the TeV energy.

It should be noticed that the spin-down energy of the Monogem pulsar is lower
than that of the Geminga pulsar with the value of $1.8\times
10^{48}/1.25\times 10^{49}erg$. It is possible that the injection power of the
Monogem SNR is much lower than that of the Geminga SNR. In fact, the
$\gamma$-ray emission at TeV energy is much lower for the Monogem compared
with that of the Geminga (Abeysekara et al., 2017). In Figure 5, the dashed
line shows the corrected injection energy resulted from the spin-down energy.
Under these circumstances, the contribution from the Monogem SNR can be
ignored.

Another reason is the propagation effect in the source region. It is well
known that the propagation coefficient in the source region is much lower
(Abeysekara et al., 2017). This means that the CRs spend more time to
disentangle from the source. The age of Monogem is at the level of 86 kyrs and
only a limited part of the CRs from Monogem arrives at our solar system.

On the other hand, the regular magnetic field could also regulate the CR flux
from the local source. The observations (Schwadron et al., 2014; McComas et
al., 2009; Zirnstein et al., 2016) have revealed that the phase of anisotropy
less than 100 TeV is consistent with the effect of the local regular magnetic
field. Ahlers (Ahlers, 2014) has pointed out that the local regular magnetic
field and the corresponding anisotropic diffusion have to be considered when
accounting for the evolution relation of the amplitude and phase of the
anisotropy with energy. When the regular magnetic field is taken into account,
the CR flux is projected into the regular magnetic field line and propagates
along the magnetic field. If the sources are far away from the regular
magnetic field, their flux to the solar system could be suppressed by the
regular magnetic field. This would be investigated in detail in the coming
work.

In a word, the contribution from the Monogem SNR is possibly limited. Other
information is still required to check our model further on this point.

Fig. 5: The calculated energy spectra (top) and anisotropy (bottom) for
protons (left) and primary electrons (right) from the Monogem SNR,
respectively. The implication of different color solid lines is the same as
that in Figure 4 and the different colored dashed lines display the corrected
injection energy in accordance with the spin-down energy of the Geminga and
the Monogem pulsars.

###  3.4 Spectra and anisotropy of the Vela SNR

The Vela SNR is a very young and nearest celestial body. It can provide very
important clue to check our model in terms of the CR spectrum and anisotropy.
Figure 6 shows the expected spectrum and anisotropy for electrons. The related
parameters are shown in Table 3. It is obvious that there exists a bump
structure in the spectrum and sharp rising in the anisotropy above the TeV
energy. In fact, the bump structure in the spectrum has been hinted by the
HESS, the VERTITAS, the DAMPE, and the CALET observations (DAMPE Collaboration
et al., 2017; Motz et al., 2021; Archer et al., 2018; Aharonian et al., 2008),
as shown in Figure 7. In addition, the electron anisotropy reveals an
remarkable hoist, which reaches $10\%$ at several TeV. What makes the fact
more convincible is that those abnormal phenomena can be tested soon. Firstly,
the space-borne experiments DAMPE and HERD will be observing the new spectral
structure in the short run. Simultaneously, the anisotropy is overwhelmingly
larger for the electrons than for the CRs, which makes it easily observed by
the ground-based experiments HAWC and LHAASO-WCDA.

Considering the importance of the Vela SNR, it is necessary to perform further
study such as on the effect of different ages. Figure 8 shows the spectra and
anisotropy at 10, 100, 300, and 500 kyrs for protons and electrons,
respectively. It is obvious that the flux reaches its maximum at 100 kyrs and
then decreases smoothly. When the source age becomes as old as 500 kyrs, the
contribution can be ignored. In a word, the Vela SNR plays the similar role of
Monogem at $\sim$ 100 kyrs, and that of Geminga at $\sim$ 300 kyrs, and
disappears above 500 kyrs.

Bckgroud | Normalization† | $\nu 1$ | $\mathcal{R}_{br1}$ | $\nu 2$ | $R_{c}$ |
---|---|---|---|---|---|---
(DAMPE) | $[({\rm m}^{2}\cdot{\rm sr}\cdot{\rm s}\cdot{\rm GeV})^{-1}]$ |  | [GV] |  | [GV] |
Electron | $2.55\times 10^{-1}$ | 1.60 | 5.10 | 2.72 | 1200 |
SNR | $r_{snr}$ | $t_{snr}$ | $\gamma_{snr}$ | $\mathcal{R}_{c}$ | $q_{0}^{elec}$ |
| [kpc] | [yrs] |  | [TV] | $[\rm GeV^{-1}]$ |
Geminga | 0.33 | $3.45\times 10^{5}$ | 2.13 | 55 | $2.30\times 10^{50}$ |
Vela | 0.30 | $1.12\times 10^{4}$ | 2.10 | 55 | $5.00\times 10^{49}$ |

†The normalization is set at kinetic energy per nucleon $E_{k}=10$ GeV/n.

Tab. 3: Electron background and local SNRs injected parameters about the DAMPE
experiment Fig. 6: Left: The calculated total positron and electron spectra of
the Vela SNR. The blue solid line is the total background flux. The yellow
solid line with shadow indicates the Geminga SNR flux. The light-red dashed
line with shadow represents the flux from the Vela SNR for different injected
energies. The corresponding solid line is the sum of the background, the
Geminga, and the Vela fluxes. The green squares are the measurements from the
DAMPE experiment (DAMPE Collaboration et al., 2017). Right: Anisotropies of
electrons with the Geminga and the Vela SNRs under the SDP model. Both the
black and red markers are the upper limits set by the Fermi-LAT experiment
(Abdollahi et al., 2017).  Fig. 7: Our model expectations of the total
electron and positron spectra compared with the experiments of DAMPE (left-
top), CALET (right-top), VERITAS (left-bottom) and HESS (right-bottom),
respectively. The blue solid lines represent the fluxes from the background
sources, the yellow lines are the fluxes from the Geminga SNR, the red lines
with shadow are the fluxes from the Vela SNR, and the black lines indicate the
corresponding sum of the above components.  Fig. 8: The calculated energy
spectra (top) and anisotropy (bottom) for protons (left) and primary electrons
(right) from the Vela SNR, respectively. For the energy spectra, the blue
lines are the fluxes from the background sources, the different colored dashed
lines are the fluxes from the Vela SNR for various ages, and the different
solid lines are their corresponding sum including the Geminga SNR flux. Note
that the Geminga SNR spectrum is not shown here. For the anisotropies, the
different colored solid lines are the model expected results.

##  4 Conclusion

With the operation of new generations of space-borne and ground-based CR
experiments, the measurements of CRs have stepped into a precise era. More and
more fine structures have been unveiled in the CR energy spectra and
anisotropy from hundreds of GeV to hundreds of TeV. To reproduce these
observational phenomena, the local sources have to be involved in the model
calculation.

In this work, we systematically study the contribution of all the local SNRs
within 1 kpc around the solar system to the spectra and anisotropy of CR
nuclei and electrons. We demonstrate that three local SNRs, i.e. Geminga,
Monogem, and Vela, could have important contribution to both proton and
electron spectra, as inferred from their distances and ages, assuming a common
injection energy. However, the expected anisotropy from Monogem is obviously
inconsistent with the observations. One of possibilities is its injection
power is lower, considering that the spin-down energy of the Monogem pulsar is
lower about an order of magnitude than that of the Geminga pulsar. When the
total injection energy lowered by one order of magnitude, the influence from
Monogem can be safely neglected. That leaves the Gminga SNR as the dominant
source responsible for the observational anomalies of the nuclear spectra and
the anisotropy.

Furthermore, we expected that there is a new bump in electron spectrum above
several TeV, which could stem from the young Vela SNR. Simultaneously, the
magnitude of anisotropy can reach as high as 10% at the corresponding energy.
Even the injected energy reuduces by an order of magnitude, the excesses from
Vela SNR in the spectrum and the anisotropy are still observable. We hope such
excesses are going to be checked soon by the space-borne DMAPE and HERD and
the ground-based HAWC and LHAASO experiments.

One has to be noted that the CR diffusion adopted in the above calculation is
approximated as isotropic, in which corresponds to the case of an ideal
turbulent magnetic field in plasma. In fact, it is well-known that there
exists a large-scale regular magnetic field in Galaxy. When the regular
component is taken into account, the diffusion coefficient is converted to a
tensor, i.e. the diffusion coefficient parallel and perpendicular to the
regular magnetic field have to be considered respectively(Cerri et al., 2017;
Liu et al., 2020). In the coming work, we would like to further study the
influence of the regular magnetic field on the energy spectra and
anisotropies.

## Acknowledgements

This work is supported by the National Key R&D Program of China under Grant
No. 2018YFA0404202, and the National Natural Science Foundation of China (Nos.
11635011, 11875264, U1831129, U1831208, and U2031110).

Note added in proof. While this paper was ready to submit, another similar
work is noticed to have appeared in https://arxiv.org/(Tang et al., 2021).

## References

  * Aartsen et al. (2013) Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, ApJ, 765, 55
  * Aartsen et al. (2016) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016, ApJ, 826, 220
  * Abbasi et al. (2010) Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2010, ApJ, 718, L194
  * Abbasi et al. (2011) —. 2011, ApJ, 740, 16
  * Abbasi et al. (2012) —. 2012, ApJ, 746, 33
  * Abdo et al. (2008) Abdo, A. A., Allen, B., Aune, T., et al. 2008, Physical Review Letters, 101, 221101
  * Abdo et al. (2009) Abdo, A. A., Allen, B. T., Aune, T., et al. 2009, ApJ, 698, 2121
  * Abdollahi et al. (2017) Abdollahi, S., Ackermann, M., Ajello, M., et al. 2017, Physical Review Letters, 118, 091103
  * Abeysekara et al. (2014) Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2014, ApJ, 796, 108
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358, 911
  * Ackermann et al. (2013) Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807
  * Adriani et al. (2011) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, Science, 332, 69
  * Aguilar et al. (2015a) Aguilar, M., Aisa, D., Alpat, B., et al. 2015a, Physical Review Letters, 115, 211101
  * Aguilar et al. (2015b) —. 2015b, Physical Review Letters, 114, 171103
  * Aguilar et al. (2017) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2017, Physical Review Letters, 119, 251101
  * Aguilar et al. (2019) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2019, Phys. Rev. Lett., 122, 041102
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, Physical Review Letters, 101, 261104
  * Ahlers (2014) Ahlers, M. 2014, Phys. Rev. Lett., 112, 021101
  * Ahn et al. (2010) Ahn, H. S., Allison, P., Bagliesi, M. G., et al. 2010, ApJ, 714, L89
  * Alemanno et al. (2021) Alemanno, F., An, Q., Azzarello, P., et al. 2021, Phys. Rev. Lett., 126, 201102
  * Alvarez et al. (2001) Alvarez, H., Aparici, J., May, J., & Reich, P. 2001, A&A, 372, 636
  * Amenomori et al. (2005) Amenomori, M., Ayabe, S., Cui, S. W., et al. 2005, ApJ, 626, L29
  * Amenomori et al. (2006) Amenomori, M., Ayabe, S., Bi, X. J., et al. 2006, Science, 314, 439
  * Amenomori et al. (2010) Amenomori, M., Bi, X. J., Chen, D., et al. 2010, ApJ, 711, 119
  * Amenomori et al. (2017) —. 2017, ApJ, 836, 153
  * An et al. (2019) An, Q., Asfandiyarov, R., Azzarello, P., et al. 2019, Science Advances, 5, eaax3793
  * Archer et al. (2018) Archer, A., Benbow, W., Bird, R., et al. 2018, Phys. Rev. D, 98, 062004
  * Aschenbach (1998) Aschenbach, B. 1998, Nature, 396, 141
  * Atkin et al. (2017) Atkin, E., Bulatov, V., Dorokhov, V., et al. 2017, J. Cosmology Astropart. Phys, 7, 020
  * Atkin et al. (2018) —. 2018, Soviet Journal of Experimental and Theoretical Physics Letters, 108, 5
  * Atoyan et al. (1995) Atoyan, A. M., Aharonian, F. A., & Völk, H. J. 1995, Phys. Rev. D, 52, 3265
  * Bartoli et al. (2013) Bartoli, B., Bernardini, P., Bi, X. J., et al. 2013, Phys. Rev. D, 88, 082001
  * Bartoli et al. (2015) —. 2015, ApJ, 809, 90
  * Bernard et al. (2012) Bernard, G., Delahaye, T., Salati, P., & Taillet, R. 2012, A&A, 544, A92
  * Blair et al. (2005) Blair, W. P., Sankrit, R., & Raymond, J. C. 2005, AJ, 129, 2268
  * Blasi (2013) Blasi, P. 2013, A&A Rev., 21, 70
  * Boschini et al. (2017) Boschini, M. J., Della Torre, S., Gervasi, M., et al. 2017, ApJ, 840, 115
  * Boulares (1989) Boulares, A. 1989, ApJ, 342, 807
  * Brisken et al. (2003) Brisken, W. F., Thorsett, S. E., Golden, A., & Goss, W. M. 2003, ApJ, 593, L89
  * Caraveo et al. (2001) Caraveo, P. A., De Luca, A., Mignani, R. P., & Bignami, G. F. 2001, ApJ, 561, 930
  * Cerri et al. (2017) Cerri, S. S., Gaggero, D., Vittino, A., Evoli, C., & Grasso, D. 2017, J. Cosmology Astropart. Phys, 2017, 019
  * Cha et al. (1999) Cha, A. N., Sembach, K. R., & Danks, A. C. 1999, ApJ, 515, L25
  * DAMPE Collaboration et al. (2017) DAMPE Collaboration, Ambrosi, G., An, Q., et al. 2017, Nature, 552, 63
  * Di Mauro et al. (2014) Di Mauro, M., Donato, F., Fornengo, N., Lineros, R., & Vittino, A. 2014, J. Cosmology Astropart. Phys, 4, 6
  * Fang et al. (2018a) Fang, K., Bi, X.-J., & Yin, P.-F. 2018a, ApJ, 854, 57
  * Fang et al. (2018b) Fang, K., Bi, X.-J., Yin, P.-F., & Yuan, Q. 2018b, ApJ, 863, 30
  * Fang et al. (2017) Fang, K., Wang, B.-B., Bi, X.-J., Lin, S.-J., & Yin, P.-F. 2017, ApJ, 836, 172
  * Feng et al. (2016) Feng, J., Tomassetti, N., & Oliva, A. 2016, Phys. Rev. D, 94, 123007
  * Gleeson & Axford (1968) Gleeson, L. J., & Axford, W. I. 1968, ApJ, 154, 1011
  * Gorham et al. (1996) Gorham, P. W., Ray, P. S., Anderson, S. B., Kulkarni, S. R., & Prince, T. A. 1996, ApJ, 458, 257
  * Green (2014) Green, D. A. 2014, Bulletin of the Astronomical Society of India, 42, 47
  * Green (2015) —. 2015, MNRAS, 454, 1517
  * Guillian et al. (2007) Guillian, G., Hosaka, J., Ishihara, K., et al. 2007, Phys. Rev. D, 75, 062003
  * Guo et al. (2016) Guo, Y.-Q., Tian, Z., & Jin, C. 2016, ApJ, 819, 54
  * Guo & Yuan (2018) Guo, Y.-Q., & Yuan, Q. 2018, Phys. Rev. D, 97, 063008
  * Iyudin et al. (1998) Iyudin, A. F., Schönfelder, V., Bennett, K., et al. 1998, Nature, 396, 142
  * Jin et al. (2016) Jin, C., Guo, Y.-Q., & Hu, H.-B. 2016, Chinese Physics C, 40, 015101
  * Joncas et al. (1989) Joncas, G., Roger, R. S., & Dewdney, P. E. 1989, A&A, 219, 303
  * Katsuda et al. (2009) Katsuda, S., Petre, R., Hwang, U., et al. 2009, PASJ, 61, S155
  * Katsuda et al. (2008) Katsuda, S., Tsunemi, H., & Mori, K. 2008, ApJ, 678, L35
  * Kothes et al. (2006) Kothes, R., Fedotov, K., Foster, T. J., & Uyanıker, B. 2006, A&A, 457, 1081
  * Kothes et al. (2008) Kothes, R., Landecker, T. L., Reich, W., Safi-Harb, S., & Arzoumanian, Z. 2008, ApJ, 687, 516
  * Lazendic et al. (2004) Lazendic, J. S., Slane, P. O., Gaensler, B. M., et al. 2004, ApJ, 602, 271
  * Leahy & Tian (2006) Leahy, D., & Tian, W. 2006, A&A, 451, 251
  * Leahy & Tian (2007) Leahy, D. A., & Tian, W. W. 2007, A&A, 461, 1013
  * Liu et al. (2017) Liu, W., Bi, X.-J., Lin, S.-J., Wang, B.-B., & Yin, P.-F. 2017, Phys. Rev. D, 96, 023006
  * Liu et al. (2019) Liu, W., Guo, Y.-Q., & Yuan, Q. 2019, J. Cosmology Astropart. Phys, 2019, 010
  * Liu et al. (2020) Liu, W., Lin, S.-j., Hu, H.-b., Guo, Y.-q., & Li, A.-f. 2020, ApJ, 892, 6
  * Liu et al. (2018) Liu, W., Yao, Y.-h., & Guo, Y.-Q. 2018, ApJ, 869, 176
  * Malkov & Moskalenko (2020) Malkov, M. A., & Moskalenko, I. V. 2020, arXiv e-prints, arXiv:2010.02826
  * Mavromatakis et al. (2002) Mavromatakis, F., Boumis, P., Papamastorakis, J., & Ventura, J. 2002, A&A, 388, 355
  * McComas et al. (2009) McComas, D. J., Allegrini, F., Bochsler, P., et al. 2009, Science, 326, 959
  * Mertsch (2011) Mertsch, P. 2011, J. Cosmology Astropart. Phys, 2, 31
  * Miceli et al. (2008) Miceli, M., Bocchino, F., & Reale, F. 2008, ApJ, 676, 1064
  * Morlino et al. (2009) Morlino, G., Amato, E., & Blasi, P. 2009, MNRAS, 392, 240
  * Motz et al. (2021) Motz, H., Link, J. T., Adriani, O., Akaike, Y., & De, N. G. A. 2021, in 37th International Cosmic Ray Conference
  * Neronov et al. (2012) Neronov, A., Semikoz, D. V., & Taylor, A. M. 2012, Phys. Rev. Lett., 108, 051105
  * Panov et al. (2007) Panov, A. D., Adams, J. H., J., Ahn, H. S., et al. 2007, Bulletin of the Russian Academy of Sciences, Physics, 71, 494
  * Panov et al. (2009) Panov, A. D., Adams, J. H., Ahn, H. S., et al. 2009, Bulletin of the Russian Academy of Sciences, Physics, 73, 564
  * Perko (1987) Perko, J. S. 1987, A&A, 184, 119
  * Plucinsky et al. (1996) Plucinsky, P. P., Snowden, S. L., Aschenbach, B., et al. 1996, ApJ, 463, 224
  * Qiao et al. (2019) Qiao, B.-Q., Liu, W., Guo, Y.-Q., & Yuan, Q. 2019, J. Cosmology Astropart. Phys, 2019, 007
  * Redman & Meaburn (2005) Redman, M. P., & Meaburn, J. 2005, MNRAS, 356, 969
  * Reich et al. (1992) Reich, W., Fuerst, E., & Arnal, E. M. 1992, A&A, 256, 214
  * Reich et al. (2014) Reich, W., Sun, X. H., Reich, P., et al. 2014, A&A, 561, A55
  * Reich et al. (2003) Reich, W., Zhang, X., & Fürst, E. 2003, A&A, 408, 961
  * Schwadron et al. (2014) Schwadron, N. A., Adams, F. C., Christian, E. R., et al. 2014, Science, 343, 988
  * Seo & Ptuskin (1994) Seo, E. S., & Ptuskin, V. S. 1994, ApJ, 431, 705
  * Serpico (2012) Serpico, P. D. 2012, Astroparticle Physics, 39, 2
  * Smith et al. (1994) Smith, V. V., Cunha, K., & Plez, B. 1994, A&A, 281, L41
  * Sun et al. (2006) Sun, X. H., Reich, W., Han, J. L., Reich, P., & Wielebinski, R. 2006, A&A, 447, 937
  * Sveshnikova et al. (2013) Sveshnikova, L. G., Strelnikova, O. N., & Ptuskin, V. S. 2013, Astroparticle Physics, 50, 33
  * Tang et al. (2021) Tang, T.-P., Xia, Z.-Q., Shen, Z.-Q., et al. 2021, Explanation of nearby SNRs for primary electron excess and proton spectral bump, arXiv:2109.12496
  * Tang & Piran (2019) Tang, X., & Piran, T. 2019, MNRAS, 484, 3491
  * Tian et al. (2020) Tian, Z., Liu, W., Yang, B., et al. 2020, Chinese Physics C, 44, 085102
  * Tomassetti (2012) Tomassetti, N. 2012, ApJ, 752, L13
  * Tomassetti (2015) —. 2015, Phys. Rev. D, 92, 081301
  * Trotta et al. (2011) Trotta, R., Jóhannesson, G., Moskalenko, I. V., et al. 2011, ApJ, 729, 106
  * Vladimirov et al. (2012) Vladimirov, A. E., Jóhannesson, G., Moskalenko, I. V., & Porter, T. A. 2012, ApJ, 752, 68
  * Xiao et al. (2009) Xiao, L., Reich, W., Fürst, E., & Han, J. L. 2009, A&A, 503, 827
  * Xu et al. (2007) Xu, J. W., Han, J. L., Sun, X. H., et al. 2007, A&A, 470, 969
  * Yamauchi et al. (2000) Yamauchi, S., Yokogawa, J., Tomida, H., Koyama, K., & Tamura, K. 2000, in Broad Band X-ray Spectra of Cosmic Sources, ed. K. Makishima, L. Piro, & T. Takahashi, 567
  * Yar-Uyaniker et al. (2004) Yar-Uyaniker, A., Uyaniker, B., & Kothes, R. 2004, ApJ, 616, 247
  * Yoon et al. (2017) Yoon, Y. S., Anderson, T., Barrau, A., et al. 2017, ApJ, 839, 5
  * Yuan & Bi (2013) Yuan, Q., & Bi, X.-J. 2013, Physics Letters B, 727, 1
  * Yuan & Feng (2018) Yuan, Q., & Feng, L. 2018, Science China Physics, Mechanics, and Astronomy, 61, 101002
  * Yuan et al. (2020a) Yuan, Q., Qiao, B.-Q., Guo, Y.-Q., Fan, Y.-Z., & Bi, X.-J. 2020a, arXiv e-prints, arXiv:2007.01768
  * Yuan et al. (2020b) Yuan, Q., Zhu, C.-R., Bi, X.-J., & Wei, D.-M. 2020b, J. Cosmology Astropart. Phys, 2020, 027
  * Zhang et al. (2021) Zhang, P.-p., Qiao, B.-q., Liu, W., et al. 2021, J. Cosmology Astropart. Phys, 2021, 012
  * Zirnstein et al. (2016) Zirnstein, E. J., Heerikhuisen, J., Funsten, H. O., et al. 2016, ApJ, 818, L18