# Spinodal decomposition in rapidly expanding quark gluon plasma

Joseph I. Kapusta<EMAIL_ADDRESS>School of Physics and Astronomy, University
of Minnesota, Minneapolis, MN 55455, USA  Mayank Singh
<EMAIL_ADDRESS>School of Physics and Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  Department of Physics and Astronomy,
Vanderbilt University, Nashville, TN 37240, USA  Thomas Welle
<EMAIL_ADDRESS>School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  Applied Research Associates, 8537 Six Forks Road,
Raleigh, NC 27615, USA

###### Abstract

Quantum Chromodynamics (QCD) is expected to have a first order phase
transition between the confined hadron gas and the deconfined quark gluon
plasma at high baryon densities. This will result in phase boundary effects in
the metastable and unstable regions. It is important to include these effects
in phenomenological models of heavy ion collisions to identify experimental
signatures of a phase transition. This requires building intuition on phase
separation in rapidly expanding fluids. In this work we present the equations
of relativistic hydrodynamics with a phase boundary, provide prescriptions to
extend the equation of state to metastable and unstable regions, and show the
effects of spinodal separation in a Bjorken flow.

††preprint: APS/123-QED

##  I Introduction

Full mapping of the QCD phase diagram has been a long standing challenge since
the discovery of QCD five decades ago. The Beam Energy Scan (BES) program at
the Relativistic Heavy Ion Collider (RHIC) recorded collisions at different
collision energies to explore the QCD phase diagram at a wide range of
temperatures and baryon chemical potentials [1, 2]. Inferring the properties
of the phase diagram from the BES data would require careful comparison with
the phenomenological models which include all the relevant physics.

Lattice QCD calculations have provided ab initio calculations on the zero
baryon chemical potential axis of the phase diagram [3, 4, 5]. Unfortunately,
lattice QCD cannot be extended to large chemical potentials due to the
infamous sign problem. Over the last decade, progress has been made in
obtaining the QCD equation of state (EOS) at small baryon chemical potentials
$\mu_{B}$ by expanding the QCD partition function in powers of $\mu_{B}/T$,
where $T$ is the temperature [6, 7, 8, 9, 10, 11, 12]. This has been
complemented by an EOS estimated from perturbative QCD which can be matched to
lattice results at $\mu_{B}=0$ [13]. Recent developments in this area can be
found in Ref. [14].

The QCD phase diagram is expected to have a first order phase transition at
low temperatures and large baryon densities [15, 16, 17]. Lattice QCD with
physical quark masses has shown that there is no true phase transition at zero
baryon chemical potential, and that the phase change from a confined hadronic
gas to a deconfined quark gluon plasma (QGP) is accomplished by a smooth
crossover [3]. Hence, the first order phase transition line is expected to end
somewhere before reaching the zero chemical potential axis. This endpoint is
expected to be a second order phase transition critical point in the same
universality class as the liquid–gas phase transition and the 3D Ising model.
While the exact location of the critical point is not known, once a position
is assumed there are different approaches to embed a critical point obeying
the right universality properties on a background EOS [18, 19, 20, 21, 22].
Pinpointing the location of critical point will have to rely on testing models
with different critical point locations by matching them to experimental data.

Successful isolation of experimental signatures of a critical point and the
first order phase transition require implementation of critical fluctuations
and the phase boundary effects into the heavy ion simulations. These have
proven to be incredibly challenging. There have been some efforts to include
the fluctuations in the simulations [23, 24, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34, 35], but a realistic model including critical fluctuations, which can
isolate the signatures of critical point has remained elusive. For a summary
of these approaches, see Refs. [36, 37].

Simulations have shown that even individual BES collisions span a wide region
within the phase diagram and the volume of matter whose trajectory is in the
neighbourhood of the critical point is not necessarily very large [38, 39,
40]. This makes it even more challenging to get clean experimental signature
of such a critical region. What may be more realistic to achieve is to detect
effects of phase boundary created by a first order phase transition. As this
is a wide region on the phase diagram, it is much more plausible that a
collision of suitable energy will have a big chunk of its volume passing
through the phase transition line. Of course this brings in its own
challenges. The phase transition can proceed by two paths: nucleation or
spinodal decomposition. Once the system supercools nucleation is the first
process to initiate the phase change [41]. Generally this will happen in the
metastable region. However, if the nucleation rate is small relative to the
expansion rate then the system will enter the unstable region and spinodal
decomposition will occur [42]. Spinodal decomposition may be initiated by
thermal fluctuations or by fluctuations already present in the system. Most
simulations assume the latter and proceed to solve the Cahn-Hilliard equation.
Obviously there are initial state fluctuations present in heavy ion
collisions. While the first order phase transition has been extensively
studied in chemical and condensed matter systems (for example, see the
textbook [43]), the relativistic collision systems bring in an additional
problem in that it is rapidly expanding. To the best of our knowledge, the
effects of phase separation in a rapidly expanding fluid have not been
explored.

As a first step towards this goal, following the pioneering work of Ref. [44,
45], we write the equations of spinodal decomposition in a rapidly expanding
fluid. We provide a prescription to extend the EOS to the metastable and
unstable regions, and show the results of numerical simulations of a 1D
rapidly expanding QCD medium undergoing the first order phase transition. The
thermodynamics of first order phase transition is discussed in Sec. II. We
provide prescriptions to extend the equation of state to metastable regions in
Sec. III and write the equations of hydrodynamics with a phase boundary in
Sec. IV. The equations for the case of Bjorken flow are given in Sec. V and we
show the effects of first order phase transition on Bjorken fluid in Sec. VI.
We discuss our findings in Sec. VII.

##  II Thermodynamics of a phase transition

Consider the Helmholtz free energy functional at temperature $T$

$F\\{n({\bf x},t)\\}=\int
d^{3}x\left[{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}+f(T,n)\right]=\int
d^{3}x\tilde{f}(T,n).$ (1)

The function $f(T,n)$ is the bulk free energy density while $\tilde{f}(T,n)$
includes the gradient term in a Ginzburg-Landau free energy functional. The
gradient term represents the cost in energy in a spatially non-uniform system.
When there are two phases, the value of $K$ is related to the surface energy
as described below. When the density changes by a small amount $\delta n$ the
free energy change to first order is

$\displaystyle F\\{n+\delta n\\}-F\\{n\\}$ $\displaystyle=$ $\displaystyle\int
d^{3}x\left[K\bm{\nabla}\delta n\cdot\bm{\nabla}n+\frac{\partial f}{\partial
n}\delta n\right]$ (2) $\displaystyle=$ $\displaystyle\int
d^{3}x\left[K\bm{\nabla}\cdot(\delta n\bm{\nabla}n)-K\delta
n\nabla^{2}n+\frac{\partial f}{\partial n}\delta n\right].$

The total derivative is converted to a surface integral which vanishes with
appropriate boundary conditions. The local chemical potential is the
functional derivative

$\tilde{\mu}=\frac{\delta F}{\delta n}=\frac{\partial f}{\partial
n}-K\nabla^{2}n=\mu-K\nabla^{2}n.$ (3)

The local isotropic pressure is

$\displaystyle\tilde{P}$ $\displaystyle=$ $\displaystyle
n\mu-[f+{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}]$ (4) $\displaystyle=$
$\displaystyle n\frac{\partial f}{\partial
n}-f-Kn\nabla^{2}n-{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}$
$\displaystyle=$ $\displaystyle
P-Kn\nabla^{2}n-{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}.$

Assuming negligible dependence of $K$ on temperature the entropy density is as
usual

$s=-\frac{\partial f}{\partial T},$ (5)

and the local energy density is

$\tilde{\epsilon}=f-T\frac{\partial f}{\partial
T}+{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}=\epsilon+{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}.$
(6)

Finally the enthalpy density is

$\displaystyle\tilde{w}$ $\displaystyle=$
$\displaystyle\tilde{\epsilon}+\tilde{P}=n\frac{\partial f}{\partial
n}-T\frac{\partial f}{\partial T}-Kn\nabla^{2}n$ (7) $\displaystyle=$
$\displaystyle\tilde{\mu}n+Ts=w-Kn\nabla^{2}n.$

This is all consistent with the work of Cahn and Hilliard more than 60 years
ago [46, 47] as well as others, notably by Yang, Fleming and Gibbs [48].

To determine the static density distribution we should minimize the Helmholtz
free energy subject to the condition that the number of particles is fixed.
This is done by introducing a Lagrange multiplier $\lambda$ and minimizing the
integral

$I=\int d^{3}x\left[\tilde{f}(T,n)-\lambda n\right].$ (8)

The resulting Euler-Lagrange equation is

$\partial_{i}\frac{\partial(\tilde{f}-\lambda
n)}{\partial(\partial_{i}n)}-\frac{\partial(\tilde{f}-\lambda n)}{\partial
n}=0,$ (9)

which gives

$\lambda=\mu-K\nabla^{2}n.$ (10)

This means setting $\tilde{\mu}=\lambda$ = constant. If the system has two
phases in equilibrium with each other, and if there is a planar interface
between them, then the surface free energy $\sigma$ can be calculated from the
well-known formula

$\sigma=K\int_{-\infty}^{\infty}dx\left(\frac{dn}{dx}\right)^{2},$ (11)

where $n(x)$ is the solution to Eq. (10). In that case we can identify
$\lambda$ with the common chemical potential of each phase,
$\lambda=\mu_{L}=\mu_{G}$, along the coexistence curve.

##  III Metastable and Unstable Regions

First principles calculations of the equation of state will usually not
describe the metastable and unstable regions of the phase diagram. Instead,
for $T<T_{c}$ and $n_{G}\leq n\leq n_{L}$ we parameterize, or interpolate, the
pressure in terms of the density as

$P_{{\rm int}}(n)=P_{X}(T)+\sum_{i=1}^{4}c_{i}(n-n_{G})^{i}.$ (12)

The subscript $X$ refers to any point along the line of phase coexistence. The
condition $P_{{\rm int}}(n_{L})=P_{X}(T)$ places one constraint on the four
$c_{i}$

$\sum_{i=1}^{4}c_{i}\Delta n^{i}=0,$ (13)

where $\Delta n\equiv n_{L}-n_{G}\geq 0$, while continuity of the first
derivative at the endpoints provide two more conditions

$\displaystyle\frac{\partial P}{\partial n}(n_{G})$ $\displaystyle=$
$\displaystyle c_{1}=\frac{n_{G}}{\chi_{G}},$ $\displaystyle\frac{\partial
P}{\partial n}(n_{L})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{4}ic_{i}\Delta n^{i-1}=\frac{n_{L}}{\chi_{L}}.$ (14)

Here $\chi_{G}$ and $\chi_{L}$ are the baryon susceptibilities
$\chi=n/(\partial P/\partial n)_{T}$ in the gas and liquid phases. The
chemical potential is

$\mu_{{\rm
int}}(n)=\mu_{X}(T)+d_{0}\ln(n/n_{G})+\sum_{i=1}^{3}d_{i}(n-n_{G})^{i},$ (15)

where

$\displaystyle d_{0}$ $\displaystyle=$ $\displaystyle
c_{1}-2c_{2}n_{G}+3c_{3}n_{G}^{2}-4c_{4}n_{G}^{3},$ $\displaystyle d_{1}$
$\displaystyle=$ $\displaystyle 2c_{2}-3c_{3}n_{G}+4c_{4}n_{G}^{2},$
$\displaystyle d_{2}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{3}{2}}}c_{3}-2c_{4}n_{G},$ $\displaystyle
d_{3}$ $\displaystyle=$ $\displaystyle{\textstyle{\frac{4}{3}}}c_{4}.$ (16)

This automatically satisfies $\mu_{{\rm int}}(n_{G})=\mu_{X}(T)$. Equivalently

$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle d_{0}+d_{1}n_{G},$
$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{1}{2}}}d_{1}+d_{2}n_{G},$ $\displaystyle
c_{3}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{2}{3}}}d_{2}+d_{3}n_{G},$ $\displaystyle
c_{4}$ $\displaystyle=$ $\displaystyle{\textstyle{\frac{3}{4}}}d_{3}.$ (17)

The final condition $\mu_{{\rm int}}(n_{L})=\mu_{X}(T)$ is satisfied when

$d_{0}\ln(n_{L}/n_{G})+\sum_{i=1}^{3}d_{i}\Delta n^{i}=0.$ (18)

Then all coefficients $c_{i}$ and $d_{i}$ are determined.

To test the veracity of this parameterization we compare to the Van der Waals
equation of state. Its dimensionless form, representing the law of
corresponding states, is

$P_{*}=\frac{8T_{*}}{3v_{*}-1}-\frac{3}{v_{*}^{2}}.$ (19)

Here $P_{*}=P/P_{c}$, $v_{*}=v/v_{c}$, $T_{*}=T/T_{c}$, and $v=1/n$. Some
results are shown in Fig. 1.

Figure 1: The parameterization of the metastable and unstable regions, as
described in the text, compared to the Van der Waals equation of state. The
comparison is done for $T_{*}$ = 0.95, 0.90, 0.85, and 0.80. The plot shows
the volume per particle even though the parameterization was done in terms of
the density $n_{*}=1/v_{*}$ to illustrate the Maxwell construction.

Alternatively, the pressure can be expanded in terms of the volume per baryon
$v=1/n$. We do not use this construction as we found that the agreement with
the Van der Waals equation of state is not nearly as good. Note that the
pressure goes negative in the metastable and unstable regions at low
temperature. This can be handled by the theory presented herein.

For illustration, in this paper we use the QCD equation of state from Ref.
[21]. The critical point is embedded on a smooth background equation of state
such that the right critical exponents from the 3D Ising model universality
class are obtained. The background equation of state is obtained by matching a
perturbative-QCD EOS in the quark gluon phase with a hadron resonance gas EOS
in the hadronic phase with a smooth matching function [13]. The critical point
is chosen to be at $T=130$ MeV and $\mu_{B}=450$ MeV. The resulting
interpolated chemical potential and pressure at a fixed temperature are shown
in Fig. 2.

Figure 2: The chemical potential (left) and pressure (right) as a function of
baryon density at temperature of 100 MeV.

It must be pointed out that the interpolation given above is relevant for an
equation of state with mean field critical exponents. The QCD equation of
state employed here uses the known critical exponents from the 3D Ising model
universality class. It should not make much difference except when one is in
the neighborhood of the critical point. In this work we stay away from that
neighborhood, focusing instead on crossing the line of first order phase
transition during a heavy ion collision. Generalization of the polynomial
interpolation is presented in the appendix.

##  IV Hydrodynamics: General

The stress tensor was derived in Ref. [48] to be

$T_{ij}=\tilde{P}\delta_{ij}+K(\partial_{i}n)(\partial_{j}n),$ (20)

which is symmetric. This is in a stationary fluid or gas. A fully relativistic
expression for the stress–energy–momentum tensor is

$T^{\mu\nu}=\tilde{P}(u^{\mu}u^{\nu}-g^{\mu\nu})+\tilde{\epsilon}u^{\mu}u^{\nu}+K(D^{\mu}n)(D^{\nu}n),$
(21)

where $u^{\mu}$ is the local flow velocity, the metric is $(+---)$, and

$D^{\mu}n\equiv\partial^{\mu}n-u^{\mu}u^{\alpha}\partial_{\alpha}n$ (22)

is a gradient orthogonal to the flow velocity. This $T^{\mu\nu}$ is obviously
a symmetric second rank tensor. In the local rest frame

$\displaystyle T^{00}$ $\displaystyle=$ $\displaystyle\tilde{\epsilon},$
$\displaystyle T^{0j}$ $\displaystyle=$ $\displaystyle 0,$ $\displaystyle
T^{ij}$ $\displaystyle=$
$\displaystyle\tilde{P}\delta_{ij}+K(\partial_{i}n)(\partial_{j}n).$ (23)

As a check, consider the equation of motion $\partial_{\mu}T^{\mu\nu}=0$ for a
planar interface between the liquid and gas phases with $n({\bf x},t)=n(x)$.
The only nontrivial equation is

$\frac{d}{dx}\left[\tilde{P}+K\left(\frac{dn}{dx}\right)^{2}\right]=0.$ (24)

Substituting for $\tilde{P}$ this can be written as

$n\tilde{\mu}-f+{\textstyle{\frac{1}{2}}}K\left(\frac{dn}{dx}\right)^{2}={\rm
constant}.$ (25)

This is easily shown [48] to be equivalent to Eq. (10).

The local thermodynanic variables can be written covariantly using $D^{\mu}$.

$\displaystyle\tilde{\mu}$ $\displaystyle=$ $\displaystyle\mu+KD^{2}n,$
$\displaystyle\tilde{P}$ $\displaystyle=$ $\displaystyle
P+KnD^{2}n+{\textstyle{\frac{1}{2}}}K(D^{\mu}n)(D_{\mu}n),$
$\displaystyle\tilde{\epsilon}$ $\displaystyle=$
$\displaystyle\epsilon-{\textstyle{\frac{1}{2}}}K(D^{\mu}n)(D_{\mu}n),$
$\displaystyle\tilde{w}$ $\displaystyle=$
$\displaystyle\tilde{P}+\tilde{\epsilon}=Ts+\tilde{\mu}n=w+KnD^{2}n.$ (26)

##  V Hydrodynamics: Bjorken

In boost-invariant hydrodynamics one expresses the time and location along the
beam direction in terms of the proper time $\tau$ and space-time rapidity
$\xi$ as

$\displaystyle t$ $\displaystyle=$ $\displaystyle\tau\cosh\xi,$ $\displaystyle
z$ $\displaystyle=$ $\displaystyle\tau\sinh\xi,$ (27)

with the inverse relations

$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\sqrt{t^{2}-z^{2}},$
$\displaystyle\xi$ $\displaystyle=$ $\displaystyle\tanh^{-1}(z/t).$ (28)

The flow velocity has the non-vanishing components

$\displaystyle u^{0}$ $\displaystyle=$ $\displaystyle\cosh\xi,$ $\displaystyle
u^{3}$ $\displaystyle=$ $\displaystyle\sinh\xi.$ (29)

It is useful to note that

$\displaystyle\frac{\partial}{\partial t}$ $\displaystyle=$
$\displaystyle\cosh\xi\frac{\partial}{\partial\tau}-\frac{\sinh\xi}{\tau}\frac{\partial}{\partial\xi},$
$\displaystyle\frac{\partial}{\partial z}$ $\displaystyle=$
$\displaystyle-\sinh\xi\frac{\partial}{\partial\tau}+\frac{\cosh\xi}{\tau}\frac{\partial}{\partial\xi},$
$\displaystyle\partial\cdot u$ $\displaystyle=$ $\displaystyle\frac{1}{\tau},$
$\displaystyle u\cdot\partial$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\tau},$ $\displaystyle D^{0}$
$\displaystyle=$
$\displaystyle-\frac{\sinh\xi}{\tau}\frac{\partial}{\partial\xi},$
$\displaystyle D^{3}$ $\displaystyle=$
$\displaystyle-\frac{\cosh\xi}{\tau}\frac{\partial}{\partial\xi}.$ (30)

Baryon number conservation $\partial_{\mu}(nu^{\mu})=0$ quickly leads to

$\frac{dn}{d\tau}+\frac{n}{\tau}=0.$ (31)

Due to boost invariance and the assumption of transverse homogeneity in the
Bjorken model, $D^{\mu}n=0$ so that energy and momentum conservation
$\partial_{\mu}T^{\mu\nu}=0$ lead to

$\frac{d\epsilon}{d\tau}+\frac{w}{\tau}=0.$ (32)

At fixed volume $d\epsilon=Tds+\mu dn$. Together with Eq. (31) this yields the
familiar equation for entropy conservation

$\frac{ds}{d\tau}+\frac{s}{\tau}=0.$ (33)

All effects of the spatial gradient terms with coefficient $K$ drop out due to
the assumption of boost invariance and transverse homogeneity.

Numerical studies of spinodal decomposition typically use a diffusion
equation. In the Landau-Lifshitz definition of flow velocity the baryon
current is

$J^{\mu}=nu^{\mu}+\sigma_{B}TD^{\mu}\left(\frac{\tilde{\mu}}{T}\right),$ (34)

where $\sigma_{B}$ is the baryon conductivity. This vanishes under the
assumption of boost invariance and transverse homogeneity. The Cattaneo
expression for the current is [49, 50, 51]

$J^{\mu}=nu^{\mu}+\sigma_{B}TD^{\mu}(1+\tau_{B}u\cdot\partial)^{-1}\left(\frac{\tilde{\mu}}{T}\right),$
(35)

where $\tau_{B}$ is the relaxation time. This expression vanishes under boost
invariance and transverse homogeneity as well. The Cattaneo expression is
causal under linear fluctuations while the diffusion equation suffers from
acausality even under small perturbations. Nevertheless, we use the diffusion
equation in this work for better correspondence with the non-relativistic
literature. The expressions can be changed to their causal version using the
Cattaneo form in a straightforward manner.

Since diffusion is probably necessary for spinodal decomposition, or phase
separation, we should allow the thermodynamic functions to depend on $\xi$ as
well as $\tau$. The energy-momentum conservation equation is

$\displaystyle\partial_{\mu}T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial\tilde{w}}{\partial\tau}+\frac{\tilde{w}}{\tau}\right)u^{\nu}-g^{\mu\nu}\partial_{\mu}\tilde{P}+K\left[D^{\nu}(\partial_{\mu}D^{\mu}n)+D^{\mu}(\partial_{\mu}D^{\nu}n)\right]=0.$
(36)

In what follows we allow the baryon density to depend on $\xi$ but the energy
density and pressure are taken to be boost invariant. As the system evolves,
the baryon density variation will break the boost invariance of the energy
density and pressure which would lead to an acceleration term. In this work we
treat this as a perturbation and neglect the feedback. As we shall see, this
is justified in our numerical example. The implications of a fully boost-
invariace breaking general system of equations is left for future work. Some
useful expressions are

$\partial_{\mu}D^{\mu}=-\frac{1}{\tau^{2}}\frac{\partial^{2}}{\partial\xi^{2}}=D_{\mu}D^{\mu}.$
(37)

Energy conservation is

$\displaystyle\left(\frac{\partial\tilde{\epsilon}}{\partial\tau}+\frac{\tilde{w}}{\tau}\right)\cosh\xi+\frac{\sinh\xi}{\tau}\frac{\partial\tilde{P}}{\partial\xi}+\frac{K}{\tau^{3}}\left[2\sinh\xi\frac{\partial^{2}n}{\partial\xi^{2}}+\cosh\xi\frac{\partial
n}{\partial\xi}\right]\frac{\partial n}{\partial\xi}$ $\displaystyle=$
$\displaystyle 0.$ (38)

and momentum conservation is

$\displaystyle\left(\frac{\partial\tilde{\epsilon}}{\partial\tau}+\frac{\tilde{w}}{\tau}\right)\sinh\xi+\frac{\cosh\xi}{\tau}\frac{\partial\tilde{P}}{\partial\xi}+\frac{K}{\tau^{3}}\left[2\cosh\xi\frac{\partial^{2}n}{\partial\xi^{2}}+\sinh\xi\frac{\partial
n}{\partial\xi}\right]\frac{\partial n}{\partial\xi}$ $\displaystyle=$
$\displaystyle 0.$

Linear combinations yield

$\frac{\partial\tilde{\epsilon}}{\partial\tau}+\frac{\tilde{w}}{\tau}+\frac{K}{\tau^{3}}\left(\frac{\partial
n}{\partial\xi}\right)^{2}=0,$ (39)

and

$\frac{\partial\tilde{P}}{\partial\xi}+\frac{2K}{\tau^{2}}\frac{\partial
n}{\partial\xi}\frac{\partial^{2}n}{\partial\xi^{2}}=0,$ (40)

where

$\tilde{P}=P-\frac{K}{\tau^{2}}n\frac{\partial^{2}n}{\partial\xi^{2}}-\frac{K}{2\tau^{2}}\left(\frac{\partial
n}{\partial\xi}\right)^{2}.$ (41)

Expressed in terms of the local thermodynamic functions they are

$\frac{\partial\epsilon}{\partial\tau}+\frac{w}{\tau}+\frac{K}{\tau^{2}}\frac{\partial
n}{\partial\xi}\frac{\partial^{2}n}{\partial\tau\partial\xi}-\frac{K}{\tau^{3}}n\frac{\partial^{2}n}{\partial\xi^{2}}=0,$
(42)

and

$\frac{\partial
P}{\partial\xi}-\frac{K}{\tau^{2}}n\frac{\partial^{3}n}{\partial\xi^{3}}=0.$
(43)

Differentiating Eq. (42) with respect to $\xi$ and subtracting Eq. (43) leads
to

$\frac{\partial}{\partial\tau}\frac{\partial}{\partial\xi}\left[\tau\epsilon+\frac{K}{2\tau}\left(\frac{\partial
n}{\partial\xi}\right)^{2}\right]=0.$ (44)

Now consider baryon number conservation $\partial_{\mu}J^{\mu}=0$. For
ordinary diffusion one finds

$\frac{\partial}{\partial\tau}(\tau
n)-\frac{\sigma_{B}T}{\tau}\frac{\partial^{2}}{\partial\xi^{2}}\left(\frac{\tilde{\mu}}{T}\right)-\frac{1}{\tau}\frac{\partial}{\partial\xi}(\sigma_{B}T)\frac{\partial}{\partial\xi}\left(\frac{\tilde{\mu}}{T}\right)=0.$
(45)

Note that

$\tilde{\mu}=\mu-\frac{K}{\tau^{2}}\frac{\partial^{2}n}{\partial\xi^{2}}.$
(46)

The $K$ is determined by the surface free energy. We take it to be a constant.
We use the simplified expression for baryon conductivity [38]

$\sigma_{B}=\frac{1}{3}\frac{C_{B}}{T}n\left[\coth\left(\frac{\mu_{B}}{T}\right)-\frac{3Tn}{w}\right],$
(47)

where $C_{B}$ is a dimensionless number. We have taken its value to be 0.5. A
more rigorous expression for $\sigma_{B}$ and its comparison with the above
form is discussed in [40, 52].

In order to avoid solving partial differential equations with derivatives
greater than two we could choose $\epsilon(\tau,\xi)$ and $n(\tau,\xi)$ as the
independent thermodynamic variables. Then $P$, $T$, $w$, $s$, and $\mu$ are
determined via the equation of state $f(n,T)$.

The equations that we solve are

$\frac{\partial\epsilon(n,T)}{\partial\tau}+\frac{w(n,T)}{\tau}+\frac{K}{\tau^{2}}\frac{\partial
n}{\partial\xi}\frac{\partial^{2}n}{\partial\tau\partial\xi}-\frac{K}{\tau^{3}}n\frac{\partial^{2}n}{\partial\xi^{2}}=0,$
(48) $\frac{\partial}{\partial\tau}(\tau
n)-\frac{\sigma_{B}T}{\tau}\frac{\partial^{2}}{\partial\xi^{2}}\left(\frac{\tilde{\mu}}{T}\right)-\frac{1}{\tau}\frac{\partial}{\partial\xi}(\sigma_{B}T)\frac{\partial}{\partial\xi}\left(\frac{\tilde{\mu}}{T}\right)=0,$
(49)
$\tilde{\mu}=\mu(n,T)-\frac{K}{\tau^{2}}\frac{\partial^{2}n}{\partial\xi^{2}}.$
(50)

Figure 3: Phase of the system as a function of rapidity and time. Figure 4:
Baryon density evolution at different times. Proper times for curves, from top
to bottom, are 10.0 fm, 10.5 fm, 11.5 fm, 12.5 fm and 13.5 fm

##  VI Spinodal separation in Bjorken hydrodynamics

We initialize the system at time $\tau_{x}=10$ fm such that the energy density
$\epsilon(\tau_{x})$ is $\xi$-independent while the density has the periodic
fluctuation $n_{fl}\sin(k\xi)$ around the average value $n(\tau_{x})$. The
fluctuation is necessary to seed spinodal decomposition. These are chosen such
that the system just reaches the coexistence curve. Specifically, the point
$\epsilon(\tau_{x})$, $n(\tau_{x})-n_{fl}$ lies on the coexistence curve. We
evolve the system following Eqs. (48-50).

A periodic boundary condition is used in the $\xi$ direction to preserve boost
invariance. The surface energy is proportional to the coefficient $K$. The
surface energy term builds up the energy in the coexistence region and causes
spontaneous separation of phases as the system approaches equilibrium. This
has been a topic of extensive study in the phase dynamics community. Our
system has a unique feature in that it is expanding rapidly. As the system
expands and cools, the phase change is guaranteed and the system which starts
with all QGP will eventually end up as all hadron gas. But different parts of
the system will undergo the phase transition at different times. Figure 3
shows the evolution of the phases for our system at different spacetime
points.

The first order phase transition will imprint its signature on the
correlations of hadrons. More specifically, it will amplify the effects of
local fluctuations which may be seeded at the initial stage or may arise
dynamically. To study the effects of the surface energy terms on the system,
we solve the system with $K=0$ and $K=5\times 10^{-5}$ MeV${}^{-4}\approx
1.5\times 10^{7}$ MeV fm5.

Figure 5: Energy density profile at different times. The curve with the
highest energy density point corresponds to $K=5\times 10^{-5}$ MeV-4 while
the other curve has $K=0$. The energy density in the QGP phase can be lower
than the energy density in the coexistence region or in the hadron gas phase
because of the gradients of baryon density.

The baryon density evolution is primarily governed by the expansion in the
longitudinal direction and by diffusion. The surface term modifies the baryon
density evolution through the modification of the the chemical potential in
Eq. (50). This modification is proportional to the fourth derivative in baryon
density and is very sensitive to local fluctuations. A small nucleus of matter
in the coexistence region can trigger the spinodal decomposition. Due to the
fourth order derivative playing a vital role here, it is incredibly difficult
to model this using finite-difference methods and a realistic heavy-ion
simulation attempting to incorporate first-order phase transition physics will
probably need to resort to finite-element methods tailored for this system.
Here, we have side-stepped this problem by choosing a smoothly varying baryon
density profile. Figure 4 shows the baryon density evolution for $K=0$. We
checked that the baryon density evolution is essentially identical for
different values of $K$ when $K\lesssim 10^{-4}$ MeV-4. This ensures that the
gradients in $n$ are governed by the initial sinusoidal profile and not the
small fluctuations making the numerics simpler. As discussed above, this will
most likely not be the case in a realistic setting where more advanced
techniques will be necessary.

The energy contribution from the surface terms is dependent on the first,
second and third order spatial derivatives of the baryon density in the
coexistence phase (see Eq. (48)). As baryon density is smoothly varying and
the surface correction to it are small, we can numerically evaluate it with
good accuracy. The surface energy contribution is not negligible here for the
chosen values of $K$. Figure 5 shown the evolution of energy density. We see
an immediate energy build up as the system enters the coexistence phase. The
energy evolution in the QGP phase is unaffected for different choices of $K$,
as it has not yet seen the phase boundary. The energy ordering is inverted in
the negative rapidity region at $\tau=13.5$ fm. This is because of the
opposite sign of the baryon density spatial derivatives when this region
reaches the coexistence phase.

We can see the variations in temperature and chemical potential as a
consequence of the phase transition in Figs. 6 and 7. It is interesting to
note that in this particular case, the hadron gas has higher temperature than
the QGP phase. This is because of our choice of the initial state where we
took a constant energy density and a spatially varying baryon density, hence
the phase is being governed by baryon density as the change in energy density
is relatively small. For a constant energy density, the region with lower
baryon density has higher temperature. Since the energy density variation is
small, temperatures are higher where baryon densities are lower and so QGP is
cooler than the hadron gas. On the other hand, the baryon chemical potential
is higher for QGP than for the hadron gas.

Figure 6: Temperature profile at different times. The curve with the highest
temperature point corresponds to $K=5\times 10^{-5}$ MeV-4 while the other
curve has $K=0$.

Figure 7: Baryon chemical potential profile at different times. The curve with
the lowest baryon chemical potential point corresponds to $K=5\times 10^{-5}$
MeV-4 while the other curve has $K=0$.

##  VII Discussion

Numerous and various model calculations have fueled speculation that QCD
exhibits a line of first order phase transition in the temperature versus
baryon chemical potential phase diagram, ending in a critical point at some
$T_{c}$ and $\mu_{c}$, which is in the same universality class as the 3D Ising
model and the liquid-gas phase transition. Due to the fermion sign problem, it
is not known (yet) whether lattice QCD confirms or denies this purported phase
structure. In any event, physics is an experimental science, and the only
conceivable way to test this hypothesis in terrestrial experiments is with
relativistic heavy ion collisions. The Standard Model of these collisions
includes a stage using second order, relativistic hydrodynamics, which
naturally needs an equation of state. If the transition is slow and adiabatic,
the Maxwell construction may be used. If the transition is faster, then
nucleation of hadronic bubbles in the quark-gluon plasma may provide a better
quantitative description. This description applies when the system super-cools
into the metastable region. When the system expands and cools very rapidly, as
is likely the case in heavy ion collisions, a more detailed description may be
necessary. This is commonly referred to as spinodal decomposition, and is
relevant even in the unstable region. Spinodal decomposition dynamics has been
extensively studied in several branches of physics and chemistry, but only for
systems with fixed volume and temperature. The dynamics in heavy ion
collisions is much more challenging.

In this manuscript we have incorporated terms into the stress-energy-momentum
tensor, which are normally included to describe spinodal decomposition in non-
relativistic systems, into covariant hydrodynamics. These terms are second-
order in thermodynamic densities, namely, the baryon density and its space-
time derivatives. Only one new parameter, $K$, enters and is associated with
the surface free energy between the two phases. It is also necessary to have a
parameterization of the equation of state when the system is in the metastable
or unstable regions. This is better accomplished by expanding in powers of the
baryon density rather than the volume per baryon. The resulting equations of
motion are fourth order in space-time derivatives. This makes solution of the
space-time evolution very numerically challenging. Therefore, in this
manuscript we only gave a simple illustration of how this works in a nearly
boost-invariant Bjorken scenario. Future work should include realistic initial
state and hydrodynamic fluctuations. Only then will comparisons of theoretical
simulations with experimental observations fulfill our goal. Given advances in
computational calculations over the past two decades we are optimistic that
the goals can be achieved.

###### Acknowledgements.

We thank U. Heinz and C. Gale for constructive comments on the manuscript.
This work was supported by the U.S. DOE Grants No. DE-FG02-87ER40328 (JIK, MS,
and TW) and DE-SC-0024347 (MS).

## appendix

Parameterization of the metastable and unstable regions in Sec. III is
consistent with mean field values of the critical exponents, specifically
$\delta=3$. For systems in the same universality class as the 3D Ising model
and the liquid-gas phase transition the more accurate value is $\delta\approx
4.79$. In that case we write [53]

$P_{{\rm
int}}(n)=P_{X}(T)+a_{1}(n-n_{G})+a_{2}(n-n_{G})^{2}+a_{3}(n-n_{G})^{\delta-1}+a_{4}(n-n_{G})^{\delta}+a_{5}(n-n_{G})^{\delta+1}$
(51)

and

$\mu_{{\rm
int}}(n)=\mu_{X}(T)+b_{0}\ln(n/n_{G})+b_{1}(n-n_{G})+b_{2}(n-n_{G})^{\delta-1}+b_{3}(n-n_{G})^{\delta}$
(52)

Note the degeneracy when $\delta=3$. At fixed $T$

$n\frac{\partial\mu_{\rm int}}{\partial n}=\frac{\partial P_{\rm
int}}{\partial n}$ (53)

This leads to

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle b_{0}+b_{1}n_{G}$
$\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{1}{2}}}b_{1}$ $\displaystyle a_{3}$
$\displaystyle=$ $\displaystyle b_{2}n_{G}$ $\displaystyle a_{4}$
$\displaystyle=$
$\displaystyle\left(\frac{\delta-1}{\delta}\right)b_{2}+b_{3}n_{G}$
$\displaystyle a_{5}$ $\displaystyle=$
$\displaystyle\left(\frac{\delta}{\delta+1}\right)b_{3}$ (54)

There are four independent parameters to be determined at each temperature $T$
as described in the text.

## References

  * [1] M. Abdallah et al. [STAR], “Measurement of the Sixth-Order Cumulant of Net-Proton Multiplicity Distributions in Au+Au Collisions at $\sqrt{s_{NN}}=$ 27, 54.4, and 200 GeV at RHIC,” Phys. Rev. Lett. 127, no.26, 262301 (2021) doi:10.1103/PhysRevLett.127.262301 [arXiv:2105.14698 [nucl-ex]].
  * [2] See the proceedings of the Workshops of the Critical Point and Onset of Deconfinement (CPOD) series.
  * [3] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, “The Order of the quantum chromodynamics transition predicted by the standard model of particle physics,” Nature 443, 675-678 (2006) doi:10.1038/nature05120 [arXiv:hep-lat/0611014 [hep-lat]].
  * [4] H. T. Ding, F. Karsch and S. Mukherjee, “Thermodynamics of strong-interaction matter from Lattice QCD,” Int. J. Mod. Phys. E 24, no.10, 1530007 (2015) doi:10.1142/S0218301315300076 [arXiv:1504.05274 [hep-lat]].
  * [5] M. Arslandok, S. A. Bass, A. A. Baty, I. Bautista, C. Beattie, F. Becattini, R. Bellwied, Y. Berdnikov, A. Berdnikov and J. Bielcik, et al. [arXiv:2303.17254 [nucl-ex]].
  * [6] P. Hegde [BNL–Bielefeld–CCNU], “The QCD equation of state to $\mathcal{O}(\mu_{B}^{4})$ from lattice QCD,” Nucl. Phys. A 931, 851-855 (2014) doi:10.1016/j.nuclphysa.2014.08.089 [arXiv:1408.6305 [hep-lat]].
  * [7] J. N. Guenther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor, C. Ratti and K. K. Szabó, “The QCD equation of state at finite density from analytical continuation,” Nucl. Phys. A 967, 720-723 (2017) doi:10.1016/j.nuclphysa.2017.05.044 [arXiv:1607.02493 [hep-lat]].
  * [8] A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa, S. Mukherjee, H. Ohno and P. Petreczky, et al. “The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD,” Phys. Rev. D 95, no.5, 054504 (2017) doi:10.1103/PhysRevD.95.054504 [arXiv:1701.04325 [hep-lat]].
  * [9] A. Monnai, B. Schenke and C. Shen, “Equation of state at finite densities for QCD matter in nuclear collisions,” Phys. Rev. C 100, no.2, 024907 (2019) doi:10.1103/PhysRevC.100.024907 [arXiv:1902.05095 [nucl-th]].
  * [10] J. Noronha-Hostler, P. Parotto, C. Ratti and J. M. Stafford, “Lattice-based equation of state at finite baryon number, electric charge and strangeness chemical potentials,” Phys. Rev. C 100, no.6, 064910 (2019) doi:10.1103/PhysRevC.100.064910 [arXiv:1902.06723 [hep-ph]].
  * [11] S. Borsányi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pásztor, C. Ratti and K. K. Szabó, “Lattice QCD equation of state at finite chemical potential from an alternative expansion scheme,” Phys. Rev. Lett. 126, no.23, 232001 (2021) doi:10.1103/PhysRevLett.126.232001 [arXiv:2102.06660 [hep-lat]].
  * [12] M. Kahangirwe, S. A. Bass, E. Bratkovskaya, J. Jahan, P. Moreau, P. Parotto, D. Price, C. Ratti, O. Soloveva and M. Stephanov, “Finite density QCD equation of state: Critical point and lattice-based T’ expansion,” Phys. Rev. D 109, no.9, 094046 (2024) doi:10.1103/PhysRevD.109.094046 [arXiv:2402.08636 [nucl-th]].
  * [13] M. Albright, J. Kapusta and C. Young, “Matching Excluded Volume Hadron Resonance Gas Models and Perturbative QCD to Lattice Calculations,” Phys. Rev. C 90, no.2, 024915 (2014) doi:10.1103/PhysRevC.90.024915 [arXiv:1404.7540 [nucl-th]].
  * [14] A. Sorensen, K. Agarwal, K. W. Brown, Z. Chajecki, P. Danielewicz, C. Drischler, S. Gandolfi, J. W. Holt, M. Kaminski and C. M. Ko, et al. “Dense nuclear matter equation of state from heavy-ion collisions,” Prog. Part. Nucl. Phys. 134, 104080 (2024) doi:10.1016/j.ppnp.2023.104080 [arXiv:2301.13253 [nucl-th]].
  * [15] K. Fukushima and T. Hatsuda, “The phase diagram of dense QCD,” Rept. Prog. Phys. 74, 014001 (2011) doi:10.1088/0034-4885/74/1/014001 [arXiv:1005.4814 [hep-ph]].
  * [16] K. Fukushima and C. Sasaki, “The phase diagram of nuclear and quark matter at high baryon density,” Prog. Part. Nucl. Phys. 72, 99-154 (2013) doi:10.1016/j.ppnp.2013.05.003 [arXiv:1301.6377 [hep-ph]].
  * [17] C. S. Fischer, “QCD at finite temperature and chemical potential from Dyson–Schwinger equations,” Prog. Part. Nucl. Phys. 105, 1-60 (2019) doi:10.1016/j.ppnp.2019.01.002 [arXiv:1810.12938 [hep-ph]].
  * [18] R. Guida and J. Zinn-Justin, “3-D Ising model: The Scaling equation of state,” Nucl. Phys. B 489, 626-652 (1997) doi:10.1016/S0550-3213(96)00704-3 [arXiv:hep-th/9610223 [hep-th]].
  * [19] C. Nonaka and M. Asakawa, “Hydrodynamical evolution near the QCD critical end point,” Phys. Rev. C 71, 044904 (2005) doi:10.1103/PhysRevC.71.044904 [arXiv:nucl-th/0410078 [nucl-th]].
  * [20] P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer and M. Stephanov, “QCD equation of state matched to lattice data and exhibiting a critical point singularity,” Phys. Rev. C 101, no.3, 034901 (2020) doi:10.1103/PhysRevC.101.034901 [arXiv:1805.05249 [hep-ph]].
  * [21] J. Kapusta, T. Welle and C. Plumberg, “Embedding a Critical Point in a Hadron to Quark–Gluon Crossover Equation of State,” Phys. Rev. C 106, no.1, 014909 (2022) doi:10.5506/APhysPolBSupp.16.1-A46 [arXiv:2112.07563 [nucl-th]].
  * [22] J. I. Kapusta and T. Welle, “Extending a scaling equation of state to QCD,” Phys. Rev. C 106, no.4, 044901 (2022) doi:10.1103/PhysRevC.106.044901 [arXiv:2205.12150 [nucl-th]].
  * [23] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, “Signatures of the tricritical point in QCD,” Phys. Rev. Lett. 81, 4816-4819 (1998) doi:10.1103/PhysRevLett.81.4816 [arXiv:hep-ph/9806219 [hep-ph]].
  * [24] M. A. Stephanov, “Non-Gaussian fluctuations near the QCD critical point,” Phys. Rev. Lett. 102, 032301 (2009) doi:10.1103/PhysRevLett.102.032301 [arXiv:0809.3450 [hep-ph]].
  * [25] Y. Akamatsu, A. Mazeliauskas and D. Teaney, “A kinetic regime of hydrodynamic fluctuations and long time tails for a Bjorken expansion,” Phys. Rev. C 95, no.1, 014909 (2017) doi:10.1103/PhysRevC.95.014909 [arXiv:1606.07742 [nucl-th]].
  * [26] J. I. Kapusta, B. Muller and M. Stephanov, “Relativistic Theory of Hydrodynamic Fluctuations with Applications to Heavy Ion Collisions,” Phys. Rev. C 85, 054906 (2012) doi:10.1103/PhysRevC.85.054906 [arXiv:1112.6405 [nucl-th]].
  * [27] C. Young, J. I. Kapusta, C. Gale, S. Jeon and B. Schenke, “Thermally Fluctuating Second-Order Viscous Hydrodynamics and Heavy-Ion Collisions,” Phys. Rev. C 91, no.4, 044901 (2015) doi:10.1103/PhysRevC.91.044901 [arXiv:1407.1077 [nucl-th]].
  * [28] M. Singh, C. Shen, S. McDonald, S. Jeon and C. Gale, “Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions,” Nucl. Phys. A 982, 319-322 (2019) doi:10.1016/j.nuclphysa.2018.10.061 [arXiv:1807.05451 [nucl-th]].
  * [29] X. An, G. Basar, M. Stephanov and H. U. Yee, “Relativistic Hydrodynamic Fluctuations,” Phys. Rev. C 100, no.2, 024910 (2019) doi:10.1103/PhysRevC.100.024910 [arXiv:1902.09517 [hep-th]].
  * [30] M. Martinez and T. Schäfer, “Stochastic hydrodynamics and long time tails of an expanding conformal charged fluid,” Phys. Rev. C 99, no.5, 054902 (2019) doi:10.1103/PhysRevC.99.054902 [arXiv:1812.05279 [hep-th]].
  * [31] A. Sakai, K. Murase and T. Hirano, “Rapidity decorrelation of anisotropic flow caused by hydrodynamic fluctuations,” Phys. Rev. C 102, no.6, 064903 (2020) doi:10.1103/PhysRevC.102.064903 [arXiv:2003.13496 [nucl-th]].
  * [32] L. Du, U. Heinz, K. Rajagopal and Y. Yin, “Fluctuation dynamics near the QCD critical point,” Phys. Rev. C 102, no.5, 054911 (2020) doi:10.1103/PhysRevC.102.054911 [arXiv:2004.02719 [nucl-th]].
  * [33] G. Pihan, M. Bluhm, M. Kitazawa, T. Sami and M. Nahrgang, “Critical net-baryon fluctuations in an expanding system,” Phys. Rev. C 107, no.1, 014908 (2023) doi:10.1103/PhysRevC.107.014908 [arXiv:2205.12834 [nucl-th]].
  * [34] C. Chattopadhyay, J. Ott, T. Schaefer and V. Skokov, “Dynamic scaling of order parameter fluctuations in model B,” Phys. Rev. D 108, no.7, 074004 (2023) doi:10.1103/PhysRevD.108.074004 [arXiv:2304.07279 [nucl-th]].
  * [35] G. Başar, J. Bhambure, R. Singh and D. Teaney, “The stochastic relativistic advection diffusion equation from the Metropolis algorithm,” [arXiv:2403.04185 [nucl-th]].
  * [36] X. An, M. Bluhm, L. Du, G. V. Dunne, H. Elfner, C. Gale, J. Grefa, U. Heinz, A. Huang and J. M. Karthein, et al. “The BEST framework for the search for the QCD critical point and the chiral magnetic effect,” Nucl. Phys. A 1017, 122343 (2022) doi:10.1016/j.nuclphysa.2021.122343 [arXiv:2108.13867 [nucl-th]].
  * [37] L. Du, A. Sorensen and M. Stephanov, “The QCD phase diagram and Beam Energy Scan physics: a theory overview,” [arXiv:2402.10183 [nucl-th]].
  * [38] C. Shen, G. Denicol, C. Gale, S. Jeon, A. Monnai and B. Schenke, “A hybrid approach to relativistic heavy-ion collisions at the RHIC BES energies,” Nucl. Phys. A 967, 796-799 (2017) doi:10.1016/j.nuclphysa.2017.06.008 [arXiv:1704.04109 [nucl-th]].
  * [39] Y. Akamatsu, M. Asakawa, T. Hirano, M. Kitazawa, K. Morita, K. Murase, Y. Nara, C. Nonaka and A. Ohnishi, “Dynamically integrated transport approach for heavy-ion collisions at high baryon density,” Phys. Rev. C 98, no.2, 024909 (2018) doi:10.1103/PhysRevC.98.024909 [arXiv:1805.09024 [nucl-th]].
  * [40] A. De, J. I. Kapusta, M. Singh and T. Welle, “Comprehensive simulation of heavy-ion collisions at nonzero baryon chemical potential,” Phys. Rev. C 106, no.5, 054906 (2022) doi:10.1103/PhysRevC.106.054906 [arXiv:2206.02655 [nucl-th]].
  * [41] L. P. Csernai and J. I. Kapusta, “Dynamics of the QCD phase transition,” Phys. Rev. Lett. 69, 737-740 (1992) doi:10.1103/PhysRevLett.69.737
  * [42] A. M. Jokisaari, P. W. Voorhees, J. E. Guyer, J. Warren, and O. G. Heinonen, “Benchmark problems for numerical implementations of phase field models,” Computational Materials Science 126, 139 (2017) doi:10.1016/j.commatsci.2016.09.022 [arXiv:1610.00622[cond-mat.mtrl-sci]].
  * [43] N. Provatas and K. Elder, “Phase-Field Methods in Material Science and Engineering,” John Wiley and Sons, Ltd, 2011, ISBN 978-3-52-763237-4 doi:10.1002/9783527631520
  * [44] J. Randrup, “Phase transition dynamics for baryon-dense matter,” Phys. Rev. C 79, 054911 (2009) doi:10.1103/PhysRevC.79.054911 [arXiv:0903.4736 [nucl-th]].
  * [45] J. Randrup, “Spinodal phase separation in relativistic nuclear collisions,” Phys. Rev. C 82, 034902 (2010) doi:10.1103/PhysRevC.82.034902 [arXiv:1007.1448 [nucl-th]].
  * [46] J. W. Cahn and J. E. Hilliard, “Free Energy of a Nonuniform System. I. Interfacial Free Energy,” J. Chem. Phys. 28, 258 (1958) doi:10.1063/1.1744102
  * [47] J. W. Cahn and J. E. Hilliard, “Free Energy of a Nonuniform System. III. Nucleation in a Two‐Component Incompressible Fluid,” J. Chem. Phys. 31, 688 (1959) doi:10.1063/1.1730447
  * [48] A. J. M. Yang, P. D. Fleming, and J. H. Gibbs, “Molecular theory of surface tension,” J. Chem. Phys. 64 3732 (1976) doi:10.1063/1.432687
  * [49] C. Cattaneo, “On the Conduction of Heat”, Atti del Semin. Mat. e Fis. Univ. Modena 3, 3 (1948).
  * [50] C. Cattaneo, “A Form of Heat-Conduction Equations Which Eliminates the Paradox of Instantaneous Propagation”, Comptes rendus de l’Académie des sciences 247, 431 (1958).
  * [51] J. I. Kapusta and C. Young, “Causal Baryon Diffusion and Colored Noise”, Phys. Rev. C 90, 044902 (2014) doi:10.1103/PhysRevC.90.044902 [arXiv:1404.4894 [nucl-th]]
  * [52] J. I. Kapusta and J. M. Torres-Rincon, “Thermal Conductivity and Chiral Critical Point in Heavy Ion Collisions,” Phys. Rev. C 86, 054911 (2012) doi:10.1103/PhysRevC.86.054911 [arXiv:1209.0675 [nucl-th]].
  * [53] J. I. Kapusta, “Equation of State and Phase Fluctuations near the Chiral Critical Point,” Phys. Rev. C 81, 055201 (2010) doi:10.1103/PhysRevC.81.055201 [arXiv:1001.3650 [nucl-th]].


††thanks: Corresponding author E-mail<EMAIL_ADDRESS>Corresponding author E-mail<EMAIL_ADDRESS>Corresponding author
E-mail<EMAIL_ADDRESS>
# A simple scheme of low phase noise microwave synthesizers based on the sub-
sampling phase lock loop

Wenbing Li  MOE Key Laboratory of Fundamental Physical Quantities Measurement,
Hubei Key Laboratory of Gravitation and Quantum Physics, School of Physics,
Huazhong University of Science and Technology,
1037 Luoyu Road, Wuhan 430074, P. R. China  Qiang Hao  Key Laboratory of Time
and Frequency Primary Standards, National Time Service Center, Chinese Academy
of Sciences, Xi’an 710600 P. R. China  Yuanbo Du  Shaoqing Huang  MOE Key
Laboratory of Fundamental Physical Quantities Measurement, Hubei Key
Laboratory of Gravitation and Quantum Physics, School of Physics, Huazhong
University of Science and Technology,
1037 Luoyu Road, Wuhan 430074, P. R. China  Peter Yun  Key Laboratory of Time
and Frequency Primary Standards, National Time Service Center, Chinese Academy
of Sciences, Xi’an 710600 P. R. China  Zehuang Lu  MOE Key Laboratory of
Fundamental Physical Quantities Measurement, Hubei Key Laboratory of
Gravitation and Quantum Physics, School of Physics, Huazhong University of
Science and Technology,
1037 Luoyu Road, Wuhan 430074, P. R. China

###### Abstract

In this paper, we demonstrate a simple scheme of 6.835 GHz microwave frequency
synthesizer based on the sub-sampling phase lock loop (PLL) technique. The
application of the sub-sampling PLL is the key to simplify the architecture of
the synthesizer in this scheme. A 100 MHz oven controlled crystal oscillator
(OCXO) with ultra-low phase noise is used as the initial signal source. Then a
dielectric resonant oscillator (DRO) of 6.8 GHz is directly phase locked to
the 100 MHz OCXO utilizing the sub-sampling PLL. Benefiting from the sub-
sampling PLL, the processes of microwave frequency multiplication and phase
lock which are necessary in the development of microwave synthesizer are
greatly simplified. Therefore, the architecture of the synthesizer is very
simple. Correspondingly, the power consumption and cost of the synthesizer are
low. The absolute phase noises of the 6.835 GHz output signal are measured to
be -47 dBc/Hz, -77 dBc/Hz, -104 dBc/Hz and -121 dBc/Hz at 1 Hz, 10 Hz, 100 Hz
and 1 kHz offset frequencies, respectively. he synthesizer can be used as the
local oscillator of the Rb atomic clocks. For the Rb atomic clocks operated in
the continuous or pulsed optically pumped (POP) mode, Tthe short-term
frequency stability limited by the absolute phase noises of the synthesizer
through the intermodulation or the Dick effect is theoretically calculated to
be better than $5.0\times 10^{-14}\tau^{1/2}$. This low phase noise microwave
frequency synthesizer can be used in other experiments of fundamental physics
measurement.

##  I INTRODUCTION

Low phase noise microwave frequency synthesizers are generally used in the
experiments of Raman atomic interferometer,Hu2013 (1) testing local position
invariance,Ashby2007 (2) precision atomic spectroscopy measurement,Fortier2007
(3) atomic spin squeezing,Chen2012 (4) and atomic clock.Du2015 (5, 6, 7, 8, 9,
10, 11) The phase noise of the microwave synthesizer is one of the key factors
in limiting the precision of those measurements.

Among the atomic clocks, the Rb atomic clocks are widely used in the fields of
telecommunications, global navigation satellite system (GNSS) and industrial
applications because of its superior characteristic of high frequency
stability, low power consumption, and small volume.Camparo2007 (12, 13, 14) In
recent years with widespread applications promoted by GNSS, the frequency
stability of the Rb atomic clocks has been greatly improved. For example, the
short-term frequency stability of the Rb atomic clocks operated in a
continuous optically pump mode has been better than $2.5\times
10^{-13}\tau^{-1/2}$.Bandi2014 (15, 16) A pulsed optically pumped (POP) Rb
atomic clock with short-term frequency stability of $1.7\times
10^{-13}\tau^{-1/2}$ has been developed.Calosso2012 (17)

However, no matter whether the Rb atomic clocks operated in the mode of
continuous optically pumping or POP, the absolute phase noise of local
oscillator is still one of the main factors to limit their short-term
frequency stability. This has been analyzed and verified by several
researches.Dick1987 (19, 18, 20) For the continuous interrogation atomic
clocks, the frequency modulation noise of the local oscillator at even
harmonics offset frequencies would affect its short-term frequency stability
via the intermodulation effect.Deng1999 (18) For the POP clocks, its short-
term frequency stability would be affected by the absolute phase noise of the
local oscillator through the Dick effect.Dick1987 (19) For both type Rb atomic
clocks, their short-term frequency stabilities would be limited by the
absolute phase noise of the local oscillator, especially at the offset
frequency range of 100 Hz to 1 kHz. Therefore, in order to improve the short-
term frequency stability of Rb atomic clocks, it is beneficial to develop a
synthesizer with low phase noise at these offset frequencies.

TThe short-term frequency stability of the atomic clocks could be considerably
improved by the ultra-low phase noise of the ultra-stable laser millo2009 (21,
22) or cryogenic sapphire oscillator (CSO) based synthesizer.Lipphardt2008
(23, 24, 25) Nevertheless, for the compact Rb clocks used in the GNSS or
industry applications, it is not realistic to use these complex, expensive and
bulky systems. Thus, it will still remain a well-adapted solution to develop
the ultra-low noise quartz crystal oscillator based microwave frequency
synthesizer, and there have been several reports on developing microwave
synthesizers for atomic clocks.Fernando2010 (7, 8, 9, 10, 11, 12) In most of
these schemes, to phase lock a dielectric resonant oscillator (DRO), the 100
MHz OCXO is usually frequency multiplied to 6.8 GHz (for the Rb atomic clocks)
or 9.2 GHz (for the Cs atomic clocks) with a non-linear transmission line
(NLTL) or a Step Recovery Diode (SRD), and then mixed with the output signal
of the DRO to generate an intermediate frequency signal which is further phase
locked to a DDS, whose referenced signal is also provided by the OCXO. In that
way the output frequency of the DRO is indirectly phase locked to the 100 MHz
of OCXO. In this process, our previous work shows that the additional phase
noise of the NLTL are greatly affected by its operation parameters such as the
input power, the input and output impedances.Li2018 (10) To suppress the
additional phase noise of the NLTL, the architecture of the schemes would be
very complex, bulky, and expensive for optimizing the operation parameters of
the NLTL. To simplify the structure, B. François et al. designed a simple
microwave synthesizer with a 1.6 GHz custom-designed frequency multiplication
module. But the subsequent frequency division and multiplication are still
needed. Therefore, the structure of the synthesizer is still quite
complicated.Francois2015 (11) In this paper, we demonstrate a simple scheme of
a 6.835 GHz low phase noise microwave synthesizer based on the sub-sampling
PLL for Rb atomic clock. With the help of the sub-sampling PLL technique,
which is usually used in the radio-frequency and radar engineeringGao2015 (26,
27), the structure of the synthesizer can be greatly simplified.T And the
absolute phase noise is a litte better than the result in ref[8]. The
microwave frequency synthesizer developed here can also be used in the
experiments of the other experiments of fundamental physics measurement.

##  II Scheme of the synthesizer

The architecture of the synthesizer is shown in Fig. 1. The 100 MHz OCXO
(Beijing Institute of Radio Metrology and Measurement, ZF550) is divided into
two arms via a power splitter (Mini-Circuits, ZX10-2-12+). The first arm 100
MHz signal is used as the reference signal for a home-made direct digital
synthesizer (DDS) that is based on AD9852 chip, which is used to achieve
functions of frequency modulation and scanning. The DDS is controlled through
serial communication with a computer. The second arm 100 MHz signal is used to
phase lock the 6.8 GHz DRO utilizing a sub-sampling PLL. With the help of the
sub-sampling PLL technique, the DRO is directly phase locked to the 100 MHz
OCXO. The microwave frequency multiplication and mixing process are greatly
simplified. The phase locked DRO is labelled as PLDRO. The 6.8 GHz output
signal of the PLDRO is passed through a microwave isolator (Fairview,
FMIR1021) to prevent the microwave signal from being reflected back into the
PLDRO. The output signal of the PLDRO is then mixed with the 34.68 MHz output
signal of the DDS with a mixer (Mini-Circuits, ZX05-83+) to generate the final
6.835 GHz output signal. Compared to the previous works on the development of
microwave synthesizer, the structure of the synthesizer of Fig. 1 is much
simpler. Correspondingly, the power consumption and the cost of the
synthesizer are considerably decreased.Fernando2010 (7, 8, 11, 10)

The detailed schematic of the sub-sampling PLL is shown in Fig. 2. The
sampling phase detector is the key component of the PLL. The sampling phase
detector mainly consisted of a step recovery diode (SRD) and a pair of
Schottky diodes. The maximum frequency of the phase detector can reach 18 GHz.
The input impedance of the SRD is usually about 4 - 10 $\Omega$, but the
output impedance of the 100 MHz OCXO is usually 50 $\Omega$. Thus, we design
an input matching circuit to match their impedances, which is realized by a
transformer. In addition for matching impedance, the transformer is also used
to convert unipolar pulses into symmetrical pulses. The loop pass filter (LPF)
of the sub-sampling PLL is realized by a positive proportional-integral (PI)
circuit. A voltage follower behind the sampling phase detector is designed to
serve functions of isolation and impedance matching between the sampling phase
detector and the LPF, and also as a low pass filter to suppress the high-
frequency noise of the PLL. The scanning circuit is designed to improve the
capture range of the PLL and to shorten phase locking time. The volume of the
developed PLDRO is only slightly larger than that of a microwave amplifier.
Thus, the volume of the synthesizer is reduced.

Figure 1: (Colors online) The architecture of the microwave synthesizer, the part enclosed in the red dashed line is the phase locked DRO, labelled as PLDRO. Figure 2: (Colors online) The detailed schematic of the sub-sampling PLL Table 1: The phase noise sensitivity floor of the FSWP26 with option B61 (Xcorr=1, dBc/Hz) Carrier frequency | 1 Hz | 10 Hz | 100 Hz | 1 kHz | 10 kHz | 100 kHz | 1 MHz
---|---|---|---|---|---|---|---
100 MHz | -92 | -115 | -140 | -166 | -170 | -173 | -175
7 GHz | -55 | -78 | -103 | -133 | -152 | -153 | -157

##  III Results and discussion

We measure the absolute phase noises of the main output signals of the
microwave synthesizer at carrier frequencies of 100 MHz, 35 MHz, 6.8 GHz and
6.835 GHz with a Phase Noise Analyzer (R$\&$S FSWP26 with option 61) using
cross-correlation technique for enhanced sensitivity. The noise floor of the
FSWP26 with option B61 is shown in Table I (taken from FSWP26 datasheet). The
numbers given in Table I is given under the condition of the cross-correlation
factor of one (Xcorr=1). The noise floor can be further lowered with an
increasing Xcorr as 5*log(Xcorr). In the actual measurement, the Xcorr and the
bandwidth are set as 50 and 3$\%$, respectively. The phase noise floor are -55
dBc/Hz, -85 dBc/Hz, -114 dBc/Hz, -138 dBc/Hz and -142 dBc/Hz at 1 Hz, 10 Hz,
100 Hz, 1 kHz and 10 kHz offset frequencies, respectively, which will not
limit the measurement accuracy of the phase noise. The measured phase noises
are shown in Fig. 3. The absolute phase noise of the OCXO 100 MHz has a slight
bump at the offset frequencies around 100 kHz, which is inherent in the
custom-made free running OCXO. The phase noise of the 6.8 GHz signal of the
PLDRO at offset frequencies smaller than 100 Hz is deteriorated by an
idealized 20*log(N). At offset frequencies larger than 100 Hz, especially at
the offset frequency of locking bandwidth, the maximum phase noise
deterioration is increased by an additional 5 dB. This is because the absolute
phase noise of free running DRO at small offset frequencies is quite poor. The
locking bandwidth of PLL is about 150 kHz. Thus, the gain within the servo
loop bandwidth is enough. As a consequence, the phase noise, especially at the
offset frequencies near the locking bandwidth, cannot be sufficiently
suppressed. Even with this limitation, the phase noise performances are still
better than the results of the scheme based on the NLTL or SRD. That is
because the impedance matching is much better than the NLTL based scheme. The
phase noise of the 6.835 GHz output signal is the same as the phase noise of
the PLDRO. The phase noise of the DDS 35 MHz signal is excellent, which does
not limit the phase noise of the 6.835 GHz output signal. Using the absolute
phase noise of the 6.835 GHz output signal. the short-term frequency stability
of the atomic clock is estimated to be better than $5.0\times
10^{-14}\tau^{1/2}$.

Figure 3: (Colors online) The absolute phase noises of the main output signals
from the microwave synthesizer. The carrier frequencies are at: (a) 100 MHz of
the OCXO, (b) 35 MHz of the DDS, (c) 6.8 GHz deteriorated by idealized
20*log(N) from 100 MHz OCXO,

(d) 6.8 GHz of the PLDRO, (e) 6.835 GHz output signal.

In conclusion, we demonstrate a simple low phase noise 6.835 GHz microwave
synthesize based on the technique of sub-sampling PLL, with the associated low
power consumption and cost. The absolute phase noises of the 6.835 GHz output
signal are measured to be -47 dBc/Hz, -77 dBc/Hz, -104 dBc/Hz and -121 dBc/Hz
at 1 Hz, 10 Hz, 100 Hz and 1 kHz offset frequencies, respectively. When it is
used as the local oscillator of the Rb atomic clocks, the short-term frequency
stability limit of the Rb atomic clocks due to the absolute phase noise via
the intermodulation or the Dick effect is theoretically estimated to be better
than $5.0\times 10^{-14}\tau^{1/2}$. While in this work we only develop a
microwave synthesizer for Rb atomic clocks, the same scheme can be easily
adapted to be used for Cs atomic clocks by simply changing the output
frequencies of the PLDRO and DDS. The microwave frequency synthesizer can also
be used in the experiments of the Raman atomic interferometers or gyroscopes
and other experiments of fundamental physics measurement. Limited by the phase
noise of the DRO and gain and bandwidth of the PLL, the phase noise
performance of the synthesizer realized via this method would not expect to
reach the performance of the ultra-stable laser or CSO based synthesizer.

This work was partially supported by the National Key R$\&$D Program of China
(Grant No. 2017YFA0304400), the National Natural Science Foundation of China
(Grant No. 91336213, 11703031, U1731132 and 11774108) and the Post-doctoral
Science Foundation of China (Grant No. 2017M612431).

## References

  * (1) Z. K. Hu, B. L. Sun, X. C. Duan, M. K. Zhou, L. L. Chen, S. Zhan,Q. Z. Zhang and J. Luo, Phys. Rev. A 88, 043601 (2013).
  * (2) N. Ashby, T. P. Heavner, S. R. Jefferts, T. E. Parker, A. G. Radnaev, and Y. O. Dudin, Phys. Rev. Lett. 98, 070802(2007)
  * (3) T. M. Fortier, N. Ashby, J. C. Bergquist, M. J. Delaney, S. A. Diddams, T. P. Heavner, L. Hollberg, W. M. Itano, S. R. Jefferts, K. Kim, F. Levi, L. Lorini, W. H. Oskay, T. E. Parker, J. Shirley, and J. E. Stalnaker, Phys. Rev. Lett. 98, 070801(2007)
  * (4) Z. L. Chen, J. G. Bohnet, J. M. Weiner and J. K. Thompson, Rev. Sci. Instrum. 83, 044701 (2012).
  * (5) Y. B. Du, R. Wei, R. C. Dong, F. Zou, and Y. Z. Wang, Chin. Phys. B 24, 070601 (2015).
  * (6) F. Levi, D. Calonico, C. E. Calosso, A. Godone, S. Micalizio and G. A. Costanzo, Metrologia 51, 270(2014).
  * (7) F. Ramírez-Martinez, M. Lours, P. Rosenbusch, F. Reinhard and J. Reichel, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 57, 88(2010).
  * (8) B. François, C. E. Calosso, J. M. Danet and R. Boudot, Rev. Sci. Instrum. 85, 094709 (2014).
  * (9) R. Boudot, S. Guerandel and E. De Clercq, IEEE Trans. Instrum. Meas. 58, 3659 (2009).
  * (10) W. B. Li, Y. B. Du, H. Li and Z. H. Lu, AIP Adv. 8, 095311 (2018).
  * (11) B. François, C. E. Calosso, M. Abdel Hafiz, S. Micalizio and R. Boudot, Rev. Sci. Instrum. 86, 094707 (2015).
  * (12) T. P. Heavner, S. R. Jefferts, E. A. Donley, T. E. Parker, F. Levi, Proceedings of the 2005 IEEE International Frequency Control Symposium and Exposition, 308(2005) 86, 094707 (2015).
  * (13) J. C. Camparo, Phys. Today 60, 33 (2007).
  * (14) F. Vannicola, R. Beard, J. White, K. Senior, M. Largay and J. A. Buisson, in proceedings of 42nd PTTI System and Applications Meeting, 181 (2014).
  * (15) S. Micalizio, F. Levi, A. Godone, C. E. Calosso and I. Nazionale, in proceedings of IFCS and EFTF, 1 (2015).
  * (16) T. Bandi, C. Affolderbach, C. Stefanucci, F. Merli, A. K. Skrivervik and G. Mileti, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 61, 1769 (2014).
  * (17) Q. Hao, W. Li, S. He, J. Lv, P. Wang and G. Mei, Rev. Sci. Instrum. 87, 123111 (2016).
  * (18) C. E. Calosso, A. Godone, F. Levi and S. Micalizio, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 59, 2646 (2012).
  * (19) J. Q. Deng, G. Mileti, R. E. Drullinger, D. A. Jennings and F. L. Walls, Phys. Rev. A 59, 773 (1999).
  * (20) G. J. Dick, in Proceedings of PTTI 19, 133 (1987).
  * (21) A. Joyet, G. Mileti, G. Dudle and P. Thomann, IEEE Trans. Instrum. Meas.50, 150 (2001).
  * (22) J. Millö, R. Boudot, M. Lours, P. Y. Bourgeois, A. N. Luiten, Y. L. Coq, Y. Kersalé and G. Santarelli, Opt. Lett. 34, 3707 (2009).
  * (23) T. M. Fortier, M. S. Kirchner, F. Quinlan, J. Taylor, J. C. Bergquist, T. Rosenband, N. Lemke, A. Ludlow, Y. Jiang and C. W. Oates, Nature Photonics 5, 425 (2011).
  * (24) B. Lipphardt, G. Grosche, U. Sterr, C. Tamm, S. Weyers and H. Schnatz, IEEE Trans. Instru. Meas. 58, 1258 (2008).
  * (25) M. Abgrall, J. Guéna, M. Lours, G. Santarelli, M. E. Tobar, S. Bize, S. Grop, B. Dubois, C. Fluhr and V. Giordano, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63, 1198 (2016).
  * (26) A. Takamizawa, S. Yanagimachi, T. Tanabe, K. Hagimoto, I. Hirano, K. Watabe, T. Ikegami and J. G. Hartnett, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 61, 1463 (2014).
  * (27) X. Gao, E. Klumperink and B. Nauta, 2015 IEEE Custom Integrated Circuits Conference (CICC), (2015).
  * (28) X. Gao, G. Eric A. M. Kpluperink, Paul F. J. Geraedts and B. Nauta, IEEE Trans. Circuits and Systems II: Express Briefs 56, 117 (2009).