# IMA-GNN: In-Memory Acceleration of Centralized and Decentralized Graph
Neural Networks at the Edge

Mehrdad Morsali, Mahmoud Nazzal, Abdallah Khreishah and Shaahin Angizi
Department of Electrical and Computer Engineering, New Jersey Institute of
Technology, Newark, NJ 07102 USA

###### Abstract.

In this paper, we propose IMA-GNN as an In-Memory Accelerator for centralized
and decentralized Graph Neural Network inference, explore its potential in
both settings, and provide a guideline for the community targeting flexible
and efficient edge computation. Leveraging IMA-GNN, we first model the
computation and communication latencies of edge devices. We then present
practical case studies on GNN-based taxi demand and supply prediction and also
adopt four large graph datasets to quantitatively compare and analyze
centralized and decentralized settings. Our cross-layer simulation results
demonstrate that on average, IMA-GNN in the centralized setting can obtain
$\sim$790$\times$ communication speed-up compared to the decentralized GNN
setting. However, the decentralized setting performs computation
$\sim$1400$\times$ faster while reducing the power consumption per device.
This further underlines the need for a hybrid semi-decentralized GNN approach.

##  1\. Introduction

Graph data structures appear naturally in many fields such as user accounts in
a social network, atoms in a chemical molecule, and vehicles in a traffic
system. Graph Neural Networks (GNN) extend deep learning to graph data by
combining graph structure and node information through message passing and
aggregation (hamilton2017representation, ; sahu2017ubiquity, ). This enables
GNNs to deliver node embeddings to serve multiple downstream graph tasks such
as node classification (inferring the class label of a node), link prediction
(estimating the possibility of a link between given nodes), and graph
classification (inferring the class label of a graph). Essentially, GNNs
obtain embeddings for the nodes in a given graph. Each node has a
computational graph composed of its $k$-hop neighboring nodes. Node embeddings
are obtained by alternating between message passing, i.e., communicating local
information across nodes, and aggregation where received messages along with
previous node information are used to obtain an updated embedding. Message
passing is done according to the graph structure, whereas aggregation is done
by the (trainable) neural network layers of the GNN model (kiningham2022grip,
; yan2020hygcn, ). As shown in Fig. 1, for a sample input graph G, first, in
the aggregation stage, each node aggregates the information from all neighbors
(nodes 2, 4, 7, 8, 9) with its own data (node 3) and creates the
$\boldsymbol{Z}$ matrix that represents the aggregated node features from node
3 and its neighbors. During the feature extraction stage, the result of the
aggregation stage ($\boldsymbol{Z}$) is fed into a Multi-Layer Perceptron
(MLP) or a Convolutional Neural Network (CNN) model to generate the output
matrix indicated by $\boldsymbol{O}$. The GNN workflow performs the same steps
for all the nodes in the graph.

Graphs in many real-world application areas are naturally enormous. For
example, social media and e-commerce graphs such as Facebook and Amazon graphs
have billions of nodes and edges. Furthermore, local node features are
typically of large dimensions ranging between hundreds to several thousands
(hu2020open, ). This creates a corresponding time and memory space demand for
training and testing on the underlying GNN models (zheng2020distdgl, ;
zhou2020graph, ). This is because nodes are mutually dependent, and therefore,
graphs can not be arbitrarily divided into smaller subgraphs. Techniques such
as neighborhood sampling (hamilton2017inductive, ) may help to some extent,
but depending on the graph structure, even a sampled computation graph and
associated features may not fit in the memory of a single GPU. Besides, off-
chip memory access is a critical issue in the Von-Neumann computing
architecture (mandal2022coin, ; li2021gcnax, ; abedin2022mr, ; angizi2022pisa,
) and to minimize the latency and power consumption, the Processing-in-Memory
(PIM) accelerators have been set forth focusing on centralized GNNs. Along
this line, HyGNN (yan2020hygcn, ) supports hybrid computing and memory access
patterns and performs GCN computations efficiently using a hybrid
architecture. Utilizing two dedicated processing engines, HyGNN tackles
irregularity with an aggregation engine and leverages regularity with a
combination engine. The irregularity arises from the topology-dependent
aggregation process which is inherently random and sparse, while the
transformation process is a neural network operation that has a static and
regular execution pattern. In the same way, GRIP (kiningham2022grip, ) divides
the GNN computations into different engines. Multiple parallel prefetch and
reduction engines have been used for aggregation to alleviate the
irregularity. For regular computation and memory access patterns, GRIP
utilizes a high-performance matrix multiply engine with dedicated memory.
Another edge-centric paradigm, EnGN (liang2020engn, ), has been also
implemented using a ring-edge-reduce dataflow. EnGN alleviates the poor
locality of sparse and random connected vertices (nodes). AWB-GCN
(geng2020awb, ) proposes a hardware-based workload distribution auto-tuning
framework consisting of three workload rebalancing techniques to alleviate the
extreme workload imbalance. PIM-GCN (challapalle2021crossbar, ) presents a
node-stationary dataflow with support for compressed sparse row and column
graph representations.

Figure 1. GNN’s aggregation and feature extraction stages for a sample input
graph. Figure 2. (a) The proposed IMA-GNN architecture with resistive CAM
traversal core, resistive MVM aggregation, and feature extraction cores, (b)
Resistive MVM crossbar, (c) Resistive CAM crossbar.

Distributed (decentralized) GNN training and inference (li2020graph, ;
tolstaya2020learning, ) is based on dividing a given graph into smaller
subgraphs that can be more easily processed with distributed devices. Despite
its ability to mitigate the computation overhead by load sharing,
decentralized GNN operation faces a major bottleneck; the excessive
communication overhead between nodes in different distributed devices
(zheng2020distdgl, ). To the best of our knowledge, this work is among the
first to explore and compare in-memory acceleration in both centralized and
decentralized GNN settings and to offer a design guideline to the community.
The main contributions of this paper are as follows: (1) We develop a PIM
architecture with RRAM arrays based on a set of innovative micro-architectural
designs that can be optimized and used for centralized and decentralized GNN
inference for efficiency and speed-up; (2) We model the latency and power
consumption of GNN accelerators implemented in centralized and decentralized
settings considering the computation and communication between edge devices;
and (3) We present a solid bottom-up evaluation framework to analyze the
performance of the whole system in real scenarios and through adopting large
graph datasets.

##  2\. Proposed IMA-GNN

###  2.1. Architecture Overview

The IMA-GNN is a high-performance and energy-efficient RRAM crossbar-based
accelerator developed to execute GNN’s pivotal operations in both centralized
and decentralized settings inspired by (challapalle2021crossbar, ). As shown
in Fig. 2(a), IMA-GNN comprises three computation cores, i.e., traversal,
aggregation, and feature extraction as well as peripherals such as a buffer
array and a controller. The traversal core consists of resistive Content
Addressable Memory (CAM) crossbars capable of search and comparison operations
(Fig. 2(c)). All resistive CAM crossbars on the bottom side are connected to a
shared vector generator & scheduler unit and then to a high-bandwidth bus to
communicate with other cores at the top. The aggregation core includes
resistive crossbars (Fig. 2(b)) to perform in-situ Matrix-Vector-
Multiplication (MVM) operations for the feature aggregation. The feature
extraction core is designed with a similar resistive crossbar but a different
size to take care of transformation in GNN inference. The resistive MVM
crossbars are connected to a shared activation unit.

###  2.2. MVM & CAM Crossbars

The RRAM crossbar memory arrays are widely explored as a potential parallel
engine to execute MVM operation and scan and search (angizi2022pisa, ;
challapalle2021crossbar, ; abedin2022mr, ). As shown in Fig. 2(b) in the MVM
crossbars, the weight parameters are first stored as resistance states in each
RRAM device in a 1-Transistor-1-RRAM (1T1R) structure, and then the input
binary bit-strings, as the inputs to the crossbar array, are converted by the
Digital-to-Analog Converter (DAC) into voltages $V_{i}$ and applied to Bit-
Lines ($BL$s) in parallel. The weighted currents generated from the RRAM cells
sharing a Source-Line ($SL$) are accumulated resulting in an intrinsic dot-
products operation. The accumulated values are then sampled by Sample & Hold
unit and then converted to binary data using Analog-to-Digital Converters
(ADC). The partial-product results from each $SL$ are further processed by the
Shift & Add unit to generate the final result. As shown in Fig. 2(c), in the
CAM crossbars, each Ternary CAM (TCAM) cell consists of 2-Transistor-2-RRAM
(2T2R) to accomplish the XNOR search operation on each pair of cells. For this
operation, $BL$ and $\overline{BL}$s are valued with the search data by the
Search Data Driver. Accordingly, the Sense Amplifier connected to Match-Lines
($ML$SA) senses whether the row is a match or mismatch with the reference
connected to Vdd. In the compare operation, $BL$s are grounded and
$\overline{BL}$s are connected to increasing calibrated voltages from the
Least Significant Bit (LSB) to the Most Significant Bits (MSB).

###  2.3. Accelerator Dataflow

Once the edge buffers shown in Fig. 2(a) on the left have been loaded with
graph data in either centralized or decentralized GNN settings, the traversal
core starts processing edges. The traversal core performs two essential CAM-
based operations, i.e., search and compare. To maximize the data reuse of
feature data in IMA-GNN, the traversal core implements an efficient node-
stationary dataflow by buffering a set of node features in the buffer array
and reusing it for the aggregation core. IMA-GNN leverages a Compressed Sparse
Row (CSR) format (lu2019redesk, ) to form the Edge weight array (E), Column
Index array (CI), and Row Pointer array (RP) and loads the graph data to
search and scan CAMs (Fig. 2(a) 1). A sample graph adjacency matrix and the
corresponding CSR format are shown in Fig. 3(a)-(b). Any destination node then
operates as an input to the search CAM as shown in the data mapping in Fig.
3(c) and rows that match the search data are activated. Matching rows are
reference inputs for comparison in the scan CAM, which determines the source
nodes with edges to the destination node by comparing the input row with RP
(Fig. 3(d)). Next, the vector generator & scheduler unit receives the result
of scan CAM and edge data to render input control vectors for the aggregation
core (Fig. 2(a) 2). This will activate particular rows of resistive
aggregation core corresponding to incoming edges. Next, the aggregation core
input buffers receive input vectors for each destination node along with the
destination node. The aggregation core starts with source node features or
feature dimensions across its own cluster 3. IMA-GNN is equipped with double
buffering for feature data and graph data. This feature enables overlapping
writing/programming phases and the traversal stage. Next, the updated
destination node features are fed to the feature extraction core’s crossbars
programmed with weights 4. Besides, similar to (mandal2022coin, ), to maximize
the crossbar utilization in the aggregation core, both the aggregation and
feature extraction cores could work in parallel.

Figure 3. IMA-GNN’s hardware mapping and acceleration in traversal core: (a)
Sample graph adjacency matrix, (b) CSR representation, (c) Search CAM
operation, (d) Scan CAM operation.

##  3\. Network Modeling

We explore both centralized and decentralized GNN landscapes to fairly model
the performance metrics in both settings. Figure 4 shows a sample graph with
$N$ ($N\in\mathbb{Z}^{+}$) nodes (edge devices). In the centralized GNN
setting, a single powerful node as the accelerator is designed with embedded
traversal, aggregation, and feature extraction cores to communicate through
fast inter-network links ($L_{n}$) (mannoni2019comparison, ) to aggregate all
edge devices’ information and handle the computation burdens of
transformation. These cores have $M_{1}$, $M_{2}$, and $M_{3}$ times larger
allocated computing hardware for traversal, aggregation, and feature
extraction operations respectively than a single node in the decentralized
mode. Thus, we assume the processing capability of the edge device in the
centralized setting is $M_{1}$, $M_{2}$, and $M_{3}$ times larger than the
processing capability of a single node in the decentralized mode in the
traversal, aggregation, and feature extraction operations, respectively. In
the decentralized GNN setting, each edge device is observed as an accelerator
with reduced traversal and aggregation cores and in addition to a copy of our
network, has an embedded feature extraction core processing $L$ layers. The
output of the feature extraction core at each edge device is only communicated
to the adjacent edge devices at a defined cluster as shown in Fig. 4(b).
Therefore, the communication between neighbors through inter-cluster links
($L_{c}$) (miya2021experimental, ) generates a communication volume as well.
The bidirectional communication volume between node-$i$ to node-$j$ is
represented as $e_{i,j}$. Therefore, the minimization of the accelerator’s
computational latency/power and communication latency/power is a pivotal need
in the research community. We estimate the centralized and decentralized GNN
accelerators’ latency as:

(1) $\small T_{Net}(N)=T_{compute}(N)+T_{communicate}(N).$

In an N-edge device graph shown in Fig. 4, denoting by $t_{1}$, $t_{2}$, and
$t_{3}$ the traversal, aggregation, and feature extraction cores’ latency,
respectively, the computation latency of a single node in the decentralized
GNNs can be estimated by:

(2) $\small T_{compute-decentralized}=t_{1}+t_{2}+t_{3},$

Figure 4. Intra- and inter-edge links in a sample (a) centralized versus (b)
decentralized GNN.

where in the centralized setting, considering the processing capability of a
single powerful edge device the computation latency is given by:

(3) $\small T_{compute-
centralized}=({{t_{1}}/M_{1}}+{{t_{2}}/M_{2}}+{{t_{3}}/M_{3}})\times(N-1).$

In the decentralized setting, the communication latency, $T_{communicate}$,
can be given by:

(4) $\small T_{communicate-centralized}=(t_{e}+(c_{s}\times t(L_{c})))\times
2,$

where $t_{e}$ is the required time for establishing a connection between two
adjacent nodes, $c_{s}$ denotes the number of adjacent nodes inside a cluster,
$t(L_{c})$ denotes the latency of the inter-cluster link, and number 2 is to
model a two-way link. We assume that data communication inside each cluster is
done in a sequential way, thus the number of adjacent nodes is multiplied by
the latency of the inter-cluster link. For the centralized setting, we assume
data transfer between the central edge device and nodes is done in a
concurrent way. Therefore, the communication latency, $T_{communicate}$, for
centralized inference can be given by:

(5) $\small T_{communicate-decentralized}=t(L_{n}),$

where $t(L_{n})$ is the latency of the inter-network link. Suppose each edge
device runs a GNN with X-layers, the number of input and output neuron
activations for a layer $x$ for 1$\leqslant x\leqslant$X can be given by
$\alpha(x)$ and $\alpha(x+1)$, respectively. The total power consumption of
GNNs implemented in the proposed accelerator can be developed as:

(6) $\small P_{Net}(N)=P_{compute}(N)+P_{communicate}(N).$

The first part accounts for the computation power and the second part
considers communication power for inter-network and inter-cluster links. In
the centralized setting, $P_{compute-centralized}$ is given by
$\frac{E_{compute-centralized}}{T_{compute-centralized}}$ and $P_{communicate-
centralized}$ can be given by $p(L_{n})\times 2$. Here $p(L_{n})$ denotes the
power consumption of the inter-network link and number 2 is to model a two-way
transfer. $E_{compute-centralized}$ is readily calculated by achieving the
energy consumption values of traversal, aggregation, and feature extraction
cores. As for the decentralized setting, $P_{compute-decentralized}$ can be
computed with respect to energy and latency parameters. We consider
$\sum_{n=1}^{c_{s}}p_{n}(c_{s}(n)(c_{s}(n)-1))$ transactions between all
accelerators inside the cluster. Considering the X-layer GNN,
$P_{communicate}$ can be expressed as follow:

(7) $\small P_{communicate-
decentralized}=\frac{1}{t(L_{c})}\times\sum_{x=1}^{X-1}\alpha(x+1)\times
E_{perBit}$

##  4\. Experiments

###  4.1. Evaluation Framework

To evaluate the performance of the proposed architecture, a comprehensive
bottom-up evaluation framework is developed as depicted in Fig. 5. At the
circuit-level, we use the SPICE model for memristors with the Ag-Si memristor
device parameters from (gao2012analog, ). We then combine the SPICE models of
CMOS transistors and memristors under NCSU 45nm CMOS PDK (NCSU_PDK, ) to fully
design and verify the IMA-GNN cores in HSPICE and to extract performance
parameters such as delay and power consumption. We use the Synopsys Design
Compiler (DC, ) to develop the controller and buffer array using a standard
industry-level 45nm technology. At the architecture-level, we modify and
configured the NVSIM-CAM (li2016nvsim, ) memory evaluation tool and MNSIM
(zhu2020mnsim, ) with our circuit-level results to extract the performance
parameters for traversal, aggregation, and feature extraction cores. The
results are then fed to an in-house MATLAB code with the graphs taken as input
to calculate the estimated latency and power consumption for various
workloads. In this study, to have a fair comparison between GNN settings, we
set up IMA-GNN’s traversal, aggregation, and feature extraction cores with
2K$\times$(512$\times$32), 1K$\times$(512$\times$512), and
256$\times$(128$\times$128), respectively, for the centralized setting and
512$\times$32, 512$\times$512, and 128$\times$128, respectively, for the
decentralized setting. It is noteworthy that in addition to the size and
number of crossbars in each core, there are several factors that determine the
total latency and power consumption of IMA-GNN in each setting, such as the
distribution of graph edges across nodes, the availability of graph data, on-
chip storage, and off-chip data accesses.

---
Figure 5. Circuit-to-application evaluation framework.

###  4.2. Traffic Demand Forecasting

As a case study on the potential of hardware-accelerated decentralized GNNs in
real-world applications, we pick a recent work on city-wide multi-relational
and spatiotemporal taxi demand and supply forecasting (nazzal2023traffic, ).
Figure 6(a) and (b) show sample taxis in a city region and their corresponding
graph representation, respectively. This graph is composed of taxi nodes
linked by three edge types; road connectivity, location proximity, and
destination similarity edges linking taxis connected by a road, being nearby,
and targeting nearby destinations, respectively. For each taxi node, the
objective is to predict the values of transportation demand and supply for a
region surrounding it. This is done based on both historical demands and
supplies in the node’s surrounding region and the corresponding messages
shared by other connected nodes (taxis). Formally, the objective is to train a
GNN operator $\mathcal{F}$ that predicts the $Q$-long future passenger demands
and supplies in an $m$ by $n$ region surrounding each taxi at a time instant
$t$, $\mathbf{X}_{t+1:t+Q}\in\mathbf{R}^{Q\times m\times n}$. This is based on
the $P$ historical values of the demands and supplies in this region, and the
corresponding $P$ historical values passed from the $k$-hop neighboring nodes.
So,
$\left[\mathbf{X}_{t-P+1:t},G\right]\stackrel{{\scriptstyle\mathcal{F}}}{{\longrightarrow}}\mathbf{X}_{t+1:t+Q},$
where $\mathbf{X}_{t-P+1:t}\in\mathbf{R}^{P\times m\times n}$, and $G$ denotes
the node’s computational graph. Due to the existence of multiple edge types
and the time-dependency, $\mathcal{F}$ is composed of a heterogeneous GNN
(hetGNN) used for message passing, followed by a Long-Short-Term Memory (LSTM)
network to incorporate time dependency, as depicted in Fig. 7.

|
---|---
(a) | (b)

Figure 6. (a) Taxi representation as a graph, and (b) Decentralized GNN
operation.

Figure 7. The architecture of the hetGNN-LSTM based prediction in
(nazzal2023traffic, ).

The enormous sizes of transportation graphs make it challenging to apply the
model in (nazzal2023traffic, ) in centralized GNN inference. This limitation
is further aggravated by the hetGNNs and the huge volumes of local node
information. As a remedy to resolve this limitation, the authors in
(nazzal2023traffic, ) propose a decentralized GNN inference approach. In this
approach, each taxi node has a copy of the model (hetGNN-LSTM), exchanges
messages with its $k$-hop neighbors, and then uses the hetGNN-LSTM model to
predict the demands and supplies in its surrounding region. A natural
advantage of this approach is handling dynamically varying graph structures.
Nevertheless, despite the promising advantages of decentralization, there is
still a demanding need for reducing the overall computation and communication
latency in the operation of the model.

In this experiment, for the centralized GNN setting, the overall latency (in
terms of transmission delay) for sending and receiving a packet of 300 Bytes
is considered 1.1 ms where the range of the network is 300 meters
(mannoni2019comparison, ). This latency is the average overall latency to
correctly receive a packet of 300 bytes. Thus, for a packet size of 864 bytes,
which is the size of our data, the overall transmission delay can be $\sim$3.3
ms.

As for the decentralized GNN setting, we assume that nodes in the graphs of
(nazzal2023traffic, ) communicate with each other using an ad-hoc wireless
network that uses channel 9 (2.452 GHz) of IEEE 802.11n, where the
transmission power is fixed to –31 dBm, and bandwidth is 20 MHz. In this
configuration, source nodes feed their messages to nearby proxy (relay) nodes
which forward the messages to the next nodes, and so on. Since source nodes
have more computation compared to proxy (relay) nodes, they incur more delay.
We leverage the proposed bottom-up evaluation framework to estimate IMA-GNN
architecture performance in both centralized and decentralized settings. In
our evaluation, the number of nodes (taxis) of the graph and the cluster size
($c_{s}$) are set as 10000 and 10, respectively. The evaluation results of the
latency and power consumption are tabulated in Table 1.

Table 1. Computation and communication latency/power of IMA-GNN accelerator.

Settings | Centralized | Decentralized |
---|---|---|---
Figure of merits | Latency | Power | Latency | Power |
Traversal | 38.43 ns | 10.8 mW | 7.68 ns | 0.21 mW |
Aggregation | 142.77 $\mu$s | 780.1 mW | 14.27 $\mu$s | 41.6 mW |
Feature extraction | 14.53 $\mu$s | 32.21 mW | 0.37 $\mu$s | 3.68 mW |
Computation (Net) | 157.34 $\mu$s | 823.11 mW | 14.6 $\mu$s | 45.49 mW |
Communication | 3.30 ms | - | 406 ms | - |

In view of the results, the decentralized setting in (nazzal2023traffic, )
improves the total computation latency by a factor of $\sim$10$\times$. This
is achieved by reducing the traversal latency, aggregation latency, and
feature extraction latency by factors of 5$\times$, 10$\times$, and
$\sim$39$\times$, respectively. Thus, a huge improvement can be observed in
terms of computation latency. However, in terms of communication, the
centralized setting acts much better than the decentralized setting by
incurring a $\sim$120$\times$ less latency. As for computation power
consumption, we observe that the decentralized setting reduces the power
budget per node by a factor of 18$\times$. The aggregation core of IMA-GNN
consumes most of the power in both centralized and decentralized settings as
well as the highest latency. Overall, each of the settings has its own
advantage from a different point of view. In the next subsection, more graph
datasets with different characteristics are studied in order to further
elucidate the pros and cons of centralized and decentralized settings.

###  4.3. Graph Datasets

In the second case study, four graph datasets, LiveJournal, Collab, Cora, and
Citeseer (yan2019alleviating, ; challapalle2021crossbar, ) are used to
evaluate the inference latency using the proposed IMA-GNN in centralized and
decentralized settings. Key graph statistics of these datasets are provided in
Table 2. A given vertex is mapped deterministically to a fixed-sized, uniform
sample of its neighbors.

Table 2. Key statistics of the graph datasets used.

Datasets | LiveJournal | Collab | Cora | Citeseer
---|---|---|---|---
Number of Nodes | 4,847,571 | 372,475 | 2708 | 3,327
Number of Edges | 68,993,773 | 24,574,995 | 5429 | 4,732
Feature Length | 1 | 496 | 1433 | 3,703
Average $C_{s}$ | 9 | 263 | 4 | 2

Figure 8 shows the latency for the four aforementioned graph datasets. Each
bar consists of two parts; computation latency and communication latency. For
each dataset, we have two bars representing the latency in the centralized
(left) and decentralized (right) settings. By close observation, it can be
realized that in all under-test datasets, the computation latency of the
decentralized setting is less than that of the centralized setting.
Especially, the difference between the computation latency of centralized and
decentralized settings is huge in LiveJournal and Collab datasets, where the
graph size is much larger than the other two datasets. In view of Fig. 8,
amongst the four datasets, LiveJournal has the largest computation latency in
the centralized settings because it owns the largest number of nodes. This is
because in the decentralized mode, as each node is responsible to do its’ own
computation, the computation latency is independent of the total number of
nodes and doesn’t increase when the number of nodes grows. Conversely, in the
centralized setting, by growing the number of nodes, the computation burden of
the edge device will increase, causing an increment in the computation
latency. On average for these four datasets, the decentralized setting
performs computations $\sim$1400$\times$ faster. However, the communication
latency which is the dominant part of the total latency is much higher in the
decentralized setting as each node is required to establish a peer-to-peer
connection and transfer data with all adjacent nodes sequentially. According
to Fig. 8, Collab has the largest communication latency amongst the four
datasets, in the decentralized settings due to its large Average $C_{s}$,
where each node is required to communicate with a large number of adjacent
nodes sequentially. However, in the centralized mode, all nodes are connected
to the edge device using a fast and mature connection and can transfer data in
a parallel way. In terms of the communication latency, for four under-test
datasets, the centralized setting is $\sim$790$\times$ faster than the
decentralized setting.

Figure 8. Breakdown of communication and computation latency of under-test
graph datasets for centralized and decentralized settings.

In addition, we observe that the performance of the IMA-GNN architecture can
increase linearly with an increase in the number of resistive CAM and MVM
crossbars in decentralized setting for various datasets and saturate once the
entire node feature data could be fitted onto the crossbars. However, it comes
at the cost of higher power consumption for each node.

##  5\. Conclusion

While the respective benefits of centralized and decentralized GNNs are known
in software implementation, there is a lack of hardware implementation
analysis to show the communication and computation loads in each setting. This
work undertakes this task by modeling and analyzing practical case studies on
GNN-based taxi demand and supply prediction and adopting large-scale graph
datasets. Our cross-layer simulation results demonstrate our proposed platform
called IMA-GNN in the centralized GNN setting can obtain $\sim$790$\times$
communication speed-up compared to the decentralized GNN setting. However, the
decentralized GNN setting performs computation $\sim$1400$\times$ faster while
reducing the power consumption per device. This study is conducted based on
certain assumptions as discussed. Nevertheless, the results extrapolate that
the decentralized GNN setting achieves gains in reducing the computation
latency. However, this comes at the expense of increasing communication
overhead and latency. This latency is more strongly pronounced with larger
graphs. This finding confirms the necessity and the potential of balancing
this communication-computation trade-off through a semi-decentralized setting
(nazzal2023traffic, ). In that setting, multiple edge devices are employed to
decentralize the operation on a graph level while each edge device region
works in a centralized fashion. Our future work will consider the hardware
acceleration of the semi-decentralized GNN setting.

## References

  * (1) W. L. Hamilton, R. Ying, and J. Leskovec, “Representation learning on graphs: Methods and applications,” _arXiv preprint arXiv:1709.05584_ , 2017.
  * (2) S. Sahu, A. Mhedhbi, S. Salihoglu, J. Lin, and M. T. Özsu, “The ubiquity of large graphs and surprising challenges of graph processing,” _Proceedings of the VLDB Endowment_ , vol. 11, no. 4, pp. 420–431, 2017.
  * (3) K. Kiningham, P. Levis, and C. Ré, “Grip: A graph neural network accelerator architecture,” _IEEE Transactions on Computers_ , 2022.
  * (4) M. Yan, L. Deng, X. Hu, L. Liang, Y. Feng, X. Ye, Z. Zhang, D. Fan, and Y. Xie, “Hygcn: A gcn accelerator with hybrid architecture,” in _HPCA_. IEEE, 2020, pp. 15–29.
  * (5) W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec, “Open graph benchmark: Datasets for machine learning on graphs,” _Advances in neural information processing systems_ , vol. 33, pp. 22 118–22 133, 2020.
  * (6) D. Zheng, C. Ma, M. Wang, J. Zhou, Q. Su, X. Song, Q. Gan, Z. Zhang, and G. Karypis, “Distdgl: distributed graph neural network training for billion-scale graphs,” in _2020 IEEE/ACM 10th Workshop on Irregular Applications: Architectures and Algorithms (IA3)_. IEEE, 2020, pp. 36–44.
  * (7) J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, and M. Sun, “Graph neural networks: A review of methods and applications,” _AI open_ , vol. 1, pp. 57–81, 2020.
  * (8) W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * (9) S. K. Mandal, G. Krishnan, A. A. Goksoy, G. R. Nair, Y. Cao, and U. Y. Ogras, “Coin: Communication-aware in-memory acceleration for graph convolutional networks,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , vol. 12, no. 2, pp. 472–485, 2022.
  * (10) J. Li, A. Louri, A. Karanth, and R. Bunescu, “Gcnax: A flexible and energy-efficient accelerator for graph convolutional neural networks,” in _HPCA_. IEEE, 2021, pp. 775–788.
  * (11) M. Abedin, A. Roohi, M. Liehr, N. Cady, and S. Angizi, “Mr-pipa: An integrated multilevel rram (hfox)-based processing-in-pixel accelerator,” _IEEE Journal on Exploratory Solid-State Computational Devices and Circuits_ , vol. 8, no. 2, pp. 59–67, 2022.
  * (12) S. Angizi, S. Tabrizchi, and A. Roohi, “Pisa: A binary-weight processing-in-sensor accelerator for edge image processing,” _arXiv preprint arXiv:2202.09035_ , 2022.
  * (13) S. Liang, Y. Wang, C. Liu, L. He, L. Huawei, D. Xu, and X. Li, “Engn: A high-throughput and energy-efficient accelerator for large graph neural networks,” _IEEE Transactions on Computers_ , vol. 70, no. 9, pp. 1511–1525, 2020.
  * (14) T. Geng, A. Li, R. Shi, C. Wu, T. Wang, Y. Li, P. Haghi, A. Tumeo, S. Che, S. Reinhardt _et al._ , “Awb-gcn: A graph convolutional network accelerator with runtime workload rebalancing,” in _MICRO_. IEEE, 2020, pp. 922–936.
  * (15) N. Challapalle, K. Swaminathan, N. Chandramoorthy, and V. Narayanan, “Crossbar based processing in memory accelerator architecture for graph convolutional networks,” in _2021 IEEE/ACM International Conference On Computer Aided Design (ICCAD)_. IEEE, 2021, pp. 1–9.
  * (16) Q. Li, F. Gama, A. Ribeiro, and A. Prorok, “Graph neural networks for decentralized multi-robot path planning,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 11 785–11 792.
  * (17) E. Tolstaya _et al._ , “Learning decentralized controllers for robot swarms with graph neural networks,” in _Conference on robot learning_. PMLR, 2020, pp. 671–682.
  * (18) K. Lu, Z. Li, L. Liu, J. Wang, S. Yin, and S. Wei, “Redesk: A reconfigurable dataflow engine for sparse kernels on heterogeneous platforms,” in _ICCAD_. IEEE, 2019, pp. 1–8.
  * (19) V. Mannoni _et al._ , “A comparison of the v2x communication systems: Its-g5 and c-v2x,” in _VTC2019-Spring_. IEEE, 2019, pp. 1–5.
  * (20) T. Miya _et al._ , “Experimental analysis of communication relaying delay in low-energy ad-hoc networks,” in _CCNC_. IEEE, 2021, pp. 1–2.
  * (21) L. Gao _et al._ , “Analog-input analog-weight dot-product operation with ag/a-si/pt memristive devices,” in _VLSI-SoC_. IEEE, 2012, pp. 88–93.
  * (22) (2011) Ncsu eda freepdk45. [Online]. Available: http://www.eda.ncsu.edu/wiki/FreePDK45:Contents
  * (23) S. D. C. P. V. . Synopsys, Inc.
  * (24) S. Li, L. Liu, P. Gu, C. Xu, and Y. Xie, “Nvsim-cam: a circuit-level simulator for emerging nonvolatile memory based content-addressable memory,” in _2016 IEEE/ACM International Conference on Computer-Aided Design (ICCAD)_. ACM, 2016, pp. 1–7.
  * (25) Z. Zhu, H. Sun, K. Qiu, L. Xia, G. Krishnan, G. Dai, D. Niu, X. Chen, X. S. Hu, Y. Cao _et al._ , “Mnsim 2.0: A behavior-level modeling tool for memristor-based neuromorphic computing systems,” in _Proceedings of the 2020 on Great Lakes Symposium on VLSI_ , 2020, pp. 83–88.
  * (26) M. Nazzal, A. Khreishah, J. Lee, and S. Angizi, “Semi-decentralized inference in heterogeneous graph neural networks for traffic demand forecasting: An edge-computing approach,” _arXiv preprint arXiv:2303.00524_ , 2023.
  * (27) M. Yan, X. Hu, S. Li, A. Basak, H. Li, X. Ma, I. Akgun, Y. Feng, P. Gu, L. Deng _et al._ , “Alleviating irregularity in graph analytics acceleration: A hardware/software co-design approach,” in _Micro_ , 2019, pp. 615–628.


# Phaseless Super-resolution in the Continuous Domain

###### Abstract

Phaseless super-resolution refers to the problem of super-resolving a signal
from only its low-frequency Fourier magnitude measurements. In this paper, we
consider the phaseless super-resolution problem of recovering a sum of sparse
Dirac delta functions which can be located anywhere in the continuous time-
domain. For such signals in the continuous domain, we propose a novel
Semidefinite Programming (SDP) based signal recovery method to achieve the
phaseless super-resolution. This work extends the recent work of Jaganathan et
al. [1], which considered phaseless super-resolution for discrete signals on
the grid.

Index Terms—  super-resolution, microscopy, phaseless, continuous domain,
atomic norm

##  1 Introduction

In engineering and science, improving the accuracy and precision of
measurement tools, such as microscopy, X-ray crystallography and MRI, is of
great interest. However, due to the physical limitations in measurement tools,
sometimes we can only indirectly or partially observe a signal of interest,
e.g., obtaining only low-frequency information, only low-resolution image, or
only the magnitude of a signal. The microscope is a good example of a
measurement tool having such physical limitations ranging from low-frequency
measurements to phaseless measurements [2, 3, 4, 5, 6].

To overcome the limitation of low-frequency measurements, researchers have
investigated recovering a signal from only its low-frequency Fourier
measurements, and referred to it as super-resolution. The authors in [7] and
[8] proposed SDP based methods for the recovery of signals in the continuous
domain under certain separation conditions, by employing Total Variation Norm
Minimization (TVNM) and Atomic Norm Minimization (ANM) respectively. Besides,
to address the issue of phaseless measurements, people have studied the phase
retrieval to obtain phase information from the magnitude measurements of a
signal [9, 10]. Recently, in [11], the authors proposed a trace-norm
minimization to solve the phase retrieval problem with the use of masks.

Super-resolving a signal from only magnitudes of low-frequency Fourier
measurements is often ill-posed due to lack of both phase information and
high-frequency information; and hence it is a challenging problem. The authors
in [12, 1] considered the phaseless super-resolution aiming at recovering
signals with only low-frequency magnitude measurements. In the noiseless
setting, the authors in [12] proposed a combinatorial algorithm for signal
recovery using only low-frequency Fourier magnitude measurements, but this
algorithm requires additional distinguishing conditions on the signal
impulses. In the noisy setting, this combinatorial algorithm suffers from
error propagation. Instead of assuming the distinguishing conditions on
signals, the authors in [1] used masks to obtain different types of magnitude
measurements. The authors provably showed that under appropriate choice of
masks, an SDP formulation can be used to recover time-domain impulse signals
on the discretized grid.

In this paper, we consider super-resolving time-domain impulse signals located
off the grid from only low-frequency Fourier magnitude measurements. To tackle
the continuous parameter domain, we propose a novel SDP formulation, employing
ANM to recover signals from Fourier magnitude measurements. For example, our
approach applies to the magnitude measurements used in [11, 1]. In numerical
experiments, we show the successful recovery of signals in the continuous
domain from only low-frequency magnitude measurements. Furthermore, we compare
our algorithm to a simple combining algorithm performing phase retrieval
followed by ANM. Our method shows better recovery performance than the simple
combining algorithm. In the future work, we will consider the noisy magnitude
measurement case.

Notations: In this paper, we denote the set of complex numbers as
$\mathbb{C}$. We reserve calligraphic uppercase letters for index sets, e.g.,
$\mathcal{N}$. We use $|\mathcal{N}|$ as the cardinality of the index set
$\mathcal{N}$. We use the superscripts $*$, $T$, and $H$ to denote conjugate,
transpose, and conjugate transpose respectively. We reserve $i$ for the
imaginary number, i.e., $i^{2}=-1$. We denote a time-domain signal as a
lowercase letter, and its frequency-domain signal as its uppercase letter. To
denote a ground true signal, we use the superscript $o$, e.g., $x^{o}$. For
the index of a vector and a matrix, we start with the index $0$; hence, we
denote the first element of the vector $X$ as $X_{0}$, and the top-left
element of a matrix $Q$ as $Q_{0,0}$.

##  2 Problem Formulation and Background

Let $x^{o}(t)$ be a sum of Dirac functions expressed as follows:

$\displaystyle x^{o}(t)=\sum_{j=1}^{k}c^{o}_{j}\delta(t-t^{o}_{j}),$ (2.1)

where $\delta(t)$ is the Dirac delta function, $c^{o}_{j}\neq 0\in\mathbb{C}$,
and $t^{o}_{j}\in[0,1)$. Its Fourier transform is given by:

$\displaystyle X^{o}_{f}=\sum_{j=1}^{k}c^{o}_{j}e^{-i2\pi
ft^{o}_{j}}=\sum_{j=1}^{k}|c^{o}_{j}|a(t^{o}_{j},\phi^{o}_{j})_{f},\;f\in\mathcal{N},$
(2.2)

where $f\in\mathcal{N}=\\{0,1,...,n-1\\}$,
$a(t^{o}_{j},\phi^{o}_{j})\in\mathbb{C}^{|\mathcal{N}|}$ is an atom vector,
with the $f$-th element given by $a(t^{o}_{j},\phi^{o}_{j})_{f}=e^{-i(2\pi
ft^{o}_{j}-\phi^{o}_{j})}$. Simply, $X^{o}=V^{o}c^{o}$, where
$X^{o}\in\mathbb{C}^{n}$, $V^{o}=[a(t^{o}_{1},0),...,a(t^{o}_{k},0)]$, and
$c^{o}=[|c^{o}_{1}|e^{i\phi^{o}_{1}},...,|c^{o}_{k}|e^{i\phi^{o}_{k}}]^{T}$.
We also define the minimum separation of $x^{o}(t)$, denoted by
$\bigtriangleup_{t}$, as the closest distance between any two different time
value $t^{o}_{j}$’s in cyclic manner [7, 8], i.e.,

$\displaystyle\Delta_{t}=\underset{t^{o}_{i},t^{o}_{j}\in[0,1),\;i\neq
j}{\min}\;\;|t^{o}_{i}-t^{o}_{j}|.$ (2.3)

The goal here is to find $x^{o}(t)$ from the low-frequency Fourier magnitude
measurements. We state the phaseless super-resolution problem with masks as
follows:

$\displaystyle\text{Find}\;\;x(t)$ $\displaystyle\text{subject
to}\;\;Z[r,l]=\bigg{|}\int_{0}^{1}D_{r}(t)\sum_{j=1}^{k}c^{o}_{j}\delta(t-t^{o}_{j})e^{-i2\pi
lt}dt\bigg{|},$ (2.4) $\displaystyle\quad\quad\quad\text{for}\;{-R\leq r\leq
R}\;\text{and}\;l\in\mathcal{N},$

where $Z[r,l]$ is the $l$-th frequency magnitude obtained by using the $r$-th
mask function $D_{r}(t)$. Depending on the mask function $D_{r}(t)$, one can
have different types of magnitude information. For example, if we choose
$1+e^{-i2\pi t}$ for $D_{r}(t)$, we have $|X^{o}_{l}+X^{o}_{l+1}|$,
$l\in\mathcal{N}$.

In [1], Jaganathan et al. consider the case when the signal $x^{o}(t)$ is
located on the grid, i.e., $t^{o}_{j}\in\\{0,1,2,....n-1\\}$. By $n$-point
DFT, (2) is equivalent to

$\displaystyle\text{Find}\;\;x$ $\displaystyle\text{subject
to}\;\;Z[r,l]=|\langle f_{l},D_{r}x\rangle|,$ (2.5)
$\displaystyle\quad\quad\quad\text{for}\;{-R\leq r\leq
R}\;\text{and}\;l\in\mathcal{N},$

where $x\in\mathbb{C}^{n}$ is a complex valued $k$-sparse vector,
$D_{r}\in\mathbb{C}^{n\times n}$ is a diagonal matrix, and $f_{l}$ is the
conjugate of the $l$-th column of the $n$ point DFT matrix. The authors in [1]
proposed the following semidefinite relaxation-based program for the phaseless
super-resolution in the discrete domain by denoting $Y=xx^{H}$ and relaxing
the rank-1 constraint on $Y$:

$\displaystyle\underset{Y}{\text{minimize}}\;\;||Y||_{1}+\lambda\text{Tr}(Y)$
$\displaystyle\text{subject
to}\;\;Z[r,l]^{2}=\text{Tr}(D_{r}^{H}f_{l}f_{l}^{H}D_{r}Y),$ (2.6)
$\displaystyle\quad\quad\quad\text{for}\;{-R\leq r\leq
R},\;l\in\mathcal{N},\;\text{and}\;Y\succeq 0,$

for some $\lambda>0$.

This paper makes no assumption of $t^{o}_{j}$ being on the grid. In the next
section, we propose an ANM based semidefinite relaxation of (2) to deal with
impulse functions off the grid.

##  3 Phaseless Super-Resolution in the Continuous Domain

We define the atomic norm of a vector $X\in\mathbb{C}^{|\mathcal{N}|}$ as
follows:

$\displaystyle||X||_{\mathcal{A}}=\inf\\{\sum_{j}|c_{j}|:X_{l}=\sum_{j}|c_{j}|a(t_{j},\phi_{j})_{l},\;\begin{subarray}{c}{t_{j}\in[0,1),}\\\
{\phi_{j}\in[0,2\pi)}\end{subarray}\\}.$ (3.1)

We have the following new proposition for the atomic norm:

######  Proposition 3.1.

For any $X\in\mathbb{C}^{|\mathcal{N}|}$, $\mathcal{N}=\\{0,1,...,n-1\\}$,

$\displaystyle||X||_{\mathcal{A}}^{2}=\inf_{u,s}\bigg{\\{}\frac{1}{|\mathcal{N}|}\text{$s$Tr(Toep($u$))}\;:\;\begin{bmatrix}\text{Toep($u$)}&X\\\
X^{H}&s\end{bmatrix}\succeq 0\bigg{\\}},$ (3.2)

where $\text{Tr}(\cdot)$ is the trace operator, and Toep($u$) is the Toeplitz
matrix whose first column is $u=[u_{0},u_{1},...,u_{n-1}]^{T}$. Moreover,
suppose after the Vandermonde decomposition [13, 14, 15],
$\text{Toep($u$)}=VDV^{H}$, where $V=[a(t_{1},0),...,a(t_{r},0)]$ and $D$ is a
positive diagonal matrix. Then, there exists a vector $c$ such that $X=Vc$ and
$\sum_{j}|c_{j}|=||X||_{\mathcal{A}}$.

Proposition 3.1 is similar to Proposition II.1 in [8]; however, Proposition
3.1 considers the trace of $s\text{Toep($u$)}$ instead of the sum of trace of
Toep($u$) and $s$. Proposition 3.1 is essential to derive our new SDP
formulation handling phaseless measurements. For readability, we place the
proof of Proposition 3.1 in Appendix.

Motivated by Proposition 3.1, we propose the following squared atomic norm
minimization for the phaseless super-resolution in the continuous domain,
simply phaseless ANM:

$\displaystyle\underset{X}{\text{minimize}}\;\;||X||_{\mathcal{A}}^{2}$
$\displaystyle\text{subject to}\;\;a_{r}(X)=b_{r},\;r=1,2,...,q,$ (3.3)

where $q$ is the total number of magnitude measurements, $a_{r}(X)$ is the
magnitude mapping function, $|\langle a_{r},X\rangle|$,
$a_{r}\in\mathbb{C}^{|\mathcal{N}|}$, and $b_{r}$’s are magnitude measurement
results.

From Proposition 3.1, we can change (3) to

$\displaystyle\underset{u,X,s}{\text{minimize}}\;\;\frac{1}{|\mathcal{N}|}\text{$s$Tr(Toep($u$))}$
$\displaystyle\text{subject
to}\;\;U\triangleq\begin{bmatrix}\text{Toep($u$)}&X\\\
X^{H}&$s$\end{bmatrix}\succeq 0,$
$\displaystyle\quad\quad\quad\quad\;\;a_{r}(X)=b_{r},\;r=1,2,...,q,$ (3.4)

where $u,X\in\mathbb{C}^{|\mathcal{N}|}$ and $s\in\mathbb{C}$. From the
positive semidefiniteness of $U$, $s\geq 0$, and $\text{Toep($u$)}\succeq 0$.
Besides, if $X_{j}\neq 0$, $j\in\mathcal{N}$, then $s\neq 0$ from the non-
negativeness of all principal minors of $U$ [16]. However, because of the
magnitude constraints, (3) is a non-convex program.

By the Schur complement lemma [17], $U\succeq 0$ implies
$\text{$s$Toep($u$)}-XX^{H}\succeq 0$. Since
$s\text{Toep($u$)}=\text{Toep($su$)}$, by defining $Q=XX^{H}$ and
$u^{\prime}=su$, and getting rid of the rank constraint on $Q$, we have the
following SDP relaxation for the phaseless ANM:

$\displaystyle\underset{Q\succeq
0,u^{\prime}}{\text{minimize}}\;\;\frac{1}{|\mathcal{N}|}\text{Tr(Toep($u^{\prime}$))}$
$\displaystyle\text{subject to}\;\;\text{Toep($u^{\prime}$)}-Q\succeq 0,$
$\displaystyle\quad\quad\quad\quad\;\;A_{r}(Q)=b^{2}_{r},\;r=1,2,...,q,$ (3.5)

where $A_{r}(Q)$ is a mapping function, $\text{Tr}(A_{r}Q)$. Here,
$A_{r}=a_{r}a_{r}^{H}$.

After solving (3), we can find the optimal $\hat{Q}$ and optimal
Toep($\hat{u}$). Our analysis of (3) in the following section shows that under
certain conditions, $\hat{Q}=X^{o}{X^{o}}^{H}$. We can recover $X^{o}$ up to
global phase by the eigenvalue decomposition of $\hat{Q}$. More importantly,
because of the structure of $\text{Toep($\hat{u}$)}=V^{o}D{V^{o}}^{H}$ for
some diagonal matrix $D$, we can apply any parameter estimation method such as
Prony’s method [18, 19, 20] or a matrix pencil method [21, 22] to find the
time location $t^{o}_{j}$’s.

##  4 Performance Analysis

We first consider the analysis of (3) given a rank-1 matrix $Q$. And then, we
provide the analysis of (3). Finally, we look at one scenario having magnitude
measurements from a set of masks, in which (3) provides the desired signal
recovery.

######  Theorem 4.1.

For a given rank-1 positive semidefinite matrix $Q=XX^{H}$,
$X\in\mathbb{C}^{|\mathcal{N}|}$, the following optimization problem provides
the squared atomic norm of $X$, i.e., $||X||_{\mathcal{A}}^{2}$:

$\displaystyle\underset{u}{\text{minimize}}\;\;\frac{1}{|\mathcal{N}|}\text{Tr(Toep($u$))}$
$\displaystyle\text{subject to}\;\;\text{Toep($u$)}-Q\succeq 0.$ (4.1)

###### Proof.

We can prove it by using Proposition 3.1. Defining $u=u^{\prime}s$, where
$s>0$ is a scalar. Then we can re-state the constraint as
$Toep(u^{\prime})-\frac{1}{s}XX^{H}\succeq 0$. By the Schur complement lemma,
we have the optimization problem in Proposition 3.1. Therefore, from
Proposition 3.1, the optimal value of (4.1) is the same as
$||X||^{2}_{\mathcal{A}}$. ∎

######  Corollary 4.2.

If (3) gives a rank-1 solution to $Q$, then (3) minimizes the squared atomic
norm of $X$ among all vectors $X$ satisfying the given constraints
$a_{r}(X)=b_{r}$, $r=1,2,...,q$.

###### Proof.

From Theorem 4.1, (3) provides the minimum squared atomic norm of $X$ among
all vectors $X$ satisfying the constraints $a_{r}(X)=b_{r}$, $r=1,2,...,q$. ∎

Let us consider the case when we have low-frequency Fourier magnitude
measurements from a set of masks. The main difference between [7, 8] and our
setting is that we have only magnitude measurements, instead of measurements
offering both phases and magnitudes.

######  Theorem 4.3.

Given the magnitude measurements $|X^{o}_{j}|$, $|X^{o}_{j}+X^{o}_{j+1}|$, and
$|X^{o}_{j}-iX^{o}_{j+1}|$, $j\in\mathcal{N}=\\{0,1,...,n-1\\}$, (3) provides
the unique solution $Q=X^{o}{X^{o}}^{H}$, and $x^{o}(t)$ is uniquely obtained
up to global phase if the following conditions hold: $X_{j}\neq 0$, $\forall
j\in\mathcal{N}$, and $\Delta_{t}\geq 4/|\mathcal{N}|$.

###### Proof.

Given magnitude data, $|X^{o}_{j}|^{2}$, $|X^{o}_{j+1}|^{2}$,
$|X^{o}_{j}+X^{o}_{j+1}|^{2}$, and $|X^{o}_{j}-iX^{o}_{j+1}|^{2}$, we can find
$Q_{j,j}$, $Q_{j+1,j+1}$, $Q_{j,j+1}$ and $Q_{j+1,j}$, which are the elements
of the diagonal, sub-diagonal, and super-diagonal of the matrix $Q$, by simply
solving linear equations on $Q$ together. From Lemma 6.1 in Appendix 6.2, we
can uniquely recover $Q=X^{o}{X^{o}}^{H}$ and $X^{o}$ up to global phase.
According to Proposition 3.1 and Theorem 4.1, (3) with $X^{o}$ is essentially
the same as the optimization problem dealing with the standard ANM [8] or TVNM
[7]. Therefore,(3) provides unique $x^{o}(t)$ up to global phase if the
separation condition holds, i.e., $\Delta_{t}\geq 4/|\mathcal{N}|$. ∎

##  5 Numerical Experiments

We compare our phaseless ANM against the standard ANM [8] using measurements
offering both phases and magnitudes, as well as against a simple algorithm
which first performs the phase retrieval [11] and then applies the standard
ANM [8] to recover the impulse functions from the recovered signal using the
phase retrieval. We use CVX [23] to solve (3).

Fig. 1 (a) and (b) show the probability of successful recovery from the
standard ANM and the phaseless ANM respectively. We conducted 50 trials for
each parameter setting and measured the success rate. At each trial, we chose
one time impulse $t^{o}_{1}$ uniformly at random in $[0,1)$, and another time
impulse $t^{o}_{2}$ by adding the separation $\Delta_{t}$ to $t^{o}_{1}$ in
the cyclic manner. We sampled the real part and imaginary part of time
coefficients $c^{o}_{j}$’s uniformly at random in (0,1). We consider low
frequencies, i.e., $\mathcal{M}=\\{0,1,...,m-1\\}$, where $m<n$,
$\mathcal{M}\subseteq\mathcal{N}$. For a set of masks in the phaseless ANM, we
use the same masks as those of Theorem 4.3 over the index set $\mathcal{M}$.
The x-axis represents the separation condition $\Delta_{t}$ varied from $1/n$
to $11/n$, and y-axis is the number of low-frequency Fourier measurements $m$,
varied from $2$ to $30$. In fact, for the phaseless ANM, the number of
magnitude measurements is $3m-2$. We evaluated the recovery performance for
the signal dimension $n=32$. We calculated the Euclidean distance between the
estimated and true time locations. If the distance is less than $10^{-3}$,
then we consider the estimation successful. Numerical experiments show that
our phaseless ANM can find the exact time locations in the continuous domain
with the same performance as the standard ANM. For large $k$, e.g., $k=10$,
our method also provides the same performance as the standard ANM. We omit the
simulation results in this paper due to the space limitation.

Fig. 1: The probability $P$ of successful recovery by varying the separation
condition $\Delta_{t}$ and the number of measurements $m$ when $n=32$. (a)
Standard ANM. (b) Phaseless ANM Fig. 2: The probability $P$ of successful
recovery by varying the number of magnitude measurements $q$ and sparsity $k$
when $n=32$. (a) Phase retrieval and then standard ANM (b) Phaseless ANM

One can think of a simple method conducting the phase retrieval first, and
then doing the standard ANM. To compare our algorithm with this simple method,
we further carried out numerical experiments by varying the number of
magnitude measurements $q$ and the number of sparsity $k$ in (2.1). In this
simulation, instead of using a set of masks used in Theorem 4.3, we randomly
chose a vector $a_{r}$ for each magnitude measurement in (3). Fig. 2 (a) and
(b) show the probability of successful recovery from the simple combining
algorithm and the phaseless ANM respectively. The x-axis is the number of
magnitude measurements $q$, and y-axis is the number of sparsity $k$. With
randomly chosen magnitude measurements, our method outperforms the simple
combining algorithm.

##  6 Appendix

###  6.1 Proof of Proposition 3.1

We follow the proof of [8, Proposition II.1].

###### Proof.

Let us denote the optimal value of the right hand side of (3.2) by SDP($X$).
In order to show $||X||_{\mathcal{A}}^{2}=\text{SDP($X$)}$, we will show that
(1) $||X||_{\mathcal{A}}^{2}\geq\text{SDP($X$)}$ and (2)
$||X||_{\mathcal{A}}^{2}\leq\text{SDP($X$)}$.

The proof of (1) is easily shown by considering a feasible solution of
SDP($X$). For $X=\sum_{j}|c_{j}|a(t_{j},\phi_{j})$, by choosing a feasible
solution,
$\text{Toep($u$)}=\sum_{j}|c_{j}|a(t_{j},\phi_{j})a(t_{j},\phi_{j})^{H}$, and
$s=\sum_{j}|c_{j}|$, we have

$\displaystyle\begin{bmatrix}\text{Toep($u$)}&X\\\
X^{H}&s\end{bmatrix}=\sum_{j}|c_{j}|\begin{bmatrix}a(t_{j},\phi_{j})\\\
1\end{bmatrix}\begin{bmatrix}a(t_{j},\phi_{j})\\\ 1\end{bmatrix}^{H}\succeq
0.$

For this feasible solution,
$\frac{1}{|\mathcal{N}|}s\text{Tr(Toep($u$))}=(\sum_{j}|c_{j}|)^{2}$, which is
$||X||_{\mathcal{A}}^{2}$. Thus, $\text{SDP($X$)}\leq||X||_{\mathcal{A}}^{2}$.

For the proof of (2), we will show that for any $u$, $s$, and $X$,
$\frac{1}{|\mathcal{N}|}s\text{Tr(Toep($u$))}\geq||X||_{\mathcal{A}}^{2}$.
Suppose for some $u$, $s\neq 0$, and $X$, the matrix $U$ in (3) is positive
semidefinite. From the positive semidefinite condition, we have
$\text{Toep($u$)}\succeq 0$ and $s>0$. From the Vandermonde decomposition [13,
14, 15], for any positive semidefinite Toep($u$), we have
$\text{Toep($u$)}=VDV^{H}$, where $V=[a(t_{1},0)\;a(t_{2},0),...a(t_{r},0)]$,
and $D$ is a diagonal matrix having $d_{j}$ as its $j$-th diagonal element.
Since $VDV^{H}=\sum_{j=1}^{r}d_{j}a(t_{j},0)a(t_{j},0)^{H}$ and
$||a(t_{j},0)||^{2}_{2}=|\mathcal{N}|$, we have
$\frac{1}{|\mathcal{N}|}\text{Tr(Toep($u$))}=\text{Tr($D$)}$. Also, from the
Vandermonde decomposition and $U\succeq 0$, $X$ is in the range space of $V$;
in fact, if $X$ is not in the range of $V$, we can always find a vector $z$
such that $z^{H}Uz<0$. Therefore, $X=Vw=\sum_{j=1}^{r}w_{j}a(t_{j},0)$, where
$w\in\mathbb{C}^{r}$. By the Schur complement lemma [17], $U$ in (3) is
expressed as follows:

$\displaystyle VDV^{H}-\frac{1}{s}Vww^{H}V^{H}\succeq 0.$ (6.1)

It is noteworthy that we can always find a vector $q$ such that
$V^{H}q=sign(w)$, where $sign(w)^{H}w=\sum_{j=1}^{r}|w_{j}|$, by choosing
$q=V(V^{H}V)^{-1}sign(w)$. This is because $V^{H}$ has full row rank. By
choosing $q$ such that $V^{H}q=sign(w)$, we have

$\displaystyle\text{Tr($D$)}=q^{H}VDV^{H}q\geq\frac{1}{s}q^{H}Vww^{H}V^{H}q=\frac{1}{s}(\sum_{j}|w_{j}|)^{2},$

where the inequality is from (6.1). Therefore, we have

$\displaystyle\frac{1}{|\mathcal{N}|}s\text{Tr(Toep($u$))}=s\text{Tr($D$)}\geq(\sum_{j}|w_{j}|)^{2}=||X||_{\mathcal{A}}^{2}.$

If $s=0$, from the sufficient and necessary condition for the positive
semidefiniteness of a Hermitian matrix, all of $U$’s principal minors need to
be non-negative [16]. Thus, $X_{j}=0$, $\forall j\in\mathcal{N}$. In this
case, Proposition 3.1 still holds. ∎

###  6.2 Lemma for the positive semidefinite matrix $Q$

######  Lemma 6.1.

Let $Q\in\mathbb{C}^{|\mathcal{N}|\times|\mathcal{N}|}$, and
$X^{o}\in\mathbb{C}^{|\mathcal{N}|}$. Suppose (1) $Q\succeq 0$, (2)
$Q_{j,j}=Q^{o}_{j,j}$, $Q_{j,j+1}=Q^{o}_{j,j+1}$, and
$Q_{j+1,j}=Q^{o}_{j+1,j}$, $j\in\mathcal{N}=\\{0,1,...,n-1\\}$, where
$Q^{o}=X^{o}{X^{o}}^{H}$, (3) $X^{o}_{j}\neq 0$, $\forall j\in\mathcal{N}$.
Then, $Q$ is uniquely determined as $Q=X^{o}{X^{o}}^{H}$.

###### Proof.

From the fact that a Hermitian matrix is positive semidefinite if and only if
all of its principal minors are non-negative [16], all of $Q$’s principal
minors are required to be non-negative. Let us prove our lemma by induction.
When $|\mathcal{N}|=3$, the determinant of $Q$ is
$-|X^{o}_{1}Q_{0,2}-X^{o}_{0}X^{o}_{1}{X_{2}^{o}}^{*}|^{2}$, where $Q_{0,2}$
is unknown. To be
$-|X^{o}_{1}Q_{0,2}-X^{o}_{0}X^{o}_{1}{X_{2}^{o}}^{*}|^{2}\geq 0$,
$X^{o}_{1}Q_{0,2}-X^{o}_{0}X^{o}_{1}{X_{2}^{o}}^{*}=0$. Since $X^{o}_{1}\neq
0$, $Q_{0,2}$ is determined uniquely as $X^{o}_{0}{X_{2}^{o}}^{*}$. When
$|\mathcal{N}|=4$, we can consider the top-left $3\times 3$ submatrix of $Q$
and the bottom-right $3\times 3$ submatrix of $Q$ to determine $Q_{0,2}$ and
$Q_{1,3}$ respectively. And then, we can deal with $3\times 3$ principal
submatrix of $Q$ having $Q_{0,4}$ to determine $Q_{0,4}$. In the similar way,
when $|\mathcal{N}|=n$, we can uniquely determine every unknown variables in
$Q$. We omit the detailed explanation due to the space limitation. ∎

## References

  * [1] K. Jaganathan, J. Saunderson, M. Fazei, Y. C. Eldar, and B. Hassibi,  “Phaseless super-resolution using masks,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2016, pp. 4039–4043.
  * [2] S. Brenner and R. W. Horne,  “A negative staining method for high resolution electron microscopy of viruses,”  Biochimica et biophysica acta, vol. 34, pp. 103–110, 1959.
  * [3] H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang,  “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,”  Nature biotechnology, vol. 24, no. 7, pp. 848–851, 2006.
  * [4] H. Benveniste and S. Blackband,  “MR microscopy and high resolution small animal MRI: applications in neuroscience research,”  Progress in neurobiology, vol. 67, no. 5, pp. 393–420, 2002.
  * [5] E. Stolyarova, K. T. Rim, S. Ryu, J. Maultzsch, P. Kim, L. E. Brus, T. F. Heinz, M. S. Hybertsen, and G. W. Flynn,  “High-resolution scanning tunneling microscopy imaging of mesoscopic graphene sheets on an insulating surface,”  Proceedings of the National Academy of Sciences, vol. 104, no. 22, pp. 9209–9212, 2007.
  * [6] B. Huang, M. Bates, and X. Zhuang,  “Super resolution fluorescence microscopy,”  Annual review of biochemistry, vol. 78, pp. 993, 2009.
  * [7] E. J. Candès and C. Fernandez-Granda,  “Towards a mathematical theory of super-resolution,”  Communications on Pure and Applied Mathematics, vol. 67, no. 6, pp. 906–956, 2014.
  * [8] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht,  “Compressed sensing off the grid,”  IEEE Transactions on Information Theory, vol. 59, no. 11, pp. 7465–7490, 2013.
  * [9] J. R. Fienup,  “Phase retrieval algorithms: a comparison,”  Applied optics, vol. 21, no. 15, pp. 2758–2769, 1982.
  * [10] R. W. Gerchberg and W. O. Saxton,  “A practical algorithm for the determination of phase from image and diffraction plane pictures,”  Optik, vol. 35, no. 2, pp. 237–246, 1972.
  * [11] E. J. Candès, Y. C. Eldar, T. Strohmer, and V. Voroninski,  “Phase retrieval via matrix completion,”  SIAM review, vol. 57, no. 2, pp. 225–251, 2015.
  * [12] Y. Chen, Y. C. Eldar, and A. J. Goldsmith,  “An algorithm for exact super-resolution and phase retrieval,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014, pp. 754–758.
  * [13] C. C. Carathéodory,  “Über ber den variabilitätsbereich der fourierschen konstanten von positiven harmonischen funktionen,”  Rendiconti del Circolo Matematico di Palermo (1884-1940), vol. 32, no. 1, pp. 193–217, 1911.
  * [14] C. C. Carathéodory and L. Fejér,  “Über ber den zusammenhang der extremen von harmonischen funktionen mit ihren koeffizienten und über den picard-landauschen satz,”  Rendiconti del Circolo Matematico di Palermo (1884-1940), vol. 32, no. 1, pp. 218–239, 1911.
  * [15] O. Toeplitz,  “Zur theorie der quadratischen und bilinearen formen von unendlichvienlen veräanderlichen,”  Mathenatische Annalen, vol. 70, no. 3, pp. 351–376, 1911.
  * [16] C. D. Meyer,  Matrix analysis and applied linear algebra, vol. 2,  Siam, 2000.
  * [17] S. Boyd and L. Vandenberghe,  Convex optimization,  Cambridge University Press, 2004.
  * [18] S. M. Kay,  Modern spectral estimation: theory and application,  Prentice Hall, 1988.
  * [19] P. Stoica and R. L. Moses,  Spectral analysis of signals,  Prentice Hall, 2005.
  * [20] T. Blu, P.-L. Dragotti, M. Vetterli, P. Marziliano, and L. Coulot,  “Sparse sampling of signal innovations,”  IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 31–40, 2008\.
  * [21] Y. Hua and T. K. Sarkar,  “Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 38, no. 5, pp. 814–824, 1990.
  * [22] T. K. Sarkar and O. Pereira,  “Using the matrix pencil method to estimate the parameters of a sum of complex exponentials,”  IEEE Antennas and Propagation Magazine, vol. 37, no. 1, pp. 48–55, 1995.
  * [23] M. Grant and S. Boyd,  “CVX: Matlab software for disciplined convex programming, version 2.0 beta,” http://cvxr.com/cvx, Sept. 2012.


††thanks: e-mail<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>
# Light Hypernuclei Production in Au+Au Collisions at $\sqrt{s_{NN}}=$ 3 GeV
within Thermodynamic Approach

M. Kozhevnikova  Veksler and Baldin Laboratory of High Energy Physics, JINR
Dubna, 141980 Dubna, Russia  Yu. B. Ivanov  Bogoliubov Laboratory of
Theoretical Physics, JINR Dubna, 141980 Dubna, Russia  National Research
Center ”Kurchatov Institute”, 123182 Moscow, Russia

###### Abstract

Simulations of the $\Lambda$-hyperon and light-hypernuclei production in Au+Au
collisions at $\sqrt{s_{NN}}=$ 3 GeV were performed within updated Three-fluid
Hydrodynamics-based Event Simulator Extended by UrQMD (Ultra-relativistic
Quantum Molecular Dynamics) final State interactions (THESEUS). The light
(hyper)nuclei are treated thermodynamically, i.e. they are considered on the
equal basis with hadrons. The only additional parameter is related to the late
freeze-out that imitates the afterburner stage for the light (hyper)nuclei
because the UrQMD is not able to dynamically treat them. The calculation of
hypernuclei production is completely similar to that of light nuclei. The
hypernuclei results are compared with recent STAR data. It is found that the
calculated midrapidity ${}_{\Lambda}^{3}$H/$\Lambda$ ratio falls within the
error bars of the experimental point. It is remarkable that large difference
between the $t/p$ and ${}_{\Lambda}^{3}$H/$\Lambda$ ratios is reproduced
without any additional parameters. Rapidity distributions of
${}_{\Lambda}^{3}$H/$\Lambda$ and ${}_{\Lambda}^{4}$He/$\Lambda$ ratios are
predicted. Midrapidity mean transverse momenta of protons, $\Lambda$s and
light (hyper)nuclei in central collisions well agree with the data. The
calculated directed flow also reasonably well reproduces of the data.

relativistic heavy-ion collisions, hydrodynamics, light nuclei

###### pacs:

25.75.-q, 25.75.Nq, 24.10.Nz

##  I Introduction

Hypernuclei are an important topic of nuclear physics. Heavy-ion reactions at
relativistic energies are an abundant source of strangeness and therefore well
suited for the production of light hypernuclei. The interest in studying
hypernuclei in collisions of heavy ions is twofold. Firstly, the heavy-ion
experiments give us information on lifetimes and binding energies of light
hypernuclei, see, e.g., [1, 2, 3, 4]. This allow us to refine out
understanding of hyperon-nucleon interactions and the role of flavor symmetry
that are relevant for nuclear structure [5, 6, 7] and astrophysics [8, 9, 10],
as well as for construction of the hadronic equation of state (EoS) for
applications to heavy-ion collisions [11]. Another aspect of studying
production of hypernuclei is directly related to diagnostics of the quark-
gluon plasma (QGP) formation in heavy-ion collisions. It was suggested that
the strangeness population factor
$S_{3}={{}_{\Lambda}^{3}}H/{{}^{3}}He\cdot{p/\Lambda}$ can serve as a probe of
the baryon number and strangeness correlation in the produced matter because
of its different behaviors in the QGP and hadronic matter [12, 13, 14].

In this paper we focus on discussing the mechanism of hypernuclei formation in
heavy-ion collisions. Similarly to the light-nuclei production, alternative
mechanisms [13, 15, 16, 17, 18, 19] are still actively debated. The
coalescence [12, 13, 17, 18] and the thermodynamic models [20, 21, 13, 18, 19,
22, 23, 24] are two of the most popular alternative approaches. The model of
Parton-Hadron Quantum Molecular Dynamics (PHQMD) [15, 16] is based on specific
procedures of recognition of light (hyper)nuclei. Light (hyper)nuclei in this
model are not dynamic objects, but rather are associated with relatively
stable cluster-like correlations. Therefore, this approach can be viewed as
kinetics of propagation of the correlations. The light nuclei act as dynamic
objects in kinetic models of Refs. [25, 26, 27] (only deuterons) [28] (all
light nuclei up to 4He). However, these dynamical treatments have not been
extended to hypernuclei so far.

As found in Refs. [13, 18], both the coalescence and thermodynamic models
agree in their predictions for the yields of the light (hyper)nuclei. The
thermodynamical approach has an important advantage. It does not need
additional parameters for the light-(hyper)nuclei treatment. It describes the
light nuclei on equal basis with hadrons, i.e. in terms of temperatures and
chemical potentials. Therefore, its predictive power is the same for light
nuclei and hadrons. This approach was first realized within the statistical
model [20]. The statistical model gave a good description of even hypernuclei
and anti-nuclei [21].

Recently data on light-hypernuclei production in Au+Au collisions at
$\sqrt{s_{NN}}=$ 3 GeV were appeared [3, 29, 30, 31], some of them [29, 30]
preliminary. These data were analyzed in Refs. [19, 18, 16]. Simulations were
performed within the PHQMD approach [16], UrQMD+coalescence and UrQMD-
hybrid+coalescence approaches were applied in Ref. [18], JAM+coalescence model
(Jet AA microscopic transport Model [32, 33]) was used in STAR papers and
presentations [3, 29, 30, 31]. These data were also analyzed within the hybrid
dynamical-statistical approach [19], which is a modification of the
Statistical Multifragmentation Model [34].

In the present paper, we extend the light-nuclei treatment of our previous
papers [22, 23, 24] to hypernuclei production at $\sqrt{s_{NN}}=$ 3 GeV. This
treatment is based on the updated THESEUS generator [22], in which the light
(hyper)nuclei production is considered within the thermodynamic treatment. Our
approach is similar to the thermal production from UrQMD hybrid model in Ref.
[13]. We calculate some bulk properties and directed flow of protons,
$\Lambda$s and light (hyper)nuclei and compare them with available STAR data.
Since the aforementioned coalescence approaches provide a quite reasonable
description of the bulk properties, we address the question whether a similar
description can be achieved by less demanding means, i.e. using
thermodynamics.

##  II The THESEUS Generator

The THESEUS event generator [35, 36] is based on the model of the three-fluid
dynamics (3FD) [38, 37] complemented by the UrQMD [39] for the afterburner
stage. The output of the 3FD model consists of fields of local flow velocities
and thermodynamic quantities defined on the freeze-out hypersurface. The
THESEUS generator transforms the 3FD output into a set of observed particles,
i.e. performs the particlization. The particlization is followed by the UrQMD
afterburning stage.

The initial version of the THESEUS [35, 36] produced the primordial nucleons
and hyperons, i.e. both observable nucleons and those bound in the light
(hyper)nuclei. These nucleons and hyperons were intended for the subsequent
use in the coalescence model, just as it was done for the light-nuclei
production in Refs. [38, 40]. After the coalescence, the nucleons and hyperons
bound in the light (hyper)nuclei should be subtracted from the primordial
ones. If both nucleons/hyperons and light (hyper)nuclei are sampled from the
3FD output, i.e. temperature and chemical potential fields obtained from the
EoS, in which the light (hyper)nuclei are not included, the bound
nucleons/hyperons become double counted. This leads to an overestimation of
the total baryon charge in the final state containing both baryons and
clusters. To avoid this double counting in the updated version of THESEUS
[22], the baryon chemical potentials are recalculated proceeding from the
local baryon number conservation in the system of hadrons extended by the
light-(hyper)nuclei species listed in Tab. 1. The list of the light
(hyper)nuclei includes stable nuclei [deuterons ($d$), tritons ($t$), helium
isotopes 3He and 4He], and low-lying 4He resonances decaying into stable
species [41]. As compared to our previous papers [22, 23, 24], here we
additionally include the stable (with respect to strong decays) hypernuclei
${}_{\Lambda}^{3}$H and ${}_{\Lambda}^{4}$He. The anti-(hyper)nuclei are also
included.

Nucleus($E$[MeV]) | $J$ | decay modes, in %
---|---|---
$d$ | $1$ | Stable
$t$ | $1/2$ | Stable
3He | $1/2$ | Stable
4He | $0$ | Stable
4He(20.21) | $0$ |  $p$ = 100
4He(21.01) | $0$ |  $n$ = 24, $p$ = 76
4He(21.84) | $2$ |  $n$ = 37, $p$ = 63
4He(23.33) | $2$ |  $n$ = 47, $p$ = 53
4He(23.64) | $1$ |  $n$ = 45, $p$ = 55
4He(24.25) | $1$ |  $n$ = 47, $p$ = 50, $d$ = 3
4He(25.28) | $0$ |  $n$ = 48, $p$ = 52
4He(25.95) | $1$ |  $n$ = 48, $p$ = 52
4He(27.42) | $2$ |  $n$ = 3, $p$ = 3, $d$ = 94
4He(28.31) | $1$ |  $n$ = 47, $p$ = 48, $d$ = 5
4He(28.37) | $1$ |  $n$ = 2, $p$ = 2, $d$ = 96
4He(28.39) | $2$ |  $n$ = 0.2, $p$ = 0.2, $d$ = 99.6
4He(28.64) | $0$ |  $d$ = 100
4He(28.67) | $2$ |  $d$ = 100
4He(29.89) | $2$ |  $n$ = 0.4, $p$ = 0.4, $d$ = 99.2
${}_{\Lambda}^{3}$H | $1/2$ | Stable
${}_{\Lambda}^{4}$He | $0$ | Stable
Table 1: Stable light (hyper)nuclei and low-lying resonances of the 4He system
(from BNL properties of nuclides [42]). $J$ denotes the total angular
momentum. The last column represents branching ratios of the decay channels,
in per cents. The $p,n,d$ correspond to the emission of protons, neutrons, or
deuterons, respectively. The hypernuclei are replicated from Refs. [13, 43].

The light (hyper)nuclei are included on equal basis with other hadrons in the
updated THESEUS [22]. They are sampled similarly to other hadrons. However,
there is an important difference between treatments of light (hyper)nuclei and
other hadrons. While the hadrons pass through the UrQMD afterburner stage
after the particlization, the light (hyper)nuclei do not, because the UrQMD is
not able to treat them. To partially overcome this shortcoming we imitate the
afterburner for light (hyper)nuclei by late freeze-out in the 3FD, following
the recipe of Ref. [22].

The calculation of hypernuclei production is completely similar to that of
light nuclei in Ref. [24]. In Ref. [24] it was found that the late freeze-out,
characterized by the freeze-out energy density $\varepsilon_{\rm frz}=$ 0.2
GeV/fm3, is preferable for deuterons, tritons, and 3He. We use precisely the
same freeze-out for the calculation of the ${}_{\Lambda}^{3}$H production. It
was also found [24] that the 4He observables are better reproduced with the
standard 3FD freeze-out, $\varepsilon_{\rm frz}=$ 0.4 GeV/fm3, which indicates
that the 4He nuclei better survive in the afterburner stage as more spatially
compact and tightly bound objects. We use this standard 3FD freeze-out for
simulations of the ${}_{\Lambda}^{4}$He production. However, the binding
energy of ${}_{\Lambda}^{4}$He ($B_{\Lambda}\simeq$ 2.4 MeV [4]) is similar to
that of 3He ($B_{N}=$ 2.6 MeV), which may imply that the late freeze-out is
more suitable for the ${}_{\Lambda}^{4}$He production. Therefore, the
calculations with the late freeze-out are also presented for
${}_{\Lambda}^{4}$He. Details of the freeze-out procedure in the 3FD are
described in Refs. [44, 45]. Similarly to Ref. [24] three different equations
of state (EoS’s) are used in the simulations: a purely hadronic EoS [46]
(hadr. EoS) and two EoS’s with deconfinement [47], i.e. an EoS with a first-
order phase transition (1PT EoS) and one with a smooth crossover transition
(crossover EoS).

##  III Bulk Observables

Rapidity distributions of ratios $t/p$, 4He$/p$,
${}_{\Lambda}^{3}$H/$\Lambda$, and ${}_{\Lambda}^{4}$He/$\Lambda$ in central
($b=$ 3 fm) Au+Au collisions at collision energy of $\sqrt{s_{NN}}=$ 3 GeV are
presented in Fig. 1.

Figure 1: Rapidity distributions of ratios $t/p$, 4He$/p$,
${}_{\Lambda}^{3}$H/$\Lambda$, and ${}_{\Lambda}^{4}$He/$\Lambda$ in central
($b=$ 3 fm) Au+Au collisions at collision energy of $\sqrt{s_{NN}}=$ 3 GeV.
Results are calculated with hadronic, 1PT and crossover EoS’s. The $t$ and
${}_{\Lambda}^{3}$H yields are calculated with the late freeze-out, while 4He
and ${}_{\Lambda}^{4}$He ones, with the conventional 3FD freeze-out (bold
lines). Results for ${}_{\Lambda}^{4}$He with the late freeze-out are also
displayed (thin lines marked as “late frz”). Protons and $\Lambda$s are
calculated with the conventional 3FD freeze-out and the subsequent UrQMD
afterburner. STAR data for protons and light nuclei (centrality 0–10%) are
from Ref. [48]. The ${}_{\Lambda}^{3}$H/$\Lambda$ midrapidity point is taken
from Refs. [30]. Full symbols display measured experimental points, whereas
the open ones are those reflected with respect to the midrapidity.

The proton and $\Lambda$ distributions are calculated within full THESEUS,
i.e. with the standard 3FD freeze-out and the UrQMD afterburner. The $t$ and
${}_{\Lambda}^{3}$H yields are calculated with the late freeze-out
($\varepsilon_{\rm frz}=$ 0.2 GeV/fm3), while 4He and ${}_{\Lambda}^{4}$He
ones, with the conventional 3FD freeze-out ($\varepsilon_{\rm frz}=$ 0.4
GeV/fm3) without the subsequent afterburner. Since the binding energy of
${}_{\Lambda}^{4}$He is similar to that of 3He, we also present calculations
for ${}_{\Lambda}^{4}$He production with the late freeze-out. Statistical
errors of the ${}_{\Lambda}^{4}$He calculations are displayed by the
respective bands.

The strange particles are rare probes at this collision energy. The canonical
ensemble with exact strangeness conservation is needed for their description.
The calculated yields of $\Lambda$s, ${}_{\Lambda}^{3}$H and
${}_{\Lambda}^{4}$He are considerably overestimated because of grand canonical
ensemble used in 3FD. Therefore, we do not demonstrate the $\Lambda$,
${}_{\Lambda}^{3}$H and ${}_{\Lambda}^{4}$He yields by themselves. This grand-
canonical overestimation is canceled in the ${}_{\Lambda}^{3}$H/$\Lambda$ and
${}_{\Lambda}^{4}$He/$\Lambda$ ratios.

As seen from Fig. 1, the calculation reasonably well reproduces the non-
strange ratios. The slight overestimation of the $t/p$ ratio is because of
underestimation of the proton yield, see Ref. [24]. The 4He is underestimated
in our calculations. Apparently, inclusion of light-nuclei resonances with
$A=5$ into the scheme, i.e. those of 5H, 5He, and 5Li [49] that decay into
4He, may correct this underestimation. The calculated result for the
midrapidity ${}_{\Lambda}^{3}$H/$\Lambda$ ratio falls within the error bars of
the experimental point [30]. It is remarkable that large difference between
the $t/p$ and ${}_{\Lambda}^{3}$H/$\Lambda$ ratios is reproduced without any
additional parameters.

The collision process develops in the hadronic phase in all considered
scenarios, i.e. hadronic, 1PT and crossover ones. The corresponding EoS’s are
very similar in the hadronic phase but not identical. Therefore, EoS-induced
differences in the ratios indicate uncertainties of the model predictions.

Figure 2: Midrapidity mean transverse momentum of protons, $\Lambda$ hyperons
and light (hyper)nuclei in central ($b=$ 3 fm) Au+Au collisions at collision
energy of $\sqrt{s_{NN}}=$ 3 GeV. Results are calculated with hadronic, 1PT
and crossover EoS’s. Protons are calculated within the conventional 3FD
freeze-out with the subsequent UrQMD afterburner. Deuterons, tritons, 3He, and
${}_{\Lambda}^{3}$H/$\Lambda$ are calculated with the late freeze-out, while
4He and ${}_{\Lambda}^{4}$He, the conventional 3FD freeze-out. STAR data are
from Refs. [29, 48].

Midrapidity mean transverse momentum of protons, $\Lambda$s and light
(hyper)nuclei in central ($b=$ 3 fm) collisions is displayed in Fig. 2. This
quantity characterizes the radial flow. Curves in Fig. 2 are displayed only
for eye guidance. As seen, these curves (for three different EoS’s) in fact
coincide. Results of the calculations are shown by dots on these curves.
Moreover, results for strange and non-strange species lie on the same curves.
The calculated points well agree with the data [29, 48]. Even slight deviation
of these curves from straight lines is reproduced.

##  IV Directed flow

The directed flow is a more delicate observable. The calculated directed flow
of protons, $\Lambda$ hyperons and light (hyper)nuclei (tritons, 4He,
${}_{\Lambda}^{3}$H and ${}_{\Lambda}^{4}$He) as function of rapidity in
semicentral ($b=$ 6 fm) Au+Au collisions at collision energy of
$\sqrt{s_{NN}}=$ 3 GeV is presented in Fig. 3.

Figure 3: Directed flow of protons, $\Lambda$ hyperons and light (hyper)nuclei
(tritons, 4He, ${}_{\Lambda}^{3}$H, and ${}_{\Lambda}^{4}$He) as function of
rapidity in semicentral ($b=$ 6 fm) Au+Au collisions at collision energy of
$\sqrt{s_{NN}}=$ 3 GeV. Results are calculated with hadronic, 1PT and
crossover EoS’s. The THESEUS simulations for light (hyper)nuclei were
performed with the late freeze-out ($\varepsilon_{\rm frz}=$ 0.2 GeV/fm3). The
flow of protons and $\Lambda$s is calculated with the conventional 3FD freeze-
out followed by the subsequent UrQMD afterburner. STAR data are from Refs.
[31, 50, 51]. Full symbols display measured experimental points, whereas the
open ones are those reflected with respect to the midrapidity.

The results are compared with STAR data [31, 50, 51]. We do not display
results for light nuclei (deuterons and 3He) because they are not directly
related to the considered hypernuclei. Results for all light nuclei can be
found in Ref. [24]. The THESEUS simulations for light (hyper)nuclei are
performed for the late freeze-out ($\varepsilon_{\rm frz}=$ 0.2 GeV/fm3) for
three EoS’s. The flow of protons and $\Lambda$ hyperons is calculated within
the full THESEUS, i.e. with the conventional 3FD freeze-out and the subsequent
UrQMD afterburner.

Figure 4: The same as in Fig. 3 but only for light hypernuclei
(${}_{\Lambda}^{3}$H and ${}_{\Lambda}^{4}$He) calculated with the
conventional 3FD freeze-out ($\varepsilon_{\rm frz}=$ 0.4 GeV/fm3).

The directed proton flow is almost independent of the used EoS [24]. The
calculated results perfectly (except for very forward and backward rapidities)
reproduce the experimental proton flow [51]. Agreement with the data [50]
becomes worse with increase of atomic number of light nucleus. If the
calculated midrapidity slope of the triton directed flow is only slightly
steeper than the experimental one, for 4He it is already noticeably steeper.

The directed $\Lambda$ flow depends on the EoS, while that of hypernuclei is
again EoS independent (up to the statistical fluctuations). Apparently the
nucleon content of the hypernuclei dominates in the $v_{1}$ formation. The
crossover scenario results in the best in reproduction of the midrapidity
slope of the $\Lambda$ flow. The ${}_{\Lambda}^{4}$He flow is reproduced to
the same extent as that of light nuclei. It is difficult to judge the degree
of agreement with the ${}_{\Lambda}^{3}$H flow data because of their large
error bars.

In Fig. 3, the directed flow of 4He and ${}_{\Lambda}^{4}$He is calculated
with the late freeze-out ($\varepsilon_{\rm frz}=$ 0.2 GeV/fm3) instead of the
conventional 3FD freeze-out that is preferable for 4He and presumably for
${}_{\Lambda}^{4}$He. The reason is that the 4He directed flow is independent
of the type (late or conventional) of the used freeze-out, as demonstrated in
Ref. [24]. Nevertheless, we additionally checked this independence for the
${}_{\Lambda}^{4}$He. The results of calculation of $v_{1}$ of
${}_{\Lambda}^{3}$H and ${}_{\Lambda}^{4}$He with the conventional 3FD freeze-
out is presented in Fig. 4. As seen, the $v_{1}$ flows with the conventional
freeze-out for both hypernuclei are just identical (up to statistical
fluctuations) to the late-freeze-out ones. The midrapidity slope of the proton
flow also remains unchanged after the afterburner [24]. All this indicates
that the baryon directed flow is formed at the early stage of the reaction.

##  V Summary

Simulations of the $\Lambda$-hyperon and light-hypernuclei production in Au+Au
collisions at $\sqrt{s_{NN}}=$ 3 GeV were performed within the updated THESEUS
event generator [22]. In the updated THESEUS, the light (hyper)nuclei are
treated thermodynamically, i.e. they are considered on the equal basis with
hadrons. The only additional parameter associated with the light (hyper)nuclei
is related to the late freeze-out that imitates the afterburner stage because
the UrQMD is not able to dynamically treat the light (hyper)nuclei. This is a
less demanding way to describe the light-(hyper)nuclei production as compared
to the coalescence.

The calculation of hypernuclei production is completely similar to that of
light nuclei in Ref. [24]. In Ref. [24] it was found that the late freeze-out
is preferable for deuterons, tritons, and 3He. We used precisely the same
freeze-out for the calculation of the ${}_{\Lambda}^{3}$H production. It was
also found [24] that the 4He observables are better reproduced with the
standard 3FD freeze-out, which indicates that the 4He nuclei better survive in
the afterburner stage as more spatially compact and tightly bound objects. We
used this standard 3FD freeze-out for simulations of the ${}_{\Lambda}^{4}$He
production. However, the binding energy of ${}_{\Lambda}^{4}$He
($B_{\Lambda}\simeq$ 2.4 MeV [4]) is similar to that of 3He ($B_{N}=$ 2.6
MeV), which may imply that the late freeze-out is more suitable for the
${}_{\Lambda}^{4}$He production. The hypernuclei results were compared with
recent STAR data [29, 30, 31], as well as with the results of light-nuclei
calculations [24].

It is found that the calculated midrapidity ${}_{\Lambda}^{3}$H/$\Lambda$
ratio falls within the error bars of the experimental point [30]. It is
remarkable that large difference between the $t/p$ and
${}_{\Lambda}^{3}$H/$\Lambda$ ratios is reproduced without any additional
parameters. Rapidity distributions of ${}_{\Lambda}^{3}$H/$\Lambda$ and
${}_{\Lambda}^{4}$He/$\Lambda$ ratios are predicted. Midrapidity mean
transverse momenta of protons, $\Lambda$s and light (hyper)nuclei in central
collisions well agree with the data [29, 30]. The calculated directed flow
also reasonably well reproduces of the data [31]. The directed flow turned out
to be independent of the type (late or conventional) of the freeze-out. This
indicates that the baryon directed flow is formed at the early stage of the
reaction.

###### Acknowledgements.

We are sincerely grateful to Iurii Karpenko and David Blaschke who made
enormous contributions on the early stage of this project. Fruitful
discussions with D.N. Voskresensky are gratefully acknowledged. This work was
carried out using computing resources of the federal collective usage center
“Complex for simulation and data processing for mega-science facilities” at
NRC ”Kurchatov Institute” [52] and computing resources of the supercomputer
”Govorun” at JINR [53].

## References

  * [1] L. Adamczyk et al. [STAR], “Measurement of the ${}^{3}_{\Lambda}$H lifetime in Au+Au collisions at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 97, no.5, 054909 (2018) doi:10.1103/PhysRevC.97.054909 [arXiv:1710.00436 [nucl-ex]].
  * [2] S. Acharya et al. [ALICE], “${}^{3}_{\Lambda}\mathrm{H}$ and ${}^{3}_{\bar{\Lambda}}\mathrm{\overline{H}}$ lifetime measurement in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 5.02 TeV via two-body decay,” Phys. Lett. B 797, 134905 (2019) [arXiv:1907.06906 [nucl-ex]].
  * [3] M. Abdallah et al. [STAR], “Measurements of $H_{\Lambda}^{3}$ and $H_{\Lambda}^{4}$ Lifetimes and Yields in Au+Au Collisions in the High Baryon Density Region,” Phys. Rev. Lett. 128, no.20, 202301 (2022) doi:10.1103/PhysRevLett.128.202301 [arXiv:2110.09513 [nucl-ex]].
  * [4] M. Abdallah et al. [STAR], “Measurement of H$\Lambda$4 and He$\Lambda$4 binding energy in Au+Au collisions at sNN = 3 GeV,” Phys. Lett. B 834, 137449 (2022) doi:10.1016/j.physletb.2022.137449 [arXiv:2207.00778 [nucl-ex]].
  * [5] A. Gal, E. V. Hungerford and D. J. Millener, “Strangeness in nuclear physics,” Rev. Mod. Phys. 88, no.3, 035004 (2016) [arXiv:1605.00557 [nucl-th]].
  * [6] M. Knöll and R. Roth, “Hyperon-nucleon interaction constrained by light hypernuclei,” Phys. Lett. B 846, 138258 (2023) [arXiv:2307.11577 [nucl-th]].
  * [7] H. Le, J. Haidenbauer, U. G. Meißner and A. Nogga, “Separation energies of light $\Lambda$ hypernuclei and their theoretical uncertainties,” [arXiv:2308.01756 [nucl-th]].
  * [8] D. Lonardoni, A. Lovato, S. Gandolfi and F. Pederiva, “Hyperon Puzzle: Hints from Quantum Monte Carlo Calculations,” Phys. Rev. Lett. 114, no.9, 092301 (2015) [arXiv:1407.4448 [nucl-th]].
  * [9] K. A. Maslov, E. E. Kolomeitsev and D. N. Voskresensky, “Solution of the Hyperon Puzzle within a Relativistic Mean-Field Model,” Phys. Lett. B 748, 369-375 (2015) [arXiv:1504.02915 [astro-ph.HE]]; “Making a soft relativistic mean-field equation of state stiffer at high density,” Phys. Rev. C 92, no.5, 052801 (2015) [arXiv:1508.03771 [astro-ph.HE]]; “Relativistic Mean-Field Models with Scaled Hadron Masses and Couplings: Hyperons and Maximum Neutron Star Mass,” Nucl. Phys. A 950, 64-109 (2016) [arXiv:1509.02538 [astro-ph.HE]].
  * [10] M. Fortin, S. S. Avancini, C. Providência and I. Vidaña, Phys. Rev. C 95, no.6, 065803 (2017) doi:10.1103/PhysRevC.95.065803 [arXiv:1701.06373 [nucl-th]].
  * [11] A. S. Khvorostukhin, V. D. Toneev and D. N. Voskresensky, “Equation of State for Hot and Dense Matter: sigma- omega- rho Model with Scaled Hadron Masses and Couplings,” Nucl. Phys. A 791, 180-221 (2007) [arXiv:nucl-th/0612058 [nucl-th]]; “Relativistic Mean-Field Model with Scaled Hadron Masses and Couplings,” Nucl. Phys. A 813, 313-346 (2008) [arXiv:0802.3999 [nucl-th]].
  * [12] S. Zhang, J. H. Chen, H. Crawford, D. Keane, Y. G. Ma and Z. B. Xu, “Searching for onset of deconfinement via hypernuclei and baryon-strangeness correlations,” Phys. Lett. B 684, 224-227 (2010) [arXiv:0908.3357 [nucl-ex]].
  * [13] J. Steinheimer, K. Gudima, A. Botvina, I. Mishustin, M. Bleicher and H. Stocker, “Hypernuclei, dibaryon and antinuclei production in high energy heavy ion collisions: Thermal production versus Coalescence,” Phys. Lett. B 714, 85-91 (2012) [arXiv:1203.2547 [nucl-th]].
  * [14] T. Shao, J. Chen, C. M. Ko, K. J. Sun and Z. Xu, “Yield ratio of hypertriton to light nuclei in heavy-ion collisions from $\rm\sqrt{s_{NN}}$ = 4.9 GeV to 2.76 TeV,” Chin. Phys. C 44, no.11, 114001 (2020) [arXiv:2004.02385 [nucl-ex]].
  * [15] J. Aichelin, E. Bratkovskaya, A. Le Fèvre, V. Kireyeu, V. Kolesnikov, Y. Leifels, V. Voronyuk and G. Coci, Parton-hadron-quantum-molecular dynamics: A novel microscopic $n$-body transport approach for heavy-ion collisions, dynamical cluster formation, and hypernuclei production, Phys. Rev. C 101, no.4, 044905 (2020) [arXiv:1907.03860 [nucl-th]].
  * [16] S. Gläßel, V. Kireyeu, V. Voronyuk, J. Aichelin, C. Blume, E. Bratkovskaya, G. Coci, V. Kolesnikov and M. Winn, “Cluster and hypercluster production in relativistic heavy-ion collisions within the parton-hadron-quantum-molecular-dynamics approach,” Phys. Rev. C 105, no.1, 014908 (2022) [arXiv:2106.14839 [nucl-th]].
  * [17] Z. Q. Feng, “Dynamics of light hypernuclei in collisions of 197Au+197Au at GeV energies,” Eur. Phys. J. A 57, no.1, 18 (2021) [arXiv:2109.01270 [nucl-th]].
  * [18] T. Reichert, J. Steinheimer, V. Vovchenko, B. Dönigus and M. Bleicher, “Energy dependence of light hypernuclei production in heavy-ion collisions from a coalescence and statistical-thermal model perspective,” Phys. Rev. C 107, no.1, 014912 (2023) [arXiv:2210.11876 [nucl-th]].
  * [19] N. Buyukcizmeci, T. Reichert, A. S. Botvina and M. Bleicher, “Nucleosynthesis of light nuclei and hypernuclei in central Au+Au collisions at sNN=3 GeV,” Phys. Rev. C 108, no.5, 054904 (2023) [arXiv:2306.17145 [nucl-th]].
  * [20] A. Andronic, P. Braun-Munzinger, J. Stachel and H. Stocker, Production of light nuclei, hypernuclei and their antiparticles in relativistic nuclear collisions, Phys. Lett. B 697, 203-207 (2011) [arXiv:1010.2995 [nucl-th]].
  * [21] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Decoding the phase structure of QCD via particle production at high energy, Nature 561, no.7723, 321-330 (2018) [arXiv:1710.09425 [nucl-th]].
  * [22] M. Kozhevnikova, Y. B. Ivanov, I. Karpenko, D. Blaschke and O. Rogachevsky, Update of the Three-fluid Hydrodynamics-based Event Simulator: light-nuclei production in heavy-ion collisions, Phys. Rev. C 103, no.4, 044905 (2021) [arXiv:2012.11438 [nucl-th]].
  * [23] M. Kozhevnikova and Y. B. Ivanov, Light-nuclei production in heavy-ion collisions within a thermodynamical approach, Phys. Rev. C 107, no.2, 024903 (2023) [arXiv:2210.07334 [nucl-th]]; Particles 6, no.1, 440-450 (2023).
  * [24] M. Kozhevnikova and Y. B. Ivanov, [arXiv:2311.08092 [nucl-th]].
  * [25] D. Oliinychenko, L. G. Pang, H. Elfner and V. Koch, Microscopic study of deuteron production in PbPb collisions at $\sqrt{s}=2.76TeV$ via hydrodynamics and a hadronic afterburner, Phys. Rev. C 99, no.4, 044907 (2019) [arXiv:1809.03071 [hep-ph]].
  * [26] J. Staudenmaier, D. Oliinychenko, J. M. Torres-Rincon and H. Elfner, Deuteron production in relativistic heavy ion collisions via stochastic multiparticle reactions, Phys. Rev. C 104, no.3, 034908 (2021) [arXiv:2106.14287 [hep-ph]].
  * [27] K. J. Sun, R. Wang, C. M. Ko, Y. G. Ma and C. Shen, [arXiv:2106.12742 [nucl-th]].
  * [28] K. J. Sun, R. Wang, C. M. Ko, Y. G. Ma and C. Shen, Unveiling the dynamics of nucleosynthesis in relativistic heavy-ion collisions, [arXiv:2207.12532 [nucl-th]].
  * [29] Y. Ji [STAR], “Measurements on the production and properties of light hypernuclei at STAR,” EPJ Web Conf. 276, 04003 (2023)
  * [30] Yuanjing Ji, talk at Quark Matter 2023, https://indico.cern.ch/event/1139644/contributions/5456392/attachments/2707583/4708403/talk_FXT_H3L_Sep08_v11.pdf
  * [31] B. Aboona et al. [STAR], “Observation of Directed Flow of Hypernuclei H$\Lambda 3andH\Lambda 4$ in $\sqrt{s_{NN}}=$ 3 GeV Au+Au Collisions at RHIC,” Phys. Rev. Lett. 130, no.21, 212301 (2023) [arXiv:2211.16981 [nucl-ex]].
  * [32] Y. Nara, N. Otuka, A. Ohnishi, K. Niita and S. Chiba, Study of relativistic nuclear collisions at AGS energies from p + Be to Au + Au with hadronic cascade model, Phys. Rev. C 61, 024901 (2000) [arXiv:nucl-th/9904059 [nucl-th]].
  * [33] M. Isse, A. Ohnishi, N. Otuka, P. K. Sahu and Y. Nara, Mean-field effects on collective flows in high-energy heavy-ion collisions from AGS to SPS energies, Phys. Rev. C 72, 064908 (2005) [arXiv:nucl-th/0502058 [nucl-th]].
  * [34] J. P. Bondorf, A. S. Botvina, A. S. Ilinov, I. N. Mishustin and K. Sneppen, “Statistical multifragmentation of nuclei,” Phys. Rept. 257, 133-221 (1995)
  * [35] P. Batyuk et al., Event simulation based on three-fluid hydrodynamics for collisions at energies available at the Dubna Nuclotron-based Ion Collider Facility and at the Facility for Antiproton and Ion Research in Darmstadt, Phys. Rev. C 94, 044917 (2016) [arXiv:1608.00965 [nucl-th]].
  * [36] P. Batyuk, D. Blaschke, M. Bleicher, Y. B. Ivanov, I. Karpenko, L. Malinina, S. Merts, M. Nahrgang, H. Petersen and O. Rogachevsky, Three-fluid Hydrodynamics-based Event Simulator Extended by UrQMD final State interactions (THESEUS) for FAIR-NICA-SPSBES/RHIC energies, EPJ Web Conf. 182, 02056 (2018) [arXiv:1711.07959 [nucl-th]].
  * [37] Y. B. Ivanov, Alternative Scenarios of Relativistic Heavy-Ion Collisions: I. Baryon Stopping, Phys. Rev. C 87, no.6, 064904 (2013) [arXiv:1302.5766 [nucl-th]].
  * [38] Y. B. Ivanov, V. N. Russkikh and V. D. Toneev, Relativistic heavy-ion collisions within 3-fluid hydrodynamics: Hadronic scenario, Phys. Rev. C 73, 044904 (2006) [nucl-th/0503088].
  * [39] S. A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland, M. Hofmann, S. Hofmann and J. Konopka, et al., Microscopic models for ultrarelativistic heavy ion collisions, Prog. Part. Nucl. Phys. 41, 255-369 (1998) [arXiv:nucl-th/9803035 [nucl-th]].
  * [40] Y. B. Ivanov and A. A. Soldatov, Light fragment production at CERN Super Proton Synchrotron, Eur. Phys. J. A 53, no. 11, 218 (2017) [arXiv:1703.05040 [nucl-th]].
  * [41] E. Shuryak and J. M. Torres-Rincon, Baryon preclustering at the freeze-out of heavy-ion collisions and light-nuclei production, Phys. Rev. C 101, no.3, 034914 (2020) [arXiv:1910.08119 [nucl-th]].
  * [42] https://www.nndc.bnl.gov/nudat3/getdataset.jsp?nucleus=4HE&unc=nds
  * [43] D. H. Davis, “50 years of hypernuclear physics. I. The early experiments,” Nucl. Phys. A 754, 3-13 (2005)
  * [44] V. N. Russkikh and Yu. B. Ivanov, Dynamical freeze-out in 3-fluid hydrodynamics, Phys. Rev. C 76, 054907 (2007) [nucl-th/0611094].
  * [45] Yu. B. Ivanov and V. N. Russkikh, On freeze-out problem in relativistic hydrodynamics, Phys. Atom. Nucl. 72, 1238 (2009) [arXiv:0810.2262 [nucl-th]].
  * [46] I. N. Mishustin, V. N. Russkikh and L. M. Satarov, Fluid dynamical model of relativistic heavy ion collision, Sov. J. Nucl. Phys. 54, 260-314 (1991).
  * [47] A. S. Khvorostukin, V. V. Skokov, V. D. Toneev and K. Redlich, Lattice QCD constraints on the nuclear equation of state, Eur. Phys. J. C 48, 531 (2006) [nucl-th/0605069].
  * [48] [STAR], “Production of Protons and Light Nuclei in Au+Au Collisions at $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV with the STAR Detector,” [arXiv:2311.11020 [nucl-ex]].
  * [49] V. Vovchenko, B. Dönigus, B. Kardan, M. Lorenz and H. Stoecker, Feeddown contributions from unstable nuclei in relativistic heavy-ion collisions, Phys. Lett. B, 135746 (2020) [arXiv:2004.04411 [nucl-th]].
  * [50] M. S. Abdallah et al. [STAR], Light nuclei collectivity from $\sqrt{s_{NN}}$ = 3 GeV Au+Au collisions at RHIC, Phys. Lett. B 827, 136941 (2022) [arXiv:2112.04066 [nucl-ex]].
  * [51] M. S. Abdallah et al. [STAR], Disappearance of partonic collectivity in sNN=3GeV Au+Au collisions at RHIC, Phys. Lett. B 827, 137003 (2022) [arXiv:2108.00908 [nucl-ex]].
  * [52] http://ckp.nrcki.ru/
  * [53] http://hlit.jinr.ru/supercomputer_govorun/