# Optimizing optical potentials with physics-inspired learning algorithms

M. Calzavara1,4,∗, Y. Kuriatnikov2,∗, A. Deutschmann-Olek3, F. Motzoi1, S.
Erne2, A. Kugi3, T. Calarco1,4, J. Schmiedmayer2, M. Prüfer2,§
1Forschungszentrum Jülich GmbH, Peter Grünberg Institute, Quantum Control
(PGI-8), 52425 Jülich, Germany
2Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien,
Stadionallee 2, 1020 Vienna, Austria
3Automation and Control Institute, TU Wien, Gußhausstraße 27-29, 1040 Vienna,
Austria
4Institute for Theoretical Physics, Universität zu Köln, 50937 Cologne,
Germany

###### Abstract

We present our new experimental and theoretical framework which combines a
broadband superluminescent diode (SLED/SLD) with fast learning algorithms to
provide speed and accuracy improvements for the optimization of 1D optical
dipole potentials, here generated with a Digital Micromirror Device (DMD). To
characterize the setup and potential speckle patterns arising from coherence,
we compare the superluminescent diode to a single-mode laser by investigating
interference properties. We employ Machine Learning (ML) tools to train a
physics-inspired model acting as a digital twin of the optical system
predicting the behavior of the optical apparatus including all its
imperfections. Implementing an iterative algorithm based on Iterative Learning
Control (ILC) we optimize optical potentials an order of magnitude faster than
heuristic optimization methods. We compare iterative model-based “offline”
optimization and experimental feedback-based “online” optimization. Our
methods provide a new route to fast optimization of optical potentials which
is relevant for the dynamical manipulation of ultracold gases.

**footnotetext: These authors contributed equally to this work §
§footnotetext<EMAIL_ADDRESS>
##  I Introduction

The precise control and manipulation of light fields are required for many
diverse areas of research ranging from microscopy [1] to quantum simulators
[2]. In particular, optical beam shaping constitutes a common task, for which
wavefront manipulating devices, such as Spatial Light Modulators (SLM) or
Digital Micromirror Devices (DMD), are especially suited. The beam shaping is
important for experiments with ultracold gases, where optical dipole
potentials have proven to be a versatile tool to provide the demanded level of
control. In combination with a DMD, almost arbitrary shaping of the optical
potential is possible, both in 1D [3, 4, 5, 6] and 2D [7, 8, 9, 10, 11]
settings. These potentials can, as an example, be used for generating
homogeneous box potentials in ultracold gases experiments [9]. In addition to
static potentials, dynamical perturbations and time-averaged potentials can
also be created, by projecting sequences of patterns onto the DMD [4, 12].

There exist two main approaches to the optimization of optical potentials:
precalculating DMD patterns based on physical assumptions [13, 14] and models
[7] or iteratively updating DMD patterns based on experimental feedback [7,
6]. The first approach avoids the need for feedback measurement but is limited
by model precision and thus requires detailed system characterization, while
the second gives the most accurate results but typically requires a large
number of experimental iterations. Here we implement a purely data-driven
approach that combines different learning algorithms to get the best of both
worlds.

Using a digital twin of the system makes it possible to shape different types
of target potentials without the need for experimental feedback. Yet because
of speckles caused by imperfections, a model featuring just a few known
experimental parameters (parametric model) can only predict its behavior up to
limited precision and might not always be suitable for precise potential
shaping. In some cases, though, such as in “clean” systems with pin-hole
filtering, simulations combined with input beam characterization deliver very
low errors on precalculated DMD patterns [7]. We improve the prediction
performance compared to parametric models by employing data-driven learning
techniques. Learning methods are already used for estimating the transfer
matrix of complex optical systems [15] and provide good results. As their main
disadvantage, they generally require a large amount of data for sufficient
training. In our approach, we develop a physics-inspired model that requires a
smaller amount of training data and thus saves experimental time.

Despite any improvement in model precision, the effect of residual error
sources can only be mitigated by using experimental feedback [7, 6].
Therefore, we introduce a feedback optimization method based on Iterative
Learning Control (ILC) [16, 17]. This method is directly applicable to various
types of experiments with wavefront manipulating devices. Since system
knowledge is directly employed in the update law that adjusts the DMD setting
based on feedback, the number of required experimental iterations is
significantly reduced compared to heuristic methods [6].

Because of the learning-based nature of both the ML model and ILC method, they
benefit from high predictability of the system behavior, which can be achieved
by using light sources with low coherence. White light was used for Bose-
Einstein condensates (BEC) trapping in order to minimize the impact of
speckles on density fluctuations [18], and the advantages of using a
superluminescent diode (SLED) in comparison with a single-frequency laser were
shown also in combination with a DMD [12]. Indeed, we find that using the SLED
improves model prediction results, while the feedback-based methods give
errors comparable to measurement errors for both light sources.

In perspective, the ability to efficiently shape a large number of potentials
using a DMD is a precious tool for the dynamical manipulation of quantum
gases. The realization of non-harmonic optimal protocols [19] and quantum
thermal machines [20] are but two examples among countless possible
applications.

The paper is organized as follows: in Sec.II we describe our experimental
setup and test the coherence properties of the light sources by measuring
interferences of a SLED in comparison with a single-frequency laser; in
Sec.III we introduce the physics-inspired learning model used for representing
the optical system including its imperfections; in Sec.IV we use iterative
learning control algorithms for optimizing optical dipole potentials; in Sec.V
we summarize our results.

##  II Experimental setup

Figure 1: (a) Schematic of the experimental setup. The light source is an
interchangeable SLED/laser connected with a fiber to the collimator and
projected onto the DMD. The DMD is an array of individually controllable
micromirrors which allows for projecting arbitrary patterns. The system is
shown as simplified schematic omitting mirrors. The first two lenses form a
4-f telescope with the DMD in its image plane, and the slit placed in the
Fourier plane. The red path corresponds to the beam and the blue to a point
source in the DMD plane. (b)-(e) Images taken with the CCD camera in the image
plane of the DMD with the same pattern projected for (b) SLED without slit (c)
laser without slit (d) SLED with slit (e) laser with slit. Figure 2: (a)
Schematics of test measurement of coherence properties. Two pairs of 3x3 DMD
superpixels a/b and 1/2 are translated in a direction orthogonal and parallel
to the slit, respectively. The resolution of the optical system is
$\sigma_{\parallel}=3.3\,\mu$m in the longitudinal (parallel to slit) and
$\sigma_{\perp}=25\,\mu$m in the transverse (orthogonal to slit) direction,
with respect to the slit orientation. (b)-(e) Intensity measured between two
superpixels as a function of distance $d$ between them; the blue solid line
shows the intensity measured for both superpixels simultaneously “on”
(electric fields from superpixels contribute in measured intensity as
$|E_{a/1}+E_{b/2}|^{2}$), the red solid line is the sum of intensities
measured with only one of two superpixels “on” (contribution is independent
$|E_{a/1}|^{2}+|E_{b/2}|^{2}$). The distance $d$ is given in units of
$\sigma_{\parallel}$ and $\sigma_{\perp}$. When the laser-illuminated
superpixels are moved orthogonal to the slit (b), positive interference is
visible in the $0^{th}$ slit diffraction maximum and negative interference in
the $1^{st}$. When SLED-illuminated superpixels are moved orthogonal to the
slit (c) a slight interference effect is visible only in the $0^{th}$ slit
diffraction maximum for very small distances between the superpixels. Panel
(d) shows laser-illuminated superpixels moving parallel to the slit, and panel
(e) shows SLED-illuminated superpixels moving parallel to the slit. There are
no noticeable interference effects in the longitudinal direction on the scales
larger than the camera pixel size for both light sources.

In this section, we describe our experimental setup and compare the optical
coherence properties of a single-frequency laser and broadband SLED. The
optical apparatus was designed and optimized for manipulating 1D optical
dipole potentials in our atom chip experiment [21]. The simplified optical
setup is shown in Fig.1 (a). In the setup, the light source is easily
interchangeable: we here use either a superluminescent 110 mW fiber-coupled
diode with a spectral width of 12.7 nm and a central wavelength of 837 nm or a
single-frequency 780 nm laser. The fiber is connected to a collimator and
illuminates the DMD with collimated light. The DMD is placed in the focus of
the first lens. The first two lenses together form a 4-f telescope. The focal
point between the two lenses is a Fourier plane with respect to the DMD’s
image plane. The slit, which is adjustable in the transverse direction and is
parallel to the observable 1D optical potential, is placed in the Fourier
plane and closed to 0.625 mm. Being placed in the Fourier plane, the slit acts
as a low-pass filter for the transverse spatial frequencies of the light
field. This means it cuts off high-frequency $k$-modes of the DMD pattern
effectively leading to a lowering of the resolution. The other 4 lenses form
an objective designed to correct chromatic aberrations for the broadband SLED.
The system acts with a demagnification of 17.5 and resolution
$\sigma_{\parallel}=3.3\,\mu$m in the direction parallel to the slit (the
resolution is given by our atom chip experiment for which the optical
apparatus was designed and optimized) and $\sigma_{\perp}=25$ $\mu m$ in the
transverse direction (orthogonal to the slit). In the end, an image of the
optical potential is acquired with a CCD camera sensor with 2.4 $\mu m$ pixel
size.

We use a $10.8$ $\mu m$ pitch near-infrared DMD, which is a 2D array of
1280x800 micromirrors. In the image plane (on the CCD camera) 5 DMD pixels
correspond to $\sigma_{\parallel}$ in the longitudinal direction and 40 DMD
pixels correspond to $\sigma_{\perp}$ in the transverse direction. Since many
pixels in the transverse direction contribute to the local intensity (close to
the center pixels contribute almost equally to the optical potential while far
off-center pixels are contributing less), we can perform very smooth gray
scaling, which is very important for high-precision optimization of 1D optical
dipole potentials via DMD. Yet, while a narrower slit allows for better gray
scaling and with that lower discretization error and higher accuracy, it leads
to significantly reduced light intensity.

The main advantage of using the SLED compared to the single-frequency laser as
a light source is its reduced temporal coherence [22, 23] which is why random
interferences (speckles) are reduced compared to the laser (see Fig.1 (b,c)).
We anticipate the speckles influencing the predictability of the system.
Therefore, we characterize coherence effects for both light sources. More
specifically, we want to understand whether we can consider the individual
pixels as coherent or incoherent sources. The imaged intensity is in general
given by the squared sum of the electric fields of the individual pixels. In
the incoherent case, it reduces to the sum of intensities from the individual
pixels such that the system behaves linearly with respect to intensity.

We prepare two series of patterns, each of them containing two superpixels
moving away from each other (by superpixel, we understand a square pattern of
3x3 micromirrors). In the image plane on the CCD camera, the size of the
superpixel is below the resolution of the optical apparatus and the superpixel
can be considered as a point source, at the same time minimizing diffraction
from individual DMD mirrors. We measure the intensity in the center between
two superpixels. First, we projected both superpixels simultaneously such that
their optical fields can interfere. Second, we only project one of the
superpixels onto the DMD at a time to infer the intensity of a single
superpixel. Then the intensity measured with the first configuration was
compared to the sum of intensities measured with the second configuration.

In the first series of patterns (purple squares in Fig.2), the pair of
superpixels is moving orthogonal to the slit orientation (transverse
direction), where the slit is closed to 0.625 mm and resolution is lowered to
$\sigma_{\perp}$=25 $\mu m$. For both SLED and laser in Fig.2 (b,c), the slit
diffraction maxima are clearly distinguishable. For the laser, positive
interference can be observed in the $0^{th}$ maximum and negative in the
$1^{st}$. For the SLED, interference is observable on a very small scale
(below $\sigma_{\perp}$) and in general, the behavior is linear.

In the second series (gray squares in Fig.2) superpixels are moving in the
longitudinal direction, where the slit is fully open to 13 mm. For both SLED
and laser, in Fig.2 (d,e), the behavior of the system can be considered
linear. It is worth mentioning that the resolution $\sigma_{\parallel}$=3.3
$\mu m$ is very close to the camera pixel size 2.4 $\mu m$ and small
interference effects are hardly accessible on this scale.

##  III Physics-inspired learning model

Figure 3: (a) Training of the physics-inspired model. The training dataset is
composed of 1D virtual inputs $\boldsymbol{\nu}$ which are transformed into 2D
DMD patterns (according to Eq. (2)) and projected in the image plane. We read
out a pixel row from the CCD camera to infer the 1D optical potential
generated by the system. After data acquisition, the 4-layer model is trained
by minimizing the loss function. The polynomial layer represents the nonlinear
relation between virtual input and local grayscale, while the other three
layers take into account physical properties, i.e. the beam shape, the
effective point spread function, and the background. The trained model (red
solid segment) is a digital twin of the physical system which converts 1D
virtual input to predicted optical potential. All the signals and layers shown
in the figure are actual experimental data. (b) Test loss dependence on the
size of the training dataset. The loss converges to $2.4\%$ for a dataset size
of $\sim 10^{2}$, indicating that the model does not exhibit overfitting even
for small training datasets. (c) Test loss dependence on the spatial frequency
content of training and test datasets (see main text for details). (d) Two
cuts from figure (c) for $\sigma_{train}=0.3$ $\mu m$ and $\sigma_{train}=2.0$
$\mu m$. Figure 4: Scheme for optimizing optical potentials using Iterative
Learning Control. The iterative scheme is initialized by choosing an initial
guess and the target potential. We use two methods that differ in the feedback
source: offline and online. The offline method uses the trained ML model to
predict the optical potential. The online method projects the pattern on the
DMD and measures the optical potential. Then the predicted (measured)
potential is subtracted from the target to get the error that ILC uses to
update the 1D virtual input. The procedure is repeated iteratively until the
error converges or matches conditions for loop break. For the best
performance, the output of the offline optimization can be used as the initial
guess for online optimization. The result of the optimization loop is an
optimized DMD pattern. All the signals in the figure are actual experimental
data.

In this section, we describe how we use Machine Learning (ML) tools and
experimental data to obtain a digital twin of the system. It should represent
precisely the experimental system, while at the same time having a small
number of parameters to reduce the required number of experimental
measurements. In the following, we outline our approach based on a physics-
inspired model.

To formulate our problem, we treat the DMD as a 2D array to which we associate
a binary configuration matrix $u_{ij}$. For each pixel, a value of 1
corresponds to the mirror position that reflects light to the optical system,
while 0 corresponds to the mirror position that reflects light away from the
system. The potential $\boldsymbol{V}$ is a vector obtained by selecting one
row of camera pixels in the imaged output. We refer to the coordinate along
this row (and therefore, parallel to the slit) as $z$. We then cast model
training in the language of regression [24]: given a set of data points
composed by $K$ couples $(u^{(k)}_{ij},V^{(k)}_{i})$, and a family of
functions $\mathcal{M}_{\boldsymbol{\alpha}}$, parametrized by
$\boldsymbol{\alpha}=[\alpha_{1},...,\alpha_{N}]$, we find the values of
$\boldsymbol{\alpha}$ that minimize a loss function
$L(u^{(k)}_{ij},V^{(k)}_{i},\boldsymbol{\alpha})$ defined below.

Regression problems are known to be affected by the bias-variance tradeoff
[25]. This implies that large models, having
$\dim(\boldsymbol{\alpha})>>K\dim(\boldsymbol{V})$, tend to approximate well
the training data (low bias) but fail to accurately generalize the prediction
on test data (high variance). This phenomenon is known as overfitting, and it
puts a challenging limitation to our program. In order to alleviate this
problem, our approach consists in developing a physics-inspired model (see
Fig.3 (a)) based on the knowledge we have about optical systems. This way, we
can reduce the number of its coefficients to the minimum necessary to
represent the system precisely enough, while avoiding overfitting.

More specifically, we employ a dimensional reduction technique which we refer
to as the “virtual input” (cf. [17]). The virtual input measures the relative
optical intensity induced by the pixels in the corresponding column. Since the
narrow slit averages along the transverse direction, we can use the pixels on
each column to create a set of gray scales at each point along the
longitudinal direction $z$. In the limit of a very narrow slit, all the pixels
in a column contribute almost equally to the final intensity, so we can define
the virtual input $\boldsymbol{\nu}$ as the vector:

$\nu_{j}=\frac{1}{N_{row}}\sum_{i}u_{ij}\,.$ (1)

This quantity represents the fraction of pixels that are turned on in each
column, with respect to the total number of rows inside the area of interest
(AOI). The AOI has a size of $(N_{row},N_{col})$, which is determined by the
beam size and position on the DMD. In the case of non-zero slit widths, the
mapping between a virtual input according to Eq. eq:DMDtoVRT and the actual
relative optical intensity due to the column pattern is non-linear and depends
on the beam and on the transverse shape of the point spread function (PSF).
Therefore, we have to take this effect into account when we design our model.

Once $u_{ij}$ is transformed into $\nu_{j}$ some information is lost since
different binary matrices map onto the same virtual input. To invert this
mapping we have to choose a subset of binary matrices over which the mapping
is one-to-one. We define the inverse map by turning on the pixels one by one
on alternating rows above and below $i_{c}$, according to the order
$i_{c},i_{c}+1,i_{c}-1,i_{c}+2,i_{c}-2,...$. More compactly:

$u_{ij}=\theta\left(\frac{N_{row}\nu_{j}}{2}-\lvert
i-i_{c}-\frac{1}{4}\rvert\right)\,,$ (2)

where $\theta(x)$ is the Heaviside step function and $i_{c}$ the central row
of the AOI. In practice, the subset of binary matrices defined by this mapping
is wide enough so that it can be used for shaping arbitrary potentials.

The model $\boldsymbol{V}=\mathcal{M}_{\boldsymbol{\alpha}}(\boldsymbol{\nu})$
that we propose for the task is depicted in Fig.3 (a), and reads as follows:

$V_{i}=\left\lvert\sum_{j=-M}^{M}g_{j}P(\nu_{i-j},q_{1},...,q_{N_{P}})p_{i-j}\right\lvert+c_{i}$
(3)

where $P(x,q_{1},...,q_{N_{P}})=\sum_{n=1}^{N_{P}}q_{n}x^{n}$ is a polynomial
of degree $N_{P}=5$ with no constant term, and the vector of parameters
$\boldsymbol{\alpha}$ is the concatenation of
$\boldsymbol{g}$,$\boldsymbol{q}$,$\boldsymbol{p}$ and $\boldsymbol{c}$ with
$\dim(\boldsymbol{\alpha})\sim 1400$. The polynomial function $P$ represents
in an abstract way the non-linear relation between virtual input and resulting
local relative intensity as discussed above. The remaining parameters are
chosen to mimic features of the physical system: the convolutional kernel
$\boldsymbol{g}$ of size $N_{g}=2M+1\sim 71$ plays the role of the
longitudinal PSF, the position-dependent term $\boldsymbol{p}$, of size
$N_{V}\sim 650$, mimics the inhomogeneous light beam and the offset term
$\boldsymbol{c}$, also of size $N_{V}$, gives the background. In the
convolution sum we employ zero padding, which means that the summand is zero
if $i-j$ exceeds the index range $[1,N_{V}]$.

For the training, we choose the loss function as:

$L(\boldsymbol{\nu}^{(k)},\mathbf{V}^{(k)},\boldsymbol{\alpha})=\sum_{k=1}^{K}\frac{\sum_{i=1}^{N_{V}}|[\mathcal{M}_{\boldsymbol{\alpha}}(\boldsymbol{\nu}^{(k)})]_{i}-V^{(k)}_{i}|}{\sum_{i=1}^{N_{V}}V^{(k)}_{i}}\,,$
(4)

which is the mean absolute error between prediction and measurement,
normalized by the average potential itself. This way, potentials with
different average values will contribute equally during the training.

To test the performance of the ML model, we create a data set of 10,000 random
virtual inputs by sampling a white noise probability distribution. We filter
the sample with Gaussian filters with 10 values of $\sigma_{data}$ ranging
from $0.2$ to $100$ DMD pixels (which corresponds to the range $[0.1,65]\mu$m
mapped to the image plane). This way, the dataset was subdivided into 10
subsets with varying upper spatial cutoff frequency. For each virtual input,
the corresponding potential was measured and stored (along with the
corresponding input).

We tested the dependence of the test loss on the training dataset size (see
Fig.3 (b)). Both the training and test datasets were assembled by mixing the
subsets with $\sigma_{data}\geq 2.0\ \mu$m. The model is then trained on data
chunks of increasing size $K_{train}$, while keeping the test dataset size
fixed to $K_{test}=300$. This sequence is repeated 4 times, choosing new data
sets at each time to compute the standard deviation of the test loss. Even
with $K_{train}=8$, the test loss is already below $3\%$ and the improvement
for $K_{train}\sim 10^{2}$ is around half a percent. Adding even more data
points does not appreciably improve the prediction quality. This result is
particularly important since we aim for methods that are readily applicable to
atomic density data measurement. In a typical experiment with trapped
ultracold atomic gases, taking more than $K\sim\mathcal{O}(10^{2})$ atom
densities pictures (with averaging over a few shots each) for potential
optimization is a time-costly procedure. The fact that the model we developed
can be trained with less than 100 potentials means that it is in principle
possible to employ a data set obtained by atom-density estimation of the
potential [6].

To further analyze the learning process, we tested the dependence of the test
loss on the cutoff frequency (see Fig.3 (c,d)) of both the training and test
data sets. The frequency subsets are not mixed. At each time, we choose
independently the cutoff frequency of the training $\sigma_{train}$ and test
$\sigma_{test}$ datasets. The dataset sizes for training $K_{train}=160$ and
for testing $K_{test}=39$ are fixed, and this sequence is repeated 5 times.
Fig.3 (c) shows qualitatively the test loss as a function of $\sigma_{train}$
and $\sigma_{test}$. The best results are obtained close to the diagonal, that
is, where $\sigma_{train}\sim\sigma_{test}$, so that training and test data
look most similar. The test loss becomes worse when $\sigma_{train}$ or
$\sigma_{test}$ are around or below the DMD pixel size (magenta line on the
plot), indicating a possible mismatch between the behavior of the system and
the ML model on the scale of the DMD pixel. For a more quantitative
interpretation, we show two different regimes in Fig.3 (d), where a
representative scenario of the training on low frequencies (dark blue) is
compared to another curve representing the high-frequency case. In the first
case, the test loss exhibits a low plateau above the DMD pixel size and an
abrupt increase below, corresponding to a breakdown of the low-frequency
trained model in the high-frequency regime. The second curve instead shows how
including the high frequencies during the training does not solve the problem,
as now the prediction performance severely deteriorates along the whole
frequency range.

##  IV Potential optimization

Figure 5: Optimized effective potentials assuming a 10 Hz harmonic trap. In
the first column, we show the optimized effective potentials $U$ with the
offline ILC, online ILC, and the heuristic online optimization (blue, orange,
and green solid lines, respectively). The black dotted line shows the
effective target potential for atoms (see main text for details). The insets
show DMD patterns for optimization with online ILC. All three methods give a
good qualitative match with the target potentials. The second column depicts
the differences $\delta U$ between the target potential and the optimized
potential for all 3 methods. The online ILC and the heuristic online
optimization give a compatible level of optimization while the offline ILC
exhibits a larger deviation from the target. In the third column, the error
histories for online methods are shown, computed at each point using
Eq.eq:experr. The blue constant line indicates the level of offline
optimization, showing that the offline prediction can save 50-90 heuristic
iterations. The online ILC saves more than 100 iterations compared to the
heuristic method to reach a similar error level. The dashed purple line shows
the shot-to-shot measurement error.

In this section, we introduce an optimization algorithm based on iterative
learning control (ILC) methods inspired by [17]. The algorithm was used to
optimize different target potentials on the experimental setup and we compared
the results with the heuristic optimization method described in [6] (see
Sec.VI for details).

###  IV.1 Online and Offline ILC

ILC methods employ measurements of the considered output trajectory to
iteratively solve a reference tracking problem, i.e., to find an input
trajectory such that the output of a system follows a desired target
trajectory as closely as possible, even in the presence of model errors and
uncertainties. The price to pay for this property is the requirement of
running in a kind of feedback loop using experimental data. Therefore, we show
how to employ ILC algorithms either using feedback from the physical
experiment referred to as “online” ILC or from its digital counterpart, the
physics-inspired model presented in Sec.III, referred to as “offline” ILC. See
Fig.4 for a schematic of the ILC algorithm. The only difference in the
algorithms is that in the first case the potential $\boldsymbol{V}$ is
measured, while in the second it is predicted by the model. This way, previous
(training) data is structured by the ML model while we can seamlessly improve
beyond the predictive capability of the model through further online
iterations.

Let us call $\boldsymbol{\nu}^{n}$ the virtual input at the $n$-th iteration,
and $\boldsymbol{e}^{n}=\boldsymbol{V}^{tar}-\boldsymbol{V}^{n}$ the deviation
from the target. Following standard ILC approaches, the correction of the
virtual input is obtained by convolution, denoted by $*$ (see Sec.VI), with an
update filter $\boldsymbol{L}_{\nu}$, and is then added to the old virtual
input to update it, i.e.,

$\boldsymbol{\nu}^{n+1}=\boldsymbol{\nu}^{n}+\boldsymbol{L}_{\nu}*\boldsymbol{e}^{n}\,.$
(5)

The process is repeated until either convergence or the desired error level is
reached. In order to choose an appropriate update filter, we approximate the
system in a linear and time-invariant form

$\boldsymbol{V}\approx\boldsymbol{g}_{z}*\boldsymbol{\nu}\,.$ (6)

In practice, $\boldsymbol{V}$ is recorded for a trial input, and then fitted
using a Gaussian guess for the longitudinal PSF
$g_{z}(z)=A\exp(-(z-z_{0})^{2}/(2\sigma^{2}))$, to obtain estimates for the
parameters $z_{0},\sigma,A$. We then use a pseudo-inversion-based update
filter [26]

$\boldsymbol{L}_{\nu}=\frac{\overline{\boldsymbol{G}}}{\gamma_{\nu}+\overline{\boldsymbol{G}}\boldsymbol{G}}$
(7)

where $\boldsymbol{G}=\mathcal{F}[\boldsymbol{g}_{z}]$ is the discrete Fourier
transform of the Gaussian PSF, $\overline{\boldsymbol{G}}$ its complex
conjugate and all the operations are to be understood element-wise. Here,
$\gamma_{\nu}>0$ is a regularization parameter of the system inversion which
effectively reduces the high-frequency content of the input updates and,
therefore, of the explored virtual inputs $\boldsymbol{\nu}^{n}$. In the
presence of measurement noise, optimal choices for $\gamma_{\nu}$ are
ultimately given by the experimental signal-to-noise ratio [27]. The results
shown in the main text are obtained using $\gamma_{\nu}=0.1\
\max\limits_{i}|G_{i}|^{2}$.

###  IV.2 Experimental results

We tested both the offline and online ILC algorithms by optimizing three
target potentials, using the heuristic method as a reference. The results are
shown in Fig.5. In order to emulate an experiment with a trapped BEC, we
represent the longitudinal magnetic trap with a harmonic potential
$V_{mag}(z)=\frac{1}{2}m\omega^{2}z^{2}$ with $\omega=2\pi\times 10$ Hz and
$m$ the mass of a 87Rb atom. The effective potential $U(z)$ as experienced by
the atoms is the sum of the magnetic trap $V_{mag}(z)$ and the optical dipole
potential $V(z)$, which is estimated from measured intensity (see Sec.VI for
details). We design targets so that the effective potential $U(z)$ will be
either constant, sinusoidal, or linear. All methods start with the DMD
completely turned off.

Qualitatively all the methods are successful, as the final results are barely
distinguishable from the target potential. For quantitative evaluation of the
optimization performance, we compute the locally normalized root mean square
error (RMS):

$\epsilon_{RMS}=\sqrt{\sum_{i=1}^{N_{V}}\left(\frac{V^{tar}_{i}-V^{n}_{i}}{V^{tar}_{i}}\right)^{2}}\,.$
(8)

We say that the normalization is local because it is computed point by point.
As a measure of the minimum error that can be achieved due to shot-to-shot
fluctuations, we acquire 100 pictures of the optimized potential and compute
the average $\epsilon_{RMS}$, substituting $V^{tar}$ by the average potential
over the sample.

The measured optimization error of the offline ILC is shown with a constant
blue line in Fig.5. This error level is already comparable to what we obtain
with $\sim 100$ iterations of the heuristic algorithm, which would be
equivalent to more than 1 hour of experimental time per potential shape for
BEC experiments (and more than 5 hours with averaging over at least 5 images).
For a quasi-1D BEC with an atomic density of $100$ $\mu m^{-1}$, atomic shot
noise is at the level of 10 % [28]. Therefore, the offline ILC optimization
alone already provides a measurement-limited error. Moreover, since any
constant offset in optimized potential does not have an effect on the atoms,
the effective error is lower than $\epsilon_{RMS}$. Neglecting constant
offsets, the energy scale of the remaining imperfections is around 10 Hz,
which is 1$\%$ of the typical chemical potential ($\sim$ 1 kHz). This means
that offline optimization is already potentially useful for practical
applications in experiments with ultracold gases.

The online ILC and heuristic algorithms both converge to the same error level.
Yet online ILC reaches convergence in around $\sim 10$ iterations while the
heuristic algorithm needs around $\sim 100$ iterations, which is a great
advantage for schemes incorporating experimental feedback. Moreover, the ILC
algorithm does not need to manually select an optimization schedule, resulting
in increased flexibility and bypassing tedious parameter tuning. Convergence
time can be decreased even more by using the result of offline optimization as
an initial guess for online ILC. In this case, we find that the online ILC
reaches convergence in only a few iterations, resulting in an even larger
speed-up compared to the heuristic method. In any case, the final error is far
below the atomic shot noise, so it is hardly accessible in static BEC
configurations.

We ran the same tests using a single-frequency laser as a light source. The
physics-inspired model performs worse in the prediction of laser-generated
potentials, therefore the output of offline ILC is noticeably worse than the
predictions for SLED. On the other hand, the online ILC algorithm optimizes
the optical potentials created with the laser to the same level as with the
SLED. Based on these findings, we can state that our method is perfectly
suited for experimental setups employing SLEDs as well as single-frequency
lasers.

##  V Conclusions

In this paper, we presented our experimental setup for generating and
efficiently optimizing 1D optical potentials. We combined a digital
micromirror device for potential shaping control together with a SLED light
source. We performed measurements estimating the quantitative difference
between SLED’s and laser’s coherence properties and showing the advantages of
using SLED due to its generally linear behavior.

We have implemented learning algorithms that enable efficient optimization of
optical potentials. We have shown how to build a physics-inspired model, which
acts as a digital twin of the experimental setup. The model is able to
recreate the main features of the optical system based on a small set of
experimental data without the need to use deep (neural) networks and large
training data sets, with the advantage of saving experimental time.

The application of the Iterative Learning Control optimization method provides
a more than ten-fold speed-up compared to heuristic approaches. The ILC
algorithms used offline with the trained models are able to optimize optical
potentials with a precision acceptable for most experiments with trapped 1D
ultracold gases. Using online ILC with experimental feedback allows us to
optimize optical potentials to error levels comparable to measurement error,
giving a ten-fold speed-up compared to more straightforward heuristic
algorithms. By combining both ILC strategies, namely using the result of an
optimized digital twin configuration as an initial guess for online ILC, we
get the optimized potential with just a few experimental iterations.

Regarding optimization performance for the SLED and laser, we find that the
SLED outperforms the laser using offline ILC. However, the online ILC performs
equally well for laser and SLED giving the same level of optimization error.
The model we developed combined with Iterative Learning Control provides a
very fast way to optimize optical potentials with a DMD which might be used in
a large variety of experimental setups. Our work offers a prospect for fast
optimization of optical dipole potentials which is extremely important for
time-costly experiments or for very large sequences of patterns in dynamic
situations.

##  VI Appendix

###  VI.1 DMD mount

Due to the specific construction of the Texas Instruments DLP650LNIR DMD’s
micromirror control mechanism [29], the DMD is mounted 45∘ rotated so all the
optical elements are placed in one plane. Any pattern getting rotated 45∘
right before projecting on the DMD. We verified that during optimization the
rotation only leads to distortions of the potential which are below the
resolution and therefore optimization is not affected.

###  VI.2 Intensity to optical dipole potential conversion

The number of pixels in the output $V_{i}$ as obtained from the camera does
not necessarily coincide with the number of pixels in the input $\nu_{i}$. In
order to use the model, we first interpolate $V_{i}$ to the input grid size,
using the interp1 function in Matlab.

We assume the relation between light intensity and optical dipole potential to
be linear $V=\alpha_{V}I$ with $\alpha_{V}$ as found in [30]. Since we work
with red-detuned light, $\alpha_{V}$ is negative. We also suppose the relation
between the CCD readout $R$ and intensity to be linear. We then compute
$\alpha_{CCD}=Ir_{pow}/R$ by measuring the light intensity with a power meter.
To not saturate the CCD we use a reduced amount of light intensity. To
calculate the finally expected dipole potentials we employ a factor
$r_{pow}=I_{full}/I_{low}$ that accounts for the source operating at low
power.

###  VI.3 Mathematical details

This paper heavily relies on the discretization of functions of a real
variable $f(z)$ in order to obtain finite size vectors. If we define a
coordinate grid $z_{i}=(i-1)\Delta-\bar{z}$ for $i=1,...,n$, then we refer to
any discretized function as $f_{i}=f(z_{i})$ and we denote with
$\boldsymbol{f}$ the $\mathbb{R}^{n}$ vector whose elements are $f_{i}$.

Let us call $\mathcal{F}$ the discrete Fourier transform acting on a vector of
size $n$, and $\mathcal{F}^{-1}$ its inverse: F[a]_k &= ∑_j=1^n a_j e^-2
πin(j-1)(k-1)
F^-1[b]_k = 1n ∑_j=1^n b_j e^2 πin(j-1)(k-1) Then, we can define the discrete
convolution of two vectors $\boldsymbol{a}*\boldsymbol{b}$ as

$\boldsymbol{a}*\boldsymbol{b}=\mathcal{F}^{-1}[\mathcal{F}[\boldsymbol{a}]\mathcal{F}[\boldsymbol{b}]]$
(9)

where the product on the right-hand side is element-wise.

###  VI.4 Heuristic algorithm

In order to assess the advantages of the ILC methods, we employed an adapted
version of the heuristic algorithm described in [6] as a reference. It is an
iterative algorithm that updates the state of each pixel based on the local
differences between measured and target potentials. The optimization happens
in two phases, the first fast but rough and the second slower but more
precise. During the rough phase (see Fig.5 first $\sim$20 iterations) in each
column DMD pixels are turned on until the difference gets lower than the
chosen threshold. During the precise phase pixels can be moved away or turned
off.

### Acknowledgments

We would like to thank Mohammadamin Tajik and João Sabino for the discussions
and technical support. This project was funded by the DFG/FWF CRC 1225
’Isoquant’, Project-ID 273811115, by the Austrian Science Fund (FWF) P
36236-N, by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – Cluster of Excellence Matter
and Light for Quantum Computing (ML4Q) EXC 2004/1 – 390534769, and from the
German Federal Ministry of Education and Research through the funding program
quantum technologies—from basic research to market under the project FermiQP,
13N15891. M.P. has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No 101032523.

## References

  * Zhao  _et al._ [2018] G. Zhao, C. Zheng, C. Kuang, R. Zhou, M. M. Kabir, K. C. Toussaint, W. Wang, L. Xu, H. Li, P. Xiu, and X. Liu, Phys. Rev. Lett. 120, 193901 (2018).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Science 357, 995 (2017).
  * Henderson  _et al._ [2009] K. Henderson, C. Ryu, C. MacCormick, and M. G. Boshier, New Journal of Physics 11, 043030 (2009).
  * Gauthier  _et al._ [2016] G. Gauthier, I. Lenton, N. M. Parry, M. Baker, M. J. Davis, H. Rubinsztein-Dunlop, and T. W. Neely, Optica 3, 1136 (2016).
  * Ville  _et al._ [2017] J. L. Ville, T. Bienaimé, R. Saint-Jalm, L. Corman, M. Aidelsburger, L. Chomaz, K. Kleinlein, D. Perconte, S. Nascimbène, J. Dalibard, and J. Beugnon, Phys. Rev. A 95, 013632 (2017).
  * Tajik  _et al._ [2019] M. Tajik, B. Rauer, T. Schweigler, F. Cataldini, J. Sabino, F. S. Møller, S.-C. Ji, I. E. Mazets, and J. Schmiedmayer, Opt. Express 27, 33474 (2019).
  * Liang  _et al._ [2010] J. Liang, J. Rudolph N. Kohn, M. F. Becker, and D. J. Heinzen, Appl. Opt. 49, 1323 (2010).
  * Gaunt  _et al._ [2013] A. L. Gaunt, T. F. Schmidutz, I. Gotlibovych, R. P. Smith, and Z. Hadzibabic, Phys. Rev. Lett. 110, 200406 (2013).
  * Navon  _et al._ [2021] N. Navon, R. P. Smith, and Z. Hadzibabic, Nature Physics 17, 1334–1341 (2021).
  * Gauthier  _et al._ [2021] G. Gauthier, T. A. Bell, A. B. Stilgoe, M. Baker, H. Rubinsztein-Dunlop, and T. W. Neely,  Advances In Atomic, Molecular, and Optical Physics, 70, 1 (2021).
  * Zupancic  _et al._ [2016] P. Zupancic, P. M. Preiss, R. Ma, A. Lukin, M. E. Tai, M. Rispoli, R. Islam, and M. Greiner, Opt. Express 24, 13881 (2016).
  * Smith  _et al._ [2021] A. Smith, T. Easton, V. Guarrera, and G. Barontini, Phys. Rev. Research 3, 033241 (2021).
  * Floyd and Steinberg [1976] R. W. Floyd and L. Steinberg, Proceedings of the Society for Information Display 17, 75 (1976).
  * Dorrer and Zuegel [2007] C. Dorrer and J. D. Zuegel, J. Opt. Soc. Am. B 24, 1268 (2007).
  * Caramazza  _et al._ [2019] P. Caramazza, O. Moran, R. Murray-Smith, and D. Faccio, Nature Communications 10, 2029 (2019).
  * Bristow  _et al._ [2006] D. Bristow, M. Tharayil, and A. Alleyne, IEEE Control Systems Magazine 26, 96 (2006).
  * Deutschmann-Olek  _et al._ [2022] A. Deutschmann-Olek, M. Tajik, M. Calzavara, J. Schmiedmayer, T. Calarco, and A. Kugi,  (2022), arXiv:2208.00706 .
  * Sackett and Deissler [2003] C. A. Sackett and B. Deissler, Journal of Optics B: Quantum and Semiclassical Optics 6, 15 (2003).
  * Campo and Boshier [2012] A. d. Campo and M. G. Boshier, Scientific Reports 2, 648 (2012).
  * Gluza  _et al._ [2021] M. Gluza, J. Sabino, N. H. Ng, G. Vitagliano, M. Pezzutto, Y. Omar, I. E. Mazets, M. Huber, J. Schmiedmayer, and J. Eisert, PRX Quantum 2, 030310 (2021).
  * Manz  _et al._ [2010] S. Manz, R. Bücker, T. Betz, C. Koller, S. Hofferberth, I. E. Mazets, A. Imambekov, E. Demler, A. Perrin, J. Schmiedmayer, and T. Schumm, Phys. Rev. A 81, 031610(R) (2010).
  * Hitzenberger  _et al._ [1999] C. K. Hitzenberger, M. Danner, W. Drexler, and A. F. Fercher, Journal of Modern Optics 46, 1763 (1999).
  * Redding  _et al._ [2012] B. Redding, M. A. Choma, and H. Cao, Nature Photonics 6, 355–359 (2012).
  * Vapnik [2006] V. N. Vapnik,  Estimation of Dependences Based on Empirical Data, 1 (2006).
  * Geman  _et al._ [1992] S. Geman, E. Bienenstock, and R. Doursat, Neural Computation 4, 1 (1992).
  * Ghosh and Paden [2002] J. Ghosh and B. Paden, IEEE Transactions on Automatic Control 47, 831 (2002).
  * Deutschmann-Olek  _et al._ [2021] A. Deutschmann-Olek, G. Stadler, and A. Kugi, IEEE Transactions on Automatic Control 66, 3856 (2021).
  * Berrada  _et al._ [2013] T. Berrada, S. van Frank, R. Bücker, T. Schumm, J.-F. Schaff, and J. Schmiedmayer, Nature Communications 4, 2077 (2013).
  * Lee [2008] B. Lee, Texas Instruments , 1 (2008).
  * Grimm  _et al._ [2000] R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov, Academic Press Advances In Atomic, Molecular, and Optical Physics, 42, 95 (2000).


overestimated by a factor of 2.

By fitting the experimental data with $M\simeq<A_{part}>^{\alpha}$ the HADES
Collaboration extracted $\alpha_{HADES}=1.45\pm 0.06$ for $K+,K^{-},K^{0}_{s}$
and $\Lambda+\Sigma^{0}$ Adamczewski-Musch:2018xwg . On the other hand in Ref.
Adamczewski-Musch:2018xwg  the comparison of the HADES results to the model
predictions (by HSD, IQMD (with and without medium effects) as well as by
UrQMD) has been shown. A similar tension in multiplicity of $K^{0}_{s}$
between all three models and experimental data has been reported. We note that
in Ref. Adamczewski-Musch:2018xwg  the comparison has been done using the HSD
model with G-matrix from Ref. laura03 . There the HSD results for $K^{0}_{s}$
and $\Lambda+\Sigma^{0}$ have been fitted with $\alpha=1.35\pm 0.02$ (cf.
Table II in Ref. Adamczewski-Musch:2018xwg ). Our new results for $K^{0}_{s}$
and $\Lambda+\Sigma^{0}$ are consistent with this finding, too. Moreover, the
$K^{+}$ and $K^{-}$ yields can be fitted with the same $\alpha$ which is in
agreement with KaoS results within errorbars.

###  VII.6 Isospin ratios of $K^{+}$ over $K^{0}$

Finally we check the isospin dependence of the $K^{+}$ to $K^{0}$ production.
As advocated in Refs. DiToro06 ; Li:2008gp  the inclusive $K^{+}/K^{0}$ ratio
at subthreshold energies for kaon production is sensitive to the symmetry
energy of the nucleon-nucleon potential and thus can be used to obtain
experimental information on this issue. On the other hand, the isospin
decomposition of strangeness production is defined to a large extend by
hadronic inelastic and charge exchange reactions with baryons and mesons.
Their cross sections are not well known experimentally and thus the transport
calculations have to rely on phenomenological models. In this respect the
experimental data on the $K^{+}$ to $K^{0}$ production in collisions of nuclei
of different sizes provide an opportunity to access the production channels,
too.

Experimentally the isospin dependence of $K^{+}$ and $K^{0}$ production at
near sub-threshold energies has been investigated by the FOPI Collaboration in
Ref. Lopez:2007rh  which measured the ratios of $K^{+}$ to $K^{0}$ mesons for
Ru+Ru to Zr+Zr collisions at the kinetic energy 1.528 A GeV (i.e. at
$\sqrt{s}=2.527$ GeV, while the threshold for kaon production in $pp$
collisions is $\sqrt{s_{0}}=2.546$ GeV).

a Ratio |  PHSD | FOPI
---|---|---
$K^{+}_{Ru}/K^{+}_{Zr}$ | 1.09 | $1.06\pm 0.07(stat.)\pm 0.09(syst.)$
a $K^{0}_{Ru}/K^{0}_{Zr}$ | 0.96 | $0.94\pm 0.12(stat.)\pm 0.06(syst.)$
a ${(K^{+}/K^{0})_{Ru}\over(K^{+}/K^{0})_{Zr}}$ | 1.14 | $1.13\pm 0.16(stat.)\pm 0.12(syst.)$
Table 3: Comparison of the inclusive ratios of $K^{+}$ and $K^{0}$ mesons for
Ru+Ru to Zr+Zr collisions at the kinetic energy 1.528 A GeV in the PHSD versus
the FOPI data Lopez:2007rh .

In Table 3 we show the comparison of the inclusive ratios of $K^{+}$ and
$K^{0}$ mesons for Ru+Ru to Zr+Zr collisions - $K^{+}_{Ru}/K^{+}_{Zr}$,
$K^{0}_{Ru}/K^{0}_{Zr}$, ${(K^{+}/K^{0})_{Ru}\over(K^{+}/K^{0})_{Zr}}$ at the
kinetic energy 1.528 A GeV from the PHSD versus the FOPI data Lopez:2007rh .
One can see that the PHSD results are in a good agreement with the measured
ratios as well as the model calculations from Ref. DiToro06 . We note,
however, that the symmetry energy is not implemented in the potential used in
the PHSD version for this study. Thus, we attribute this agreement to the
consequence of the isospin dependence of strangeness production and
rescattering.

##  VIII Summary

In this study we have investigated strangeness production in heavy ion
collisions at(sub-)threshold energies between 1.0 and 2.0 A GeV within the
microscopic Parton-Hadron-String Dynamics (PHSD 4.5) transport approach,
extended for the incorporation of the in-medium effects for strangeness
production in terms of the state-of-the-art coupled-channel G-matrix approach
for the modification of antikaon properties in a dense and hot medium.

The PHSD is an off-shell transport approach based on the first order gradient
expansion (in phase-space representation) of the Kadanoff-Baym equations which
allow for a propagation of particles with dynamical spectral functions, i.e.
for a microscopic description of the strongly interacting system. It includes
also the possibility to implement detailed balance on the level of
$m\leftrightarrow n$ reactions (although its technical implementation is
restricted to selected channels of interest due to the limitation in
computational resources). For the present study we have incorporated the
detailed balance on $2\leftrightarrow 3$ level for the dominant channels for
strangeness production/absorbtion by baryons ($B=N,\Delta$) and pions:
$B+B\leftrightarrow N+Y+K$ and $B+\pi\leftrightarrow N+N+\bar{K}$, as well as
the non-resonant reactions $N+N\leftrightarrow N+N+\pi$ and
$\pi+N\leftrightarrow N+\pi+\pi$ to have a better control on pion dynamics.

The modification of the strangeness properties in the hot and dens medium,
which is created in the heavy-ion collisions, has been incorporated in the
PHSD in the following way:

The medium effects for antikaons are determined by a self-consistent coupled-
channel G-matrix approach Cabrera:2014lca , based on a chiral Lagrangian,
which provides the spectral function of antikaons as well as the scattering
amplitudes (i.e. in-medium cross sections) as a function of baryon density,
temperature and the three-momentum of the antikaon with respect to the matter
at rest. This novel G-matrix approach Cabrera:2014lca  incorporates the full
self-consistency in $s$\- and $p$-waves in the $\bar{K}N$ interaction at
finite density and temperature based on a SU(3) chiral Lagrangian. It differs
therefore from the previous work of Ref. laura03  where a self-consistent
unitary approach based on the Jülich meson-exchange model Tolos:2000fj ;
Tolos:2002ud  was used.

We point out that the novel G-matrix approach provides a more shallow
dependence of the real part of the complex antikaon self-energy $Re\Sigma$
(and, respectively, $\bar{K}N$ potential) on baryon density than that from
Refs. Tolos:2000fj ; Tolos:2002ud , which has been used in the past in the
off-shell HSD laura03  for antikaon dynamics at SIS energies.

Nuclear medium effects on the production and interactions of kaons are
implemented by a repulsive kaon potential of $V_{K}=+25$ MeV
$(\rho/\rho_{0})$.

We have studied the effects of the (anti)kaon medium modifications on
different observables in heavy-ion collisions at SIS energies measured by the
KaoS, FOPI and HADES Collaborations. Our study confirms the sensitivity of
antikaon observables to the details of the in-medium models: the in-medium
effects on observables obtained with the novel G-matrix are slightly less
pronounced as compared to the previous G-matrix calculations which provided a
much stronger attraction of antikaons in nuclear matter.

Our findings can be summarized as follows:
$\bullet$ The kaon-nuclear potential increases the threshold energy for the
kaon production in $BB$ and $mB$ reactions in a hot and dense medium.
Consequently, the $K^{+},K^{0}$ production in A+A collisions is suppressed
when accounting for the potential compared to the case without potential. We
note, that in A+A collisions at (sub-)threshold energies the kaons are
dominantly produced by the reactions with $\Delta$’s (as $N+\Delta\to Y+K+N$)
and mesons (as $\pi+N\to Y+K$). That enhances kaon production substantially
relative to the production by $N+N$ reactions only. Moreover, production of
kaons below the $N+N$ threshold is possible only in A+A collisions due to the
Fermi-motion of nucleons and by secondary nucleon reactions with pions.

$\bullet$ The self-consistent coupled-channel unitarized scheme based on a
SU(3) chiral Lagrangian leads to a broadening of the $\bar{K}$ spectral
function without drastic changes of the pole mass. Moreover, the strength of
the attractive $\bar{K}N$ potential decreases with increasing three-momentum
of $\bar{K}$ mesons in nuclear matter. The $\bar{K}$ production in A+A
collisions is enhanced by the broad spectral width (which increases with
density) which leads to a reduction of the production threshold. The antikaons
are dominantly produced by the secondary hyperon-meson reactions $Y+m\to
B+\bar{K}$ which are delayed in time. Consequently, the baryon density at the
antikaon production point is lower compared to that at the kaon production
point, that reduces the in-medium effects.

$\bullet$ The kaon-nuclear potential hardens the $p_{T}$-spectra of
$K^{+},K^{0}$ mesons while the in-medium effects in terms of self-energies
within the G-matrix approach soften the $K^{-}$ spectra. These effects are
stronger in larger systems such as Au+Au collisions, where a higher density is
reached and therefore the density gradient (which is proportional to the
force) is larger as compared to smaller systems like C+C collisions.

$\bullet$ For kaons, both collisions with baryons and mesons as well as the
repulsive interaction with the nuclear medium lead to an increase of the
effective temperature $T_{eff}$ of the $p_{T}$-spectrum compared to that at
the kaon production. For antikaons the hardening of the $p_{T}$-spectra due to
collisions and the softening due to the in-medium effects counterbalance
partially.

$\bullet$ The repulsive kaon potential widens the rapidity distribution of
$K^{+}$ mesons, while an attractive antikaon interaction shrinks that of
$K^{-}$ mesons. It is also visible in the enhancement of the polar angular
distribution of $K^{+}$ mesons in forward and backward directions, although
the effect of $K^{+}$ rescattering in the medium plays an even larger role
than the $KN$ potential, especially for heavy systems such as Au+Au. The polar
angular distribution of $K^{-}$ mesons stays rather flat due to the
competition of two effects: the scattering pushes $K^{-}$ mesons backward and
forward similar to $K^{+}$, however, the antikaon potential attracts them to
the middle zone of the colliding system; this push is stronger for central
collisions and heavy systems due to the larger density.

$\bullet$ The collective flow coefficients for $K^{+}$ and $K^{-}$ are
strongly affected by the flow of nucleons, since (anti)kaons are produced by
nucleon scattering and secondary interactions with nucleons. Considering only
scattering and neglecting the kaon-nuclear potential, the directed flow
$v_{1}$ of $K^{+}$ mesons follows that of nucleons such that the $v_{1}$ of
$K^{+}$ mesons is positive at forward rapidity and negative at backward
rapidity. However, the repulsive kaon potential pushes the $K^{+}$ mesons in
the opposite direction to the nucleons and, as a result, the $v_{1}$ of
$K^{+}$ mesons for low transverse momenta increases strongly and becomes
positive, while the $v_{1}$ of $K^{+}$ mesons with large transverse momenta
stays negative and finally changes only slightly.
The elliptic flow $v_{2}$ of $K^{+}$ and $K^{-}$ mesons are also affected by
the flow of nucleons. In non-central heavy-ion collisions, the nucleons gain a
negative elliptic flow due to the geometric shape of the medium and attractive
Skyrme potential. The $v_{2}$ of the produced $K^{+}$ mesons is initially
positive, following that of nucleons, and becomes negative due to the
scattering and repulsive interactions with nucleons. This is also seen in the
azimuthal angular distribution of kaons. On the other hand, the antikaon-
nuclear potential acts in opposite direction: while the $v_{2}$ of the
initially produced antikaons is also positive (as for $K^{+}$) and becomes
also negative due to the rescattering in the matter, the attractive antikaon
potential pushes the $v_{2}$ towards less negative values such that the final
$v_{2}$ of antikaons is close to zero in the present PHSD calculations.

$\bullet$ Moreover, we have investigated the sensitivity of strangeness
production to the nuclear equation-of-state. In line with the previous
findings Aichelin:1986ss ; Fuchs ; Hartnack:2011cn  we have observed an
enhancement of (anti)kaon production by a softening of the EoS. However, we
stress that solid conclusions on the EoS can be obtained only by considering a
variety of observables (including the non-strange hadrons) since the dynamics
of strangeness is tightly bound to the dynamics of the ’bulk’ matter in heavy-
ion collisions and the density of pions.

$\bullet$ Finally, our study with the improved description for the
modifications of the (anti)kaon properties in the hot and dense medium (within
the novel G-matrix approach) as well as an improved general dynamics of the
medium itself (within the PHSD) confirms the previous findings by different
groups (summarized in the review Hartnack:2011cn ) on the observation of the
visible medium effects in heavy-ion collisions at SIS energies.

## Acknowledgements

The authors acknowledge valuable discussions with Wolfgang Cassing and his
suggestions. We are grateful to Daniel Cabrera for providing us the G-matrix
calculations. We are also thankful for inspiring discussions with Christoph
Blume, Laura Fabbietti, Yvonne Leifels, Manuel Lorenz, Pierre Moreau and Iori
Vassiliev. Furthermore, we acknowledge support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation): grant BR 4000/7-1,
by the Russian Science Foundation grant 19-42-04101 and by the GSI-IN2P3
agreement under contract number 13-70. L.T. acknowledges support from the
Deutsche Forschungsgemeinschaft (DFG, German research Foundation) under the
Project Nr. 383452331 (Heisenberg Programme), Nr. 411563442 (Hot Heavy
Mesons). Her research has been also supported by the Spanish Ministerio de
Economía y Competitividad under contract FPA2016-81114-P and Ministerio de
Ciencia e Innovación under contract PID2019-110165GB-I00. This project has,
furthermore, received funding from the European Union’s Horizon 2020 research
and innovation program under grant agreement No 824093 (STRONG-2020). Also we
thank the COST Action THOR, CA15213 and CRC-TR 211 ’Strong-interaction matter
under extreme conditions’- project Nr.315477589 - TRR 211. The computational
resources have been provided by the LOEWE-Center for Scientific Computing and
the ”Green Cube” at GSI, Darmstadt.

##  Appendix A Cross sections for (anti)kaon production

In this Appendix we collect the parametrizations of the cross sections for
$K,\bar{K}$ production by $BB,mB$ and $mm$ reactions in the vacuum used in the
PHSD. The theoretical origin of these parametrizations as well as the
comparison to the available experimental data have been discussed in the
original publications (and the references therein) Cass97 ; brat97 ; CB99 ;
laura03 ; Mishra:2004te ; Hartnack:2011cn ; Kolomeitsev:2004np  and are not
repeated here. The in-medium modifications of (anti)kaon production cross
sections are explained in Section V.

###  A.1 $N+N\rightarrow N+Y+K$

The reaction cross sections for the channels $N+N\rightarrow N+Y+K$ can be
approximated by Sibirtsev:1995xb ,

$\displaystyle\sigma_{pp\rightarrow p\Lambda K^{+}}$ $\displaystyle=$
$\displaystyle
0.732\bigg{(}1-\frac{s_{01}}{s}\bigg{)}^{1.8}\bigg{(}\frac{s_{01}}{s}\bigg{)}^{1.5}~{}[\rm
mb],$ $\displaystyle\sigma_{pp\rightarrow p\Sigma^{+}K^{0}}$ $\displaystyle=$
$\displaystyle
0.339\bigg{(}1-\frac{s_{02}}{s}\bigg{)}^{2.25}\bigg{(}\frac{s_{02}}{s}\bigg{)}^{1.35}~{}[\rm
mb],$ $\displaystyle\sigma_{pp\rightarrow p\Sigma^{0}K^{+}}$ $\displaystyle=$
$\displaystyle
0.275\bigg{(}1-\frac{s_{02}}{s}\bigg{)}^{1.98}\bigg{(}\frac{s_{02}}{s}\bigg{)}~{}[\rm
mb],$ (54)

where $\sqrt{s_{01}}=m_{\Lambda}+m_{p}+m_{K}$ and
$\sqrt{s_{02}}=m_{\Sigma}+m_{p}+m_{K}$. For simplicity, an isospin-averaged
cross section is introduced by multiplying with a factor of 1.5 the cross
sections of Eq. (54):

$\displaystyle\sigma_{NN\rightarrow N\Lambda K}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\sigma_{pp\rightarrow p\Lambda K^{+}},$ (55)
$\displaystyle\sigma_{NN\rightarrow N\Sigma K}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\bigg{(}\sigma_{pp\rightarrow
p\Sigma^{+}K^{0}}+\sigma_{pp\rightarrow p\Sigma^{0}K^{+}}\bigg{)}.$

The different isospin channels are then evenly weighted, for example,

$\displaystyle\sigma_{pn\rightarrow n\Lambda K^{+}}=\sigma_{pn\rightarrow
p\Lambda K^{0}}=\frac{1}{2}\sigma_{NN\rightarrow N\Lambda K},$
$\displaystyle\sigma_{pn\rightarrow n\Sigma^{0}K^{+}}=\sigma_{pn\rightarrow
n\Sigma^{+}K^{0}}=\sigma_{pn\rightarrow p\Sigma^{-}K^{+}}$
$\displaystyle=\sigma_{pn\rightarrow
p\Sigma^{0}K^{0}}=\frac{1}{4}\sigma_{NN\rightarrow N\Sigma K}.$ (56)

The $\Delta$ baryon is assumed to have the same cross section for strangeness
production as a nucleon:

$\displaystyle\sigma_{N\Delta\rightarrow N\Lambda
K}=\sigma_{\Delta\Delta\rightarrow N\Lambda K}=\sigma_{NN\rightarrow N\Lambda
K},$ (57)

where the production thresholds are modified accordingly depending on $\Delta$
mass.

###  A.2 $N+N\rightarrow N+N+K+\bar{K}$

As the center of mass energy $\sqrt{s}$ increases, a strange meson pair can be
produced in the channel $N+N\rightarrow N+N+K+\bar{K}$ with a cross section
given by Sibirtsev:1996rh :

$\displaystyle\sigma_{NN\rightarrow
NNK\bar{K}}=1.5\bigg{(}1-\frac{s_{0}}{s}\bigg{)}^{3.17}\bigg{(}\frac{s_{0}}{s}\bigg{)}^{1.96}~{}[\rm
mb],$ (58)

where $\sqrt{s_{0}}=2m_{N}+m_{K}+m_{\bar{K}}$. We assume the same form of the
cross section if the initial nucleon is replaced by a $\Delta$ baryons, while
baryons in the final state are always $n$ or $p$.

###  A.3 $m+N\rightarrow K+Y,~{}K+\bar{K}+N$

The cross section for $\Lambda$ production ($\pi N\rightarrow K\Lambda$) is
given by Huang:1994xq

$\displaystyle\sigma_{\pi^{-}p\rightarrow\Lambda
K^{0}}=\frac{0.007665(\sqrt{s}-\sqrt{s_{o}})^{0.1341}}{(\sqrt{s}-1.72)^{2}+0.007826}~{}[\rm
mb],$ (59)

where $\sqrt{s_{0}}$ is the threshold energy. Using isospin relations
Huang:1994xq  we find

$\displaystyle\sigma_{\pi^{0}p\rightarrow\Lambda K^{+}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{0}n\rightarrow\Lambda
K^{0}}=\frac{1}{2}\sigma_{\pi^{-}p\rightarrow\Lambda K^{0}},$ (60)
$\displaystyle\sigma_{\pi^{+}n\rightarrow\Lambda K^{+}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{-}p\rightarrow\Lambda K^{0}}.$ (61)

In heavy-ion collisions many baryons are in resonance states. The contribution
of resonances to the $\Lambda$ production through $\pi+\Delta$ scattering is
included by substituting $p$ by a $\Delta^{+}$ and $n$ by a $\Delta^{0}$ and
assuming

$\displaystyle\sigma_{\pi^{-}\Delta^{++}\rightarrow\Lambda
K^{+}}=\sigma_{\pi^{+}\Delta^{-}\rightarrow\Lambda
K^{0}}=\sigma_{\pi^{-}p\rightarrow\Lambda K^{0}}.$ (62)

The cross sections for $\Sigma$ baryon production ($\pi N\rightarrow K\Sigma$)
are given by Tsushima:1994rj ; Tsushima:1994pv

$\displaystyle\sigma_{\pi^{+}p\rightarrow\Sigma^{+}K^{+}}$ $\displaystyle=$
$\displaystyle\frac{0.03591(\sqrt{s}-\sqrt{s_{o}})^{0.9541}}{(\sqrt{s}-1.89)^{2}+0.01548}$
$\displaystyle+$
$\displaystyle\frac{0.1594(\sqrt{s}-\sqrt{s_{o}})^{0.01056}}{(\sqrt{s}-3.0)^{2}+0.9412}~{}[\rm
mb],$ $\displaystyle\sigma_{\pi^{-}p\rightarrow\Sigma^{-}K^{+}}$
$\displaystyle=$
$\displaystyle\frac{0.009803(\sqrt{s}-\sqrt{s_{o}})^{0.6021}}{(\sqrt{s}-1.742)^{2}+0.006583}$
$\displaystyle+$
$\displaystyle\frac{0.006521(\sqrt{s}-\sqrt{s_{o}})^{1.4728}}{(\sqrt{s}-1.94)^{2}+0.006248}~{}[\rm
mb],$ $\displaystyle\sigma_{\pi^{+}n\rightarrow\Sigma^{0}K^{+}}$
$\displaystyle=$
$\displaystyle\frac{0.05014(\sqrt{s}-\sqrt{s_{o}})^{1.2878}}{(\sqrt{s}-1.73)^{2}+0.006455}~{}[\rm
mb],\,\,\,~{}~{}~{}~{}~{}$ (65)
$\displaystyle\sigma_{\pi^{0}n\rightarrow\Sigma^{-}K^{+}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{+}n\rightarrow\Sigma^{0}K^{+}},$ (66)
$\displaystyle\sigma_{\pi^{0}p\rightarrow\Sigma^{0}K^{+}}$ $\displaystyle=$
$\displaystyle\frac{0.003978(\sqrt{s}-\sqrt{s_{o}})^{0.5848}}{(\sqrt{s}-1.74)^{2}+0.00667}$
$\displaystyle+$
$\displaystyle\frac{0.04709(\sqrt{s}-\sqrt{s_{o}})^{2.165}}{(\sqrt{s}-1.905)^{2}+0.006358}~{}[\rm
mb],$

and isospin symmetry is assumed. Consequently we find

$\displaystyle\sigma_{\pi^{-}n\rightarrow\Sigma^{0}K^{0}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{+}p\rightarrow\Sigma^{+}K^{+}},$
$\displaystyle\sigma_{\pi^{+}n\rightarrow\Sigma^{+}K^{0}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{-}p\rightarrow\Sigma^{-}K^{+}},$
$\displaystyle\sigma_{\pi^{-}p\rightarrow\Sigma^{0}K^{0}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{0}p\rightarrow\Sigma^{+}K^{0}}=\sigma_{\pi^{+}n\rightarrow\Sigma^{0}K^{+}},$
$\displaystyle\sigma_{\pi^{0}n\rightarrow\Sigma^{0}K^{0}}$ $\displaystyle=$
$\displaystyle\sigma_{\pi^{0}p\rightarrow\Sigma^{0}K^{+}},$ (68)

As for the $\Lambda$ production, an initial $p$ can be substituted by a
$\Delta^{+}$ and an initial $n$ by a $\Delta^{0}$. $\Delta^{++}$ and
$\Delta^{-}$ are excluded for $\Sigma$ production.

If the collision energy $\sqrt{s}$ exceeds 1.7 GeV, the cross section for
$\Sigma+K$ production in Eqs. (A.3) to (68) is modified to

$\displaystyle\sigma_{\pi N\rightarrow\Sigma
K}(\sqrt{s})=A\sqrt{s}-\sigma_{0}-\sigma_{\pi N\rightarrow\Lambda K}$
$\displaystyle-\sigma_{\pi N\rightarrow NK\bar{K}}~{}[\rm mb],$ (69)

where $A=1.0$ (${\rm mb/GeV}$) and

$\displaystyle\sigma_{0}=1.7~{}{\rm mb}-\sigma_{\pi N\rightarrow\Sigma
K}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle-\sigma_{\pi N\rightarrow\Lambda K}-\sigma_{\pi N\rightarrow
NK\bar{K}}\bigg{|}_{\sqrt{s}=1.7~{}{\rm GeV}}.$ (70)

This means that the total cross section for strangeness production in
baryon+meson scattering is a linear function of $\sqrt{s}$ above 1.7 GeV and
smoothly connected at $\sqrt{s}=1.7$ GeV.

The cross section $\sigma_{\pi N\rightarrow NK\bar{K}}$ in Eq. (69) is given
by Sibirtsev:1996rh

$\displaystyle\sigma_{\pi^{-}p\rightarrow
pK^{0}K^{-}}=1.121\bigg{(}1-\frac{s_{0}}{s}\bigg{)}^{1.86}\bigg{(}\frac{s_{01}}{s}\bigg{)}^{2}~{}[\rm
mb],$ (71)

where $\sqrt{s_{0}}=m_{N}+m_{K}+m_{\bar{K}}$, and the other channels with
various isospin combinations are related to $\sigma_{0}$ of Eq. (71) as
follows Sibirtsev:1996rh :

$\displaystyle 2\sigma(\pi^{-}p\rightarrow
pK^{0}K^{-})=2\sigma(\pi^{-}n\rightarrow nK^{0}K^{-})$
$\displaystyle=2\sigma(\pi^{+}p\rightarrow
pK^{+}\bar{K}^{0})=2\sigma(\pi^{+}n\rightarrow nK^{+}\bar{K}^{0})$
$\displaystyle=\sigma(\pi^{+}n\rightarrow
pK^{+}K^{-})=\sigma(\pi^{+}n\rightarrow pK^{0}\bar{K}^{0})$
$\displaystyle=\sigma(\pi^{0}p\rightarrow
nK^{+}\bar{K}^{0})=4\sigma(\pi^{0}p\rightarrow pK^{+}K^{-})$
$\displaystyle=4\sigma(\pi^{0}p\rightarrow
pK^{0}\bar{K}^{0})=\sigma(\pi^{0}n\rightarrow pK^{0}K^{-})$
$\displaystyle=4\sigma(\pi^{0}n\rightarrow
nK^{+}K^{-})=4\sigma(\pi^{0}n\rightarrow pK^{0}\bar{K}^{0})$
$\displaystyle=\sigma(\pi^{-}p\rightarrow
nK^{+}K^{-})=\sigma(\pi^{-}p\rightarrow pK^{0}\bar{K}^{0}).$ (72)

Again an initial $p$ can substituted by a $\Delta^{+}$ and an initial $n$ by a
$\Delta^{0}$ assuming the same form of the cross section, while in the final
state only $p$ or $n$ are admitted.

The inverse reaction $K+Y\rightarrow\pi+N$ is realized by detailed balance and
$N+K+\bar{K}\rightarrow N+\pi$ by using an equation similar to Eq. (90),
including a $\Delta$ baryon in the final state. For example,

$\displaystyle\sigma_{\Lambda
K^{+}\rightarrow\pi^{+}\Delta^{0}}=2\sigma_{\Lambda
K^{+}\rightarrow\pi^{0}\Delta^{+}}=\sigma_{\Lambda
K^{+}\rightarrow\pi^{-}\Delta^{++}}$
$\displaystyle=\bigg{(}\frac{p_{\Delta}}{p_{\Lambda}}\bigg{)}^{2}\sigma_{\pi^{+}\Delta^{0}\rightarrow\Lambda
K^{+}}=\bigg{(}\frac{p_{N}}{p_{\Lambda}}\bigg{)}^{2}\sigma_{\pi^{+}n\rightarrow\Lambda
K^{+}},$ (73)

where $p_{i}$ is the three-momentum of particle $i$ in the c.m. frame.

###  A.4 $\bar{K}+B\rightarrow\pi+Y$

The reactions $\bar{K}+N\rightarrow\pi+Y$ are modeled (also in the vacuum)
according to the G-matrix as indicated by Eqs. (6)-(7). However, the present
coupled channels of G-matrix do not incorporate reactions with baryonic
resonances which are abundant at SIS energies. For such reactions (with
$B=\Delta,N(1440),N(1535)$) we use the parametrizations presented below.

The cross section for $\bar{K}$ absorption by baryons is equal to the total
cross section of $\bar{K}B$ scattering subtracted by the $\bar{K}N$ elastic
scattering cross sectionEfremov:1994ie . The total cross section of $\bar{K}B$
scattering Efremov:1994ie  is

$\displaystyle\sigma_{\bar{K}+B\to\pi+Y}^{tot}=$ $\displaystyle
22.6p_{K}^{-1.14}~{}[\rm mb],$ $\displaystyle p_{K}<0.7~{}{\rm GeV/c},$
$\displaystyle 47.54p_{K}^{-0.94}~{}[\rm mb],$ $\displaystyle 0.7~{}{\rm
GeV/c}\leq p_{K}<1.1~{}{\rm GeV/c},$ $\displaystyle 69.87p_{K}^{-3.1}~{}[\rm
mb],$ $\displaystyle 1.1~{}{\rm GeV/c}\leq p_{K}<1.3~{}{\rm GeV/c},$
$\displaystyle 32.4p_{K}^{-0.17}~{}[\rm mb],$ $\displaystyle 1.3~{}{\rm
GeV/c}\leq p_{K}<10~{}{\rm GeV/c},$ $\displaystyle 22~{}~{}[\rm mb],$
$\displaystyle 10~{}{\rm GeV/c}\leq p_{K},$ (74)

while the $\bar{K}N$ elastic scattering cross section is

$\displaystyle\sigma_{\bar{K}+B\to\pi+Y}^{elast}=$ $\displaystyle
10.58p_{K}^{-0.98}~{}[\rm mb],$ $\displaystyle 0.03~{}{\rm GeV/c}\leq
p_{K}<0.7~{}{\rm GeV/c},$ $\displaystyle 23p_{K}^{1.2}~{}[\rm mb],$
$\displaystyle 0.7~{}{\rm GeV/c}\leq p_{K}<1~{}{\rm GeV/c},$ $\displaystyle
23p_{K}^{-2.6}~{}[\rm mb],$ $\displaystyle 1~{}{\rm GeV/c}\leq
p_{K}<1.5~{}{\rm GeV/c},$ $\displaystyle 9.56p_{K}^{-0.44}~{}[\rm mb],$
$\displaystyle 1.5~{}{\rm GeV/c}\leq p_{K}<20~{}{\rm GeV/c},$ $\displaystyle
2.56~{}~{}[\rm mb],$ $\displaystyle 20~{}{\rm GeV/c}\leq p_{K},$ (75)

where $p_{K}$ is the momentum of the $\bar{K}$ in the $B$ rest frame. The
$\bar{K}$ absorption cross section is weighted by 0.3 for $Y=\Lambda$ and by
0.2 for $Y=\Sigma$.

###  A.5 Elastic and isospin-exchange scattering of kaons

Finally $K+N$ elastic cross section and the cross section for isospin exchange
are parameterized in the units of mb as below

$\displaystyle\sigma_{elastic}^{K^{+}p}$ $\displaystyle=$
$\displaystyle\sigma_{elastic}^{K^{0}n}$ $\displaystyle=$ $\displaystyle
12.3+1.88p_{K}-2.32p_{K}^{2}~{}~{}({\rm p_{K}<2.3~{}GeV/c}),$ $\displaystyle=$
$\displaystyle 5~{}[\rm mb]~{}~{}~{}({\rm p_{K}>2.3~{}GeV/c}).$

Here $p_{K}$ is the three momentum of kaon (GeV/c) in center-of-mass frame,

$\displaystyle\sigma_{elastic}^{K^{+}n}$ $\displaystyle=$
$\displaystyle\sigma_{elastic}^{K^{0}p}$ $\displaystyle=$ $\displaystyle
10.51p_{K}^{0.82}~{}~{}~{}({\rm p_{K}<0.6~{}GeV/c}),$ $\displaystyle=$
$\displaystyle-18.93+76.68p_{K}-66.65p_{K}^{2}+17.74p_{K}^{3}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}({\rm
0.6~{}GeV/c<p_{K}<1.5~{}GeV/c}),$ $\displaystyle=$ $\displaystyle 6~{}[\rm
mb]~{}~{}~{}({\rm 1.5~{}GeV/c<p_{K}}),$

and

$\displaystyle\sigma^{K^{+}n\leftrightarrow K^{0}p}$ $\displaystyle=$
$\displaystyle 5.63+4.996p_{K}-4.519p_{K}^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}({\rm
p_{K}<1.5~{}GeV/c}),$ $\displaystyle=$ $\displaystyle
6.502/p_{K}^{1.932}~{}~{}~{}({\rm 1.5~{}GeV/c<p_{K}}).$

The above cross sections are important for the generation of collective flows
and isospin diffusion of $K^{+}$ and $K^{0}$ in heavy-ion collisions.

##  Appendix B 3-to-2 interactions

The scattering cross section for 2-to-3 processes can be expressed as

$\displaystyle\sigma_{2\rightarrow 3}$ $\displaystyle=$
$\displaystyle\frac{1}{4E_{1}E_{2}v_{\rm
rel}}\int\frac{d^{3}p_{1}^{\prime}}{(2\pi)^{3}2E_{1}^{\prime}}\int\frac{d^{3}p_{2}^{\prime}}{(2\pi)^{3}2E_{2}^{\prime}}\int\frac{d^{3}p_{3}^{\prime}}{(2\pi)^{3}2E_{3}^{\prime}}$
(79) $\displaystyle\times$
$\displaystyle\overline{|M|}^{2}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{1}^{\prime}-p_{2}^{\prime}-p_{3}^{\prime}),$

where $v_{\rm rel}$ is the relative velocity of the incident particles and
$\overline{|M|^{2}}$ is the transition amplitude squared averaged over initial
states:

$\displaystyle\overline{|M|}^{2}=\frac{|M|^{2}}{D_{1}D_{2}}$ (80)

with $D_{1}$, $D_{2}$ denoting the spin-flavor degeneracy of the initial 2
particle states. Assuming that the transition amplitude does not depend on
scattering angle but only on the collision energy $\sqrt{s}$, the cross
section can be rewritten as

$\displaystyle\sigma_{2\rightarrow 3}$ $\displaystyle=$
$\displaystyle\frac{1}{4E_{1}E_{2}v_{\rm
rel}}\frac{|M|^{2}}{D_{1}D_{2}}(PS)_{3}$ (81) $\displaystyle=$
$\displaystyle\frac{1}{4p_{1}^{c.m.}\sqrt{s}}\frac{|M|^{2}}{D_{1}D_{2}}(PS)_{3},$

where $p_{1}^{c.m.}=|\vec{p_{1}}|$ in the c.m. frame and

$\displaystyle(PS)_{3}=\int\frac{d^{3}p_{1}^{\prime}}{(2\pi)^{3}2E_{1}^{\prime}}\int\frac{d^{3}p_{2}^{\prime}}{(2\pi)^{3}2E_{2}^{\prime}}\int\frac{d^{3}p_{3}^{\prime}}{(2\pi)^{3}2E_{3}^{\prime}}$
$\displaystyle\times(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{1}^{\prime}-p_{2}^{\prime}-p_{3}^{\prime}).$
(82)

The three-body phase-space integral $(PS)_{3}$ is simplified in the c.m. frame
of $p_{2}^{\prime}+p_{3}^{\prime}$ as follows:

$\displaystyle(PS)_{3}=\int\frac{d^{3}p_{1}^{\prime}}{(2\pi)^{3}2E_{1}^{\prime}}\int\frac{d^{3}p_{2}^{\prime}}{(2\pi)^{2}4E_{2}^{\prime}E_{3}^{\prime}}$
$\displaystyle\times\delta(E_{1}+E_{2}-E_{1}^{\prime}-E_{2}^{\prime}-E_{3}^{\prime})$
$\displaystyle=\frac{1}{16\pi^{3}\sqrt{s}}\int^{\sqrt{s}-m_{1}^{\prime}}_{m_{2}^{\prime}+m_{3}^{\prime}}dM_{23}p_{1}^{\prime~{}c.m.}p_{2}^{*}\
,$ (83)

where $M_{23}=\sqrt{(p_{2}+p_{3})^{2}}$ and
$p_{1}^{\prime~{}c.m.}=|p_{1}^{\prime}|$ in the c.m. frame of $p_{1}+p_{2}$
and $p_{2}^{*}=|p_{2}^{\prime}|$ in the c.m. frame of
$p_{2}^{\prime}+p_{3}^{\prime}$, i.e.

$\displaystyle
p_{1}^{\prime~{}c.m.}=\frac{\sqrt{\\{s-(m_{1}^{\prime}+M_{23})^{2}\\}\\{s-(m_{1}^{\prime}-M_{23})^{2}\\}}}{2\sqrt{s}},$
$\displaystyle
p_{2}^{*}=\frac{\sqrt{\\{M_{23}^{2}-(m_{2}^{\prime}+m_{3}^{\prime})^{2}\\}\\{M_{23}^{2}-(m_{2}^{\prime}-m_{3}^{\prime})^{2}\\}}}{2M_{23}}.$
(84)

A cross section for 3-to-2 processes cannot be defined but a Lorentz invariant
interaction rate for 3-to-2 processes is given by Cassing:2001ds ; Xu:2004mz ;
Song:2012at

$\displaystyle\frac{\Delta N^{3\rightarrow 2}}{\Delta t\Delta
V}=\frac{1}{D_{1}^{\prime}D_{2}^{\prime}D_{3}^{\prime}}\int\frac{d^{3}p_{1}}{(2\pi)^{3}2E_{1}}$
$\displaystyle\times\int\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\int\frac{d^{3}p_{1}^{\prime}}{(2\pi)^{3}2E_{1}^{\prime}}f_{1}(p_{1}^{\prime})$
$\displaystyle\times\int\frac{d^{3}p_{2}^{\prime}}{(2\pi)^{3}2E_{2}^{\prime}}f_{2}(p_{2}^{\prime})\int\frac{d^{3}p_{3}^{\prime}}{(2\pi)^{3}2E_{3}^{\prime}}f_{3}(p_{3}^{\prime})$
$\displaystyle\times|M|^{2}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{1}^{\prime}-p_{2}^{\prime}-p_{3}^{\prime}),$
(85)

where $f_{i}(p_{i}^{\prime})$ is a particle distribution function including
its degeneracy factor $D_{i}^{\prime}$. Keeping in mind that $|M|^{2}$ does
not depend on the scattering angle in the c.m. frame of $p_{1}+p_{2}$, the
rate may be written as

$\displaystyle\frac{\Delta N^{3\rightarrow 2}}{\Delta t\Delta
V}=\frac{1}{D_{1}^{\prime}D_{2}^{\prime}D_{3}^{\prime}}\frac{p_{1}^{c.m.}}{4\pi\sqrt{s}}|M|^{2}\int\frac{d^{3}p_{1}^{\prime}}{(2\pi)^{3}2E_{1}^{\prime}}f_{1}(p_{1}^{\prime})$
$\displaystyle\times\int\frac{d^{3}p_{2}^{\prime}}{(2\pi)^{3}2E_{2}^{\prime}}f_{2}(p_{2}^{\prime})\int\frac{d^{3}p_{3}^{\prime}}{(2\pi)^{3}2E_{3}^{\prime}}f_{3}(p_{3}^{\prime}).$
(86)

Expressing the distribution functions by Xu:2004mz

$\displaystyle f_{i}=(2\pi)^{3}\frac{\Delta N_{i}}{\Delta p_{i}^{3}\Delta V},$
(87)

the interaction rate turns to

$\displaystyle\frac{\Delta N^{3\rightarrow 2}}{\Delta t\Delta
V}=\frac{1}{D_{1}^{\prime}D_{2}^{\prime}D_{3}^{\prime}}\frac{p_{1}^{c.m.}}{4\pi\sqrt{s}}|M|^{2}\frac{\Delta
N_{1}^{\prime}\Delta N_{2}^{\prime}\Delta
N_{3}^{\prime}}{8E_{1}^{\prime}E_{2}^{\prime}E_{3}^{\prime}(\Delta V)^{3}}.$
(88)

Now substituting $|M|^{2}$ from Eq. (81) one obtains,

$\displaystyle\frac{\Delta N^{3\rightarrow 2}}{\Delta
t}=\frac{D_{1}D_{2}}{D_{1}^{\prime}D_{2}^{\prime}D_{3}^{\prime}}\frac{(p_{1}^{c.m.})^{2}}{\pi(PS)_{3}}\sigma_{2\rightarrow
3}$ $\displaystyle\times\frac{\Delta N_{1}^{\prime}\Delta N_{2}^{\prime}\Delta
N_{3}^{\prime}}{8E_{1}^{\prime}E_{2}^{\prime}E_{3}^{\prime}(\Delta V)^{2}},$
(89)

where $\Delta t$, $\Delta V$ and $\Delta N$ correspond to the simulation time
interval, volume of the grid cell, and particle number in the grid cell,
respectively. We note that Eq. (89) is equivalent to the expression in Ref.
Seifert:2017oyb .

The transition probability for the transition $NYK\rightarrow
N^{\prime}N^{{}^{\prime\prime}}$ in the volume $\Delta V$ during the time
interval $\Delta t$ is given by

$\displaystyle P_{NYK\rightarrow N^{\prime}N^{{}^{\prime\prime}}}$
$\displaystyle=$
$\displaystyle\frac{D_{N^{\prime}}D_{N^{{}^{\prime\prime}}}}{D_{N}D_{Y}D_{K}}\frac{(p_{N^{\prime}}^{c.m.})^{2}}{\pi(PS)_{3}}\sigma_{N^{\prime}N^{{}^{\prime\prime}}\rightarrow
NYK}$ (90) $\displaystyle\times$ $\displaystyle\frac{\Delta N_{N}\Delta
N_{Y}\Delta N_{K}}{8E_{N}E_{Y}E_{K}(\Delta V)^{2}}\Delta t,$

with the three-body phase-space integral

$\displaystyle(PS)_{3}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{N}}{(2\pi)^{3}2E_{N}}\int\frac{d^{3}p_{Y}}{(2\pi)^{3}2E_{Y}}\int\frac{d^{3}p_{K}}{(2\pi)^{3}2E_{K}}\hfill$
(91) $\displaystyle\times$
$\displaystyle(2\pi)^{4}\delta^{(4)}(p_{N^{\prime}}+p_{N^{\prime\prime}}-p_{N}-p_{Y}-p_{K}).\hfill\,\,\,~{}~{}~{}~{}~{}~{}$

The degeneracy factors need additional specification: for example,
($\Sigma^{0}~{}K^{+}~{}n$), ($\Sigma^{+}~{}K^{0}~{}n$),
($\Sigma^{0}~{}K^{0}~{}p$), or ($\Sigma^{-}~{}K^{+}~{}p$) can produce
($p~{}n$), ($p~{}\Delta^{0}$), ($n~{}\Delta^{+}$),
($\Delta^{+}~{}\Delta^{-}$), or ($\Delta^{++}~{}\Delta^{-}$), and the ratio of
the degeneracy factors in Eq. (90) is 1/4, 1/2, and 1/2 for the final states
of $(N~{}N)$, $(N~{}\Delta)$ and $(\Delta~{}\Delta)$, respectively, which are
distinguished by $p_{N^{\prime}}^{c.m.}$. As another example,
$(\Lambda~{}K^{+}~{}p)$ to $(p~{}p)$ has an additional degeneracy factor of
1/2 because two protons in the final state are indistinguishable. The mass of
the $\Delta$ baryon is randomly sampled from its spectral function, and if the
total mass of the final states is above $\sqrt{s}$, the transition does not
take place. This Monte Carlo sampling is effectively equivalent to the
integral of the $\Delta$ spectral function over invariant mass.

##  Appendix C cross section for $\bar{K}$ production in the medium

The phase space for the scattering of on-shell particles is given by

$\displaystyle\int\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}\frac{d^{3}p_{4}}{(2\pi)^{3}2E_{4}}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4}).~{}~{}~{}$
(92)

The phase space for $p_{3}$ can be changed into a more explicitly covariant
form

$\displaystyle\int\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}\rightarrow\int\frac{d^{4}p_{3}}{(2\pi)^{4}}2\pi\delta^{+}(p_{3}^{2}-m_{3}^{3})$
(93)

where $\delta^{+}(x-x_{0})$ is nonzero only for positive $x_{0}$. Comparing
the normalization condition of Eq. (28) and Eq. (93), one can find that the
spectral function $A(\omega,{\bf p}_{3})$ corresponds to
$2\pi\delta^{+}(p_{3}^{2}-m_{3}^{2})$ with the same normalization condition,

$\displaystyle 2\int_{0}^{\infty}\frac{d\omega}{2\pi}\omega
2\pi\delta^{+}(p_{3}^{2}-m_{3}^{2})=\int_{0}^{\infty}d\omega
2\omega\delta^{+}(\omega^{2}-E_{3}^{2})$
$\displaystyle=\int_{0}^{\infty}d\omega\frac{\omega}{E_{3}}\delta(\omega-
E_{3})=1,~{}~{}~{}~{}~{}$ (94)

where $E=\sqrt{m_{3}^{2}+p_{3}^{2}}$. Therefore, the covariant form of Eq.
(93) can be modified for off-shell particles into

$\displaystyle\int\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}\rightarrow\int\frac{d^{4}p_{3}}{(2\pi)^{4}}A(p_{3}^{0},{\bf
p}_{3})$ $\displaystyle=\int\frac{dp_{3}^{0}}{2\pi}A(p_{3}^{0},{\bf
p}_{3})\int\frac{d^{3}p_{3}}{(2\pi)^{3}}$
$\displaystyle=\int_{0}^{(\sqrt{s}-m_{4})^{2}}\frac{dm_{3}^{2}}{2\pi}A(m_{3}^{2})\int\frac{d^{3}p_{3}}{(2\pi)^{3}2p_{3}^{0}},$
(95)

where $m_{3}^{2}=(p_{3}^{0})^{2}-{\bf p}_{3}^{2}$ and $m_{4}$ is the mass of
the particle other than particle $3$ in the final state. The cross section
turns to

$\displaystyle\sigma(\sqrt{s};M_{3})\rightarrow\int_{0}^{(\sqrt{s}-m_{4})^{2}}\frac{dm_{3}^{2}}{2\pi}A(m_{3}^{2})~{}\sigma(\sqrt{s};m_{3}),$
(96)

where $M_{3}$ and $m_{3}$ are on-shell and off-shell masses, respectively.
Applying the same assumption of Eq. (43) to Eq. (96), the cross section for
off-shell particle production reads

$\displaystyle\int_{0}^{(\sqrt{s}-m_{4})^{2}}\frac{dm_{3}^{2}}{2\pi}A(m_{3}^{2})~{}\sigma(\sqrt{s}-m_{3}+M_{3}).$
(97)

The same spectral function can be applied to the decay of a $\phi$ or $K^{*}$
meson. Assuming that the transition amplitude does not depend on the mass of
the daughter particles, the decay width is calculated as

$\displaystyle\Gamma\sim\frac{1}{M}\int_{0}^{(M-m_{2})^{2}}\frac{dm_{1}^{2}}{2\pi}A(m_{1}^{2})\int\frac{d^{3}p_{1}}{(2\pi)^{3}2E_{1}}$
$\displaystyle\times\int\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}(2\pi)^{4}\delta^{(4)}(p-p_{1}-p_{2})$
$\displaystyle=\frac{1}{4\pi
M^{2}}\int_{0}^{(M-m_{2})^{2}}\frac{dm_{1}^{2}}{2\pi}A(m_{1}^{2})p_{1}^{c.m.}(m_{1}^{2}),$
(98)

where $M$ is the mass of the mother particle and $p_{1}^{c.m.}$ is the three-
momentum of the off-shell particle in the c.m. frame. Therefore $m_{1}$ \-
produced through the decay - has the mass differential distribution

$\displaystyle\frac{dP(m_{1}^{2})}{dm_{1}^{2}}=\frac{A(m_{1}^{2})p_{1}^{c.m.}(m_{1}^{2})}{8\pi^{2}M^{2}\Gamma}=A(m_{1}^{2})~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\times\frac{\sqrt{\\{M^{2}-(m_{1}+m_{2})^{2}\\}\\{M^{2}-(m_{1}-m_{2})^{3}\\}}}{(4\pi)^{2}M^{3}\Gamma}.$
(99)

The probability is proportional to $p_{1}^{c.m.}$ because the phase space
increases with $p_{1}^{c.m.}$.

##  Appendix D Time-evolution of the off-shell particle mass

As shown in Eq. (37) the energy of the off-shell particle is updated as

$\displaystyle\frac{dE}{dt}=\frac{1}{2E}\bigg{[}\partial_{t}{\rm
Re}\Sigma+\frac{M^{2}-M_{0}^{2}}{{\rm Im}\Sigma}\partial_{t}{\rm
Im}\Sigma\bigg{]},$ (100)

where $C$ is neglected for simplicity. Since the partial derivatives are
expressed in terms of the total derivatives as

$\displaystyle\frac{\partial{\rm Re}\Sigma}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{d{\rm Re}\Sigma}{dt}-{\bf v}\cdot\nabla{\rm Re}\Sigma,$
(101) $\displaystyle\frac{\partial{\rm Im}\Sigma}{\partial t}$
$\displaystyle=$ $\displaystyle\frac{d\Gamma}{dt}-{\bf v}\cdot\nabla{\rm
Im}\Sigma,$ (102)

Eq. (100) turns to

$\displaystyle\frac{dE}{dt}=\frac{1}{2E}\bigg{\\{}\frac{d{\rm
Re}\Sigma}{dt}+\frac{M^{2}-M_{0}^{2}}{{\rm Im}\Sigma}\frac{d{\rm
Im}\Sigma}{dt}\bigg{\\}}$ $\displaystyle-\frac{{\bf
v}}{2E}\cdot\bigg{\\{}\nabla{\rm Re}\Sigma+\frac{M^{2}-M_{0}^{2}}{{\rm
Im}\Sigma}\nabla{\rm Im}\Sigma\bigg{\\}}$
$\displaystyle=\frac{1}{2E}\bigg{\\{}\frac{d{\rm
Re}\Sigma}{dt}+\frac{M^{2}-M_{0}^{2}}{{\rm Im}\Sigma}\frac{d{\rm
Im}\Sigma}{dt}\bigg{\\}}+{\bf v}\cdot\frac{d{\bf p}}{dt},$ (103)

where Eq. (36) has been substituted in the second equation. Multiplying by
$2E$ on both sides,

$\displaystyle 2E\frac{dE}{dt}-2{\bf p}\cdot\frac{d{\bf p}}{dt}-\frac{d{\rm
Re}\Sigma}{dt}=\frac{M^{2}-M_{0}^{2}}{{\rm Im}\Sigma}\frac{d{\rm
Im}\Sigma}{dt},$ (104)

and using

$\displaystyle 2E\frac{dE}{dt}=\frac{dE^{2}}{dt},~{}~{}~{}~{}~{}2{\bf
p}\cdot\frac{d{\bf p}}{dt}=\frac{d{\bf p}^{2}}{dt},$ (105)

the left-hand-side of Eq. (104) is simplified to

$\displaystyle\frac{d(E^{2}-{\bf p}^{2}-{\rm
Re}\Sigma)}{dt}\equiv\frac{dM^{2}}{dt}=\frac{M^{2}-M_{0}^{2}}{{\rm
Im}\Sigma}\frac{d{\rm Im}\Sigma}{dt},$ (106)

where $M^{2}=E^{2}-{\bf p}^{2}-{\rm Re}\Sigma$, which is equivalent to Eq.
(39).

## References

  * (1) D. B. Kaplan and A. E. Nelson, Phys. Lett. B 175, 57 (1986); A. E. Nelson and D. B. Kaplan, ibid, 192, 193 (1987).
  * (2) D. Best et al., FOPI Collaboration, Nucl. Phys. A 625, 307 (1997).
  * (3) F. Laue et al., KaoS Collaboration, Phys. Rev. Lett. 82, 1640 (1999).
  * (4) M. Menzel et al. [KaoS], Phys. Lett. B 495, 26 (2000).
  * (5) C. Sturm et al., KaoS Collaboration, Phys. Rev. Lett. 86, 39 (2001).
  * (6) A. Förster et al., KaoS Collaboration, J. Phys. G 28, 2011 (2002).
  * (7) V. Zinyuk et al. [FOPI], Phys. Rev. C 90, 025210 (2014).
  * (8) G. Agakishiev et al. [HADES], Eur. Phys. J. A 47, 21 (2011).
  * (9) J. Adamczewski-Musch et al. [HADES], Phys. Lett. B 793, 457 (2019).
  * (10) J. Aichelin and C. M. Ko, Phys. Rev. Lett. 55, 2661 (1985).
  * (11) X. S. Fang, C. M. Ko, G. Q. Li, and Y. M. Zheng, Phys. Rev. C 49, R608 (1994); G. Q. Li, C. M. Ko, and X. S. Fang, Phys. Lett. B 329, 149 (1994).
  * (12) C. M. Ko and G. Q. Li, J. Phys. G 22, 1673 (1996); C. M. Ko, V. Koch, and G. Q. Li, Ann. Rev. Nucl. Part. Sci. 47, 505 (1997).
  * (13) G. Q. Li, C.-H. Lee, and G. E. Brown, Nucl. Phys. A 625, 372 (1997).
  * (14) C. M. Ko, J. Phys. G 27, 327 (2001).
  * (15) S. Pal, C. M. Ko, and Z.-W. Lin, Phys. Rev. C 64, 042201 (2001).
  * (16) W. Cassing, E.L. Bratkovskaya, U. Mosel, S. Teis, and A. Sibirtsev, Nucl. Phys. A 614, 415 (1997).
  * (17) E. L. Bratkovskaya, W. Cassing, and U. Mosel, Nucl. Phys. A 622, 593 (1997).
  * (18) W. Cassing and E. L. Bratkovskaya, Phys. Rep. 308, 65 (1999).
  * (19) W. Cassing, L. Tolós, E. L. Bratkovskaya, and A.Ramos, Nucl. Phys. A 727, 59 (2003).
  * (20) J. Schaffner-Bielich, V. Koch, and M. Effenberger, Nucl. Phys. A 669, 153 (2000).
  * (21) C. Hartnack, H. Oeschler, and J. Aichelin, Phys. Rev. Lett. 90, 102302 (2003); C. Hartnack and J. Aichelin, J. Phys. G 28, 1649 (2002).
  * (22) C. Fuchs, Amand Faessler, E. Zabrodin, and Yu-Ming Zheng, Phys. Rev. Lett. 86, 1974 (2001); C. Fuchs, Prog. Part. Nucl. Phys. 56, 1 (2006).
  * (23) A. Mishra, E. L. Bratkovskaya, J. Schaffner-Bielich, S. Schramm and H. Stoecker, Phys. Rev. C 70, 044904 (2004).
  * (24) C. Hartnack, H. Oeschler, Y. Leifels, E. L. Bratkovskaya and J. Aichelin, Phys. Rept. 510, 119 (2012).
  * (25) E. E. Kolomeitsev, C. Hartnack, H. W. Barz, M. Bleicher, E. Bratkovskaya, W. Cassing, L. W. Chen, P. Danielewicz, C. Fuch, T. Gaitanos, C. M. Ko, A. Larionov, M. Reiter, G. Wolf and J. Aichelin, J. Phys. G 31, S741 (2005).
  * (26) G. Q. Li, C. M. Ko, and B. A. Li, Phys. Rev. Lett. 74, 235 (1995); Nucl. Phys. A 594, 460 (1995); G. Q. Li and C. M. Ko, Phys. Rev. C 54, 1897 (1996); G. Q. Li, C. M. Ko, and G. E. Brown, Phys. Lett. B 381, 17 (1996).
  * (27) G. E. Brown and M. Rho, Phys. Rev. Lett. 66, 2720 (1991).
  * (28) E. Friedman, A. Gal, and C.J. Batty, Nucl. Phys. A 579, 518 (1994).
  * (29) A. Gal, Nucl. Phys. A 691, 268 (2001).
  * (30) E. Friedman and A. Gal, Phys. Rept. 452, 89 (2007).
  * (31) M. Lutz, Phys. Lett. B 426, 12 (1998).
  * (32) T. Waas, N. Kaiser and W. Weise, Phys. Lett. B 365, 12 (1996).
  * (33) T. Waas, N. Kaiser and W. Weise, Phys. Lett. B 379, 34 (1996).
  * (34) M. Lutz and E. E. Kolomeitsev, Nucl. Phys. A 700, 193 (2002).
  * (35) M. Lutz and C. L. Korpa, Nucl. Phys. A 700, 309 (2002).
  * (36) A. Ramos and E. Oset, Nucl. Phys. A 671, 481 (2000).
  * (37) A. Mishra, K. Balazs, D. Zschiesche, S. Schramm, H. Stoecker and W. Greiner, Phys. Rev. C 69, 024903 (2004).
  * (38) L. Tolos, A. Ramos, A. Polls and T. T. S. Kuo, Nucl. Phys. A 690, 547 (2001).
  * (39) L. Tolos, A. Ramos and A. Polls, Phys. Rev. C 65, 054907 (2002).
  * (40) V. Koch, Phys. Lett. B 337, 7 (1994).
  * (41) L. Tolos and L. Fabbietti, Prog. Part. Nucl. Phys. 112, 103770 (2020).
  * (42) W. Cassing and S. Juchem, Nucl. Phys. A 672, 417 (2000).
  * (43) W. Cassing and S. Juchem, Nucl. Phys. A 665, 377 (2000).
  * (44) W. Cassing, Eur. Phys. J. ST 168, 3 (2009).
  * (45) W. Ehehalt and W. Cassing, Nucl. Phys. A 602, 449 (1996).
  * (46) W. Cassing and E. L. Bratkovskaya, Phys. Rev. C 78, 034919 (2008).
  * (47) W. Cassing and E. L. Bratkovskaya, Nucl. Phys. A 831, 215 (2009)
  * (48) E. L. Bratkovskaya, W. Cassing, V. P. Konchakovski and O. Linnyk, Nucl. Phys. A 856, 162 (2011)
  * (49) O. Linnyk, E. L. Bratkovskaya and W. Cassing, Prog. Part. Nucl. Phys. 87, 50 (2016)
  * (50) P. Moreau, O. Soloveva, L. Oliva, T. Song, W. Cassing and E. Bratkovskaya, Phys. Rev. C 100, 014911 (2019).
  * (51) A. Peshier and W. Cassing, Phys. Rev. Lett. 94, 172301 (2005).
  * (52) W. Cassing, Nucl. Phys. A 795, 70 (2007).
  * (53) W. Cassing, Nucl. Phys. A 791, 365 (2007).
  * (54) L. Tolos, A. Ramos and E. Oset, Phys. Rev. C 74, 015203 (2006).
  * (55) L. Tolos, D. Cabrera and A. Ramos, Phys. Rev. C 78, 045205 (2008).
  * (56) D. Cabrera, L. Tolos, J. Aichelin and E. Bratkovskaya, Phys. Rev. C 90, 055207 (2014).
  * (57) S. A. Bass, M. Belkacem, M. Bleicher, et al., Prog. Part. Nucl. Phys. 41, 255 (1998).
  * (58) V. Steinberg, J. Staudenmaier, D. Oliinychenko, F. Li, Ö. Erkiner and H. Elfner, Phys. Rev. C 99, 064908 (2019).
  * (59) E. Oset and A. Ramos, Nucl. Phys. A 679, 616 (2001).
  * (60) E. Oset and A. Ramos, Nucl. Phys. A 635, 99 (1998).
  * (61) J. A. Oller and U. G. Meissner, Phys. Lett. B 500, 263 (2001).
  * (62) D. Jido, E. Oset and A. Ramos, Phys. Rev. C 66, 055203 (2002).
  * (63) H. Schopper (Ed.), Landolt-Börnstein, New Series, Vol. I/12, Springer-Verlag, 1988.
  * (64) B. Nilsson-Almqvist and E. Stenlund, Comp. Phys. Comm. 43, 387 (1987); B. Andersson, G. Gustafson, and H. Pi, Z. Phys. C 57, 485 (1993).
  * (65) B. Nilsson-Almqvist and E. Stenlund, Comp. Phys. Comm. 43, 387 (1987); B. Andersson, G. Gustafson, and H. Pi, Z. Phys. C 57, 485 (1993).
  * (66) T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605, 026 (2006).
  * (67) V. Kireyeu, I. Grishmanovskii, V. Kolesnikov, V. Voronyuk and E. Bratkovskaya, Eur. Phys. J. A 56, 223 (2020).
  * (68) J. Schwinger, Phys. Rev. 83, 664 (1951).
  * (69) W. Cassing, A. Palmese, P. Moreau and E.L. Bratkovskaya, Phys. Rev. C 93, 014902 (2016).
  * (70) A. Palmese, W. Cassing, E. Seifert, T. Steinert, P. Moreau, and E.L. Bratkovskaya Phys. Rev. C 94, 044912 (2016).
  * (71) S. Borsanyi et al., Phys. Lett B 730, 99 (2014); Phys. Rev. D 92, 014505 (2015).
  * (72) V. Ozvenchuk, O. Linnyk, M. I. Gorenstein, E. L. Bratkovskaya and W. Cassing, Phys. Rev. C 87, 024901 (2013).
  * (73) W. Cassing, Nucl. Phys. A 700, 618 (2002).
  * (74) E. Seifert and W. Cassing, Phys. Rev. C 97, 024913 (2018).
  * (75) E. L. Bratkovskaya and W. Cassing, Nucl. Phys. A 807, 214 (2008).
  * (76) A. Ilner, D. Cabrera, C. Markert and E. Bratkovskaya, Phys. Rev. C 95, 014903 (2017); A. Ilner, J. Blair, D. Cabrera, C. Markert and E. Bratkovskaya, Phys. Rev. C 99, 024914 (2019).
  * (77) H. W. Barz and H. Iwe, Phys. Lett. B 153, 217 (1985).
  * (78) S. Protopopescu, M. Alston-Garnjost, A. Barbaro-Galtieri, S. M. Flatte, J. Friedman, T. Lasinski, G. Lynch, M. Rabin and F. Solmitz, Phys. Rev. D 7, 1279 (1973).
  * (79) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, 030001 (2018).
  * (80) E. L. Bratkovskaya, W. Cassing, C. Greiner, M. Effenberger, U. Mosel and A. Sibirtsev, Nucl. Phys. A 675, 661 (2000).
  * (81) J. Cleymans and H. Satz, Z. Phys. C 57, 135 (1993).
  * (82) N. Bastid et al. [FOPI], Phys. Rev. C 76, 024906 (2007).
  * (83) J. Adamczewski-Musch et al. [HADES], Phys. Lett. B 778, 403 (2018).
  * (84) A. Forster et al. [KaoS], Phys. Rev. C 75, 024906 (2007).
  * (85) A. Feijoo, V. Magas and A. Ramos, Phys. Rev. C 99, 035211 (2019).
  * (86) H. Joos, Phys. Lett. B 24, 103 (1967).
  * (87) G. Agakishiev et al. [HADES], Phys. Rev. C 80, 025209 (2009).
  * (88) G. Agakishiev, A. Balanda, B. Bannier, R. Bassini, D. Belver, A. V. Belyaev, A. Blanco, M. Bohmer, J. L. Boyard and P. Cabanelas, et al. Phys. Rev. C 82, 044907 (2010).
  * (89) J. L. Ritman et al. [FOPI], Z. Phys. A 352, 355 (1995).
  * (90) P. Crochet et al. [FOPI], Phys. Lett. B 486, 6 (2000).
  * (91) M. Ploskon [KaoS], Nucl. Phys. A 749, 170 (2005).
  * (92) F. Uhlig et al. [KaoS], Phys. Rev. Lett. 95, 012301 (2005).
  * (93) X. Lopez et al. [FOPI], Phys. Rev. C 75, 011901 (2007).
  * (94) B. A. Li, L. W. Chen and C. M. Ko, Phys. Rept. 464, 113 (2008).
  * (95) T.H.R. Skyrme, Nucl. Phys. 9, 615 (1958–1959).
  * (96) G. Ferini, M. Colonna, T. Gaitanos, M. Di Toro, H. Wolter, Phys. Rev. Lett. 97, 202301 (2006); M. Di Toro, et al., Nucl. Phys. A 782, 267 (2007).
  * (97) T. Song and C. M. Ko, Phys. Rev. C 91, 014901 (2015).
  * (98) A. Sibirtsev, Phys. Lett. B 359, 29 (1995).
  * (99) A. Sibirtsev, W. Cassing and C. Ko, Z. Phys. A 358, 101 (1997).
  * (100) S. Huang, K. Tsushima, A. Faessler, E. Lehmann and R. K. Puri, [arXiv:nucl-th/9412011 [nucl-th]].
  * (101) K. Tsushima, S. Huang and A. Faessler, Phys. Lett. B 337, 245 (1994).
  * (102) K. Tsushima, S. Huang and A. Faessler, J. Phys. G 21, 33 (1995).
  * (103) S. Efremov and E. Parev, Phys. Atom. Nucl. 57, 532 (1994).
  * (104) T. Song, K. C. Han and C. M. Ko, Phys. Rev. C 85, 054905 (2012).
  * (105) Z. Xu and C. Greiner, Phys. Rev. C 71, 064901 (2005).