# Resilient Fleet Management for Energy-Aware Intra-Factory Logistics

Mithun Goutham and Stephanie Stockar Mithun Goutham and Stephanie Stockar are
with the Department of Mechanical and Aerospace Engineering, The Ohio State
University, Columbus, OH 43210, USA<EMAIL_ADDRESS>
###### Abstract

This paper presents a novel fleet management strategy for battery-powered
robot fleets tasked with intra-factory logistics in an autonomous
manufacturing facility. In this environment, repetitive material handling
operations are subject to real-world uncertainties such as blocked passages,
and equipment or robot malfunctions. In such cases, centralized approaches
enhance resilience by immediately adjusting the task allocation between the
robots. To overcome the computational expense, a two-step methodology is
proposed where the nominal problem is solved a priori using a Monte Carlo Tree
Search algorithm for task allocation, resulting in a nominal search tree. When
a disruption occurs, the nominal search tree is rapidly updated a posteriori
with costs to the new problem while simultaneously generating feasible
solutions. Computational experiments prove the real-time capability of the
proposed algorithm for various scenarios and compare it with the case where
the search tree is not used and the decentralized approach that does not
attempt task reassignment.

##  I INTRODUCTION

Advancements in material handling have identified fleets of autonomous mobile
robots and automated guided vehicles as key components of autonomous
operations within flexible manufacturing systems (FMS) [1]. These battery-
powered robots handle the repetitive material handling tasks (MHT) integral to
manufacturing a product [2]. In this context, an energy-aware fleet management
policy aims to assign tasks and route robots in a manner that minimizes energy
expenses while adhering to the constraints imposed by battery charging
policies [3]. Furthermore, the operational requirements of an FMS introduce
additional constraints that ensure that routes visit pick up locations before
the corresponding delivery locations, and also respect the robot’s payload
limits.

This paper examines the scenario where a robot fleet, while actively executing
a nominal MHT based on a pre-determined fleet policy, experiences disruptions
due to real-world uncertainties. These perturbations in the definition of the
MHT arise from various factors such as machine failures, robot malfunctions,
battery degradation, fluctuating charge power or blocked passageways due to
fallen objects [4]. These disruptions can drastically affect the optimal
solution and the nominal policy may no longer be optimal or even feasible.
Resilience in this context refers to the ability of the fleet management
system to immediately adapt the fleet policy to guarantee uninterrupted
operations [5].

The $\mathcal{NP}$-hard nature of the task assignment problem has resulted in
the development of numerous decentralized approaches for fleet management that
focus on reassigning tasks only to the affected robots to rapidly recover from
a perturbation [6]. However, this approach typically leads to suboptimal
solutions as it overlooks the possible assistance that other unaffected AMRs
in the fleet could provide to enhance resilience [7]. In contrast, centralized
fleet management (CFM) considers all the available AMRs for task reassignment
and rerouting, thereby harnessing the fleet’s collective resilience to achieve
an optimal policy [7]. However, using exact methods for CFM is intractable for
real-time execution due to the $\mathcal{NP}$-hardness of the problem [8]. For
this reason, metaheuristic algorithms such as simulated annealing, genetic
algorithms, and tabu search are typically used to quickly improve multiple
trial solutions [9, 10]. However, these methods do not provide algorithmic
guarantees on the convergence or optimality of the resulting solution [11].
Another approach that enables real-time computation uses supervised machine
learning to map expert-identified problem perturbations to pre-computed
solutions [12, 13]. However, its performance is adversely affected when the
disruption differs significantly from the training data.

A gap in the CFM literature in the context of repetitive MHTs is the under-
utilization of prior knowledge of the nominal search space when recomputing
the policy for a perturbed problem. This stems from the intractable memory
requirements needed to store information about the task assignment and routing
search space when using the typical approach of using a single decision
variable to define both the task assignment and the robot routes[14].
Consequently, an entirely new problem is solved each time a small change to
the nominal problem is realized, restricting real-time applicability to small
problems [8].

In this paper, the routing problem is solved using a heuristic and the task
assignment search space is explored using a Monte Carlo Tree Search (MCTS)
algorithm. For task assignment problems, MCTS algorithms build a search tree
that stores cost estimates of task assigning decisions while exploring the
search space and producing solutions with low optimality gaps [15, 16]. The
contribution of this paper lies in the re-utilization of the search tree
topology and cost estimates as prior knowledge when the problem is perturbed.
This prior knowledge is used in a transfer learning framework to rapidly
update cost estimates while also generating solutions, after which the MCTS
algorithm is re-initialized. Computational experiments are performed on a
modified TSPLIB instance [17] to capture realistic FMS operational constraints
of charging policies, payload and battery constraints, and pickup - delivery
requirements of unique items. Results show that the solutions obtained using
prior knowledge have a lower optimality gap than when the perturbed problem is
approached as an entirely new problem.

##  II PROBLEM DEFINITION

Consider that $n$ material handling tasks are to be completed, and the
different commodities are represented by the set
$\mathcal{H}:=\\{h_{1},h_{2},...,h_{n}\\}$. The set of paired pickup and
delivery locations are defined by $\mathcal{V}^{P}:=\\{1,2,...,n\\}$ and
$\mathcal{V}^{D}:=\\{n+1,n+2,...,2n\\}$ respectively. Define
$\mathcal{V}:=\mathcal{V}^{P}\cup\mathcal{V}^{D}$, and let each location
$i\in\mathcal{V}$ be associated with a cargo mass $q_{im}\in\mathbb{R},\forall
m\in\mathcal{H}$. A commodity picked up at $i\in\mathcal{V}^{P}$ is paired
with a delivery location $n+i\in\mathcal{V}^{D}$, such that
$q_{im}+q_{i+n,m}=0$. The start and end locations of each robot are at the
depot defined by the nodes $\\{0,2n+1\\}$, and also serve as the charger
location. Define $\overline{\mathcal{V}}:=\mathcal{V}\cup\\{0,2n+1\\}$ so that
the graph representation is given by
$\mathcal{G}:=(\overline{\mathcal{V}},\mathcal{E})$, where
$\mathcal{E}:=\\{(i,j)\in\overline{\mathcal{V}}\times\overline{\mathcal{V}}:i\neq
j\\}$ denotes the set of edges. The set $\mathcal{T}:=\\{1,2,...,t_{max}\\}$
denotes the types of robots available for performing the defined material
handling tasks. For each type of robot $t\in\mathcal{T}$, the battery size is
defined by $B^{t}$ and its payload capacity is $Q^{t}$. The set
$\mathcal{R}^{t}:=\\{1,2,...,r^{t}_{max}\\}$ represents robots of type $t$
present in the fleet. For each type $t\in\mathcal{T}$, the energy to travel
between each node pair $(i,j)\in\overline{\mathcal{V}}$ is defined as
parameter $\delta e_{ij}^{t}\in\mathbb{R}^{+}$, normalized to be a fraction of
the battery capacity $B^{t}$. The problem formulation is defined in Eq.1
below:

$\displaystyle J=\min_{\begin{subarray}{c}x_{ij}^{a_{t}}\\\
\end{subarray}}~{}\sum_{(ij)\in\mathcal{E}}~{}\sum_{t\in\mathcal{T}}~{}\sum_{a_{t}\in\mathcal{R}^{t}}E_{ij}^{a_{t}}x_{ij}^{a_{t}}$
(1a) s.t. $\displaystyle
x_{ij}^{a_{t}}\in\\{0,1\\}\quad\quad\forall(i,j)\in\mathcal{E},~{}a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$
(1b) $\displaystyle\sum_{j\in\mathcal{V}^{p}}x_{0j}^{a_{t}}\leq
1~{}~{}\quad\forall a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1c)
$\displaystyle\sum_{i\in\mathcal{V}^{D}}x_{i,2n+1}^{a_{t}}\leq 1\quad\forall
a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1d)
$\displaystyle\sum_{(ij)\in\mathcal{E}}x_{ij}^{a_{t}}\leq
1~{}\quad\quad\quad\forall a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1e)
$\displaystyle\sum_{i\in\mathcal{V}}x_{ij}^{a_{t}}=\sum_{k\in\mathcal{V}}x_{jk}^{a_{t}}\quad\forall
j\in\mathcal{V},~{}a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1f)
$\displaystyle y_{0m}^{a_{t}}=0\quad\forall
m\in\mathcal{H},a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1g)
$\displaystyle\begin{split}&y_{jm}^{a_{t}}=y_{im}^{a_{t}}+\sum_{i\in\mathcal{V}}q_{jm}x_{ij}^{a_{t}}\\\
&\quad\quad\quad\quad\quad\forall
m\in\mathcal{H},j\in\mathcal{V},~{}a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}\end{split}$
(1h)
$\displaystyle\begin{split}&\sum_{i\in\mathcal{V}}y_{im}^{a_{t}}x_{ij}^{a_{t}}=-q_{jm}\\\
&\quad\quad\quad\quad\quad\forall
m\in\mathcal{H},j\in\mathcal{V}^{D},a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}\end{split}$
(1i) $\displaystyle\sum_{m\in\mathcal{H}}y_{im}^{a_{t}}\leq
Q^{t}\quad\quad\forall
i\in\overline{\mathcal{V}},a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1j)
$\displaystyle z_{j}^{a_{t}}=\begin{cases}z_{i}^{a_{t}}-\delta
e_{ij}^{t}&\text{if }x_{ij}^{a_{t}}=1\wedge z_{i}^{a_{t}}-\delta
e_{ij}^{t}>0\\\ 1-\delta e_{0j}^{t}&\text{if }x_{ij}^{a_{t}}=1\wedge
z_{i}^{a_{t}}-\delta e_{ij}^{t}\leq 0\\\ \end{cases}$
$\displaystyle\quad\quad\quad\quad\quad~{}~{}\quad\quad\forall(i,j)\in\mathcal{E},~{}a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$
(1k) $\displaystyle z_{0}^{a_{t}}=1;~{}0\leq z_{i}^{a_{t}}\leq 1~{}~{}\forall
i\in\overline{\mathcal{V}},~{}a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$ (1l)
$\displaystyle E_{ij}^{a_{t}}$ $\displaystyle=\begin{cases}B^{t}\delta
e_{ij}^{t}&\text{if }x_{ij}^{a_{t}}=1\wedge z_{i}^{a_{t}}-\delta
e_{ij}^{t}>0\\\ B^{t}(\delta e_{i0}^{t}+\delta e_{0j}^{t})&\text{if
}x_{ij}^{a_{t}}=1\wedge z_{i}^{a_{t}}-\delta e_{ij}^{t}\leq 0\\\ \end{cases}$
$\displaystyle\quad\quad\quad\quad\quad~{}~{}\quad\quad\forall(i,j)\in\mathcal{E},~{}a_{t}\in\mathcal{R}^{t},~{}t\in\mathcal{T}$
(1m)

The goal of minimizing the total energy traveled by all robots in the fleet is
captured in Eq. (1a) of the MHT problem formulation. Here,
$E_{ij}^{a_{t}}\in\mathbb{R}^{+}$ accounts for charge events and is the energy
expense of a robot $a_{t}\in\mathcal{R}^{t}$ between a pair of nodes
$(i,j)\in\mathcal{E}$. Binary variables $x_{ij}^{a_{t}}$ are used to indicate
whether robot $a_{t}$ of type $t\in\mathcal{T}$ uses edge
$(i,j)\in\mathcal{E}$. If a robot is assigned a task, it must start and end at
the depot, as specified by Eq. (1c) and (1d) respectively. Additionally, the
robot is permitted to visit each location at most once, as enforced by Eq.
(1e), and must leave the location after completing the visit, as defined in
Eq. (1f). Payload variables $y_{im}^{a_{t}}$ are used to define the mass of
commodity $m\in\mathcal{H}$ being carried by robot $a_{t}\in\mathcal{R}^{t}$
as it leaves node $i\in\overline{\mathcal{V}}$. All robots start their tour
with no payload at the depot, as defined in Eq. (1g). The evolution of the
commodity-wise payload is defined in Eq. (LABEL:eq:_mass_evolution) as the
robot visits locations in its tour. Precedence constraints for each commodity
are defined in Eq. (LABEL:eq:cargoPrecedence), meaning that a robot can visit
a delivery location if and only if the corresponding commodity has been
previously picked up. Payload limitations are captured in Eq. (1j). The state
of charge (SOC) of robot $a_{t}$ as it arrives at location $j$ is given by
$z_{j}^{a_{t}}$ and Eq. (1) defines the charging policy that requires a robot
to head to the depot for a recharge if required. As described in Eq. (1), the
energy expense $E_{ij}^{a_{t}}$ between locations $i$ and $j$ is dependent on
whether a recharge event occurs between the two locations.

Figure 1: Illustration of perturbations to the problem formulation

The uncertainties of real-world deployment cause disruptions that manifest as
changes in the parameters of Eq. (1). In Fig. 1, an illustrative plant layout
is shown, with green arrows marking 4 material handling tasks and blue dotted
lines indicating the centralized communication between the FMS and the robots.
A blocked aisle or a change in pickup position results in a change in
parameter $\delta e_{ij}$ in Eq. (1). Similarly, a degraded robot battery
capacity changes parameter $B^{t}$. The FMS objective is to quickly reassign
tasks and routes to the robots to adapt to the updated MHT definition.

##  III METHODOLOGY

The fleet policy both assigns tasks to the robots and also routes each robot,
and is defined by the value of the binary variable $x$. Finding the optimal
solution to Eq. (1) is computationally expensive because of the
$\mathcal{NP}$-hardness of the problem and the nonlinear constraints. The
proposed framework uses an offline MCTS algorithm to first compute a near-
optimal solution to the nominal problem. This utilizes a sufficiently high
computational time budget since the MHT parameters are known well in advance.
This produces a richly populated MCTS search tree with cost estimates for task
assignment decisions. When a perturbation is realized, an online algorithm
uses these cost estimates to rapidly obtain feasible solutions to the updated
problem.

###  III-A Solving the nominal problem offline

The MCTS algorithm explores the task assignment search space, generating a
search tree whose nodes represent decisions related to assigning a robot to a
task. The root node represents the start of the decision-making process where
no tasks have been assigned. The terminal node of the tree represents the
final outcome of the task-assigning process, signifying that all the tasks
have been assigned to the available robots. As the tree is traversed from the
root node to a terminal node, tasks are assigned to robots based on their
order in the defined task list. The parent of a node is the node that precedes
it in the decision-making process, that is, the robot assigned the previous
task in the task list. Similarly, its child nodes are the nodes that
immediately follow it, and represent the robots available for selection at the
next task in the task list. Each child node is connected to its parent by a
branch that represents the decision of assigning the next task to that robot
while fixing the previous decisions from the parent node to the root node. A
leaf node does not have any child nodes, and if non-terminal, indicates that
some task-assignment decisions have not yet been made. The MCTS algorithm is
as enumerated in Fig. 2:

1) Selection: Starting from the root node, the algorithm traverses the search
tree by selecting child nodes based on a selection policy. For the cost
minimization objective, the Lower Confidence Bound (LCB) selection policy is
used:

$\text{LCB}(s)=\operatorname*{arg\,min}_{s^{\prime}\in\text{children of
}s}\frac{J(s^{\prime})}{N(s^{\prime})J_{max}}-\gamma\sqrt{\frac{\ln
N(s)}{N(s^{\prime})}}$ (2)

where $N(s)$ is the number of cost explorations at node $s\in S$ and $S$ is
the set of nodes that constitute the search tree. $J(s)$ is the sum of costs
from all the previous visits to node $s$. The constant $\gamma$ balances the
exploitation of promising nodes with the exploration of unfavorable nodes that
are visited less often. This ensures that the entire search space is
systematically explored when given sufficient computation time. During the
conducted explorations, the maximum cost found is denoted by $J_{max}$, and is
used as a normalization factor that is continually updated as the search
proceeds.

Figure 2: Flow diagram of the Monte Carlo Tree Search Algorithm

The process of selecting child nodes by applying Eq. (2) starts from the root
node and continues until a leaf or terminal node is reached. If the selected
node is a previously visited child node, an expansion of the tree is conducted
as defined in step 2) while if it was unvisited, a rollout is conducted per
step 3). On the other hand, if it is a terminal node, the costs of routing are
computed for each robot and the costs are backpropagated according to step 4).

2) Expansion: After selecting a node based on the LCB policy, the tree is
expanded by adding child nodes to the selected node to represent possible task
assignments. Once new child nodes have been added, one is selected for
rollout.

3) Rollout: If the selected node is a leaf node, Monte Carlo sampling randomly
assigns the remaining tasks to fleet robots until the terminal node is
reached. Here, a robot has been assigned to each task and the total cost
associated with the assignment is obtained by solving the routing problem for
each robot. While numerous approaches exist to obtain the routing cost for a
single robot, in this paper the recursive B&B Alg. 1 derived from [18] is used
to account for the nonlinear constraints associated with precedence and
charging policies. To limit computation time, the B&B of Alg. 1 is terminated
after a 0.1 second time cap, since reasonably good routes are expected due to
the best first order of exploration in the recursive algorithm. The total
routing cost for the entire fleet then provides an estimate of the cost of
selecting that node in Step 1.

4) Backpropagation: To update the tree based on the outcome of the conducted
rollouts, the algorithm traverses the search tree from the selected leaf or
terminal node $s_{l}$ up to the root node. For each parent node $s$ whose
selection by Eq. (2) resulted in the evaluation at $s_{l}$, the number of
visits is updated as $N(s)\leftarrow N(s)+r$, where $r$ is the number of
rollouts conducted. The accumulated costs for these nodes are also updated as
$J(s)\leftarrow J(s)+\sum_{i=1}^{r}J_{r}(s_{l})$, where $J_{r}(s_{l})$
represents costs obtained from the rollout at $s_{l}$.

The four steps are repeated until the pre-defined and problem-specific
computational budget of time or number of iterations is exhausted. Throughout
the MCTS exploration, the task assignment that resulted in the minimum cost
$J_{min}$ is referred to as its incumbent solution. Like $J_{max}$, the
incumbent solution is also continually updated as the search progresses and is
the output of the MCTS algorithm when terminated. The resulting search tree
topology and the average cost $J(s)/N(s)$ at each node $s$ are saved as the
prior knowledge of the nominal problem which will be utilized when a
perturbation occurs.

1: sequenceCost = _B &B_(robotState, taskList, _location_) 2: Find feasible
next locations based on payload, cargo, SOC  3: Sort locations by operational
cost of branching to that location  4: for  _i_ in feasible locations do 5:
_branchCost_ = tourCost + operational cost(i)  6: if  _branchCost_ $\geq$
robotState.bestCost then 7: continue { skip to next _location i+_}  8: else if
_branchCost_ $<$ robotState.bestCost then 9:  _State+_ = Update robotState:
SOC, position, remaining locations  10: if number of remaining locations $>0$
then 11: Cost = _B &B_(robotState+, taskList, _location_(i)) Recursive Alg. 1
12: else 13: State.bestCost = Cost  14: end if 15: end if 16: end for 17:
Return robotState

Algorithm 1 Routing B&B

###  III-B Solving the perturbed problem online

Figure 3: Schematic of the proposed algorithm

Consider the case where the nominal Eq (1) has been addressed using the
offline MCTS algorithm, generating a task assignment search tree, referred to
as the nominal tree. In real-world operation, when perturbations affect the
problem definition, the proposed method makes use of the nominal tree topology
and cost information as follows:

  1. 1.

The leaf nodes $s_{l}$ of the nominal search tree are first ordered in
increasing average costs to the nominal problem given by $J(s_{l})/N(s_{l})$.

  2. 2.

A search tree that replicates the topology of nodes and branches in the
nominal tree is initialized. However, at each node of this perturbed search
tree, the number of visits $N(s)$ and accumulated costs $J(s)$ are set to
zero.

  3. 3.

For a predefined parameter $k$, select the $k^{th}$ percentile of ordered leaf
nodes. Rollouts are then conducted for these nodes in a nominally cheapest-
first order. During the rollout process, leaf node costs are re-evaluated
under the perturbed problem parameters. The updated costs and number of visits
are then backpropagated through the perturbed tree as defined in step 4) of
Section III-A. This simultaneously updates the search tree while exploring
promising leaf nodes. The incumbent solution to the perturbed problem is
continuously updated and available for a policy update if necessary.

  4. 4.

Once all the selected leaf nodes have undergone re-evaluation, the MCTS
algorithm, as described in Section III-A, is re-initialized on the updated
perturbed tree. This creates new nodes and utilizes the balance of
exploitation and exploration to further reduce the cost for the duration of
the remaining computation time.

The proposed method is illustrated in Fig. 3, demonstrating how the topology
of the nominal search tree is replicated, and the costs associated with
promising leaf nodes are re-evaluated with backpropagation. The figure also
shows the generation of new nodes within the perturbed tree once the online
MCTS algorithm is initialized in Step 4). An inherent assumption in this
approach is that the perturbation does not change the topology of the existing
search tree, that is, the number of robots is not changed, since this would
not permit the topology of the nominal tree to be reused. When perturbations
are bounded, which is a reasonable assumption for the controlled environment
of an FMS, the online approach is expected to yield solutions with a reduced
optimality gap compared to not utilizing the nominal tree.

##  IV Computational Experiments

To test the effectiveness of the proposed algorithm across a variety of
perturbation types, the nominal MHT problem was first defined using a TSPLIB
benchmark instance [17]:

  1. Step 1:

Load the TSPLIB $eil51$ point cloud to obtain a set of $n$ points with defined
Cartesian coordinates.

  2. Step 2:

Find the centroid of the point cloud.

  3. Step 3:

Sort and assign indices $1,2,...,n$ to the points by order of increasing
distance from the centroid.

  4. Step 4:

Designate the point with index 1 as the depot.

  5. Step 5:

Define precedence constraints between points with pairs of indices as $(2\prec
n),(3\prec n-1)$, and so on.

The resulting point cloud defines a depot and 25 MHTs, each with paired pickup
and delivery locations. In the nominal case being studied, the MHTs are to be
completed by two robots with payload capacities $Q^{1}=Q^{2}=10$ commodities.
The energy expense associated with traveling a Euclidean distance of $d$ units
is defined to be $d$ kJ, and the battery capacity of each robot is nominally
$B^{1}=B^{2}=20$ kJ.

Extensive computational experiments were conducted in a Matlab R2022a
environment on an Intel Xeon E5-2680 v4 CPU clocked at 2.4 GHz at the Ohio
Super Computer [19]. Three fleet management strategies were compared for
perturbations associated with battery degradation, payload capacity
variations, and shifts in pickup and delivery locations. The optimal task
assignment solution was first found for the nominal problem and each
perturbation, requiring 370 hours of processing time each. This involved an
exhaustive search that first listed the $2^{25}$ possible task assignments,
and then obtained the cost of each task assignment using Alg. 1, after which
the lowest-cost task assignment was found.

Each experiment was repeated 25 times to account for the stochastic nature of
MCTS algorithms. In each repetition, the nominal problem was first solved
using the offline algorithm of Section III-A to populate a nominal search tree
over a computational time budget of 12 hours. The 25 resulting task assignment
solutions varied slightly due to the stochasticity but were within 5% of the
optimal solution.

For benchmarking, a decentralized approach was evaluated that continued to use
one of the 25 nominal task assignment solutions that were computed offline.
The problem perturbation was only addressed by rerouting the affected robots
by using Alg. 1, providing near-instantaneous recovery from the perturbation.
Another comparison was made with a centralized approach which executed the
offline algorithm of Section III-A ab initio when a perturbation is realized,
without utilizing nominal search tree information. This was also repeated 25
times to account for stochasticity. Finally, to evaluate the proposed online
CFM algorithm, each of the 25 nominal trees was used as prior knowledge, thus
producing 25 updated search trees for each perturbation.

For every MCTS algorithm used in these experiments, the parameter $\gamma$ in
Eq. (2) is set to $\sqrt{0.5}$, and the number of rollouts $r$ is set to 20.
When the nominal tree is utilized, the parameter $k$ of the online algorithm
is set to 0.05, implying that the 5th percentile of low-cost nominal leaf
nodes is first explored to acquire updated costs for the new problem.

(a) Incumbent solutions obtained by the algorithms

(b) 1 second  (c) 1 minute

(d) 1 hour

Figure 4: Battery capacity of one robot changed from 20 to 16 kJ

###  IV-A Variation in battery capacity

The battery capacity of one of the robots was perturbed from its original
capacity of $20$ kJ to $16$ kJ, without altering the capacity of the other
robot. Since the decentralized method does not optimize the task assignment,
the 25 task assignments obtained from the repeated solving of the nominal
problem produced 25 heuristic solutions, many of which overlap, as seen in
Fig. 4a. These solutions are found near-instantaneously because only the time-
capped heuristic Alg. 1 is used, but they are significantly outperformed
within 10 seconds by every repetition of the centralized methods. In the case
of these centralized algorithms, it is evident that when the perturbed problem
is solved without utilizing the nominal tree, the incumbent solution costs at
any instant are typically higher than when the nominal tree is utilized.
Histograms of incumbent solutions at one second, minute, and hour of
computation time are shown in Fig. 4b, 4c, and 4d respectively, indicating a
significant advantage to using the nominal tree, especially when computation
time is limited. Given sufficient computation time, it is seen that both the
centralized algorithms converge to the optimum solution in each of their 25
repetitions. Similar results are seen in Fig. 5 for the case when $B^{2}$ is
changed to $12$ kJ.

(a) Incumbent solutions obtained by the algorithms

(b) 1 second  (c) 1 minute

(d) 1 hour

Figure 5: Battery capacity of one robot changed from 20 to 12 kJ

(a) Nominal problem($\xi=0\%$)  (b) Perturbed problem: $\xi=4\%$

(c) Incumbent solutions obtained for $\xi$ = 4%

Figure 6: Spatial variations in pickup and delivery locations

###  IV-B Spatial variations of locations

The $eil51$ point cloud is shown in Fig. 6a, where grey line segments show the
precedence constraints. Let $\overline{x}$ and and $\overline{y}$ denote the
range of $x$ and $y$ coordinates respectively. The box uncertainty parameter,
denoted by $\xi$, causes deviations within the range of $\pm\xi\overline{x}$
and $\pm\xi\overline{y}$ in the $x$ and $y$ coordinates respectively, and
results in a change in the problem parameter $\delta e_{ij}$ in Eq. (1). A box
uncertainty of $\xi=4\%$ affects every location, as shown in Fig. 6b. The
performance comparisons shown in Fig. 6c show that both centralized methods
have a lower optimality gap than the decentralized approach, and the nominal
tree provides an advantage to the CFM.

(a) Payload capacity of one robot changed to 8

(b) Payload capacity of one robot changed to 6

Figure 7: Performance comparison for payload capacity variations

###  IV-C Variation in payload capacity

Cases studied in Fig. 7 relate to changes in the payload capacity of one of
the robots. It is seen that the centralized approach that uses the nominal
tree is able to improve upon the decentralized solutions within 10 seconds,
unlike the approach that does not use the nominal tree.

###  IV-D Discussion

Nominal task assignment solutions that were obtained using the offline MCTS
algorithm were within a 5% optimality gap for the 25 repetitions. However,
when the decentralized algorithm used these nominal solutions to adapt to a
perturbation, there was significantly higher variation in the solutions found.
This was in addition to a higher optimality gap as compared to the centralized
approaches that fully utilize the other robots of the fleet. When the payload
capacity was changed, the proposed centralized method found significantly
improved solutions as compared to when the nominal tree was not used. When a
change occurs in battery capacity or pickup-delivery locations, this
improvement is not as significant. This performance variation is as expected
because small changes to $\mathcal{NP}$-hard problems can result in drastic
changes in the optimal solution. In all cases and at any instance, solutions
obtained from the proposed method had a lower optimality gap than when the
nominal tree was not used.

##  V CONCLUSIONS

This paper presented a centralized fleet management strategy that utilizes
prior knowledge of the search space when there is a change in the nominal task
definition. The nominal material handling problem is first solved offline
using an MCTS algorithm for the task assignment problem, and using a heuristic
for the routing sub-problem. When the problem is perturbed, the proposed
online method evaluates the lowest-cost leaf nodes of the search tree first,
rapidly producing feasible low-cost solutions. The approach is verified to be
real-time capable and is shown to perform better than computing without using
the search tree and also the decentralized approach that does not attempt task
reassignment. Future work will seek to define the magnitude of perturbations
that can be handled by the developed algorithm, and also its capability for
larger fleets and other combinatorial problems such as the vehicle routing
problem with time windows.

## References

  * [1] Z. Ghelichi and S. Kilaru, “Analytical models for collaborative autonomous mobile robot solutions in fulfillment centers,” _Applied Mathematical Modelling_ , vol. 91, pp. 438–457, 2021.
  * [2] G. Fragapane, D. Ivanov, M. Peron, F. Sgarbossa, and J. O. Strandhagen, “Increasing flexibility and productivity in industry 4.0 production networks with autonomous mobile robots and smart intralogistics,” _Annals of operations research_ , vol. 308, no. 1-2, pp. 125–143, 2022.
  * [3] S. Scholz, “Decentral decision-making for energy-aware charging of intralogistics equipment,” _Logistics Research_ , vol. 16, no. 1, 2023.
  * [4] M. Daub, F. Duddeck, and M. Zimmermann, “Optimizing component solution spaces for systems design,” _Structural and Multidisciplinary Optimization_ , vol. 61, pp. 2097–2109, 2020.
  * [5] H.-G. Beyer and B. Sendhoff, “Robust optimization–a comprehensive survey,” _Computer methods in applied mechanics and engineering_ , vol. 196, no. 33-34, pp. 3190–3218, 2007.
  * [6] M. Sauer, A. Dachsberger, L. Giglhuber, and L. Zalewski, “Decentralized deadlock prevention for self-organizing industrial mobile robot fleets,” in _2022 IEEE International Conference on Omni-layer Intelligent Systems (COINS)_. IEEE, 2022, pp. 1–6.
  * [7] G. Fragapane, R. De Koster, F. Sgarbossa, and J. O. Strandhagen, “Planning and control of autonomous mobile robots for intralogistics: Literature review and research agenda,” _European Journal of Operational Research_ , vol. 294, no. 2, pp. 405–426, 2021.
  * [8] M. De Ryck, M. Versteyhe, and F. Debrouwere, “Automated guided vehicle systems, state-of-the-art control algorithms and techniques,” _Journal of Manufacturing Systems_ , vol. 54, pp. 152–173, 2020.
  * [9] B. H. O. Rios, E. C. Xavier, F. K. Miyazawa, P. Amorim, E. Curcio, and M. J. Santos, “Recent dynamic vehicle routing problems: A survey,” _Computers & Industrial Engineering_, vol. 160, p. 107604, 2021.
  * [10] R. Elshaer and H. Awad, “A taxonomic review of metaheuristic algorithms for solving the vehicle routing problem and its variants,” _Computers and Industrial Engineering_ , vol. 140, 2 2020.
  * [11] B. Fahimnia, H. Davarzani, and A. Eshragh, “Planning of complex supply chains: A performance comparison of three meta-heuristic algorithms,” _Computers & Operations Research_, vol. 89, pp. 241–252, 2018.
  * [12] A. Malus, D. Kozjek _et al._ , “Real-time order dispatching for a fleet of autonomous mobile robots using multi-agent reinforcement learning,” _CIRP annals_ , vol. 69, no. 1, pp. 397–400, 2020.
  * [13] A. Lodi, L. Mossina, and E. Rachelson, “Learning to handle parameter perturbations in combinatorial optimization: an application to facility location,” _EURO Journal on Transportation and Logistics_ , vol. 9, no. 4, p. 100023, 2020.
  * [14] R. F. Fachini and V. A. Armentano, “Logic-based benders decomposition for the heterogeneous fixed fleet vehicle routing problem with time windows,” _Computers & Industrial Engineering_, vol. 148, p. 106641, 2020.
  * [15] S. Edelkamp, M. Gath, C. Greulich, M. Humann, O. Herzog, and M. Lawo, “Monte-Carlo Tree Search for Logistics,” in _Lecture Notes in Logistics_. Springer Cham, 2015, pp. 427–440.
  * [16] C. Barletta, W. Garn, C. Turner, and S. Fallah, “Hybrid fleet capacitated vehicle routing problem with flexible Monte–Carlo Tree search,” _International Journal of Systems Science: Operations and Logistics_ , 2022.
  * [17] G. Reinhelt, “$\\{$TSPLIB$\\}$: a library of sample instances for the tsp (and related problems) from various sources and of various types,” _URL: http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/_ , 2014.
  * [18] T. Baltussen, M. Goutham, M. Menon, S. Garrow, M. Santillo, and S. Stockar, “A parallel monte-carlo tree search-based metaheuristic for optimal fleet composition considering vehicle routing using branch & bound,” _arXiv preprint arXiv:2303.03156_ , 2023.
  * [19] O. S. Center, “Ohio supercomputer center,” 1987. [Online]. Available: http://osc.edu/ark:/19495/f5s1ph73


# Crossover Equation of State Based on Color-Molecular-Dynamics

Nobutoshi Yasutake<EMAIL_ADDRESS>Department of Physics,
Chiba Institute of Technology, 2-1-1 Shibazono, Narashino, Chiba 275-0023,
Japan  Advanced Science Research Center, Japan Atomic Energy Agency, Tokai,
Ibaraki 319-1195, Japan  Toshiki Maruyama<EMAIL_ADDRESS>Advanced Science Research Center, Japan Atomic Energy Agency, Tokai, Ibaraki
319-1195, Japan

###### Abstract

Equation of State for dense matter is studied with color molecular dynamics,
in which hadron matter and quark matter are automatically distinguished only
from quark color state. The quark-quark interactions are optimized to be
consistent with saturation properties: symmetric energy, $L-$parameter, and
incompressibility around nuclear density. In the calculations, the degrees of
freedom of colors are solved at each numerical step, although the flavors are
fixed as up or down quarks. The resultant mass-radius relations also satisfy
the observational constraints such as the gravitational wave observations,
NICER, and “the two-solar mass observations” of neutron stars. In this model
with the allowed parameter range, deconfined quark matter appears in the core
of neutron stars via crossover. Although the current constraints from the
observations are not enough to conclude whether quark matter appears at high-
density region, our method would help to understand high-density material
properties inside neutron stars in the future.

###### pacs:

26.60.+c, 24.10.Cn, 97.60.Jd, 12.39.-x

##  I Introduction

The equation of state (EOS) is a key to understand neutron star (NS) physics.
In particular, it is an interesting and big question how hadron matter
composed of baryons deconfine to quark matter at high density. However, since
the lattice QCD simulations are not yet feasible for finite density and low-
temperature $T\simeq 0$, there remains significant theoretical uncertainty
regarding the behavior of matter inside NSs. It is a fundamental issue to
clarify deconfinement behaviors.

Instead of solely relying on theoretical and direct understanding of the
confinement-decontainment phase changes, the investigation of astronomical
observations and nuclear experiments has emerged as a promising approach. In
this context, it is the first step to reproduce the mass-radius (MR) relations
of NSs. One of the well known constraints on the MR relations is “two solar
mass constraints”, which was derived from the observation of PSR J0348+432
[1]. Gravitational wave observations also provide strong constraints on EOS.
The observation GW170817 has yielded an upper limit of the radius of neutron
stars assuming a binary system with canonical masses [2]. The electromagnetic
counterpart observations of the gravitational wave event have further
constrained the MR relations of NSs [48]. On the other hand, Neutron star
Interior Composition Explorer (NICER) has gradually narrowed the MR regions of
NSs, e.g. PSR J0030+045 [4, 5] and PSR J0740+6620 [6, 7]. Consequently, it is
expected that astrophysical observations will significantly constrain the EOSs
for high-density regions of NSs.

Thanks to these developments in astronomical observation techniques, it would
be nice to consider the possibility of the quark-hadron phase change in NSs.
This possibility has been pointed out by past studies, and many of them have
assumed a first-order phase transition [8, 9, 10, 11, 12, 13]. The phase
transitions, in this case, apply to that of multi-component systems, where two
phases are in equilibrium with different chemical compositions. In such phase
transitions nonuniform structures, namely “pasta structures” may appear, after
the balance between the charge interactions and surface tensions. Not only the
first-order phase transitions but also the possibility of the crossover are
discussed [14, 15, 16, 17, 18].

One of the points that we should keep in mind is that the model should
describe nuclear matter at lower densities as a hadron phase and at higher
densities as a quark phase, on the latter phase we don’t have enough
information. Furthermore, big uncertainties lie in the mechanisms of
confinement and deconfinement, and the EOS in between. At lower densities the
property of hadronic matter is rather well-known: the binding energy per
baryon of symmetric nuclear matter, $S_{0}$, has the minimum value of -16 MeV
at the saturation density $n_{0}$, around 0.16 fm-3. The range of symmetric
energy $J$, which is the energy difference per baryon between the symmetric
nuclear matter and the pure neutron matter, is relatively narrowly restricted:
a plausible range for $J$ of 29-33 MeV [19, 20]. However, the values of the
$L-$parameter (the slope of energy for neutron matter), and the
incompressibility $K$ (second derivative of symmetric matter energy by the
density) are strongly model dependent.

The fiducial range of the $L$-parameter is reported as $L=60\pm 20$ MeV, which
is constrained by analyses of terrestrial experiments, as summarized well in
Refs.[21, 22]. Compared with them, there are also the possibilities for larger
$L$-values: Danielewicz et al. have derived the range as $70<L<101$ MeV from
the isobaric analog states and iso-vector skin results [23]. Radioactive
Isotope Beam Factory (RIBF) at RIKEN in Japan has reported as $42<L<117$ MeV
[24]. The Pb radius experiment (PREX)-II projects have suggested the range of
$L$ as $L=106\pm 37$ MeV [25], while the analyses after the Ca radius
experiment (CREX) result have suggested the ordinal $L$-values recently.
Hence, the current terrestrial experiments did not solely conclude the value
of $L$.

It is well known that the $L-$parameter described above is strongly correlated
with neutron skin thickness. It was theoretically predicted that the neutron
skin thickness would change by considering alpha-clustering at the surface of
neutron-rich heavy nuclei. This point is mentioned in the paper on EOS by
Typel et al. which is often used in numerical simulations in theoretical
astronomy and widely known as “DD2” [26]. Recently, this prediction was
confirmed experimentally by Tanaka et al.[27]. In this study, DD2 EOS is
compared to our EOS as one benchmark.

As for the incompressibility $K$, Danielewicz et al.[28] pointed out 167 MeV
$<K<$ 300 MeV, where they analyzed the flow of matter to extract pressures in
nuclear collisions. Piekarewicz has given the range of $K=248\pm 8$ MeV from
the iso-scalar giant monopole resonance (ISGMR) in 208Pb adopted the
relativistic mean-field model with a random-phase-approximation [29]. G. Colò
et al. have also analyzed the measurements of the ISGMR in medium-heavy nuclei
introducing some types of Skyrme forces, then predicted the $K$ around 230 MeV
[30]. As for the prediction of a small value of $K$, Sturm et al. and Hartnack
et al. have concluded $K$ is around 200 MeV from a comparison of the results
of transport theories with the experimental data of heavy-ion collisions with
the production of $K^{+}$ mesons [31, 32].

The EOS of the deconfinement quark phase should be consistent with all
constraints as described above. For this purpose, Kojo et al. constructed EOSs
describing quark deconfinement with crossover [18]. These EOSs are based on
hadron EOSs derived from the chiral effective field theory nuclear equation of
state [33, 34, 35] at low densities. By interpolating this EOS with the quark
EOS at high density, they construct a crossover EOS. Such an interpolative
method is practical as a first step to obtain EOS, but one should be careful
to obtain thermodynamic quantities such as thermal conductivity since one
cannot clearly distinguish hadron matter and quark matter in the crossover
case. One possible approach to obtaining such physical quantities is to
discuss the deconfinement-confinement mechanism only within the framework of
the quark model.

Keeping this in mind, we employ a color molecular dynamics (CMD) simulation
which deals with constituent quarks with color degrees of freedom [36]. In
this method we solve the time evolutions of positions, momenta, and color
coordinates of quarks, which are governed by the Hamiltonian with the
potential term consisting of the color confining potential, the perturbative
gluon-exchange potential, and the meson-exchange potentials.

At low density, quarks are clusterized into baryons, in which the color-
singlet state is favored. As density increases, the baryons start to overlap
with each other, then they start deconfined into quark matter. Our CMD
simulations show such percolative behavior during the deconfinement. Compared
with our previous studies [37], we improved the scheme to include relativistic
kinetic energy and solve the time-dependence of the color [36]. We also
include color-independent non-linear quark-quark repulsions to keep the
consistency with mass-radius relations of compact stars constrained by the
astrophysical observations. The non-linear interaction can be understood as
the quark many-body effects.

This paper is organized as follows. In Sec. II, we outline our framework for
CMD. Sec. III contains numerical results consistent with the constraints of
mass-radius relations from astronomical observations, and with saturation
properties around nuclear density. Sec. IV is devoted to the conclusion and
the discussion of our results.

##  II Color Molecular Dynamics

Our formulation is based on our previous papers [36, 37]. Throughout this
paper, we assume quark isospin symmetry with the constituent quark mass, $m$.

We start with the total wave function defined by a direct product of single-
particle wave packets of quarks, the position and the momentum of which are
centered at time-dependent variables ${\bf R}_{i}$ and ${\bf P}_{i}$,
respectively, and $\chi_{i}$ is the internal degree of freedom given by a
direct production of the fixed flavor ${\chi_{i}}^{f}$, the time-dependent
color ${\chi_{i}}^{c}$, and the fixed spin orientation ${\chi_{i}}^{s}$,

$\Psi=\prod_{i=1}^{N}\frac{1}{(\pi{L_{q}}^{2})^{3/4}}\exp\left[-\frac{({\bf
r}_{i}-{\bf R}_{i})^{2}}{2{L_{q}}^{2}}+\frac{i}{\hbar}{\bf P}_{i}{\bf
r}_{i}\right]\chi_{i},$ (1)

here $N$ and $L_{q}$ denote the total number of quarks and the fixed width of
wave packets, respectively. We employ the width $L_{q}=0.37$ fm in this work.
In this paper, we fix the flavors and spins, hence ${\chi_{i}}^{f}$ and
${\chi_{i}}^{s}$ are treated as constants. The explicit form of the time-
dependent degree of freedom for color is shown as

$\displaystyle{\chi_{i}}^{c}=\left(\begin{array}[]{c}\cos\alpha_{i}\;e^{-i\beta_{i}}\;\cos\theta_{i}\\\
\sin\alpha_{i}\;e^{+i\beta_{i}}\;\cos\theta_{i}\\\
\sin\theta_{i}\;e^{i\varphi_{i}}\\\ \end{array}\right),$ (5)

hence $\alpha_{i}$, $\beta_{i}$, $\theta_{i}$, $\varphi_{i}$ are the variables
for the color of each particle.

The system follows the Hamiltonian,

$H=H_{0}+V_{\rm Pauli}-T_{\rm spur},$ (6)

where $H_{0}$ is the conventional Hamiltonian expressed as

$H_{0}=\left<\Psi\left|\sum_{i=1}^{N}\sqrt{m+\hat{\textbf{p}}^{2}_{i}}+\frac{1}{2}\sum_{i,j\neq
i}^{N}\hat{V}_{ij}\right|\Psi\right>.$ (7)

The first term is the kinetic term with relativistic kinematics. The second
term is the interaction term expressed as

$\displaystyle\hat{V}_{i,j}$ $\displaystyle=$
$\displaystyle-\sum_{a=1}^{8}\tau_{i}^{a}\tau_{j}^{a}\hat{V}_{\rm
C}+\hat{V}_{\rm M},$ (8) $\displaystyle\hat{V}_{\rm C}$ $\displaystyle=$
$\displaystyle\kappa\hat{r}_{ij}-\frac{\alpha_{s}}{\hat{r}_{ij}},$ (9)
$\displaystyle\hat{V}_{\rm M}$ $\displaystyle=$
$\displaystyle\frac{g_{\omega}^{2}C_{\omega}}{4\pi}\left(\sum^{n}_{j\neq
i}\frac{e^{-\mu_{\omega}\hat{r}_{ij}}}{\hat{r}_{ij}}\right)^{1+\epsilon_{\omega}}-\frac{g_{\sigma}^{2}C_{\sigma}}{4\pi}\left(\sum^{n}_{j\neq
i}\frac{e^{-\mu_{\sigma}\hat{r}_{ij}}}{\hat{r}_{ij}}\right)^{1+\epsilon_{\sigma}}$
(10) $\displaystyle+$
$\displaystyle\frac{\sigma_{i}^{3}\sigma_{j}^{3}}{4}\frac{g_{\rho}^{2}}{4\pi}\frac{e^{-\mu_{\rho}\hat{r}_{ij}}}{\hat{r}_{ij}}$

where $\hat{r}_{ij}\equiv|\hat{\bf r}_{i}-\hat{\bf r}_{j}|$ is the distance
between $i$-th and $j$-th quarks, $\tau^{a}_{i}=\lambda^{c}_{i}$/2 with
$\lambda^{c}_{i}$ being the Gell-Mann matrices. The color-dependent
interaction $V_{C}$ consists of the linear confining potential (the first
term) and the one-gluon exchange potential (the second term) as shown in
Eq.(9). The string tension of confinement $\kappa$ and the QCD fine structure
constant $\alpha_{s}$ are set as $\kappa=0.75$ GeV fm-1, and
$\alpha_{s}=1.25$, as shown in Ref. [36]. Note that these two interactions are
the main components to determine the nucleon mass $M$, and these parameter
sets are typical values [38]. In this method, we do not conduct the anti-
symmetrization of the total wave function, in which the numerical cost is
$N^{4}$ [39, 40]. Hence, the interaction is underestimated by a factor of 4
when we take the matrix element of $\tau^{a}_{i}\tau^{a}_{j}$, and then we
introduce the effective coupling constants as $\kappa_{\rm eff}=4\kappa$ and
$\alpha_{s}^{\rm eff}=4\alpha_{s}$ to get consistency with color SU(3)
algebra. The constituent quark mass is then obtained as $m=361.8$ MeV, which
is determined to reproduce the nucleon mass, $M=938$ MeV, with the above
parameter sets for the confinement and the one-gluon exchange.

As for the nuclear force, the non-perturbative gluon exchanges in the color
singlet channels would be the essential part. Hence, we introduce the
$\sigma+\omega+\rho$ quark-meson couplings acting between quarks [41]. These
coupling constants are set as $g_{\omega}=5.46$, $g_{\sigma}=3.23$, and
$g_{\rho}=8.19$. The meson-quark coupling constants $g_{\omega}$,
$g_{\sigma}$, and $g_{\rho}$, are estimated from the meson-nucleon couplings
as $g_{\omega}=g_{\omega N}/3=4.98$, $g_{\sigma}=g_{\sigma N}/3=3.09$, and
$g_{\rho}=g_{\rho N}/3=9$, in previous works [37]. These values of
$g_{\omega}$, $g_{\sigma}$, and $g_{\rho}$ in this paper are different by
$\sim$10 % from our previous works. We also introduce the small non-linearity
parameters $\epsilon_{\omega}$, and $\epsilon_{\sigma}$ for the $\omega$-, and
$\sigma$-exchange potentials. This non-linearity provides the density-
dependence to the two-body interaction, which enables us to control the
stiffness of matter in the present framework. We also introduce $C_{\omega}$,
$C_{\sigma}$to make the coupling constant $g_{\omega}$, $g_{\sigma}$
dimensionless, set as $C_{\omega}=1/(1+{\epsilon_{\omega}})$,
$C_{\sigma}=1/(1+{\epsilon_{\sigma}})$, which are same as Ref.[37]. We have
chosen $\epsilon_{\omega}=0.20$, $\epsilon_{\sigma}=-0.13$ to make our EOS
consistent with the constraints from astrophysical observations and
experimental nuclear physics described later. We have introduced the effective
widths $L_{\omega}$, $L_{\sigma}$, and $L_{\rho}$ as $L_{\omega}=$ 0.75 fm,
$L_{\sigma}=1.35$ fm, and $L_{\rho}=$ 1.30 fm for each meson exchange term.
Also, these values are almost the same as the previous work: $L_{\omega}=$
0.70 fm, and $L_{\sigma}=L_{\rho}=$ 1.30 fm [37]. The values $L_{\omega}$,
$L_{\sigma}$, and $L_{\rho}$ are set to larger values than $L_{q}$ so that the
quark-meson interactions replicate the nucleon-meson interactions.

Instead of the anti-symmetrization of the wave function, we introduce an
effective Pauli-potential, $V_{\rm Pauli}$. It acts as a repulsive force
between quarks with the same intrinsic degree of freedom such as flavor,
color, and spin, collectively denoted as $\chi_{i}$. According to our previous
studies, we employ the following form of Pauli potential,

$\displaystyle V_{\rm Pauli}=\sum_{i=1}^{N}\frac{C_{p}}{(q_{0}p_{0})^{3}}$
$\displaystyle\exp{\left[-\frac{({\bf R}_{i}-{\bf
R}_{j})^{2}}{2q_{0}^{2}}\right]}$ (11)
$\displaystyle\times\exp{\left[-\frac{({\bf P}_{i}-{\bf
P}_{j})^{2}}{2p_{0}^{2}}\right]}\delta_{\chi_{i},\chi_{j}}$

The parameters are set as $q_{0}=2.46$ fm, $p_{0}=240$ MeV, and $C_{p}=$ 131
MeV, to reproduce relativistic kinetic energy for free fermions at zero
temperature.

In this paper, we do not take into account the spurious zero-point energy,
$T_{\rm spur}$, which comes from the center-of-mass motion of clusters as
shown in Ref.[37]. These effects remain future problems.

Time evolution of the system is given by the Euler-Lagrange equation for
{${\bf R}_{i}$, ${\bf P}_{i}$, $\alpha_{i}$, $\beta_{i}$, $\theta_{i}$,
$\varphi_{i}$} with the classical Lagrangian [36]. Accordingly, the explicit
form of the equations of motion reads:

$\displaystyle\dot{\bf R}_{i}$ $\displaystyle=$ $\displaystyle{\partial
H\over\partial{\bf P}_{i}},\ \ \ \ \dot{\bf P}_{i}=-{\partial
H\over\partial{\bf R}_{i}},$ (12) $\displaystyle\dot{\beta}_{i}$
$\displaystyle=$ $\displaystyle-{1\over 2\hbar\sin
2\alpha_{i}\;\cos^{2}\theta_{i}}{\partial H\over\partial\alpha_{i}},$ (13)
$\displaystyle\dot{\theta}_{i}$ $\displaystyle=$
$\displaystyle{1\over\hbar\sin{2\theta_{i}}}{\partial
H\over\partial\varphi_{i}},$ (14) $\displaystyle\dot{\alpha}_{i}$
$\displaystyle=$ $\displaystyle{1\over 2\hbar\sin
2\alpha_{i}\;\cos^{2}\theta_{i}}{\partial H\over\partial\beta_{i}}-{\cos
2\alpha_{i}\over 2\hbar\sin 2\alpha_{i}\;\cos^{2}\theta_{i}}{\partial
H\over\partial\varphi_{i}},$ (15) $\displaystyle\dot{\varphi}_{i}$
$\displaystyle=$ $\displaystyle-{1\over\hbar\sin{2\theta_{i}}}{\partial
H\over\partial\theta_{i}}+{\cos 2\alpha_{i}\over 2\hbar\sin
2\alpha_{i}\;\cos^{2}\theta_{i}}{\partial H\over\partial\alpha_{i}}.$ (16)

In the calculations, all quarks are, initially, distributed randomly without
momenta in a box with the periodic boundary condition. The ground state
(matter at zero temperature) is obtained by the frictional cooling [36, 37].
For this purpose, we solve a damping equation of motion instead of Eq.(12),

$\displaystyle\dot{\bf R}_{i}={\partial H\over\partial{\bf
P}_{i}}+\mu_{R}{\partial H\over\partial{\bf R}_{i}},$ (17)
$\displaystyle\dot{\bf P}_{i}=-{\partial H\over\partial{\bf
R}_{i}}+\mu_{P}{\partial H\over\partial{\bf P}_{i}},$ (18)

where $\mu_{R}$ and $\mu_{P}$ are damping coefficients set as
$\mu_{R}=-0.00002$ and $\mu_{P}=-0.02$.

Based on the above formulation, we search for parameter sets consistent with
the constraints from astronomical observations and nuclear experiments. Once
the parameters are given, the energy per nucleon for symmetric nuclear matter
and neutron matter are obtained. Then, we have made fits of the ground-state
energies by regressions. Using the regressions, we calculate the EOS of matter
in the charge neutral and beta-equilibrium including the contribution of
electrons. Thanks to GPU parallel computing, it takes just one day to check
one parameter set. In molecular dynamics, the heaviest numerical cost is in
the particle correlation part. There are some well-known techniques to speed
up calculations in molecular dynamics, such as the tree method and Fast
Multipole Method (FMM) [42], but, to our knowledge, there does not exist such
methods involving confining potential, which does no decrease at a long
distance, hence we do not use such coarse graining techniques.

##  III RESULTS

###  III.1 EQUATION OF STATE
BY COLOR MOLECULAR DYNAMICS

Figure 1:  Energy components per nucleon for $ud$ matter (upper panel) and
$udd$ matter (lower panel). Each dot represents the numerical results by CMD.
The energies, which come from quark-meson couplings of $\omega$, $\sigma$,
$\rho$, are shown as $E_{\omega}$, $E_{\sigma}$, $E_{\rho}$. The captions,
$E_{\rm kin}$ and $E_{\rm Pauli}$, are the quark kinetic energy and the quark
Pauli energy. The energies originated from the confinement and the one-gluon
exchange are expressed as $E_{\rm conf}$ and $E_{\rm OGE}$.

Let us first show the energy components per nucleon for $ud$ and $udd$ matter
in FIG.1. In symmetric nuclear matter, $u$ and $d$quarks are equally present,
hence we call it $ud$ matter. In the same way, neutron matter is called $udd$
matter in this paper. The total energy consists of the potential energies of
quark-meson coupling such as $q$-$\omega$, $q$-$\sigma$, and $q$-$\rho$, the
confining potential, the one-gluon exchange potential, the Pauli potential,
and the kinetic energy, expressed as $E_{\omega}$, $E_{\sigma}$, $E_{\rho}$,
$E_{\rm conf}$, $E_{\rm OGE}$, $E_{\rm Pauli}$, and $E_{\rm kin}$,
respectively. As the density increases, the main component of the energy
changes to $q$-$\omega$ and $q$-$\sigma$ couplings which are non-linearly
dependent on the density, while color-dependent potentials are rather moderate
since they have roughly linear dependence on the density.

Figure 2:  Density dependence of energy per baryon $E/A$ for each matter $ud$
matter and $udd$ matter. Each dot represents the numerical results by CMD.
Thin dot lines are the regression curves of them. The thick line, which is
labeled by “$\beta$-eq.+e”, is obtained $E/A$ under the charge-neutral and
beta-equilibrium conditions including contribution of electrons. The dashed
curve represents $E/A$ for $ud$ matter with electrons under the charge-neutral
condition.

The total energies of matter are shown in FIG.2. Note that these quantities
are obtained based on the quark models, namely CMD only. In this figure, the
regression curves are obtained assuming a sigmoid function with regularization
terms as

$\displaystyle(E/A)_{ud}(x)=3\Big{(}a_{1}x+a_{2}x^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle+\frac{a_{3}}{1+\exp(-a_{4}x+a_{5})}-\frac{a_{3}}{1+\exp(a_{5})}\Big{)},$
(19)
$\displaystyle(E/A)_{udd}(x)=3\Big{(}b_{1}x^{2}+b_{2}x^{4}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle+\frac{b_{3}}{1+\exp(-b_{4}x+b_{5})}-\frac{b_{3}}{1+\exp(b_{5})}\Big{)},$
(20)

where $x$ is normalized baryon density as $x=n/{\bar{n}}$. Note that $\bar{n}$
is not necessarily the saturation density, but just a normalization factor set
as ${\bar{n}}=0.16$ fm-3. Factor 3 at the beginning represents “three quarks”,
i.e. per nucleon. In other words, by removing this factor, we can obtain the
energy per quark, $E/Q$. The first two terms in both equations are introduced
for the regularization of the regression. The last terms are introduced to
make the regression curves on the origin. The obtained parameters are shown in
Table 1. The resultant fitting curves of energy per baron for $ud$ matter and
$udd$ matter are also shown in Fig.2. From these curves, we can obtain the
characteristic values for EOS around saturation density $n_{0}$: the
saturation energy $S_{0}$, the symmetry energy $J$, the density gradient of
neutron-matter energy $L$, and the curvature of symmetric-matter energy $K$.
They are summarized in Tab.2.

On the other hand, when we focus on high-density matter inside NSs, the
charge-neutral condition and the beta-equilibrium are considered to be
realized. However, even for hadron matter, it is described by quark dynamics
in this model.

Hence, the charge-neutral condition is evaluated as,

$\displaystyle 0=\sum_{i=u,d,e}Q_{i}n_{i},$ (21)

where $Q_{i}$, $n_{i}$ are the particle charge and density. We do not take
into account muons for simplicity. In this study, we also neglect anti-
particles, since we obtain ground state at zero temperature conducting the
frictional coolings shown in Eq.(17) and Eq.(18). In Fig.2, we show both cases
of $ud$ matter with/without the contribution of Coulomb energy. The former one
is obtained under the charge-neutral condition, and shown as the dashed curve
in Fig.2. As for the latter one, the role of electrons is necessarily
important, hence it does not include electrons.

In addition, the beta-equilibrium must be also fulfilled between $u,d$ quarks
and electrons in any phases:

$\displaystyle\mu_{u}+\mu_{e}=\mu_{d},$ (22)

where $\mu_{u}$, $\mu_{d}$, and $\mu_{e}$ are the chemical potentials for $u$,
$d$ quarks, and electrons, respectively.

Once we obtain $(E/A)_{ud}$ and $(E/A)_{udd}$ curves, the energy per nucleon
$E/A$ under the charge-neutral condition and beta-equilibrium is calculated as
the mixture of $ud$ matter and $udd$ matter.

$\displaystyle E/A$ $\displaystyle=$
$\displaystyle(1-\beta^{2})(E/A)_{ud}+\beta^{2}(E/A)_{udd},$ (23)

where $\beta$ is defined to denote the asymmetry of $u,d$ quark density as

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle
3\cfrac{n_{d}-n_{u}}{n_{d}+n_{u}}.$ (24)

In our calculation, we have to find the optimal value of $\beta$ consistent
with both conditions at each density. This result is also shown as the thick
line in FIG.2, including the energy contribution from electrons.

$a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$
---|---|---|---|---
$b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$
$\rm{[MeV]}$ | [MeV] | [MeV] |  |
-92.13 | 4.52325 | 1636.17 | 0.202968 | 0.263397
2.70704 | -0.00203469 | 38.0054 | 0.728765 | 2.83639
Table 1: The optimized parameters for Eq.(19) and Eq.(20). $n_{0}$ | $J$ | $L$ | $K$ | $S_{0}$
---|---|---|---|---
$[\rm fm^{-3}]$ | [MeV] | [MeV] | [MeV] | [MeV]
0.167 | 31.0 | 74.2 | 260 | -15.8
Table 2: The characteristic physical values around the saturation density
obtained by CMD calculations.

In the upper panel of FIG.3, we show the relationship between pressure and
baryon density corresponding $E/A$ under charge-neutral and beta-equilibrium
conditions shown in Eq.(23). However, below the subnuclear density
$n_{B}<0.06$ fm-3, we use the EOS by Baym et al.(BPS) [43]. There is no
density jump which is a characteristic of the first-order phase transition.
This suggests that the deconfinement of quarks appears as the crossover. The
EOSs by Akmal et al.(APR) [44], Kojo et al.(QHC21A-D) [18] and S. Typel et
al.(DD2) [26], available on the CompOSE archive [45], are included in the
figure for comparison. The EOSs by Akmal et al.(APR) and S. Typel et al.(DD2)
have often been adopted as benchmarks, but the EOSs by Kojo et al. are also
chosen for comparison in this paper since they are the state-of-the-art
crossover EOSs.

The corresponding sound speeds normalized by light speed are shown in the
lower panel of FIG.3. The band (painted region) shows the constraint from
“PSR+GW +J0030+J0740” traced Ref. [18], which has been originally deduced by
Legred et al.[46]. Since the EOSs except for APR are designed to be within
this range, it is obvious that they satisfy the constraint. It can be seen
that the sound velocity of our result is monotonically increasing with
density, but not for QHC A-D models. This difference is due to the quark
deconfinement occurring over a wide density range in our crossover EOS,
compared with QHC A-D models. Our results suggest that the pure quark matter
does not appear even in the core of neutron stars with the maximum mass.
However, we confirm that the peak of sound velocity appears in our CMD
calculations at high-density region, which is over the maximum density inside
of stable NSs: $n_{B}=1.08$ fm-3. Additionally, our model does not violate the
causality law in the stable interior of neutron stars with respect to the
speed of sound, but it does at higher density, $n_{B}=$ 1.42 fm-3. When we
extend our code to be fully relativistic in the future, it might be possible
that it may no longer violate the causality law at such high densities.

Figure 3:  The upper panel shows the pressure as a function of baryon density
under charge-neutral and beta-equilibrium conditions, along with the results
from Akmal et al.(APR) [44], Kojo et al.(QHC21A-D) [18], and S. Typel et
al.(DD2) [26] for comparison. The lower panel shows the corresponding sound
speeds normalized by light speed. The band (painted region) shows the
constraint from “PSR+GW +J0030+J0740”. See the text for the details.

The conditions for confinement/deconfinement are set as follows: If three
quarks are within a certain distance $d_{\rm cluster}$ and the total color of
three quarks is white with an accuracy $\varepsilon$, these quarks are
considered confined. It is formulated as

$\left\\{\begin{array}[]{l}\left|{\bf R}_{i}-{\bf R}_{j}\right|<d_{\rm
cluster}~{}~{}~{}~{}~{}(i,j=1,2,3),\\\
\displaystyle{\sum_{a=1}^{8}\left[\sum_{i=1}^{3}\langle\chi_{i}|\lambda^{a}|\chi_{i}\rangle\right]^{2}<\varepsilon}.\end{array}\right.$

The values of the criteria are set as $d_{\rm cluster}=0.33$ fm and
$\varepsilon=0.01$ in this study. The perspectives of $ud$ matter and $udd$
matter depended on the density are shown in FIG.4 and FIG.5. In these figures,
the colors of the quarks are represented using the discontinuous colors red,
blue, and green (RGB) for visibility, however in the actual calculations, the
quarks are in a mixed state of RGB reflecting the internal degrees of freedom
in Eq. (5). Namely, the color is selected and depicted among RGB, which is
closest to the color mixed state for each particle. Particles with white
colors indicate quarks in the confined state, while red, blue, and green ones
show in the deconfined quark matter. In these figures, it can be seen that the
fraction of deconfined quarks gradually increases with density. This behavior
is similar to the percolative depiction [47].

0.833 fm-3 1.166 fm-3 1.583 fm-3 1.916 fm-3

Figure 4:  The perspectives for $ud$ matter depended on density. Each color
corresponds to the color’s internal degree of freedom for each quark: the
white color balls represent the quarks in the baryon state, while red, blue,
and green ones do in the deconfined quark matter.

0.833 fm-3 1.166 fm-3 1.583 fm-3 1.916 fm-3

Figure 5:  Same as Fig. 4, but for $udd$ matter.

The fractions of confined quarks to total quarks in $ud$ and $udd$ matter are
shown in FIG.6. For both matters, deconfined quarks appear around 0.6 fm-3. As
we will show later, the maximum density in stable neutron stars is 1.08 fm-3,
hence our calculations suggest that the interior of neutron stars is mostly
filled with confined quark matter, i.e., baryons. The fraction of deconfined
quark matter is roughly estimated to be less than 0.3 in neutron stars.

Figure 6:  Fractions of quark number within baryon state against total quark
number for $ud$ matter and $udd$ matter, shown in Fig. 4 and Fig. 5. The
painted area indicates the region below the maximum density 1.08 fm-3, i.e.,
at the density allowed in the interior of stable NSs.

###  III.2 MASS-RADIUS RELATION

Finally, the mass-radius (MR) relation for neutron stars obtained by CMD
calculation is shown FIG. 7. The results by Akmal et al.(APR) [44], Kojo et
al.(QHC21A-D) [18], and S. Typel et al.(DD2) [26] are also included for
comparison. The radii obtained from CMD calculations are larger than those
from the two models, but still consistent with the constraints imposed by
astronomical observations, as shown below.

The constraints obtained from the astronomical observations with NICER, PSR
J0030+0451 [5], and MSP J0740+6620 [6, 7], are also shown in the figure. The
area bounded by two jagged circles represents 1$\sigma$ and 2$\sigma$ ranges
of the observational result based on the analysis assuming three uniform-
temperature oval spots [5]. Note that these ranges agree with the other result
assuming two uniform-temperature oval spots. As for PSR J0030+0451, see also
Ref. [4]. The meshed area on the constraint of MSP J0740+6620,
$R={11.5}_{-1.2}^{+1.8}$ km (68%), has been obtained from the NICER-only
analysis, and the painted are, $R={13.7}_{-1.5}^{+2.6}$ km (68%), were given
based on NICER+XMM-Newton data sets [7]. The mass of MSP J0740+6620 has been
suggested as ${2.08\pm 0.07}M_{\odot}$. Our CMD calculations give the maximum
mass 2.19 $M_{\odot}$, which is over the lower limit of the mass suggested
from MSP J0740+6620. The maximum density in stable neutron stars is $n_{B}$ =
1.08 fm-3 in our CMD calculations, which indicates that the interior of
neutron stars is occupied by mostly hadrons, as shown in FIG. 6. Note that our
result lies outside the NICER+XMM-Newton range (painted area). We believe that
it would be satisfied after taking into account the other physics, which
should be improved in the future, such as color-magnetic interactions.

Besides the above two constraints from NICER, we also take into account the
observations related to event GW170817, and show them in FIG. 7.

First, the gravitational wave observation, GW170817, itself gave us valuable
information on the tidal deformability of the neutron stars, which could be a
constraint for the radii [2]. The observation suggests that the radii should
be less than 13.6 km for neutron stars with a canonical mass, 1.4 $M_{\odot}$.
The dimensionless tidal deformability from our CMD calculation is
$\Lambda_{1.4M\odot}=458$, and the corresponding radius is less than the
constraint by GW170817.

The electromagnetic (EM) observation accompanied by GW170817 has also provided
complementary information: the merger was suggested not to be in a prompt
collapse to a black hole (BH) because of the large quantity of ejecta and its
high electron fraction. Given that the threshold for prompt collapse depended
on the NS compactness, Bauswein et al. had placed a lower limit of NS radius
as $R_{1.6M_{\odot}}$=10.3-10.7 km with the mass of $M=1.6M_{\odot}$ [48].
Moreover, comparing numerical relativistic simulations as for a supermassive
neutron star remnant with the measured binary mass $M_{\rm
tot}=2.74+0.04M_{\odot}$ [2], Shibata et al. have given an upper limit on the
NS mass of $M<2.3M_{\odot}$ [49]. For more details on EM counterparts, see
also the review by Metzger [50] and the reference therein.

Although not included in FIG 7 to avoid complications, it should be mentioned
that Margalit and Metzger gave a more strict constraint, $M<2.17M_{\odot}$
[51]. On the other hand, Keck-telescope spectrophotometry and imaging of the
companion of the “black widow” pulsar PSR J09520607, suggests 2.35 $\pm
0.17M_{\odot}$ [52]. Therefore, careful consideration is still needed
regarding the maximum mass of neutron stars. With this background in mind,
this paper adopts the more conservative constraint by Shibata et al. [49]
rather than the one by Margalit and Metzger [51].

Figure 7:  The mass-radius relation obtained from the numerical results of
CMD, where the maximum mass and the central baryon number density of NSs are
$M_{\rm max}=2.19M_{\odot}$ and $n_{B}=$ 1.08 fm-3, respectively. The circle
represents the maximum mass point, and each number close to the filled circles
indicates the central density normalized by the factor $\bar{n}$. For
comparison, the results by Akmal et al.(APR) [44], Kojo et al.(QHC21A-D) [18],
and S. Typel et al.(DD2) [26] are also included. The constraints obtained from
the astronomical observations with NICER, PSR J0030+0451 [4, 5], and MSP
J0740+6620 [6, 7], are also shown. The meshed area on the constraint of MSP
J0740+6620 has been obtained from the NICER-only analysis, and the painted
area was given based on NICER+XMM-Newton data sets. We also show the
constraints inferred from event GW170817, labeled by “GW170817”, “Bauswein et
al.”, and “Shibata et al.”. See the text for the details.

##  IV DISCUSSION

In this study, CMD calculations were performed to obtain the EOS, which is
consistent with the constraints from the astrophysical observations. In our
current framework, the deconfinement of quarks is crossover, and the insides
of NSs are occupied mainly by confined hadron matter. At low densities, our
results deviate slightly from the values suggested by various nuclear
experiments: our model satisfies just barely the range of $K$ by Danielewicz
et al.[28] but not for the others [29, 30, 31, 32], while it does the fiducial
values of $n_{0}$, $J$, and $L$ [19, 20].

However, we still miss a lot of physics: we adopt simple Newtonian
interactions without strangeness, quarks spin-spin interactions (color-
magnetic interactions), etc. One can easily anticipate that the EOS would be
softened when we further take into account strangeness, but it is unclear how
relativistic interactions and color magnetic interactions would change the EOS
on the other side. Note that our model is based on the constituent quark
models, and then we do not take into account the chiral condensations. It is a
future challenge to incorporate such physics into our CMD code.

###### Acknowledgements.

We are grateful to N. Hoshi, Y. Mukobara, H. Tanihisa, S. Ozaki, C. J. Xia, P.
Gubler, K. Kyutoku, T. Muto, A. Park, S-H. Lee, M. Oka, and T. Hatsuda for
fruitful discussions. This work was supported by JSPS KAKENHI Grant Numbers
20H04742, 20K03951. The calculations were performed by the supercomputing
system HPE SGI 8600 at the Japan Atomic Energy Agency.

## References

  * [1] J. Antoniadis et al., Science 340, 1233232 (2013).
  * [2] B. P. Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * [3] A. Bauswein, O. Just, H-T. Janka, and N. Stergioulas, Astrophys. J. Lett. 850, L34 (2017).
  * [4] T. E. Riley et al., Astrophys. J. 887, L21, (2019).
  * [5] M. C. Miller et al., Astrophys. J. 887, L24, (2019).
  * [6] T. E. Riley et al., Astrophys. J. 918, L27, (2021).
  * [7] M. C. Miller et al., Astrophys. J. 918, L28, (2021).
  * [8] N. Yasutake, T. Maruyama, and T. Tatsumi, Phys. Rev. D 80, 123009 (2009).
  * [9] H. Chen, G. F. Burgio, H.-J. Schulze, and N. Yasutake, Astro. & Astrophys. 551, A13 (2013).
  * [10] N. Yasutake et al., Phys. Rev. C 89, 065803 (2014).
  * [11] K. Maslov et al., Phys. Rev. C 100, 025802 (2019).
  * [12] C-J. Xia, T. Maruyama, N. Yasutake, and T. Tatsumi, Phys. Rev. D 99, 103017 (2019).
  * [13] C-J. Xia, T. Maruyama, N. Yasutake, T. Tatsumi, H. Shen, and H. Togashi, Phys. Rev. D 102, 023031 (2020).
  * [14] K. Masuda, T. Hatsuda & T. Takatsuka, Eur. Phys. Jour. A 52, 65 (2016).
  * [15] G. Baym, S. Furusawa, T. Hatsuda, T. Kojo, and H. Togashi, Astrophys. J. 885, 42 (2019).
  * [16] T. Minamikawa, T. Kojo, M. Harada, Phys. Rev. C 103, 045205 (2021).
  * [17] D. Blaschke, E. -O. Hanu, S. Liebing, Phys. Rev. C 105, 035804 (2022).
  * [18] T. Kojo, G. Baym, and T. Hatsuda, Astrophys. J. 934, 46 (2022).
  * [19] M. B. Tsang et al., Phys. Rev. C 86, 015803 (2012).
  * [20] J. M. Lattimer and Y. Lim, Astrophys. J. 771, 51 (2013).
  * [21] X. Viñas, M. Centelles, X. Roca-Maza, and M. Warda, Eur. Phys. J. A 50, 27 (2014).
  * [22] B.-A. Li, P. G. Krastev, D.-H. Wen, and N.-B. Zhang, Eur. Phys. J. A 55, 117 (2019).
  * [23] P. Danielewicz, P. Singh and J. Lee, Nucl. Phys. A958, 147 (2017).
  * [24] J. Estee et al. , Phys. Rev. Lett. 126, 162701 (2021).
  * [25] B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021)
  * [26] S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010)
  * [27] J. Tanaka et al., Science 371, 260 (2021)
  * [28] P. Danielewicz, R. Lacey & W. G. Lynch, Science, 298, 1592 (2002).
  * [29] J. Piekarewicz, Phys. Rev. C 69, 041301 (2004).
  * [30] G. Colò et al., Phys. Rev. C 70, 024307 (2004).
  * [31] C. Sturm et al., Phys. Rev. L 86, 39 (2001).
  * [32] Ch. Hartnack et al., Phys. Rev. L 96, 012302 (2006).
  * [33] D. Lonardoni, I. Tews, I., S. Gandolfi, & J. Carlson, PhRvR, 2, 022033 (2020).
  * [34] C. Drischler, S. Han, J. M. Lattimer, et al. Phys. Rev. C 103, 045808 (2021).
  * [35] C. Drischler, J. W. Holt, & C. Wellenhofer, ARNPS, 71, 403 (2021).
  * [36] T. Maruyama, T. Hatsuda, Phys. Rev. C 61, 062201 (2000).
  * [37] Y. Akimura, T. Maruyama, N. Yoshinaga, and S. Chiba, Eur. Phys. Jour., A 25, 405 (2005).
  * [38] T. Yoshimoto, T. Sato, M. Arima, T.S.H. Lee, Phys. Rev. C 61, 065203 (2000).
  * [39] H. Feldmeier, Nucl. Phys. A 515, 147 (1990).
  * [40] A. Ono, H. Horiuchi, T. Maruyama, A. Ohnishi, Prog. Theor. Phys. 87, 1185 (1992).
  * [41] P. A. M. Guichon, Phys. Lett. B 200, 235(1988).
  * [42] H. Cheng, L. Greengard, V. Rokhlin, J. Comput. Phys. 155 468 (1999).
  * [43] G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * [44] A. Akmal, V. R. Pandharipande and D. G. Ravenhall, Phys. Rev. C 58 1804 (1998).
  * [45] https://compose.obspm.fr/eos/
  * [46] I. Legred, K. Chatziioannou, R. Essick, S. Han, & P. Landry, Phys. Rev. D104, 063003 (2021).
  * [47] T. Kojo, P. D. Powell, Y. Song, and G. Baym, Phys. Rev. D 91, 045003 (2015).
  * [48] A. Bauswein, O. Just, H-T. Janka, N. Stergioulas, Astrophy. J. Lett., 850, L34 (2017).
  * [49] M. Shibata, E. Zhou, K. Kiuchi, S. Fujibayashi, Phys. Rev., D100, 023015 (2019).
  * [50] B. D. Metzger, Living Reviews in Relativity, 23, 1 (2020).
  * [51] B. Margalit, B. D. Metzger, Astrophys. J. Lett., 850, L19 (2017).
  * [52] R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, Astrophys. J. Lett., 934 L17 (2022).