# Characteristics of Power Loss in SMC a Key for Desining the Best Values of
Technological Parameters

Barbara Ślusarek<EMAIL_ADDRESS>Tele and Radio Research
Institute, ul. Ratuszowa 11, 03-450 Warszawa, Poland  Bartosz Jankowski
<EMAIL_ADDRESS>Tele and Radio Research Institute, ul. Ratuszowa
11, 03-450 Warszawa, Poland  Krzysztof Sokalski<EMAIL_ADDRESS>Institute of Computer Science, Czȩstochowa University of Technology, Al. Armii
Krajowej 17, 42-200 Czȩstochowa, Poland  Jan Szczygłowski
<EMAIL_ADDRESS>Institute of Power Engineering, Czȩstochowa University
of Technology, Al. Armii Krajowej 17, 42-200 Czȩstochowa, Poland
Eletrotechnical Institute, Ul. Pozaryskiego 28, 04-703 Warszawa, Poland

###### Abstract

Optimization of power loss in soft magnetic components basis on the choice of
the best technological parameters values. Therefore, the power losses have
been measured in Somaloy 500 samples for a wide range of frequency and
magnetic induction. These samples have been prepared under a wide range of the
hardening temperatures and pressures. The power loss characteristics have been
derived by assuming that investigated samples obeyed the scaling law.
Agreement obtained between experimental data and the scaling theory has
confirmed this assumption. Moreover, the experimental data of the given sample
have been collapsed to a single curve which represented measurements for all
values of frequency an magnetic induction pick. Therefore, the scaling
transforms the losses characteristics from the two dimensional surfaces to the
one dimensional curves. The samples were produced according two methods: for
different pressures with constant temperature and at different temperatures
with constant pressure. In both cases the power losses decrease with
increasing pressure and with increasing temperature. These trends in
decreasing losses stopped for certain critical values of pressure and
temperature, respectively. Above these values the power losses increase
suddenly. Therefore, the mentioned above the critical pressure and the
critical temperature are sought after solutions for optimal values. In order
to reduce the parameters values set the limit curve in the pressure-
temperature plane has been derive. This curve constitutes a separation curve
between the parameters values corresponding to high and low losses.

###### pacs:

44.05.+e,75.90.+w

##  I Introduction

The expanding interest in the introduction of Fe-based composite materials,
also called the Soft Magnetic Composite (SMC), in electrical devices, such as
electrical motors, is in obvious ways connected with their properties. The
principle feature of these materials is that iron particles are insulated by a
thin organic or inorganic coating. SMC offers several advantages over
laminated steel sheet, e.g. isotropic magnetic properties, the opportunity to
tailor their physical properties, very low eddy current loss, relatively high
resistivity, high magnetic permeability (bib:Shokrollahi, ). Fe-based
composite materials are used in the manufacture of electrical devices owing to
their physical properties. The most common use of these materials is in
various electrical motors, e.g. a high speed permanent magnet motor, a
synchronous electric motor, an axial-flux permanent magnet synchronous
machine, a claw pole permanent magnet motor (bib:Wojciechowski, ),
(bib:Hamlera, ), (bib:Stefano, ), (bib:Guo, ). Furthermore, soft magnetic
composites are used in an electromagnetic actuator (bib:Jack, ) and in
electromagnetic shielding (bib:Kim, ). The manufacturing process of magnetic
composites can be divided into a few stages. Every stage is characterized by a
set of technological parameters, which may ultimately affect the outcome
parameters of magnetic composite materials. Three of the most important
processing considerations of compression moulding techniques are pressure,
temperature and time of curing. The value of compaction pressure has a strong
influence on the initial, ’green’, and final compacted densities, which is
reflected in some of the physical properties of magnetic composites. On the
other hand, changes in temperature or time of curing lead to variations in the
dielectric parameters of insulation layer. It is all the more relevant because
inappropriate selection of these parameters can lead to an increase in total
energy losses. Development of electrical devices, such as electric motors or
actuators, is closely connected with development of magnetic materials,
obviously. However, the mere application of modern magnetic materials is no
guarantee of success. This implies the need to look for new ways to design the
best value of technological parameters. In the last decade some methods for
determining these values have been elaborated by a few laboratories
(bib:Dobrz, ),(bib:Dobrz2, ), (bib:Pang, ),(bib:Gilbert, ),(bib:Lefe,
),(bib:Lemi, ),(bib:Gelin, ). Although none of them have formulated their
method in the algorithmic form, that mainly one can distinguish the following
common steps in the applied procedures:

  * •

The series of samples are prepared of different powders compositions and
processing methods.

  * •

The core loss versus magnetic induction and frequency is measured. Usually the
results of these measurements are presented as energy loss characteristics
versus $B_{m}$ for the constant values of $f$.

  * •

The particles content and the values of technological parameters corresponding
to the lowest energy loss are reported for desired values of $B_{m}$ and $f$.

The aim of this paper is to work out systematic procedures for designing the
best values of technological parameters and optimizing the power losses.
Realization of this venture has been done on the basis of the scaling method
applied to modeling of the power loss in the soft magnetic materials
(bib:Sokal, )\- (bib:Sokal4, ). Additionally, two advantages have been
achieved. The first one consists of algorithmic structure of the applied
procedure which will enable creation of a computer implementation. The second
one is the case with which SMC’s properties can be tailored. As an example the
technological parameters for production of a material with less excess losses
has been designed and such a sample has been prepared in accordance with that.

##  II SOFT MAGNETIC COMPOSITE

Soft magnetic composite (SMC) consist generally of ferromagnetic particles
distributed in a matrix of binding agent. Metal powder is produced by
atomization techniques; a stream of molten metal is atomized by inert gas or
water under high pressure, leading to formation of fine liquid metal droplets
which then solidify resulting in fine powder particles. The nature of impact
media determines shape of powder particles. A jet impacting on water produces
irregular particles (Fig.1). In turn, use of the gas as the media produces
spherical particles. Grains of iron powder for soft magnetic composites are
covered chemically by thin inorganic insulation layer-phosphate glass with
thickness of a few nanometers.

Figure 1: Irregular particles of soft magnetic composite

Magnetic composites are prepared by compression moulding; a die is filled with
magnetic powder and it is compacted. Then the green compact is cured in a
furnace at melting temperature of binding agent. In the case of compression
moulding techniques an air atmosphere is used during the curing, frequently. A
commercially available iron powder, produced by water atomization, with a
binding agent has been used in the experiment as the soft magnetic composite.
It is called Somaloy 500 and is produced by H gan s AB Company.

##  III Power Loss’ Measurements

The samples investigated in this study were produced from commercially
available pure iron powders produced by water atomization techniques, with a
special surface coating on each and every particle. The powder used in the
investigation was Somaloy 500 produced by H gan s AB Company. The samples used
were ring-shaped with a square cross-section. Samples had the following
dimensions, external diameter of 55 mm, internal diameter of 45 mm and
thickness of 5 mm. The power losses Ptot were obtained from measurements of
the AC hysteresis cycle according to the IEC Standards 60404-6 using
hysteresisgraph AMH-20K-HS produced by Laboratorio Elettrofisico Walker LDJ
Scientific. Total power losses Ptot were measured at maximum flux density
$B_{m}$ = 0,1…1,3 T, over a frequency range from 10 to 5000 Hz. During
measurements of the total power losses $P_{tot}$ the shape factor of the
secondary voltage was equal to $1.111\pm 1.5\%$. Measurement error of the
total energy losses was equal to $\pm 3\%$. Measurements of power losses were
conducted by recording individual points and the integration area of the
hysteresis loop. Optimization of power loss has been performed with respect to
the compacting pressure $p$ and temperature $T$, while the other technological
parameters (powder compositions, volume fraction) were constant. This
procedure was based on the set of $16$ samples which has been produced
according to the measurement matrix of the technological parameters. This
matrix is readable from the two first columns of the TABLE 1. The particle
morphology observations were carried out using a scanning electron microscope
JSM-7600F produced by Joel.

##  IV Power Loss’ Modelling and Technological Parameters Values’ Desining

Soft Magnetic Composites are assumed to be complex system where the function
of energy losses obeys the scaling law. This assumption leads to the
expression of total energy loss in form of general homogeneous function
bib:Sokal :

$P_{tot}(f,B_{m})=B_{m}^{\beta}F\left(\frac{f}{B_{m}^{\alpha}}\right),$ (1)

where $F(\cdot)$ is an arbitrary function, $\alpha$ and $\beta$ are scaling
exponents. This function depends on features of phenomena to be describe. In
order to much $F(\cdot)$ with the eddgy current and hysteresis losses we
choose for $F(\cdot)$ the power series as a rough description of the energy
losses $P_{tot}$ (bib:Sokal, ):

$P_{tot}=B_{m}^{\beta}\,(\Gamma_{1}\,\frac{f}{B_{m}^{\alpha}}+\Gamma_{2}\,\left(\frac{f}{B_{m}^{\alpha}}\right)^{2}+\Gamma_{3}\,\left(\frac{f}{B_{m}^{\alpha}}\right)^{3}+\Gamma_{4}\,\left(\frac{f}{B_{m}^{\alpha}}\right)^{4}),$
(2)

where $f$ \- frequency, $B_{m}$ \- amplitude of magnetic field’s induction.
Values of $\alpha$, $\beta$ and amplitudes $\Gamma_{n}$ have been estimated
for different values of pressure and temperature. It is easy to recognize by
the powers of $f$ that the two first terms in (2) correspond to hysteresis
losses $P_{h}$ and eddy current losses $P_{c}$, whereas higher terms
correspond to excess losses $P_{ex}$, respectively:

$P_{tot}=P_{h}+P_{c}+P_{ex}.$ (3)

For more details concerning this separation formula and scaling see to
(bib:Sokal4, ).

###  IV.1 Energy losses’ data of SMC and the scaling law

The parameters’ values of (2) has been estimated for each sample’s measurement
data by minimization of $\chi^{2}$

$\chi^{2}=\Sigma_{i=1}^{n}\left(\frac{P_{tot\,i}-P_{exp\,i}}{B_{m\,i}^{\beta}}\right)^{2}$
(4)

using the simplex method of Nelder and Mead (bib:Nelder, ).

Table 1:  Somaloy 500. Values of scaling exponents and coefficients of (2) v.s. hardening temperature and pressing pressure  Temp.[${}^{o}{C}$] | Pressure [MPa] | $\alpha$ | $\beta$ | $\Gamma_{1}[\frac{m^{2}}{s^{2}}T^{\alpha-\beta}]$ | $\Gamma_{2}[\frac{m^{2}}{s}T^{2\alpha-\beta}]$ | $\Gamma_{3}[m^{2}T^{3\alpha-\beta}]$ | $\Gamma_{4}[{m^{2}}{s}\,T^{4\alpha-\beta}]$
---|---|---|---|---|---|---|---
500 | 500 | -1,312223 | -0,01052897 | 0,170825 | 3,60591E-05 | 1,95309E-08 | -2,25448E-12
500 | 600 | -1,383117 | -0,1249742 | 0,1525481 | 3,32813E-05 | 9,25383E-09 | -1,17698E-12
500 | 700 | -1,73502 | -0,5173269 | 0,1558437 | 2,39335E-05 | 2,30877E-09 | -8,07475E-14
500 | 900 | -1,394503 | -0,082346 | 0,1012103 | 6,06463E-05 | -8,03122E-09 | 7,87711E-13
400 | 800 | -1,47298 | -0,028389 | 0,1831606 | 1,34725E-05 | 3,68924E-09 | -1,18518E-13
450 | 800 | -1,596151 | -0,1231805 | 0,1493641 | 2,481657E-05 | -1,217578E-09 | 6,120261E-14
550 | 800 | -2,034396 | -1,325493 | 0,1055319 | 1,40738E-04 | -1,06631E-08 | 4,54085E-13
600 | 800 | -1,608202 | -0,2323374 | 1,219838 | 0,000894116 | -5,30152E-08 | -1,66434E-11
438 | 764 | -2,059196 | -1,014382 | 0,1395258 | 7,941049E-05 | -1,446713E-08 | 1,482247E-12
469 | 764 | -1,402494 | -0,0137938 | 0,1340093 | 5,373893E-05 | -6,098702E-09 | 6,883649E-13
531 | 764 | -3,09792 | -2,318892 | 0,1293812 | 6,054176E-05 | -3,967688E-09 | 2,019266E-13
480 | 780 | -1,634466 | -0,3416979 | 0,1242102 | 6,852739E-05 | -9,746355E-09 | 1,001089E-12

For results see TABLE 1. Using these values and (2) and the estimators of
energy losses in considered samples has been obtained. Each continuous line
represents the corresponding model predictions. Comparisons of estimations and
the measurement results were depicted in Fig.3 and Fig. 2 for $T=constant$ and
in Fig.4 and Fig. 5 for p= constant. Basing on these results one can derive
the following conclusions:

  * •

The energy losses’ experimental data obey the scaling law (1) which means that
$P_{tot}/B_{m}^{\beta}$ depends only on the one effective variable
$\tilde{f}=\frac{f}{B_{m}^{\alpha}}$. This property enables to present all
measurement results for $P_{tot}$ v.s. $f$ and $B_{m}$ on a one-dimensional
curve instead of the two-dimensional surface, see Fig.2-Fig.5.

  * •

The power loss characteristics $T=constant$ for different compacting pressures
do not cross each other except the origin point
$(\tilde{P}_{tot}=0,\frac{f}{B_{m}^{\alpha}}=0)$ , and similarly the power
loss characteristics $p=constant$ for different compacting temperatures do not
cross each other except the mentioned origin, see Fig.2-Fig.5.

  * •

Let
$\tilde{P}_{tot,T_{j},p_{k}}(\tilde{f})=P_{tot,T_{j},p_{k}}(\tilde{f})/B_{m}^{\beta}$
be scaled total energy loss in sample material which has been compacted in
$T_{j}$ temperature and under $p_{k}$ pressure. Let $T_{1}=T_{2}$, then
$\exists p_{c}(T)$ such that $\forall\tilde{f}\geq 0$ and $\forall
p_{c}(T_{1})\geq p_{1}>p_{2}$.the following relation holds:
$\tilde{P}_{tot,T_{1},p_{1}}(\tilde{f})<\tilde{P}_{tot,T_{1},p_{2}}(\tilde{f})$,
see Fig.2-Fig.3.

  * •

Let
$\tilde{P}_{tot,T_{j},p_{k}}(\tilde{f})=P_{tot,T_{j},p_{k}}(\tilde{f})/B_{m}^{\beta}$
be scaled total energy loss in sample material which has been compacted in
$T_{j}$ temperature and under $p_{k}$ pressure. Let $p_{1}=p_{2}$. Then
$\exists T_{c}(p)$ such that $\forall\tilde{f}\geq 0$ and $\forall
T_{c}(p)\geq T_{1}>T_{2}$ the following relation holds:
$\tilde{P}_{tot,T_{1},p_{1}}(\tilde{f})<\tilde{P}_{tot,T_{2},p_{1}}(\tilde{f})$,
see Fig.4-Fig.5.

The data presented here of measured properties and the mathematical model (2)
constitute algorithmic structure. Moreover, they enable one to optimize the
power losses globally for every values of independent variables $f$ and
$B_{m}$, therefore for the optimization process their values need not to be
assigned. Graph presenting working algorithm is presented in Fig.6 (bib:Jank,
).

Figure 2: Somaloy 500. $P_{tot}/B_{m}^{\beta}$ vs. $f/B_{m}^{\alpha}$ for
$T=450^{o}C$, $p=500,600,700$ and $900MPa.$. The continuous lines correspond
to (2). Figure 3: Somaloy 500. $P_{tot}/B_{m}^{\beta}$ vs. $f/B_{m}^{\alpha}$
for $T=500^{o}C$, $p=500,600,700,800$ and $900MPa$. The continuous lines
correspond to (2). Figure 4: Somaloy 500. $P_{tot}/B_{m}^{\beta}$ vs.
$f/B_{m}^{\alpha}$ for $p=800MPa$, $T=400,450,500$ and $550^{o}C$. The
continuous lines correspond to (2). Figure 5: Somaloy 500.
$P_{tot}/B_{m}^{\beta}$ vs. $f/B_{m}^{\alpha}$ for $p=800MPa$,
$T=400,450,500,550$ and $600^{o}C.$ The continuous lines correspond to (2).
Figure 6: Algorithm for designing the best values of technological parameters.

##  V Desining excess power loss-less material

Based on the results of (bib:Sokal4, ) and on the subsection IV of this paper
it is possible to design and produce a SMC’s material which does not exhibit
the excess contribution to the total power loss, e.g.: $P_{tot}=P_{h}+P_{c}$.
Let us call the values of technical parameters leading to such a material
technical critical parameters: $T_{tc}$ and $p_{tc}$. Let us truncate formula
(2) above the second order term. The obtained form described materials with
less excess power loss:

$P_{tot}=B_{m}^{\beta}\,\left(\Gamma_{1}\,\frac{f}{B_{m}^{\alpha}}+\Gamma_{2}\,\left(\frac{f}{B_{m}^{\alpha}}\right)^{2}\right).$
(5)

Comparing (2) and (5) one can conclude that materials with less excess power
loss are characterized by $\Gamma_{3}=0$ and $\Gamma_{4}=0$. Therefore, in
order to satisfy the formulated task one should produce a sample under the
technological conditions leading to these values of $\Gamma_{3}$ and
$\Gamma_{4}$, e.g. $T_{tc}$ and $p_{tc}$. How to do it? Having at least two
samples which have been produced in different points $(T,p_{1}),\,(T,p_{2})$
or $(T_{1},p),\,(T_{2},p)$ one can approximate the required values of $T_{tc}$
and $p_{tc}$ by inter or extrapolation. Fig.7 \- Fig.9 presents three examples
of such calculations. The best accuracy has been achieved for the example
presented in Fig.9. Note that both curves
$\Gamma_{3}=\Gamma_{3}(T=constant,p)$ and
$\Gamma_{4}=\Gamma_{4}(T=constant,p)$ get $0$ for the same values of $p$. By
this manner the technological critical parameters have been determined for
three materials with less excess power loss (see TABLE 2). These points are
sufficient to draw a curve in the $(T,p)$ plane (see Fig.10) which is
approximation representation of the whole set of technological parameters
values leading to materials with less excess power loss. Note that this curve
divides $T,p$ space into two parts. The upper part consists of better values
of $T$ and $p$ which lead to smaller values of the power loss in comparison to
the lower part. However, this space is bounded from above, both from $p$ and
$T$ sides (see Fig.3-Fig.5). This revelation helps to limit $T,p$ space to a
small subset whose boundaries are determined by the set of the best values of
technological parameters.

Figure 7: Somaloy 500, Interpolation of the technical critical point $T_{tc}=438^{o}C,p_{tc}=800MPa$. Figure 8: Somaloy 500. Extrapolation of the technical critical point $T_{tc}=495^{o}C,p_{tc}=764MPa$ Figure 9: Somaloy 500. Interpolation of the technical critical point $T_{tc}=500^{o}C,p_{tc}=728MPa$ Table 2:  Values of hardening temperature and pressing pressure corresponding to excess power loss less materials  $T_{tc}[{}^{o}{C}]$ | 438 | 480 | 495 | 500
---|---|---|---|---
$p_{tc}[MPa]$ | 800 | 780 | 764 | 728

Figure 10: Somaloy 500. Technical critical curve in the $T,p$ plane.
Markers correspond to the points determined in this paper.

##  VI Conclusions

Modern engineering of SMM designs materials which posses required physical and
technological properties. Particularly, technology of SMC enables one to
optimize both, the power loss and the magnetic properties. The results
presented here concern only the power losses and they are part of the Complex
researching programme (No. N N507 249940) which includes investigations and
optimization also that ones. However, in the case of this paper’s topic the
following results have been achieved. Power losses have been measured in
Somaloy 500 samples for a wide range of frequency and magnetic induction.
These samples have been prepared under a wide range of the hardening
temperatures and pressures. The power loss characteristics have been derived
by assuming that investigated samples obeyed the scaling law. Agreement
obtained between experimental data and the scaling theory has confirmed that
the assumptions of this paper are true. Moreover, the experimental data of the
given sample have been collapsed to a single curve which represented
measurements for all values of $f$ and $B_{m}$. Therefore, the the scaling
transforms the losses characteristics from the two dimensional surfaces to the
one dimensional curves. The samples were produced according two methods: for
different pressures with $T=constant$ and at different temperatures with
$p=constant$. In both cases the power losses decrease with increasing pressure
and with increasing temperature. These trends in decreasing losses stopped for
certain critical values $p_{c}(T)$ and $T_{c}(p)$, respectively. Above these
values the power losses increase suddenly. Therefore, $p_{c}(T)$ and
$T_{c}(p)$ are sought after solutions for optimal values of pressure and
temperature. The formula of presented power losses consists of four terms.
According to (bib:Sokal4, ) the two first terms in (2) correspond to $P_{h}$
and $P_{c}$. However, the sum of the third and fourth terms describes
$P_{ex}$. By the appropriate choice of $p$ and $T$ one can reduce $P_{ex}$ to
zero. In this way it is possible to design materials with less excess losses.
Three such points have been determined experimentally and a curve was drawn
representing all Somaloy 500 materials with less excess losses in the $(p,T)$
plane. The obtained curve divides $(p,T)$ plane into subsets of sample high
and low losses.

Finally we hope to achieve compact algorithm enabling the design of both
magnetic and power loss properties.

## References

  * (1) H. Shokrollahi, K. Janghorban, Soft magnetic composite materials (SMCs), Journal Of Materials Processing Technology, 189, 1-12(2007)
  * (2) R. M. Wojciechowski, C. Jedryczka, P. Lukaszewicz, D. Kapelski, Analysis of high speed permanent magnet motor with powder core material, COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 31(2012)
  * (3) A. Hamlera, V. Goricana, B. Sustarsicb, A. Sircc, The use of soft magnetic composite materials in synchronous electric motor, Journal of Magnetism and Magnetic Materials, 304,816 819(2006)
  * (4) R. Di Stefano, F. Marignetti, Electromagnetic Analysis of Axial-Flux Permanent Magnet Synchronous Machines With Fractional Windings With Experimental Validation, IEEE Transactions on Industrial Electronics, 59, 2573-2582(2012)
  * (5) Y. Guo, J. Zhu, D. G. Dorrell, Design and Analysis of a Claw Pole Permanent Magnet Motor With Molded Soft Magnetic Composite Core, IEEE Transactions on Industrial Electronics 45,4582-4585(2009)
  * (6) A. Jack, B. Mecrow, G. Atkinson, Current Developments in Soft Magnetic Composites for Electromagnetic Actuators, Sensor Letters 7,369-373(2009)
  * (7) S. W. Kim, Y.W. Yoon, S.J. Lee, G.Y. Kim, Y.B. Kim, Yun Yeo Chun, K.S. Lee, Electromagnetic shielding properties of soft magnetic powder polymer composite films for the application to suppress noise in the radio frequency range, Journal of Magnetism and Magnetic Materials 316,472 474(2007)
  * (8) P. Gramatyka, R. Nowosielski, P. Sakiewicz, T. Raszka, Soft magnetic composite based on nanocrystalline $Fe_{73.5}Cu_{1}Nb_{3}Si_{13.5}B_{9}$ and $Fe$ powders, Journal of Achievements in Materials and Manufacturing Engineering, 15, 28-31(2006).
  * (9) P. Gramatyka, A. Kolano-Burian , R. Kolano , M. Polak Nanocrystalline iron based powder cores for high frequency applications, Journal of Achievements in Materials and Manufacturing Engineering, 18, 99-102(2006).
  * (10) Y. X. Pang, S. N. B. Hodgson, J. Koniarek, B. Weglinski, The influence of the dielectric on the properties of dielectromagnetic soft magnetic composities. Investigations with silica and silica hybrid sol-gel derived model dielectric, J. Mag. Mag. Mater. 310, 83-91 (2007).
  * (11) I. Gilbert, S. Bull, T. Evans, A. Jack, D. Stephenson, and A. De Sa, Effects of processing upon the properties of soft magnetic composities for low loss application, J. Matter Sci., 39, 457 (2004).
  * (12) L. P. Lefebvre, S. Pelletier and C. Gélinas, Effect of electrical resistivity on core losses in soft magnetic iron powder materials, J. Mag. Mag. Mater. 176, L93 (1997).
  * (13) P. Lemieux, R. Guthrie, M. Isac, Optimizing Soft Magnetic Composites for Power Frequency Applications and Power-Trains, JOM, 64, 374-387 (2012).
  * (14) C. Gelinas, S. Pelletier, P. Lemieux, L. Azzi, Properties and Processing of Improved SMC Materials, Advances in Powder Metallurgy and Particulate Materials (Princeton, NJ, USA: MPIF), 9-108 ( 2005).
  * (15) K. Sokalski, J. Szczygłowski, M. Najgebauer and W. Wilczyński, COMPEL: Int. J. Comput. Math. Electr. Electron. Eng.,26, 640-649 ( 2007).
  * (16) K. Sokalski, J. Szczygłowski, Acta. Phys. Pol. A, 115, 920-924 (2009).
  * (17) K. Sokalski, J. Szczygłowski, Acta. Phys. Pol. A, 117, 497-499 (2010).
  * (18) K. Sokalski, J. Szczyg owski and W. Wilczy ski , Scaling Conception of Energy Loss’ Separation in Soft Magnetic Materials , arXiv:1111.0939 (2011).
  * (19) J. A. Nelder and R. Mead, A simpleks method for a function minimization, The Computer Journal, 7, 308-313 (1965).
  * (20) B. Jankowski, PhD Thesis (in Polish), Warsaw 2013.


# The privacy protection effectiveness of the video conference platforms’
virtual background and the privacy concerns from the end-users

Shijing He<EMAIL_ADDRESS>Cisco Systems, Inc.Shanghai CityChina and
Yaxiong Lei  University of St.AndrewsSt.AndrewsUK<EMAIL_ADDRESS>
###### Abstract.

Due to the abrupt arise of pandemic worldwide, the video conferencing
platforms are becoming ubiquitously available and being embedded into either
various digital devices or the collaborative daily work. Even though the
service provider has designed many security functions to protect individual’s
privacy, such as virtual background (VB), it still remains to be explored that
how the instability of VB leaks users’ privacy or impacts their mentality and
behaviours. In order to understand and locate implications for the contextual
of the end-users’ privacy awareness and its mental model, we will conduct
survey and interviews for users as the first stage research. We will raise
conceptual challenges in terms of the designing safety and stable VB, as well
as provide design suggestions.

video conference platforms, virtual background, privacy, mental model,
contextual

††ccs: Security and privacy Usability in security and privacy

##  1\. Introduction

Due to the global circumstances caused by the COVID-19 pandemic, an
unprecedented increase in working from home (WFH) or remotely working happens
worldwide. Everything has become ”remotely”, and WFH is getting a facelift and
becoming more popular. Video conference platforms (e.g. Zoom, Webex and
Microsoft Teams) have played an alternative role, and have replaced the
office-based work tradition. Nowadays, people are getting used to the WFH
mode. Bloom et al. (Bloom et al., 2015) noticed WFH could increase
productivity by 13%. The reason why the performance increase lies in the quiet
and convenient working environment that could lead to more calls per minute.

In addition, the aforementioned platforms and companies continue to actively
optimize global customer’s feedback and usage worldwide, and provide several
new functions to protect the meeting security. Virtual background (VB), as one
of a privacy protection design, has been widely applied on different video
conference platforms. Users could apply VB to limit distractions and maintain
the surrounding environment to avoid the privacy leaking, as well as satisfy
their aesthetic needs. VB also allows users to blur background, change
pictures or customize background accordingly. In particular, Yao et al. (Yao
et al., 2017) has pointed out that blurring is an effective method to protect
people’s sensitive data.

Although the VB has provided a good chance to show user’s personality
(aesthetic taste or personal preference), or as a business requirement (sales
meet with customers on behalf of the company), education demands (academic
conference, teaching needs), people started to concern about whether the VB
will leak their privacy or not. The instability of the segmentation algorithm
also has its limitations on leaking people’s environment data or potential
privacy. Our future research will aim to investigate this question and explore
the privacy mental model of participants, as well as their privacy awareness
in different scenarios.

Figure 1. 2020 Cisco. Webex Virtual Background. (https://blog.webex.com/wp-
content/uploads/2020/08/Virtual-Background-in-Meeting-
Spheres-e1596570058514.png)

##  2\. Related works

There has been a wide investigation about using video conference platforms as
effective remotely virtual tools in different fields. For instance, in
educational field (Barbosa and Barbosa, 2019; Serhan, 2020; Lowenthal et al.,
2020; Kohnke and Moorhouse, 2020),qualitative data collection and research
method (Gray et al., 2020; Archibald et al., 2019), Bailenson (Bailenson,
2021) has provided 4 explanations to furtherly illustrate the causes of the
Zoom fatigue: ”excessive amounts of close-up eye gaze, cognitive load,
increased self-evaluation from staring at video of oneself, and constraints on
physical mobility.” Due to the amount of time in video conferencing, nearly
half users prefer to work from the living room or bedroom (Nulab, 2020). Conti
et al. (Conti et al., 2021) also noticed the similar privacy concerns about
VB, that is, users could conceal where they are with the VB for the purpose of
cheating other participants in the meeting.

Online meeting or live streaming platforms, similar to aforementioned video
conferencing platforms, have provided VB features and offered some room for
customization. For example, users can upload their own preferred VB images. By
semantic segmentation of each frame in the video in real time, the algorithm
separates the user from the surrounding environment and replaces it with a VB.
Video meeting as a high-frequency scenario in daily work, along with the
popularity of WFH, the protection of user privacy has put forward a higher
demand.

Compared with the high-performance servers in the cloud, the video
conferencing carriers of personal scenarios are mainly laptops of various
colors. There are three main factors—various models of laptops’ performances,
highly real-time requirements of video conferencing, various online-meeting
contexts, that have put more demanding demands on the performance of the end-
side (or say user-side) algorithm.

The implementation of the real-time feature needs the end-side portrait
segmentation model to be light enough, while such a small model processing
requires more sensitive data and is weak in some difficult scenarios (such as
portrait edge and background similarity, etc.). These will easily lead to
background misclassification into portraits, blurred portrait edges, etc (Lin
et al., 2021). These problems have directly resulted in the performance and
accuracy of the algorithm, the drawback of the segmentation algorithm would
cause the virtual background to be unstable. There is also a certain privacy
exposure problem, along with the aliasing and other segmentation problems,
which will lead to the risk of potential leakage to the environment behind the
participants.

##  3\. Research Aim and Question

According to the internal user experience feedback and user testing data when
using Cisco Webex VB from Webex Design Team, we have found that the
circumstance of WFH could bring a certain risk of an individual’s privacy
leakage, especially if users apply the VB or blur background. We will not only
study how privacy concerns proceed practically in a range of circumstances or
scenarios with classic qualitative research approaches, but also be grounded
in actual design practices. We believe that the experimental result will be
more likely to be adopted. Our main research questions are:

RQ1: What is the possibility of causing panic and discomfort for end-users
when there are sudden privacy leaks caused by VBs in meetings?

RO2: Is the power differential (Bernd et al., 2020) or identity (e.g.,
interviewer and interviewee) in meetings as a key factor that can impact the
mental model of privacy awareness when the sudden privacy leaks caused by VBs
in meetings?

To answer these questions, we will conduct an empirical investigation from the
end-user’s perspective, which contains scenario-based surveys (see Methodology
and Experimental Design). We will conduct semi-structured interviews based on
survey results, in order to provide personal information and insights in the
context of using video conferencing platforms.

##  4\. Methodology and Experimental Design

The early stage of research would be based on the Grounded theory that was
applied to the existing privacy qualitative studies (Abdi et al., 2021; Westin
and Chiasson, 2021). In particular, the interviewing and analysis will be
conducted simultaneously and iteratively by using the Grounded theory. We will
conduct online surveys for end-users, and set semi-structured interviews for
the respondents who joined the survey. This aims to investigate the privacy
perception and personal awareness in different scenarios and contextual
environment factors. The participants will be the experienced end-users.
Moreover, the bystanders (The bystander types are basically include the
cohabitants, acquaintances or strangers (e.g. delivery man or someone suddenly
flashed ), or other potential factors maybe impact the meeting status, such as
pets, or sudden accidents.) will also be evaluated because they may appear in
the camera randomly, due to the instability of the recognition algorithm. This
might be a good start to understand the end-user’s mental model in which the
respondents face the emergency situation in the various scenarios; as well as
to investigate the differences of the power differentials and identity. To sum
up, we will focus on end-users and their mental model depending on contexts
and scenarios.

Additionally, including aforementioned cases, we also have extended the
privacy factors and dimensions that containing: emergencies (e.g. external
influence, technical limitation), the number of participants, social relations
(stranger-intimacy), power differential, meeting attributes (official or
unofficial), benefit, network types, duration time, recurrence, environment
and number of bystanders. Considering different privacy awareness factors, we
will create six main scenarios—- 1) Business requirements, 2) Education
demands, 3) Family meeting, 4) Job interview, 5) Online conference or event
and 6) Enterprise internal meeting.

The ethics forms, survey and remuneration mechanism will be reviewed and
provided by the Institutional Review Board (IRB) of the universities. The
researchers who also participate in the related training will get the
certificates. Once complete the pilot test and final version data collection
from Qualtrics, we will conduct interviews for participants who joined the
survey investigation. The series of surveys consisted of multiple types of
questions in both open-ended and closed-ended form.

##  5\. Future Work

The future work will basically conduct in two mainly direction based
contextual scenarios: the mental model of privacy awareness and algorithm
optimization. To meet up with the expected result, we will continue to
evaluate the mental model and user’s privacy awareness so as to face future
emergency situations. However, it is a bit difficult to observe the
participants’ behaviours directly and immersively as the lab-based interview
could cause bias due to the Hawthorne effect. Hence, we will apply an
ethnographic approach (Gan et al., 2020) and use quantitative data analysis
(e.g., mental workload data or gaze data) to conduct a field study to get
participants’ behaviours and micro-expression data from closer observation. In
order to improve the existing user experience and the effectiveness of privacy
protection, we will expect to deliver the end-users’ mental model of privacy
awareness. Moreover, due to the special contextual integrity (e.g.,
environment) of WFH, it should also be noticed whether these WFH circumstances
are affected by personal smart home devices. For instance, the smart home
device could potentially track voice data or even sensitive data when a user
is having a confidential meeting.

In addition, the algorithm design strategy has a noticeable impact on the VB’s
privacy improvement, improving the segmentation quality and effectiveness of
algorithms is also urgently needed. Nowadays, due to the processing
performance of end devices, lightweight models are usually used for portrait
segmentation. Some work (Lin et al., 2021) applied the idea of PointRend
(Kirillov et al., 2020) to optimize the edge blurring problem and the
performance problem under high resolution image conditions. It has to be
noticed that only a few areas need fine segmentation and most areas need
coarse segmentation. The problem of information exposed from coarse
segmentation should also be brought to our attention.

## References

  * (1)
  * Abdi et al. (2021) Noura Abdi, Xiao Zhan, Kopo M Ramokapane, and Jose Such. 2021\.  Privacy Norms for Smart Home Personal Assistants. In _Proceedings of the 2021 CHI Conference on Human Factors in Computing Systems_. 1–14.
  * Archibald et al. (2019) Mandy M Archibald, Rachel C Ambagtsheer, Mavourneen G Casey, and Michael Lawless. 2019.  Using zoom videoconferencing for qualitative data collection: perceptions and experiences of researchers and participants.  _International Journal of Qualitative Methods_ 18 (2019), 1609406919874596\.
  * Bailenson (2021) Jeremy N Bailenson. 2021\.  Nonverbal overload: A theoretical argument for the causes of Zoom fatigue.  _Technology, Mind, and Behavior_ 2, 1 (2021).
  * Barbosa and Barbosa (2019) Tamara JG Barbosa and Mary Jo Barbosa. 2019.  Zoom: An innovative solution for the live-online virtual classroom.  _HETS Online Journal_ 9, 2 (2019).
  * Bernd et al. (2020) Julia Bernd, Ruba Abu-Salma, and Alisa Frik. 2020\.  Bystanders’ Privacy: The Perspectives of Nannies on Smart Home Surveillance. In _10th $\\{$USENIX$\\}$ Workshop on Free and Open Communications on the Internet ($\\{$FOCI$\\}$ 20)_.
  * Bloom et al. (2015) Nicholas Bloom, James Liang, John Roberts, and Zhichun Jenny Ying. 2015.  Does working from home work? Evidence from a Chinese experiment.  _The Quarterly Journal of Economics_ 130, 1 (2015), 165–218.
  * Conti et al. (2021) Mauro Conti, Simone Milani, Ehsan Nowroozi, and Gabriele Orazi. 2021.  Do Not Deceive Your Employer with a Virtual Background: A Video Conferencing Manipulation-Detection System.  _arXiv preprint arXiv:2106.15130_ (2021).
  * Gan et al. (2020) Yumei Gan, Christian Greiffenhagen, and Stuart Reeves. 2020\.  Connecting distributed families: Camera work for three-party mobile video calls. In _Proceedings of the 2020 CHI Conference on Human Factors in Computing Systems_. 1–12.
  * Gray et al. (2020) Lia M Gray, Gina Wong-Wylie, Gwen R Rempel, and Karen Cook. 2020.  Expanding qualitative research interviewing strategies: Zoom video communications.  _The Qualitative Report_ 25, 5 (2020), 1292–1301.
  * Kirillov et al. (2020) Alexander Kirillov, Yuxin Wu, Kaiming He, and Ross Girshick. 2020\.  Pointrend: Image segmentation as rendering. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 9799–9808.
  * Kohnke and Moorhouse (2020) Lucas Kohnke and Benjamin Luke Moorhouse. 2020.  Facilitating synchronous online language learning through Zoom.  _RELC Journal_ (2020), 0033688220937235\.
  * Lin et al. (2021) Shanchuan Lin, Andrey Ryabtsev, Soumyadip Sengupta, Brian L Curless, Steven M Seitz, and Ira Kemelmacher-Shlizerman. 2021.  Real-time high-resolution background matting. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 8762–8771.
  * Lowenthal et al. (2020) Patrick Lowenthal, Jered Borup, Richard West, and Leanna Archambault. 2020.  Thinking beyond Zoom: Using asynchronous video to maintain connection and engagement during the COVID-19 pandemic.  _Journal of Technology and Teacher Education_ 28, 2 (2020), 383–391.
  * Nulab (2020) Nulab. 2020.  _Adjusting to remote work_.  Retrieved August 26, 2021 from https://nulab.com/blog/collaboration/adjusting-to-remote-work/
  * Serhan (2020) Derar Serhan. 2020\.  Transitioning from face-to-face to remote learning: Students’ attitudes and perceptions of using Zoom during COVID-19 pandemic.  _International Journal of Technology in Education and Science_ 4, 4 (2020), 335–342.
  * Westin and Chiasson (2021) Fiona Westin and Sonia Chiasson. 2021.  “It’s So Difficult to Sever that Connection”: The Role of FoMO in Users’ Reluctant Privacy Behaviours. In _Proceedings of the 2021 CHI Conference on Human Factors in Computing Systems_. 1–15.
  * Yao et al. (2017) Yaxing Yao, Huichuan Xia, Yun Huang, and Yang Wang. 2017\.  Privacy mechanisms for drones: Perceptions of drone controllers and bystanders. In _Proceedings of the 2017 CHI Conference on Human Factors in Computing Systems_. 6777–6788.


# Dynamics of the conserved net-baryon density near QCD critical point within
QGP profile

Shanjin Wu1<EMAIL_ADDRESS>1School of Nuclear Science and Technology,
Lanzhou University, Lanzhou 730000, China

###### Abstract

This paper investigates the dynamics of the net-baryon multiplicity
fluctuations near the QCD critical point, within the inhomogeneous temperature
and baryon chemical potential profile of quark-gluon plasma from the
hydrodynamics. We numerically solve the Langevin dynamics of conserved net-
baryon density and borrow the temperature and chemical potential profile from
hydrodynamic simulation. It is found that the local systems at different
rapidities reach the critical point at different proper times, owing to the
inhomogeneous temperature and chemical potential. As a result, we find the
pronounced enhancement of the magnitude for the net-baryon multiplicity
fluctuations with large rapidity acceptance at the freeze-out surface, which
is the consequence of the combined effect of critical slowing down and
inhomogeneous profile.

##  I Introduction

Searching the Quantum Chromodynamics (QCD) critical point is one of the most
important goals of the relativistic heavy-ion collisions. The transition from
the quark-gluon plasma (QGP) phase to the hadron phase is revealed as a
crossover at the vanishing baryon chemical potential ($\mu_{B}\simeq 0$) by
lattice QCD [1, 2, 3, 4]. While the first-order phase transition together with
the critical point is predicted by the effective theories of QCD at finite
chemical potential [5, 6, 7, 8]. The most important property of the critical
point is the long-range correlation and large fluctuations. As a consequence,
the non-monotonic behavior is conjectured as the characteristic signature of
the QCD phase transition [9, 10].

The Beam Energy Scan program at RHIC has been dedicated to exploring the QCD
phase structure. The preliminary non-monotonic behavior of net-proton
multiplicity fluctuations as a function of colliding energy has been observed
[11, 12]. However, the statistics of the measurement are insufficient to
conclude the observation of this non-monotonicity so far, which is looking
forward to the higher statistics in the coming Bean Energy Scan phase two
program.

On the other hand, the final confirmation of the existence of the QCD critical
point requires the comparison between the experimental measurement and
theoretical prediction. After decades of efforts, remarkable progress has been
made in the theoretical modeling of the dynamics near the QCD critical point
within the complex evolution of the relativistic heavy-ion collisions. Please
see e.g., Refs.[13, 14, 15, 16, 17, 18] for recent review. Due to the fast-
expanding QGP fireball, the non-equilibrium critical fluctuations become non-
trivial compared with the equilibrium ones. For example, the magnitude of the
fluctuations is suppressed [19, 20], the sign could be reversed [21], and the
largest fluctuations not necessarily correspond to the trajectory closest to
the critical point [22]. Subsequently, several models have been developed to
incorporate the dynamics of the critical effects. For instance, the dynamics
of conserved charge were constructed [23, 24, 25, 26]and non-monotonic
behavior of multiplicity fluctuations with respect to acceptance was predicted
[25, 26]. To investigate the dynamics of the QGP fireball near the QCD
critical point, the dynamics of the critical fluctuations were also coupled to
the dynamics of the fireball as an additional degree of freedom, such as the
stochastic dynamics of order parameter field in non-equilibrium chiral
hydrodynamics [27, 28, 29], the stochastic noise terms in the fluctuating
hydrodynamics [30, 31, 32, 17] and the slow mode in the hydro+ [33].

The dynamics of the conserved quantity are of particular interest because the
corresponding diffusion process requires time. The correlation between
particles far away from each other preserves the early evolution of diffusion
[25], and the non-monotonic behavior with respect to acceptance could be
regarded as the imprint of the evolving trajectory passing through the
critical region. In addition, multi-particle correlation becomes indispensable
as the system is extremely close to the critical point. Thus, the dynamics of
conserved net-baryon density with higher-order correlation was constructed in
the longitudinal Bjorken expansion system with uniform temperature and
chemical potential profile [26]. It was found that the pronounced minimum of
kurtosis presents at the intermediate rapidity due to the existence of the
critical point.

In the realistic context, the temperature and chemical point are not uniform
across the QGP profile because of the baryon stopping effects [34, 35]. As a
result, different regions of the QGP profile experience distinct trajectories
across the critical region [36], and the rapidity dependence is expected to
deform non-trivially [37]. In the dynamics of conserved variables, the
realistic QGP profile is particularly important because its rapidity
dependence on correlation preserves the history along the evolving trajectory.
Therefore, it is essential to investigate the rapidity dependence of the
multiplicity fluctuations within the more realistic QGP profile. This work
studies the dynamics of the conserved net-baryon density near the critical
point within the inhomogeneous temperature and chemical potential QGP profile.
With the profile obtained from the hydrodynamic simulation, we can study the
multiplicity fluctuations at freeze-out hypersurface. It is found that the
magnitude of the second-order cumulants and kurtosis enhance significantly at
large rapidity, which is the consequence of combined effects of critical
slowing down and inhomogeneous QGP profile.

##  II Model and Setups

To study the dynamics of QCD phase transition, one of the natural and simple
dynamical models is just focusing on the dynamics of the order parameter
field. However, it was pointed out that the dynamics of conserved baryon
density is the slowest mode of the system near the critical point when
considering the coupling between baryon density and order parameter field
[38]. Therefore, the simplest description of the dynamics near the QCD
critical point is only considering the dynamics of conserved baryon density
and treating the dynamics of other degrees of freedom as the heat bath.

This study focuses on the $1+1$-dimensional conserved dynamics of net-baryon
density $n_{B}$ near the QCD critical point along the longitudinal direction
in Milne frame, i.e. the proper-time $\tau=\sqrt{t^{2}-z^{2}}$ and space-time
rapidity $\eta=\frac{1}{2}\ln[(t+z)/(t-z)]$:

$\displaystyle\partial_{\tau}n_{B}(\tau,\eta)=\frac{D\chi_{2}T}{\tau}\frac{\partial^{2}}{\partial\eta^{2}}\biggl{(}\frac{1}{T}\frac{\delta
F}{\delta n_{B}}\biggr{)}-\partial_{\eta}\zeta(\tau,\eta),$ (1)

with the white noise $\zeta(\tau,\eta)$ has only one non-zero correlation:

$\displaystyle\langle\zeta(\tau,\eta)\zeta(\tau^{\prime},\eta^{\prime})\rangle=\frac{2D\chi_{2}T}{A\tau}\delta(\tau-\tau^{\prime})\delta(\eta-\eta^{\prime}).$
(2)

Here, $D$ is the diffusion coefficient, $T$ is the temperature, and the
transverse area is set as $A=1\mbox{fm}^{2}$ following Ref. [24]. The equation
of $n_{B}$ can be derived from the conservation of net-baryon number
$\partial_{\mu}N^{\mu}=0$ within the Bjorken flow
$u^{\mu}=(\cosh\eta,0,0,\sinh\eta)$. Please see Refs. [26, 39] for detail and
e.g., Refs. [40, 41, 42] for the further extension.

The effective potential $F$ near the critical point can be parameterized as
[26]:

$\displaystyle F[n_{B}]=\int d^{2}x_{\perp}d\eta\biggl{[}$
$\displaystyle\frac{1}{2}\frac{1}{\chi_{2}\tau}n_{B}(\tau,\eta)^{2}+\frac{1}{2}\frac{K}{\tau^{3}}\biggl{(}\frac{\partial}{\partial\eta}n_{B}(\tau,\eta)\biggr{)}^{2}$
$\displaystyle\quad+\frac{1}{3}\frac{\lambda_{3}}{\tau^{2}}n_{B}(\tau,\eta)^{3}+\frac{1}{4}\frac{\lambda_{4}}{\tau^{3}}n_{B}(\tau,\eta)^{4}\biggr{]},$
(3)

where $K$ represents the surface tension coefficient and is treated as a
constant in this study. The second-order baryon susceptibility $\chi_{2}$,
third- and fourth-order coupling coefficients $\lambda_{3},\lambda_{4}$
contain two parts: regular part estimated from lattice QCD simulation and
singular part
($\chi_{2}^{\mbox{cri}}(T,\mu),\lambda^{\mbox{cri}}_{3}(T,\mu),\lambda^{\mbox{cri}}_{4}(T,\mu)$)
mapped from three-dimensional Ising model
($\chi_{2}^{\mbox{cri}}(r,h),\lambda^{\mbox{cri}}_{3}(r,h),\lambda^{\mbox{cri}}_{4}(r,h)$)
(here $r,h$ are Ising variables). A complete description of the
parameterization of the coefficients can be found in Appendix A.

The parametrization of the baryon susceptibility and coupling coefficients on
the QCD phase diagram requires the time evolution profiles for the QGP
fireball: temperature $T(\tau,\eta)$ and chemical potential $\mu(\tau,\eta)$.
The time evolution of the QGP profiles is obtained from the 3+1-d hydrodynamic
MUSIC simulation [43], as in Ref.[22]. The initial profiles constructed from
the transport model AMPT [44] and the Equation of State are input with the
lattice simulation results, incorporating with a critical point [45]. To
roughly fit the experimental measurement of Au+Au collisions at 19.6GeV, such
as the multiplicity, spectra, and flow, the parameters for the hydrodynamic
model are set as [22]: the starting proper time for the AMPT initial condition
$\tau_{\tiny 0I}=0.4$fm, the starting time for the hydrodynamic evolution
$\tau_{\tiny 0h}=1.2$fm, the specific shear and bulk viscosity
$\eta/s=0.08,\zeta/s=0.01$, and the switching energy density from QGP to
hadronic phase $e_{\tiny sw}=0.1\mbox{GeV/fm}^{3}$. We fix the position of the
critical point as $(T_{c},\mu_{c})=(0.16\mbox{GeV},0.155\mbox{GeV})$. To
illustrate the effects of the realistic QGP profile, we perform the simulation
with two scenarios:

  1. 1.

Scenario I: We extract the temperature and chemical potential along the
longitudinal direction by averaging the temperature and chemical potential
profiles over the whole QGP fireball with the energy density as the weight:
$T(\tau,\eta)=\langle T(\tau,\bm{x})\rangle$ and
$\mu(\tau,\eta)=\langle\mu(\tau,\bm{x})\rangle$. $T(\tau,\eta)$ and
$\mu(\tau,\eta)$ profile from hydrodynamic simulation at 19.6GeV with $0-5\%$
centrality is presented in Figs. 1

  2. 2.

Scenario II: For comparison, we also study the net-baryon density dynamics
within the uniformed background of temperature and chemical potential along
the longitudinal direction. The temperature and chemical potential along the
longitudinal direction are fixed the ones at $\eta=0$: $T(\tau,\eta)=\langle
T(\tau,x_{\perp},\eta=0)\rangle$ and
$\mu(\tau,\eta)=\langle\mu(\tau,x_{\perp},\eta=0)\rangle$.

As depicted in the upper panel of Figs. 1, the temperature $T(\tau,\eta)$ is
highest at the center of the fireball and decreases with increasing rapidity
$\eta$, and this behavior becomes smoother as a function of proper time
$\tau$. Conversely, the chemical potential $\mu_{(}\tau,\eta)$ reaches its
minimum at the fireball’s center and exhibits a larger value for higher
rapidity. This behavior illustrates the baryon-stopping effects along the
longitudinal direction. This study only focuses on the dynamics with the QGP
fireball, there is no value outside the fireball as the system has freeze-out
and turned into the hadronic phase.

Figure 1: Time evolution of the temperature $T$ (upper panel) and chemical
potential $\mu$ (lower panel) profiles in the rapidity $\eta$ space, obtained
from the hydrodynamic simulation of Au+Au collisions at 19.6GeV with $0-5\%$
centrality. The $T$ and $\mu$ profiles are depicted from $\tau_{0}=3$fm, the
starting time of Eq.1.

It requires the solution of the Eq.(1) to obtain the time evolution of the
net-baryon density along the longitudinal axis. In the case of the linear
stochastic diffusion equation, without cubic and quartic terms of Eq.(1), the
dynamics of net-baryon density can be solved analytically. However, it
requires numerical simulation of the non-linear equation Eq.(1) when
considering cubic and quartic terms. One of the most popular algorithms for
this diffusive stochastic equation is the explicit Forward or Backward Euler
Method. However, the Euler Method is conditional stable and the numerical
simulation is extremely inefficient near the critical point. In this work, we
perform, the numerical simulation of Eq.(1) with the Saul’yev scheme, which is
an explicit unconditional stable scheme. In 1957, V.K.Saul’yev proposed the
so-called asymmetric methods in the simulation of diffusion equation [46, 47]
and has been successfully implemented in the study of Cahn-Hilliard equation
[48]. For the detailed numerical implementation of Eq.(1), please see
Appendix.B.

With the configuration of the net-baryon density $n_{B}(\tau,\eta)$ after
solving Eq.(1), we are able to study the dynamics of the multiplicity
fluctuations, which can be calculated as follows:

$\displaystyle C_{1}$ $\displaystyle=\langle N_{B}\rangle,\quad
C_{2}=\langle(\delta N_{B})^{2}\rangle,$ $\displaystyle C_{4}$
$\displaystyle=\langle(\delta N_{B})^{4}\rangle-3\langle(\delta
N_{B})^{2}\rangle^{2},$ (4)

and $\sigma^{2}=C_{2},\kappa\sigma^{2}=C_{4}/C_{2}$. Here
$N_{B}=\int^{\Delta\eta}_{-\Delta\eta}d\eta n_{B}(\tau,\eta)$ denotes the
multiplicity of net-baryon, $\langle\cdots\rangle$ represents averaging over
the events and the multiplicity event-by-event fluctuations are defined as
$\delta N_{B}=N_{B}-\langle N_{B}\rangle$.

After the evolution of the net-baryon density with the QGP profile, systems
freeze out and turn into the hadronic phase when the energy density is below
the switching energy density $e_{sw}$. The freeze-out hyper-surface from the
hydrodynamic simulation is employed and shown as the edge of the QGP profile
in Figs.1. To calculate the multiplicity fluctuations after freeze-out, we
employ the Cooper-Frye formula [49]:

$\displaystyle N_{B}=g\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{p^{0}}\int
d\sigma_{\mu}p^{\mu}f(\bm{x},\bm{p}),$ (5)

where $g$ is the spin degeneracy and we take the Boltzmann approximation:

$\displaystyle f(\bm{x},\bm{p})=\exp[-(p_{\mu}u^{\mu}-\mu)/T].$ (6)

Following Refs. [30, 39, 50], we employ the Bjorken limit:
$d\sigma_{\mu}p^{\mu}=dx_{\perp}d\eta\tau_{f}m_{\perp}\cosh(\eta-y)$ and
$p_{\mu}u^{\mu}=m_{\perp}\cosh(\eta-y)$, with
$m_{\perp}=\sqrt{\bm{p}_{\perp}^{2}+m^{2}}$ and $\bm{p}_{\perp}$ is transverse
momentum. Therefore, the expression of $N_{B}$ reads:

$\displaystyle
N_{B}=\frac{2gA}{(2\pi)^{2}}\int^{\Delta\eta}_{-\Delta\eta}d\eta\tau_{f}\exp\biggl{(}\frac{\mu}{T}\biggr{)}T^{3}\biggl{(}\frac{m}{T}\biggr{)}^{2}K_{2}\biggl{(}\frac{m}{T}\biggr{)},$
(7)

where $K_{n}(x)$ is the irregular modified Bessel function of order-$n$. The
integral over the transverse space $dx_{\perp}$ is represented by
$A=1\mbox{fm}^{2}$ [24]. As there is only one fluctuating variable $n_{B}$ in
this framework, following Ref.[39], the multiplicity fluctuations is obtained
as

$\displaystyle\delta
N_{B}=\frac{2gAm^{2}}{(2\pi)^{2}}\int^{\Delta\eta}_{-\Delta\eta}d\eta\exp\biggl{(}\frac{\mu}{T}\biggr{)}\frac{\delta
n_{B}(\tau_{f},\eta)}{\chi_{2}}K_{2}\biggl{(}\frac{m}{T}\biggr{)},$ (8)

by considering the definition of the susceptibility $\delta\mu=\delta
n_{B}(\tau_{f},\eta)/(\tau_{f}\chi_{2})$ in Milne frame. Here, the integration
over $d\eta$ is performed at the edge of the hydrodynamic profile, the freeze-
out surface. The fluctuations of the net-baryon density
$n_{B}(\tau_{f},\eta)$, as well as the temperature $T(\tau_{f},\eta)$,
chemical potential $\mu(\tau_{f},\eta)$ and susceptibility
$\chi_{2}(\tau_{f},\eta)$ are also extracted from the hypersurface. Various
orders of the cumulants of the multiplicity fluctuations are calculated as
Eq.(II) accordingly.

Figure 2: Scenario I:Time evolution of the second order net-baryon
susceptibility $\chi_{2}$ (upper panel) and fourth order coupling coefficient
$\lambda_{4}$ (lower panel) profiles in the rapidity $\eta$ space.

Figure 3: Scenario II: Similar with Fig.2, but within the uniform temperature
$T$ and chemical potential $\mu_{B}$ profile across rapidity $\eta$.

##  III RESULT AND DISCUSSION

This work focuses on the inhomogeneous QGP background effects on the dynamics
of the conserved net-baryon fluctuations near the QCD critical point. The
inhomogeneous temperature and chemical potential profile are extracted from
the hydrodynamic simulation. This section will show the net-baryon
susceptibility and coupling constants within the inhomogeneous profile, as
well as its impact on the dynamics of the net-baryon density near the QCD
critical point.

Figure 4: Rapidity acceptance $\Delta\eta$ dependence of the second order
multiplicity fluctuations $\sigma^{2}/\Delta\eta$ (upper panel) and kurotsis
$\kappa\sigma^{2}$ (lower panel). Curves with different colors correspond to
the multiplicity fluctuations at different proper times $\tau$. Dashed for the
inhomogeneous profile (Scenario I) and solid for the uniform ones (Scenario
II).  Figure 5: Time evolution of the second order net-baryon susceptibility
$\sigma^{2}/\Delta\eta$ (upper panel) and kurtosis $\kappa\sigma^{2}$ (lower
panel). Curves with different colors correspond to the multiplicity
fluctuations with different rapidity regions $\eta$. Dashed for the
inhomogeneous profile (Scenario I) and solid for the uniform ones (Scenario
II).

Figs.2 plot the time evolution of second-order net-baryon susceptibility
$\chi_{2}(\tau,\eta)$ and fourth-order coupling coefficient
$\lambda_{4}(\tau,\eta)$ across the QGP profile. Due to the inhomogeneous
temperature and chemical potential across rapidity space 1, the local regions
at different rapidities reach the critical point $(T_{c},\mu_{c})$ at varying
times. Consequently, the strength of critical effects varies along the
rapidity axis for a given proper time. As shown in Fig.2, the local regions at
finite rapidities (around $\eta=\pm 2$) are the points closest to the QCD
critical point, owing to the combined influence of the inhomogeneous $T(\eta)$
and $\mu_{B}(\eta)$. Therefore, it is natural to expect that multiplicity
fluctuations would exhibit nontrivial behavior in response. For the
comparison, Scenario II, the time evolution of the coefficients for the
uniform temperature and chemical potential has also been shown in Fig.3, where
the proper time evolution of $\chi_{2}$ and $\lambda_{4}$ are uniform across
the $\eta$ space.

With the temperature and chemical potential profile obtained from hydrodynamic
simulation, it is possible to study the net-baryon number fluctuations during
the QGP evolution. Let’s start with the rapidity dependence of the net-baryon
number fluctuations. Figs.4 show the rapidity dependence of the net-baryon
number fluctuations with fixed proper time $\tau=3.9,4.8$ and $7.5$fm,
respectively. The second order net-baryon number cumulant
$\sigma^{2}/\Delta\eta$ (kurtosis $\kappa\sigma^{2}$) shows non-monotonic
behavior with increasing rapidity and reaches a maximum (minimum) at small
$\Delta\eta$, which agrees with previous studies [25, 26]. One of the most
important properties of the dynamics for the conserved variable is that the
diffusion process consumes time (Please see Ref. [25] for detail). As a
result, the correlation between particles with small $\Delta\eta$ encodes the
late-stage dynamics of diffusion, while the one with large $\Delta\eta$
preserves the early evolution of diffusion. As the system scans the critical
regime, the susceptibility exhibits a peak with increasing proper time $\tau$.
Therefore, the correlation behaves non-monotonically with increasing rapidity
interval $\Delta\eta$, which can be regarded as the imprints from the critical
fluctuations. We also compare the $\sigma^{2}/\Delta\eta$ and
$\kappa\sigma^{2}$ in uniform (Scenario II, solid curves) with inhomogeneous
(Scenario I, dashed curves) profile. In the case of second-order cumulant
$\sigma^{2}/\Delta\eta$, the difference between Scenario II and I is
negligible. For kurtosis, the discrepancy is also small at small rapidity but
becomes significant for large rapidity. This can be understood that the
difference of the transport coefficients ($\chi_{2}$ and $\lambda_{4}$) in
Eq.(1) between Scenario I and II becomes pronounced at large $\eta$, by
comparing Figs.2 and 3.

Since the realistic inhomogeneous QGP profile has an impact on the number
fluctuations at large rapidity, let’s focus on the time evolution of the
system fluctuations at different $\eta$. With the number fluctuations obtained
from Eq.(II), Figs.5 present the time evolution of $\sigma^{2}/\Delta\eta$ and
$\kappa\sigma^{2}$ with different $\eta$ intervals both for uniform (Scenario
II) and inhomogeneous (Scenario I) profile. In these figures, solid curves
represent the time evolution of the fluctuations with realistic inhomogeneous
QGP profile (Scenario I), while the dashed curves correspond to the uniform
profile case (Scenario II). As shown in Figs.5, the discrepancy between these
two cases is negligible for midrapidity interval ($-0.7<\eta<0.7$), while
becomes significant at large rapidity intervals ($1.0<\eta<2.4$ and
$1.6<\eta<3.0$) 111Noted that the integral interval of $\sigma^{2}$ in Figs.4
is from $-\Delta\eta$ to $\Delta\eta$. The difference between Scenario I and
II is small when the integration is also performed from $-\Delta\eta$ to
$\Delta\eta$ for the time evolution of $\sigma^{2}/\Delta\eta$ in Figs.5. This
is consistent with Figs.4 and the difference between Scenario I and II in
different rapidity is smeared out from $-\Delta\eta$ to $\Delta\eta$. .

Figure 6: The second order multiplicity fluctuations $\sigma^{2}/\Delta\eta$
(upper panel) and kurtosis $\kappa\sigma^{2}$ (lower panel) along the freeze-
out surface as shown in Fig.1. Dashed for the inhomogeneous profile (Scenario
I) and solid for the uniform ones (Scenario II).

The hydrodynamic simulation provides us with the temperature and chemical
potential profile. After the energy density drops and below the switching
energy density $e_{sw}$, QGP evolution terminates and freezes out into the
hadronic phase. This procedure defines the freeze-out hypersurface, shown as
the edge of the QGP fireball of the ($\tau,\eta$) plane in Figs.1, and enables
us to study the multiplicity fluctuations at the freeze-out hyper-surface.
Following Eq.(8), net-baryon multiplicity fluctuations at the freeze-out
surface are calculated and shown in Figs.6. For illustrative purposes, we
present the multiplicity fluctuations obtained with the uniform $T$ and $\mu$
profile (Scenario II) but still at the edge of the profile 1. Note that the
freeze-out surface for Scenario II is not the edge of the profile 1, and the
solid curves are presented for comparison. In the scenario of a uniform
temperature and chemical potential profile, second order fluctuations
$\sigma^{2}/\Delta\eta$ (kurtosis $\kappa\sigma^{2}$) display an increase
(decrease) first followed by a rapid decline (increase) with increasing
rapidity interval $\Delta\eta$, exhibiting a maximum (minimum) at small values
of rapidity interval $\Delta\eta$. This behavior is consistent with studies on
diffusive dynamics near the critical point [25, 26], suggesting that the
correlation with large $\Delta\eta$ preserves the early dynamics of diffusion.
On the contrary, $\sigma^{2}/\Delta\eta$ and $\kappa\sigma^{2}$ in the case of
the realistic QGP profile deviate significantly with uniform temperature and
chemical potential case at large rapidity interval $\Delta\eta$. As shown in
Figs.1, the fluctuations with different rapidity at the freeze-out surface are
determined by the dynamics with different time $\tau_{f}(\eta)$. The
significant enhancement (decrease) of $\sigma^{2}/\Delta\eta$ and
$\kappa\sigma^{2}$ is a result of the dynamics occurring at small $\tau_{f}$.

To see this, it is instructive to study the rapidity dependence of the
multiplicity fluctuations before freeze-out. This can be achieved by moving
the freeze surface backward to the smaller $(\tau,\eta)$ region. In Fig.7, the
curves in QGP fireball are depicted, along which the variables at smaller
$(\tau,\eta)$ region are extracted. Figs.8 show the susceptibility $\chi_{2}$
and $C_{2}/\Delta\eta$ along the curves in Fig.7. Here
$C_{2}=\langle[\int^{\Delta\eta}_{-\Delta\eta}d\eta\delta n_{B}]^{2}\rangle$,
where the integral is performed along the curves in Fig.7. With increasing
$\tau$ in Scenario I, one could see the peak of susceptibility $\chi_{2}$ only
exhibits along the curves with small $\tau$ and $\eta$, and been significantly
suppressed near the freeze-out surface. This is the impact of the
inhomogeneous QGP profile, in which the system passes through the critical
regime and a pulse of the susceptibility appears and vanishes rapidly. On the
other hand, during the process of diffusion, the pulse of the susceptibility
near the QCD critical point results in the large fluctuations $C_{2}$ at small
$\Delta\eta$, as shown in the dashed-magenta curve in the lower panel of
Figs.8. The large fluctuations $C_{2}$ smears and diffuses with increasing
proper time and rapidity, but can not catches up with the rapid decreasing
$\chi_{2}$, which is known as the critical slowing down effects. In short, the
evolution of the $T$ and $\mu$ profiles drives the susceptibility passing
through the critical point, while the critical slowing down effects preserve
the memory of the large net-baryon density fluctuations at the critical point.
Consequently, the fluctuations in Eq.(8) behave with pronounced enhancement at
large rapidity intervals. On the contrary, one could see the susceptibility
$\chi_{2}$ always has a peak and moves to larger $\eta$ in Scenario II, which
is expected in Fig.3. Therefore, the fluctuations at Eq.(8) drop down rapidly
with increasing $\Delta\eta$ even with the diffusion of net-baryon density
(lower panel of Figs.8).

Figure 7: A sketch of the curves before freeze-out. Along these curves, the
susceptibility $\chi_{2}$ and second order density cumulants
$C_{2}/\Delta\eta$ are extracted (See Fig.8) to study the diffusion of the
fluctuations before freeze-out. Figure 8: Dashed: The susceptibility
$\chi_{2}$ (upper panel) and second order density cumulants $C_{2}/\Delta\eta$
(lower panel) along the curves in Fig.7. $\chi_{2}$ and $C_{2}/\Delta\eta$
with different colors are extracted along the corresponding curves in
Fig.7(Scenario I). Solid: The susceptibility $\chi_{2}$ (upper panel) and
second order density cumulants $C_{2}/\Delta\eta$ (lower panel) along the
curves in Fig.7 with $\chi_{2}$ profile replaced by the Fig.3 (Scenario II).

##  IV Conclusion and Outlook

This work studies the conserved dynamics of the net-baryon multiplicity
fluctuations near the QCD critical point with the inhomogeneous temperature
and chemical potential profile. The dynamics of conserved net-baryon density
are based on the numerical simulation of the 1+1-dimensional stochastic
diffusion equation. The realistic QGP profile is obtained by the hydrodynamic
simulation and the susceptibilities in the diffusion equation are constructed
based on the temperature and chemical potential profile. In this context, the
QGP profile is inhomogeneous and the susceptibilities are non-trivial across
the proper time $\tau$ and rapidity $\eta$ plane. As discovered by the early
studies, the fluctuations of the conserved net-baryon density behave non-
monotonically with the increasing rapidity interval, owing to the fact that
the correlation function in the diffusion process preserves the evolution
history at large rapidity interval. It is found that the influence of the
inhomogeneous profile is negligible at small rapidity but relatively
pronounced at large rapidity, as the susceptibility in the realistic case
deviates from the case of the uniform profile at large rapidity. Furthermore,
the fluctuations on the freeze-out hypersurface have also been investigated.
Comparing the case with a uniform QGP profile, the magnitudes of second-order
cumulants as well as kurtosis in inhomogeneous profile present significant
enhancement at large rapidity. This is the result of the combined effect of
critical slowing down and inhomogeneous profile.

Finally, it deserves to be pointed out that this study of the critical
fluctuations is based on a simplified model, where only the 1+1-dimensional
conserved net-baryon density is considered. Other degrees of freedom in
hydrodynamic evolution have not been considered comprehensively, and are only
regarded as the background. As shown in this paper, the dynamical simulation
of the relevant quantity with a realistic setup is essential for comparison
with experimental measurement. For example, the hydrodynamics coupling with
the additional slow modes has been developed [33], and further incorporating
the higher order slow modes, their extension at freeze-out surface [52], and
corresponding phenomenological study are required for the persuasive
prediction.

## Acknowledgments

The author would like to thank Shian Tang for his contribution at the
beginning of this study and discussion with Huichao Song. This work is
supported by the NSFC under grant No. 12305143 and the China Postdoctoral
Science Foundation under Grant No. 2023M731467.

##  Appendix A Parametrization

The coefficients of the Eq.(1) require to be specified for the simulation of
net-baryon density near the QCD critical point. According to the universal
analysis, one can obtain the diffusion coefficient $D$, second-order baryon
susceptibility $\chi_{2}$, and third- and fourth-order coupling coefficients
$\lambda_{3},\lambda_{4}$ from the mapping three-dimensional Ising model.
Following Refs. [25, 26], these coefficients include regular and singular
parts:

$\displaystyle\chi_{2}$
$\displaystyle=\chi^{\mbox{cri}}_{2}+\chi^{\mbox{reg}}_{2},$
$\displaystyle\lambda_{3}$
$\displaystyle=\lambda^{\mbox{cri}}_{3}+\lambda^{\mbox{reg}}_{3},$
$\displaystyle\lambda_{4}$
$\displaystyle=\lambda^{\mbox{cri}}_{4}+\lambda^{\mbox{reg}}_{4},$ (9)

where the regular coefficients read

$\displaystyle\lambda^{\mbox{reg}}_{3}=-\frac{1}{2}\chi^{\mbox{reg}}_{3}\chi^{-3}_{2},\quad\lambda^{\mbox{reg}}_{4}=\frac{1}{2}(\chi^{\mbox{reg}}_{3})^{2}\chi^{-5}_{2}-\frac{1}{6}\chi^{\mbox{reg}}_{4}\chi^{-4}_{2}.$

With the susceptibility, one has the diffusion coefficient $D=D_{c}/\chi_{2}$
near the critical point. $D_{c}$ and surface tension coefficient are treated
as constants: $D_{c}=16.5$ and $K=0.005$.

The regular susceptibilities are obtained by interpolating interpolate between
the regular susceptibilities at hadronic phase $\chi^{\mbox{\tiny H}}_{n}$ and
QGP phase $\chi^{\mbox{\tiny QGP}}_{n}$

$\displaystyle\chi^{\mbox{reg}}_{n}=\chi^{\mbox{\tiny
H}}_{n}+(\chi^{\mbox{\tiny QGP}}_{n}-\chi^{\mbox{\tiny H}}_{n})S(T),$

where the interpolating function is $S(T)=\frac{1}{2}(1+\tanh((T-T_{c})/\delta
T))$ and the width of the transition region is set as $\delta T=0.01$GeV.
Following Ref. [39], the second-order susceptibility $\chi^{\mbox{\tiny
H,latt}}_{2}/T^{2}_{\mbox{\tiny H}}=1/3$ for hadronic phase and
$\chi^{\mbox{\tiny QGP,latt}}_{2}/T^{2}_{\mbox{\tiny QGP}}=2/3$ for QGP phase,
where we set the low-temperature limit $T_{\mbox{\tiny H}}=0.1$GeV and the
high-temperature limit $T_{\mbox{\tiny QGP}}=0.25$GeV. For the fourth order
susceptibility, $\chi^{\mbox{\tiny H,latt}}_{4}/\chi^{\mbox{\tiny
H,latt}}_{2}=1$ at hadronic phase and $\chi^{\mbox{\tiny
QGP,latt}}_{4}/\chi^{\mbox{\tiny QGP,latt}}_{2}=2/(3\pi^{2})$ from the lattice
simulation [53, 54]. We assume the third-order susceptibility has the relation
$\chi^{\mbox{\tiny latt}}_{3}/\chi^{\mbox{\tiny latt}}_{2}=\chi^{\mbox{\tiny
latt}}_{4}/\chi^{\mbox{\tiny latt}}_{2}$ from Refs. [55, 56]. Finally, we have
the regular susceptibilities both for the hadronic and QGP phases:

$\displaystyle\chi^{\mbox{\tiny H}}_{2}$ $\displaystyle=T^{2}_{\mbox{\tiny
H}}\chi^{\mbox{\tiny H,latt}}_{2},\quad\chi^{\mbox{\tiny
QGP}}_{2}=T^{2}_{\mbox{\tiny QGP}}\chi^{\mbox{\tiny QGP,latt}}_{2},$
$\displaystyle\chi^{\mbox{\tiny H}}_{3}$ $\displaystyle=T_{\mbox{\tiny
H}}\chi^{\mbox{\tiny H,latt}}_{3},\quad\chi^{\mbox{\tiny
QGP}}_{3}=T_{\mbox{\tiny QGP}}\chi^{\mbox{\tiny QGP,latt}}_{3},$
$\displaystyle\chi^{\mbox{\tiny H}}_{4}$ $\displaystyle=\chi^{\mbox{\tiny
H,latt}}_{4},\qquad\chi^{\mbox{\tiny QGP}}_{4}=\chi^{\mbox{\tiny
QGP,latt}}_{4}.$

The critical contributions $\chi^{\mbox{\tiny cri}}_{2},\lambda^{\mbox{\tiny
cri}}_{3}$ and $\lambda_{4}^{\mbox{\tiny cri}}$ are constructed from the
cumulants of Ising model $\kappa_{2}$

$\displaystyle\chi^{\mbox{\tiny cri}}_{2}$
$\displaystyle=T^{2}_{A}C_{c}\kappa_{2}\equiv\xi^{2},\quad\lambda^{\mbox{\tiny
cri}}_{3}=\tilde{\lambda}_{3}T(T\xi)^{-3/2}/T^{3}_{A},$
$\displaystyle\lambda_{4}^{\mbox{\tiny cri}}$
$\displaystyle=\tilde{\lambda}_{4}(T\xi)^{-1}/T^{6}_{A},$ (10)

where $C_{c}$ is mapping constant and chosen as $C_{c}=4$ in this work.
$T_{A}$ is the constant for the dimensional consistency and is set as
$T_{A}=0.5$GeV. The dimensionaless coupling constants $\tilde{\lambda}_{3}$
and $\tilde{\lambda}_{4}$ have values range (0,8) and (4,20), respectively
[57]. $\tilde{\lambda}_{3}=1$ and $\tilde{\lambda}_{4}=4$ are used in this
study. The second-order cumulant of the Ising model reads [20, 25]

$\displaystyle\kappa_{2}=\frac{M_{0}}{H_{0}}\frac{1}{R^{4/3}(3+2\theta^{2})},$
(11)

where the normalization constants are $M_{0}\simeq 0.605,h_{0}\simeq 0.394$.
On the phase diagram near the critical point, the distance $R$ and angle
$\theta$ to the critical point are calculated with the equations

$\displaystyle r(R,\theta)=R(1-\theta^{2}),\quad
h(R,\theta)=R^{5/2}(3\theta-2\theta^{3}).$ (12)

The Ising variables ($r,h$) are connected to the temperature and chemical
potential of the QCD system by the mapping:

$\displaystyle\frac{T-T_{c}}{\Delta T}=\frac{h}{\Delta
h},\quad\frac{\mu-\mu_{c}}{\Delta\mu}=-\frac{r}{\Delta r},$ (13)

where $T_{c}$ and $\mu_{c}$ are the critical temperature and chemical
potential, respectively. $\Delta T$ and $\Delta\mu$ represent the widths of
the critical region of the QCD phase diagram, $\Delta h$ and $\Delta r$ are
the ones in the Ising model. These are non-universal parameters and we use
$\Delta T=T_{c}/8,\Delta\mu=0.1$GeV,$\Delta r=(5/3)^{3/4},\Delta h=1$ in this
work. The temperature $T(\tau,\eta)$ and chemical potential $\mu(\tau,\eta)$
of the QGP profile are borrowed from hydrodynamic simulation, which is
addressed in Sec. II

##  Appendix B Numerical details of Eq.(1)

This appendix presents the details of the numerical simulation of Eq.(1) and
its verification by comparison with the analytical calculation in the linear
limit of Eq.(1).

The explicit form of Eq.(1) reads

$\displaystyle\partial_{\tau}n_{B}$
$\displaystyle=D_{\eta}\chi_{\eta}\frac{\partial^{2}}{\partial\eta^{2}}\biggl{[}\frac{n_{B}}{\chi_{\eta}}-\frac{K}{\tau^{3}}\frac{\partial^{2}}{\partial\eta^{2}}n_{B}+\frac{\lambda_{3}}{\tau^{2}}n_{B}^{2}+\frac{\lambda_{4}}{\tau^{3}}n_{B}^{3}\biggr{]}$
$\displaystyle\qquad-\partial_{\eta}\zeta(\tau,\eta),$ (14)

with the noise correlator

$\displaystyle\langle\zeta(\tau,\eta)\zeta(\tau^{\prime},\eta^{\prime})\rangle=2D_{\eta}\chi_{\eta}\delta(\tau-\tau^{\prime})\delta(\eta-\eta^{\prime}),$
(15)

where $D_{\eta}=D/\tau^{2},\chi_{\eta}=\chi_{2}\tau T$.

Figure 9: Time evolution of second-order net-baryon multiplicity fluctuations
$\sigma^{2}$ as a solution of the dynamics of the net-baryon density without
higher order terms (17). The green dashed curve of $\sigma^{2}$ is obtained
from the simulation with the Saul’yev scheme (B) and the red one corresponds
to the analytical calculation of Eq. (B). $\Delta\eta=0.25$ here.

One of the most popular algorithms for this diffusive stochastic equation is
the explicit Forward or Backward Euler Method. In this numerical algorithm,
Eq.(1) can be discretized in an explicit form:
$(n^{i+1}_{B,j}-n^{i}_{B,j})/\Delta\tau=f(n^{i}_{B,j})$, where the next time
step $n^{i+1}_{B,j}$ only shown in the left-hand side and can be obtained
explicitly. However, this scheme is conditionally stable because it is only an
approximation to Eq.(1), error will gradually accumulate and eventually lead
to instability. The situation is exacerbated as the system approaches the
critical point. To attain stable solutions, the temporal step should be
significantly smaller than the spatial step, such as
$\Delta\tau\leq\Delta\eta^{2}/(2D)$ for diffusion equation $\partial\tau
n_{B}=D\partial^{2}_{\eta}n_{B}$, namely the conditional stable. This makes
the numerical simulation extremely inefficient.

This work implements the numerical simulation with Saul’yev scheme [46, 47],
where this equation can be discretized as

$\displaystyle n^{i+1}_{j}$
$\displaystyle\approx\frac{1}{r}\biggl{\\{}\biggl{(}1-D_{\eta}\frac{\delta\tau}{\delta\eta^{2}}\biggr{)}n^{i}_{j}+D_{\eta}\frac{\delta\tau}{\delta\eta^{2}}\biggl{(}n^{i}_{j+1}+n^{i+1}_{j-1}\biggr{)}$
$\displaystyle\quad+D_{\eta}\chi_{\eta}\frac{\delta\tau}{\delta\eta^{2}}\biggl{(}S[n^{i}_{j+1}]+S[n^{i}_{j-1}]-2S[n^{i}_{j}]\biggr{)}$
$\displaystyle\quad-
D_{\eta}\chi_{\eta}\frac{K}{\tau_{i}^{3}}\frac{\delta\tau}{\delta\eta^{4}}\biggl{(}n^{i}_{j+2}-4n^{i}_{j+1}-4n^{i+1}_{j-1}+n^{i+1}_{j-2}\biggr{)}$
$\displaystyle\quad+\sqrt{\frac{2D_{\eta}\chi_{\eta}}{A\delta\tau\delta\eta}}\frac{\delta\tau}{\delta\eta}\biggl{(}W^{i}_{j+1}-W^{i}_{j}\biggr{)}\biggr{\\}},$
(16)

where
$r=(1+D_{\eta}\frac{\delta\tau}{\delta\eta^{2}}+6D_{\eta}\chi_{\eta}\frac{K}{\tau^{3}}\frac{\delta\tau}{\delta\eta^{4}})$
and $S[n^{i}_{j}]$ denotes the discretization of the higher order terms
$\frac{\lambda_{3}}{\tau^{2}}n_{B}^{2}+\frac{\lambda_{4}}{\tau^{3}}n_{B}^{3}$.
The noise term is discretized into $W^{i}_{j}$, corresponding to the Gaussian
white noise with unit variance. $n^{i}_{j}$ is the net-baryon density at
proper time $\tau_{i}=\tau_{0}+i\cdot\delta\tau$ and rapidity
$\eta_{j}=-L/2+j\cdot\delta\eta$. In this simulation, the increment in the
proper time is chosen as $\delta\tau=0.0018$fm and the spacing of the rapidity
is $\delta\eta=0.115$. The grid size of the simulation is $L=\delta\eta\cdot
N$ with $N=128$. The initial proper time is $\tau_{0}=3$fm with the initial
condition: $\langle n_{B}(\tau_{0},\eta)\rangle=0,\,\langle
n_{B}(\tau_{0},\eta)n_{B}(\tau_{0},\eta^{\prime})\rangle=\chi_{\eta}(\tau_{0})\delta(\eta-\eta^{\prime})$.

As shown in Eq.(B), the net-baryon density for the next temporal step $i+1$ at
grid site $j$ only depends on the ones of time step $i$, expect
$n^{i+1}_{j-1}$ and $n^{i+1}_{j-2}$ at right hand side of Eq.(B). Therefore,
$n^{i+1}_{j}$ can be computed explicitly, with the boundary condition of the
cell $j=1$ given, and the simulation starts the cell from the left to the
right $j=1,\cdots,N$. On the other hand, the freeze-out hyper-surface of the
QGP profile naturally provides the boundary of the simulation. At the edge of
this hypersurface, we employ the boundary condition as follows: evolving the
net-baryon density with Eq.(B), but replacing the cell outside
(e.g.,$n^{i}_{j+1}$) with the closest boundary (e.g.,$n^{i}_{j}$) if the cell
$i$ fall at the right side of the profile.

To verify the implemented numerical scheme Eq.(B), the dynamics of the net-
baryon density without the higher order terms are investigated:

$\displaystyle\partial_{\tau}n_{B}$
$\displaystyle=D_{\eta}\frac{\partial^{2}}{\partial\eta^{2}}n_{B}-\partial_{\eta}\zeta(\tau,\eta),$
(17)

with noise (15). This is the stochastic diffusion equation employed in the
previous works [25, 58] and the corresponding analytical solution of second-
order baryon multiplicity fluctuations reads

$\displaystyle\sigma^{2}$
$\displaystyle=\Delta\eta\chi_{\eta}(\tau)-\Delta\eta\int^{\tau}_{\tau_{0}}d\tau_{1}\chi^{\prime}_{\eta}(\tau_{1})\biggl{\\{}\mbox{erf}\biggl{[}\frac{1}{\bar{D}(\tau_{1},\tau,\Delta\eta)}\biggr{]}$
$\displaystyle\qquad+\frac{\bar{D}(\tau_{1},\tau,\Delta\eta)}{\sqrt{\pi}}\biggl{[}\exp\biggl{(}-\frac{1}{\bar{D}(\tau_{1},\tau,\Delta\eta)^{2}}\biggr{)}-1\biggr{]}\biggr{\\}},$
(18)

where
$\bar{D}(\tau_{1},\tau,\Delta\eta)=\Delta\eta/\sqrt{8\int^{\tau}_{\tau_{1}}d\tau^{\prime}D_{\eta}(\tau^{\prime})}$
and $\chi^{\prime}_{\eta}(\tau_{1})$ is the derivative of $\chi_{\eta}$ over
the proper time.

It requires the time evolution of $\chi_{\eta}$ and $D_{\eta}$ for the
comparison of analytical second order multiplicity Eq. (B) with the numerical
results of Eq. (17) with the algorithm Eq.(B). This can be constructed
following the method of Eq. (A) in Append A. Here, the time evolution of
temperature $T(\tau,\eta)$ is not obtained from QGP profile, but we employ the
Hubble-like expansion [25, 58], $T(\tau)=T_{I}(\tau/\tau_{0})^{-c^{2}_{s}}$
and fixed chemical potential $\mu=0.13$GeV, with $T_{I}=0.325$GeV and
$c^{2}_{s}=1$. In this calculation, temperature and chemical potential are
constant across the rapidity space. Fig.9 shows the comparison of the time
evolution of $\sigma^{2}$ obtained both from numerical simulation and
analytical calculation. The agreement between these two methods verifies the
reliability of the Saul’yev scheme in the simulation of Eq. (1).

## References

  * Aoki  _et al._ [2006] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature 443, 675 (2006), arXiv:hep-lat/0611014 .
  * Ding  _et al._ [2015] H.-T. Ding, F. Karsch, and S. Mukherjee, Int. J. Mod. Phys. E 24, 1530007 (2015), arXiv:1504.05274 [hep-lat] .
  * Bazavov  _et al._ [2019] A. Bazavov, F. Karsch, S. Mukherjee, and P. Petreczky (USQCD), Eur. Phys. J. A 55, 194 (2019), arXiv:1904.09951 [hep-lat] .
  * Ratti [2018] C. Ratti, Rept. Prog. Phys. 81, 084301 (2018), arXiv:1804.07810 [hep-lat] .
  * Fischer [2019] C. S. Fischer, Prog. Part. Nucl. Phys. 105, 1 (2019), arXiv:1810.12938 [hep-ph] .
  * Fukushima and Hatsuda [2011] K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Fukushima and Sasaki [2013] K. Fukushima and C. Sasaki, Prog. Part. Nucl. Phys. 72, 99 (2013), arXiv:1301.6377 [hep-ph] .
  * Fu [2022] W.-j. Fu, Commun. Theor. Phys. 74, 097304 (2022), arXiv:2205.00468 [hep-ph] .
  * Stephanov [2011] M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
  * Athanasiou  _et al._ [2010] C. Athanasiou, K. Rajagopal, and M. Stephanov, Phys. Rev. D 82, 074008 (2010), arXiv:1006.4636 [hep-ph] .
  * Adam  _et al._ [2021] J. Adam _et al._ (STAR), Phys. Rev. Lett. 126, 092301 (2021), arXiv:2001.02852 [nucl-ex] .
  * Abdallah  _et al._ [2021] M. Abdallah _et al._ (STAR), Phys. Rev. C 104, 024902 (2021), arXiv:2101.12413 [nucl-ex] .
  * Asakawa and Kitazawa [2016] M. Asakawa and M. Kitazawa, Prog. Part. Nucl. Phys. 90, 299 (2016), arXiv:1512.05038 [nucl-th] .
  * Bzdak  _et al._ [2020] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * Bluhm  _et al._ [2020] M. Bluhm _et al._ , Nucl. Phys. A 1003, 122016 (2020), arXiv:2001.08831 [nucl-th] .
  * Wu  _et al._ [2021] S. Wu, C. Shen, and H. Song, Chin. Phys. Lett. 38, 081201 (2021), arXiv:2104.13250 [nucl-th] .
  * An  _et al._ [2022] X. An _et al._ , Nucl. Phys. A 1017, 122343 (2022), arXiv:2108.13867 [nucl-th] .
  * Du  _et al._ [2024a] L. Du, A. Sorensen, and M. Stephanov (2024) arXiv:2402.10183 [nucl-th] .
  * Berdnikov and Rajagopal [2000] B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000), arXiv:hep-ph/9912274 .
  * Nonaka and Asakawa [2005] C. Nonaka and M. Asakawa, Phys. Rev. C 71, 044904 (2005), arXiv:nucl-th/0410078 .
  * Mukherjee  _et al._ [2015] S. Mukherjee, R. Venugopalan, and Y. Yin, Phys. Rev. C 92, 034912 (2015), arXiv:1506.00645 [hep-ph] .
  * Tang  _et al._ [2023] S. Tang, S. Wu, and H. Song, Phys. Rev. C 108, 034901 (2023), arXiv:2303.15017 [nucl-th] .
  * Nahrgang  _et al._ [2019] M. Nahrgang, M. Bluhm, T. Schaefer, and S. A. Bass, Phys. Rev. D 99, 116015 (2019), arXiv:1804.05728 [nucl-th] .
  * Nahrgang and Bluhm [2020] M. Nahrgang and M. Bluhm, Phys. Rev. D 102, 094017 (2020), arXiv:2007.10371 [nucl-th] .
  * Sakaida  _et al._ [2017] M. Sakaida, M. Asakawa, H. Fujii, and M. Kitazawa, Phys. Rev. C 95, 064905 (2017), arXiv:1703.08008 [nucl-th] .
  * Pihan  _et al._ [2023] G. Pihan, M. Bluhm, M. Kitazawa, T. Sami, and M. Nahrgang, Phys. Rev. C 107, 014908 (2023), arXiv:2205.12834 [nucl-th] .
  * Nahrgang  _et al._ [2011] M. Nahrgang, S. Leupold, C. Herold, and M. Bleicher, Phys. Rev. C 84, 024912 (2011), arXiv:1105.0622 [nucl-th] .
  * Nahrgang  _et al._ [2013] M. Nahrgang, C. Herold, S. Leupold, I. Mishustin, and M. Bleicher, J. Phys. G 40, 055108 (2013), arXiv:1105.1962 [nucl-th] .
  * Herold  _et al._ [2014] C. Herold, M. Nahrgang, Y. Yan, and C. Kobdaj, J. Phys. G 41, 115106 (2014), arXiv:1407.8277 [hep-ph] .
  * Kapusta  _et al._ [2012] J. I. Kapusta, B. Muller, and M. Stephanov, Phys. Rev. C 85, 054906 (2012), arXiv:1112.6405 [nucl-th] .
  * An  _et al._ [2020] X. An, G. Başar, M. Stephanov, and H.-U. Yee, Phys. Rev. C 102, 034901 (2020), arXiv:1912.13456 [hep-th] .
  * An  _et al._ [2021] X. An, G. Başar, M. Stephanov, and H.-U. Yee, Phys. Rev. Lett. 127, 072301 (2021), arXiv:2009.10742 [hep-th] .
  * Stephanov and Yin [2018] M. Stephanov and Y. Yin, Phys. Rev. D 98, 036006 (2018), arXiv:1712.10305 [nucl-th] .
  * Shen and Schenke [2018] C. Shen and B. Schenke, Phys. Rev. C 97, 024907 (2018), arXiv:1710.00881 [nucl-th] .
  * Du  _et al._ [2024b] L. Du, H. Gao, S. Jeon, and C. Gale, Phys. Rev. C 109, 014907 (2024b), arXiv:2302.13852 [nucl-th] .
  * Du  _et al._ [2021] L. Du, X. An, and U. Heinz, Phys. Rev. C 104, 064904 (2021), arXiv:2107.02302 [hep-ph] .
  * Brewer  _et al._ [2018] J. Brewer, S. Mukherjee, K. Rajagopal, and Y. Yin, Phys. Rev. C 98, 061901 (2018), arXiv:1804.10215 [hep-ph] .
  * Son and Stephanov [2004] D. T. Son and M. A. Stephanov, Phys. Rev. D 70, 056001 (2004), arXiv:hep-ph/0401052 .
  * Ling  _et al._ [2014] B. Ling, T. Springer, and M. Stephanov, Phys. Rev. C 89, 064901 (2014), arXiv:1310.6036 [nucl-th] .
  * Chao and Schaefer [2021] J. Chao and T. Schaefer, JHEP 01, 071, arXiv:2008.01269 [hep-th] .
  * Chao and Schaefer [2023] J. Chao and T. Schaefer, JHEP 06, 057, arXiv:2302.00720 [hep-ph] .
  * Hu [2024] J. Hu, arXiv e-prints , arXiv:2403.15825 (2024), arXiv:2403.15825 [cond-mat.stat-mech] .
  * Shen  _et al._ [2016] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, Comput. Phys. Commun. 199, 61 (2016), arXiv:1409.8164 [nucl-th] .
  * Lin  _et al._ [2005] Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005), arXiv:nucl-th/0411110 .
  * Parotto  _et al._ [2020] P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, and M. Stephanov, Phys. Rev. C 101, 034901 (2020), arXiv:1805.05249 [hep-ph] .
  * Saul’yev [1957] V. K. Saul’yev, Dokl. Akad. Nauk SSSR 115, 1077 (1957).
  * Saul’yev [1964] V. K. Saul’yev,  _Integration of Equations of Parabolic Type by the Method of Nets_ (Pergamon Press, Oxford, 1964).
  * Yang  _et al._ [2022] J. Yang, Y. Li, C. Lee, H. G. Lee, S. Kwak, Y. Hwang, X. Xin, and J. Kim, International Journal of Mechanical Sciences 217, 106985 (2022).
  * Cooper and Frye [1974] F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
  * Kapusta and Torres-Rincon [2012] J. I. Kapusta and J. M. Torres-Rincon, Phys. Rev. C 86, 054911 (2012), arXiv:1209.0675 [nucl-th] .
  * Note [1] Noted that the integral interval of $\sigma^{2}$ in Figs.4 is from $-\Delta\eta$ to $\Delta\eta$. The difference between Scenario I and II is small when the integration is also performed from $-\Delta\eta$ to $\Delta\eta$ for the time evolution of $\sigma^{2}/\Delta\eta$ in Figs.5. This is consistent with Figs.4 and the difference between Scenario I and II in different rapidity is smeared out from $-\Delta\eta$ to $\Delta\eta$.
  * Pradeep and Stephanov [2023] M. S. Pradeep and M. Stephanov, Phys. Rev. Lett. 130, 162301 (2023), arXiv:2211.09142 [hep-ph] .
  * Cheng  _et al._ [2009] M. Cheng _et al._ , Phys. Rev. D 79, 074505 (2009), arXiv:0811.1006 [hep-lat] .
  * Bazavov  _et al._ [2017] A. Bazavov _et al._ , Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
  * Motornenko  _et al._ [2020] A. Motornenko, J. Steinheimer, V. Vovchenko, S. Schramm, and H. Stoecker, Phys. Rev. C 101, 034904 (2020), arXiv:1905.00866 [hep-ph] .
  * Mukherjee  _et al._ [2017] A. Mukherjee, J. Steinheimer, and S. Schramm, Phys. Rev. C 96, 025205 (2017), arXiv:1611.10144 [nucl-th] .
  * Tsypin [1994] M. M. Tsypin, Phys. Rev. Lett. 73, 2015 (1994).
  * Wu and Song [2019] S. Wu and H. Song, Chin. Phys. C 43, 084103 (2019), arXiv:1903.06075 [nucl-th] .