Hydrodynamic Stability of Cosmological
Quark-Hadron Phase Transitions

Abstract

Recent results from linear perturbation theory suggest that
first-order cosmological quark-hadron phase transitions occurring
as deflagrations may be “borderline” unstable, and
those occurring as detonations may give rise to growing modes behind the
interface boundary. However, since nonlinear effects can play important
roles in the development of perturbations,
unstable behavior cannot be asserted entirely by linear analysis, and
the uncertainty of these recent studies is compounded further by
nonlinearities in the hydrodynamics and self-interaction
fields. In this paper we investigate the growth
of perturbations and the stability of
quark-hadron phase transitions
in the early Universe by solving numerically
the fully nonlinear relativistic hydrodynamics equations
coupled to a scalar field with a quartic self-interaction
potential regulating the transitions. We consider single, perturbed,
phase transitions propagating either by detonation or deflagration,
as well as multiple phase and shock front interactions in
1+2 dimensional spacetimes.

I Introduction

First-order phase transitions ocurring at either the electroweak
or QCD symmetry breaking epochs in the early Universe, as predicted
by the standard model of cosmology,
can have important consequences for the history of
our Universe. In particular, the formation of
co-existing bubbles and droplets of different phases during the QCD transition
affects the production and distribution of
hadrons, and may lead to significant baryon number density
fluctuations.
Assuming these fluctuations survive to the epoch of
Big Bang nucleosynthesis, they can strongly influence the predicted
light element abundances Fuller et al. (1988).
This, in turn, may alter our current understanding
and interpretation of homogeneous Big Bang nucleosynthesis,
as well as the dynamical and chemical
evolution of matter structures, including the origins of
primordial magnetic fields and structure
perturbations that seed the subsequent formation of
stellar, galactic, and cluster scale systems.

Baryon density fluctuations evolve hydrodynamically by the competing
effects of local phase mixing and phase separation that
may arise during the transition period. A potential mixing mechanism
that we consider here
is the hydrodynamic instability of the interface surfaces
separating regions of different phase. Unstable modes may
distort the phase boundary and
induce mixing and diffusive homogenization
through hydrodynamic turbulence.
Whether unstable modes can exist in
cosmological phase fronts has been discussed
recently by several authors with conflicting conclusions
based on first-order perturbation theory.
For example, Kamionkowski and Freese Kamionkowski and Freese (1992)
suggest that subsonic deflagration fronts in electroweak transitions
may be accelerated until they become supersonic and proceed
as detonations through an effective increase
in the front velocity due to surface distortion effects as
the transition becomes turbulent.
Link Link (1992) indicated that the phase boundary in slow
deflagrations from quark-hadron transitions may be unstable to
long wavelength perturbations, but may be stabilized by
surface tension at shorter wavelengths.
In contrast, Huet et al. Huet et al. (1993) find that for electroweak
transitions, the shape of the phase boundary is stable under
perturbations, and quark-hadron transitions are at the “borderline”
between stable and unstable. They also argue that unstable
modes are not possible in detonation fronts. However,
Abney Abney (1994) suggests that this is the case
for the quark phase ahead of the supersonic interface boundary,
but that the cold
phase hadron regions behind the bubbles might be unstable,
at least for the Chapman-Jouget type detonations
investigated.

These results are certainly not conclusive,
and in some cases are in apparent contradiction due
to the level of approximation and coupling
assumed between thermal and dynamical variables.
Also, the full consequences or even existence
of instabilities cannot be determined entirely by
linear analysis. Nonlinear effects, including
higher-order coupling between hydrodynamic, microphysical,
scalar field, and interaction potentials, as well as
surface tension, entropy flow,
and baryon transport, may all play important
roles in the stabilization (or de-stabilization) of
phase boundaries, and remain to be investigated in detail.
Hence, we undertake this study to more fully investigate the hydrodynamic
stability of cosmological phase transitions ocurring during the QCD epoch.
We take a numerical approach and thus generalize previous
results by solving the multi-dimensional
relativistic hydrodynamics equations, presented
in §II, coupled together with a model
scalar wave equation and an interaction potential
regulating the phase transitions. Using results from
perturbation theory to define initial data (as discussed
in §III), we solve numerically for the evolutions
of single distorted phase fronts as well as interactions and
collisions of multiple front systems propagating as either
deflagrations or detonations. Our numerical results are
presented in §IV, and summarized in §V.

Ii Basic Equations

The equations formulated with internal fluid energy are derived
from the 4-velocity normalization uμuμ=−1,
the parallel component of the stress–energy conservation equation
uν∇μTμν=0 for internal energy,
the transverse component of the stress–energy equation
(gαν+uαuν)∇μTμν=0 for momentum,
and an equation of state
for the fluid pressure and temperature, using traditional
high energy units in which c=ℏ=kB=1.
We consider the following special relativistic stress-energy tensor

Tμν=ρhuμuν+prgμν+∂μϕ∂νϕ−gμν(12∂αϕ∂αϕ+V(ϕ,T)),

(1)

which includes both fluid and scalar field (ϕ) contributions.
ρh is the relativistic enthalpy,
pr is the radiation pressure,
uμ is the contravariant 4-velocity, and
gμν=ημν is the flat space 4-metric.

For the
scalar field self-interaction potential V(ϕ,T), we adopt the form
Linde (1983)

The potential (2) is plotted at four
different temperatures T/Tc=(2,1,0.9,0.5) as a function of
scalar field ϕ in Figure 1. At high temperatures (T>Tc)
the potential exhibits a single minimum (ϕ=V(ϕ,T)=0),
where the cosmological fluid is essentially a
quark plasma of unbroken symmetry.
A second minimum (ϕ=ϕmin(T)), corresponding to a low temperature
hadron phase, occurs at temperatures between T0 and Tc, where

ϕmin(T)=αT2λ⎛⎝1+√1−4λβ(T2−T20)α2T2⎞⎠.

(4)

Below T0 the potential has a local maximum at ϕ=0 and a single
minimum at ϕ=ϕmin(T).

Figure 1: Interaction potential V(ϕ,T)
(equation (2) in the text) as a function of
scalar field at four different temperatures
T/Tc=(2,1,0.9,0.5). The potential parameters
used here are:
surface tension σ=0.1Tc,
correlation length ℓc=1.0/Tc,
latent heat L=2.0T4c,
and critical temperature Tc=1.

The resulting differential equations for momentum and energy, neglecting baryon
conservation, can be written in flux conservative form as

∂tSj+∂i(Sjvi)

=

−∂j[pr−V(ϕ,T)]−[κW∂tϕ+κWvi∂iϕ+∂ϕV(ϕ,T)]∂jϕ,

(5)

∂tE+∂i(Evi)

=

−[pr−V(ϕ,T)][∂tW+∂i(Wvi)]

(6)

+κ(W∂tϕ+Wvi∂iϕ)2+∂ϕV(ϕ,T)(W∂tϕ+Wvi∂iϕ),

which, together with the scalar field wave equation,

∂2tϕ−∂i∂iϕ=−∂ϕV(ϕ,T)−κW(∂tϕ+vi∂iϕ),

(7)

and a prescription for the equation of state described below,
represent the complete system of equations we solve in this paper.
In deriving the above equations, we
have implicitly assumed a flat space metric,
and defined variables so that
W=u0=1/√1−vivi is the relativistic boost factor,
vi=ui/u0 is the transport velocity,
Si=Wρhui is the covariant momentum density,
E is the generalized internal energy density

EW=3arT4+V(ϕ,T)−T∂TV(ϕ,T),

(8)

ar=g∗(π2/90) is a radiation constant defining the particle
content and degrees of freedom g∗,
and κ is a constant coefficient regulating entropy
production and kinematic and scalar field energy dissipation.
The terms associated with κ act
as an effective friction force at the phase
transition surface boundary Ignatius et al. (1994); Kurki-Suonio and Laine (1996).

The hydrodynamic equations are solved using the Cosmos code
Anninos and Fragile (2003); Anninos et al. (2002)
with time-explicit operator split methods,
second order spatial finite differencing, and artificial
viscosity to capture shocks Wilson (1979).
The scalar field wave equation is solved using a
second order (in space and time) predictor-corrector
integration, sub-cycling within a single
hydrodynamic cycle since equation (7) generally
has a tighter restriction on the time step for stability.

A second formulation of the dynamical equations considered
here is based on a simpler conservative hyperbolic
description of the hydrodynamics equations. In this case,
the equations are derived directly by
expanding the conservation of stress-energy
Tμν,ν=0 into time and space explicit parts.
In flat space, this yields the same equations
(5) and (7) for momentum conservation
and scalar field evolution, but with

The Cosmos code uses a non-oscillatory central
difference numerical scheme Jiang et al. (1998)
with second order spatial differencing and
predictor-corrector time integration
with dimensional splitting
to solve the hydrodynamics equations in this formalism.
Artificial viscosity is not needed in this approach. Instead,
shocks are captured with MUSCL-type piece-wise linear interpolants
and high order limiters for flux reconstruction, but without the
complexity of Riemann solvers.

A final ingredient needed is an equation of state
relating pr to energy E and temperature T.
Here we assume pr=arT4, treating the
plasma as a gas of relativistic particles with enthalpy density

ρh=γprγ−1−T∂TV(ϕ,T)=4arT4−T∂TV(ϕ,T)=EW+pr−V(ϕ,T),

(11)

and adiabatic index γ=4/3.
The equation of state is then written in terms of computed quantities

peff=pr−V(ϕ,T)=(γ−1)(EW+T∂TV(ϕ,T)−V(ϕ,T))−V(ϕ,T),

(12)

or

peff=pr−V(ϕ,T)=(γ−1)(E−V(ϕ,T)+TW2∂TV(ϕ,T)γW2−(γ−1))−V(ϕ,T),

(13)

depending on which energy variable is evolved.
The quantity peff is introduced as an effective pressure
for numerical convenience in solving equations (5),
(6) and/or (9).
The fluid temperature is computed from either

EW−3arT4−V(ϕ,T)+T∂TV(ϕ,T)=0,

(14)

or

E−V(ϕ,T)+TW2∂TV(ϕ,T)−3arT4[γW2−(γ−1)]=0.

(15)

Equations (14) and (15) are solved
using a Newton-Raphson method to iterate and converge on the
temperature, assuming an initial guess of
T0=(E/3War)1/4.

Iii Perturbation Analysis

In this section we carry out a first-order perturbation expansion
of the hydrodynamics equations, and review briefly the expected
dynamical behavior and stability of phase transitions in the
context of linear theory. The general aproach follows that
presented in references Link (1992); Kamionkowski and Freese (1992); Huet et al. (1993); Abney (1994).
In addition to elucidating
any unstable behavior, the results of this section are also
used to set up inhomogeneous perturbations as initial data
that is evolved numerically in §IV.

It is convenient to start with the hydrodynamics equations in
conservative hyperbolic form (5) and
(9), which we rewrite here as

∂tE+∂i(Evi+vip)

=

Σ0,

(1)

∂tSj+∂i(Sjvi+δijp)

=

Σj,

(2)

assuming a constant scalar field potential on either side
of the phase transition front, and
Σ0 and Σj represent any source or dissipative terms.
Equations (1) and (2) are expanded
out to first perturbative order assuming
two dimensional perturbations off a homogeneous background
flow of the form

p

=

p0+δp(t,x,y),

(3)

v

=

v0ˆi+[δvx(t,x,y)ˆi+δvy(t,x,y)ˆj],

(4)

for the fluid pressure and velocity.
The corresponding boost factor and conserved variables
are, to first order

W

=

1√1−v2=1√1−v20+v0(1−v20)3/2δvx,

(5)

E

=

p(γW2−γ+1)γ−1=−p0+γp0(γ−1)(1−v20)

(6)

+2γp0v0(γ−1)(1−v20)2δvx+γ(γ−1)(1−v20)δp−δp,

Sx

=

ρhvx1−v2=γ(γ−1)(1−v20)(p0v0+p0δvx+v0δp+2p0v201−v20δvx),

(7)

Sy

=

ρhvy1−v2=γp0(γ−1)(1−v20)δvy,

(8)

where ρh=(e+p)=γp/(γ−1), and v0
is the unperturbed fluid velocity.

Substituting (5) - (8)
into (1) and (2) yields
for momentum conservation along the y-axis

∂y(δp)+γp0(γ−1)(1−v20)[∂t(δvy)+v0∂x(δvy)]=Σy,

(9)

for momentum conservation along the x-axis

(γ−1)(1−v20)Σx

=

γv0∂t(δp)+γp0(1+v20)(1−v20)−1∂t(δvx)

(10)

+

(γ−1+v20)∂x(δp)+2γp0v0(1−v20)−1∂x(δvx)+γp0v0∂y(δvy),

and for energy conservation

(γ−1)(1−v20)Σ0

=

(1−v20+γv20)∂t(δp)+2γp0v0(1−v20)−1∂t(δvx)

(11)

+

γv0∂x(δp)+γp0(1+v20)(1−v20)−1∂x(δvx)+γp0∂y(δvy),

all to first perturbative order.

Equations (10) and (11)
can be simplified further by first eliminating
∂y(δvy) from the x-momentum equation

0

=

(γ−1)(1−v20)(Σx−v0Σ0)

(12)

=

v0∂t(δp)+ρh0W20∂t(δvx)+∂x(δp)+W20ρh0v0∂x(δvx),

then eliminating ∂t(δvx) from the energy equation

0

=

c2s(γ−1)(1−v20)[Σ0−(2v01+v20)Σx]

(13)

=

(1−v20c2s)∂t(δp)+v0(1−c2s)∂x(δp)+ρh0[∂x(δvx)+∂y(δvy)].

In writing (12) and (13) we
introduced the notation c2s=γ−1, W20=1/(1−v20)ρh0=γp0/(γ−1), and
explicitly ignored the source terms
Σα by setting them to zero.

The final equations (9), (12), and
(13) can be written conveniently in compact form as

At∂tW+Ax∂xW+Ay∂yW=0,

(14)

where W=(δp,δvx,δvy), and

At

=

⎛⎜
⎜⎝1−v20c2s00v0ρh0W20000ρh0W20⎞⎟
⎟⎠,

(15)

Ax

=

⎛⎜
⎜⎝v0(1−c2s)ρh0c2s01W20ρh0v0000ρh0W20v0⎞⎟
⎟⎠,

(16)

Ay

=

⎛⎜⎝00ρh0c2s000100⎞⎟⎠.

(17)

Next, following the general procedure outlined
in reference Abney (1994), we assume a solution of the form
W(t,x,y)=w(x)e−i(ωt+ky), with
w(x)=∑aje−iλ∗jxRj
and the following characteristic equation or dispersion relation
for (14):

(ω+λ∗v0)[k2(1−v20)−c−2s(ω+λ∗v0)2+(λ∗+ωv0)2]=0.

(18)

Equation (18) has either three distinct roots
which agree with those found by Huet et al. (1993); Abney (1994)

where we define Ω=[ω2+k2(v20−c2s)/(1−v20)]1/2
in equations (24) and (25).

To determine whether any unstable modes of W exist
and are consistent with imposed boundedness conditions
at x→±∞, it is convenient to
construct a coordinate system centered and moving with
the frame of the interface between quark and hadron phases
at x=0. Assuming the high temperature quark phase is to the left
of the interface (x<0), and the hadron phase to
the right (x>0), we can easily determine
if any unstable modes (defined by Im(ω)>0)
obey the conditions W→0 as
x→±∞. As concluded in references
Huet et al. (1993); Abney (1994), if Im(ω)>0 in the quark region
of a detonation front,
then Im(λ∗j)<0 and therefore
aj=0 for all j
in order for the solutions
to be bounded at x→−∞. Hence
unstable modes do not exist in regions ahead of nucleating
bubbles separated from the quark phase by a detonation front.
However, in the hadron region to the right (x>0) of
a detonation front,
the boundedness conditions at x→+∞
do not rule out unstable modes to linear order. Also, in the case
of deflagration fronts, unstable modes are
realizable on both the quark and hadron sides, at least for
cases in which surface tension, dissipation,
and strong thermal coupling effects
are negligible Link (1992); Huet et al. (1993); Kamionkowski and Freese (1992).
It is not entirely clear what role these effects
play in stabilizing the transition, and certainly
higher order nonlinear effects are not accounted for in this analysis.
Table 1 summarizes which of the regions and hydrodynamic
states may potentially give rise to unstable modes with
Im(ω)>0, at least to linear order.

Quark region (x<0)

Hadron region (x>0)

detonations:

detonations:

vq>cs

vh<vq

a1=0

λ∗1=−ω/v0

a−=0

a−=0

a+=0

Im(λ∗+)<0

deflagrations:

deflagrations:

vq<cs

vh>vq

a1=0

λ∗1=−ω/v0

Im(λ∗−)>0

a−=0

a+=0

Im(λ∗+)<0

Table 1:
Summary of parameters yielding bound solutions that are
susceptible to unstable modes. The entries represent a general
set of conditions that allow unstable modes (Im(ω)>0)
to first perturbative order for each of the phases and front types.
These conditions are used to construct initial
data for the multi-dimensional numerical stability studies presented
in section §IV.

Iv Numerical Results

iv.1 Initial Data

The initial data has just
one free dimension, energy, which we use to specify the
critical transition temperature Tc and normalize it to unity.
It is also convenient to define a characteristic
length scale as the radius of a spherical nucleated
hadron bubble in approximate equilibrium
with an exterior quark plasma. This is estimated by balancing
pressure forces acting from the quark and hadron sides and including
the effect of surface tension, resulting in
Req=2σ/(ph−pq)
for the equilibrium radius,
where σ is the surface tension, and
ph and pq are the pressure fields from the hadron and
quark sides respectively. For the preliminary one-dimensional
calculations in §IV.2, the computational
box sizes are set to a multiple of the equilibrium
bubble radius L=xReq, with
x generally in the range 20 - 2000.
Typical length scales or box sizes
in these calculations vary from a few hundred to a few thousand fermi.
By comparison, the proper horizon size at the time of
the QCD phase transition is

dH=2tage=R(t)∫t0dt′/R(t′)=16√g∗/51.25(TMeV/150MeV)2km,

(1)

or 16 km for
g∗=51.25 and TMeV=150MeV. In writing (1)
we have assumed
R(t) is the cosmological scale factor in the isotropic FLRW
background model dominated by radiation energy density of the
form e=arT4=e0R−2 such that

tage=√3c232πGarT−2=2.7×10−5√g∗/51.25(TMeV/150MeV)2s

(2)

is the age of the Universe as a function of temperature.

We consider two background temperatures
for the super-cooled initial state triggering the shock
and phase front propagation:
T=0.9Tc and T=0.9943Tc.
The latter temperature is chosen to match
the 1+1 dimensional cases considered in reference Ignatius et al. (1994),
and provides a useful benchmark
against which some of our results can be compared.

To keep the number of numerical simulations down to
a reasonable level, we fix the following additional parameters
for all calculations and for both temperature cases:
critical transition temperature Tc=150 MeV corresponding
to the QCD symmetry breaking energy occurring
when the Universe was approximately
tage=2.7×10−5 s old (see equation
(2));
particle degrees of freedom g∗=51.25 to be consistent
with previous studies Kajantie and Kurki-Suonio (1986);
surface tension σ=0.1Tc;
correlation length ℓc=1.0/Tc; and
latent heat L=2T4c.
Note that for both temperatures
considered here, the interaction potential (equation (2)
and Figure 1) has two local minima
(at ϕ=0 and ϕ=ϕmin(T)>0, which allows
for the quark (cold phase with ϕ=0) and hadron
(hot phase with ϕ=ϕmin(T)) states to co-exist.

The remaining parameter, the dissipation constant
κ, is varied over the range
0.1≤κ≤10 in the one-dimensional calculations
to determine the regimes in which detonations
and deflagrations are triggered, and to compute the
stability parameters from linear theory.

¯¯¯λ=9λβ(T2−T20)/(2α2T2), and
x=L is the right-most edge along the x-axis representing
the computational box size.
These expressions define a hadron phase with ϕ>0 at L−x→0,
and quark phase with ϕ→0 as L−x→∞,
with an exponentially steep interface boundary of width
set by the model parameters.
The numerical constants ϕc and
Δ=0.05 are used to perturb the solution
out of dynamical equilibrium by varying both the amplitude and
width of the field in the hadron phase.
For most of the calculations we set ϕc=1. However, for the
larger κ cases we found it useful to
specify ϕc such that ϕ(x=L)≈ϕmin(T), in order
to reduce the transient interval between the initial configuration
and the time that the phase front fully develops.

iv.2 One-Dimensional Results

The primary motivations for the one-dimensional calculations
presented here were to narrow the parameter space and help
illuminate the two-dimensional studies presented in the following
section for the different background temperatures considered.
In cases where the temperature was initialized to T=0.9Tc, we used anywhere
from 2×103 to 3.2×104 zones to
resolve a length scale of L=2000Req=2190 fm.
For runs in which T=0.9943Tc we
used from 5×102 to 4×103 zones to
resolve L=20Req=452 fm. In these 1D parameter studies,
we also considered an initial background temperature of T=0.95Tc.
These runs were useful for characterizing some of the intermediate states
in the 2D runs below.
For these T=0.95Tc runs we used from 1×103 to 4×103
zones to
resolve L=250Req=1157 fm.
The higher resolutions in each temperature case
were generally required for the
larger values of κ in order to maintain reasonable zone coverage
in the hadron phase, since the bubble growth velocity was sometimes
substantially less than the shock front velocity, or near the transition
from a detonation front to a deflagration front, as detailed below.
For this first
group of 1D runs, the code was allowed to
run for approximately one sound-crossing time ts=L/cs,
where cs=1/√3.

Figure 2 shows the phase front or bubble growth velocity
normalized to the sound speed vf/cs,
and the quantity η=−Tcvfdvf/dTq
as a function of the dissipation parameter κ for
temperatures T/Tc=0.9, 0.95, and 0.9943.
According to reference Huet et al. (1993), η determines the
stability of bubbles and we compute it using the
approximation η=v2f/v2c, where
v2c=(1/2)(1−T2q/T2c).
Notice the T=0.9Tc curve switches from a detonation to a
deflagration at κ≈1 where
vf/cs crosses unity, and the T=0.95Tc curve switches at
κ≈0.3. However, the precise transition points are
highly sensitive to the resolution and duration of
the simulations, and what might appear
as a deflagration at one resolution may
turn out to be a detonation at another.
This is attributed to the solutions approaching
either a Jouget or temporary strong detonation state
(which eventually decays into a weak deflagration) as the
parameter κ is increased. In particular,
the post-shock velocity for the T=0.9Tc case
computed over short time intervals
in the frame of the front decreases from being supersonic
at low values of κ, to about sonic speed at κ=0.875,
then subsonic at higher values but less than unity (e.g., κ=0.9),
corresponding to weak, nearly Jouget, and strong detonations, respectively.
Solutions for which κ is close to unity thus represent a regime
in parameter space that does not allow stable detonation modes since strong
detonations are forbidden modes of propagation
(see, for example, Kurki-Suonio and Laine (1995)).
We note that the conservative hyperbolic form of the hydrodynamics
equations and
the non-oscillatory central difference methods for solving the
equations generally out-perform artificial viscosity methods in resolving
and maintaining stable evolutions near the critical Jouget state.

Figure 2: Bubble growth velocity normalized to the sound speed
vf/cs, and stability parameter η as a function
of the dissipation constant κ for temperatures
T/Tc=0.9, 0.95, and 0.9943. The T=0.9Tc (T=0.95Tc)
curve switches from a detonation to a deflagration at κ≈1
(κ≈0.3).

Table 2 summarizes our results for the potentially unstable,
η<1, deflagration runs.
Th and Tq in Table 2 are the
average temperatures of the hadron and quark phases, respectively;
and vh and vq are the hadron and quark velocities
on either side of the phase front
as measured in the rest frame of the moving front.
The critical wave number, kc=(μ−1)ρhqv2q/σ, defines the
limit above which the system
is stabilized by surface tension Link (1992), where μ=vh/vq and
ρhq is the enthalpy density of the quark phase.
Modes with wavelengths λ<λc=k−1c
are thus expected to be stable.
The instability growth timescale, τ, is estimated as Link (1992)

Table 3 summarizes the results computed for the
detonation cases with temperature T=0.9Tc.
Here vf represents both the front and shock velocities,
and ρhh is the enthalpy density
of the hadrons immediately behind the phase front.
The quantity σ/ρh has dimensions of length,
and provides a natural scale that gauges the relative
importance of surface tension.
It is also applied as a dimensional scaling variable in the calculations
of reference Abney (1994), which we use as a guide to
estimate perturbation growth rates and physical run
times for the two-dimensional calculations
presented in §IV.3.

T/Tc

κ

|vf/cs|

(ρhh)/σ

0.9

0.1

1.85

191.7

0.3

1.71

195.7

0.5

1.57

201.1

0.75

1.36

216.2

0.8

1.33

220.4

Table 3:
Front velocity vf and enthalpy to surface tension
ratio (ρhh)/σ on the hadron side
as a function of dissipation constant κ for
a selected set of 1D detonations
with initial temperature T=0.9Tc

In Figures 3-5, we present
more detailed results of three particular 1D runs that will help illuminate
the 2D calculations with similar parameters presented in the next section.
Figure 3
shows a series of outputs of the scalar field and temperature for a
weak (subsonic post-front velocity relative to the phase front)
deflagration case with T=0.9Tc and κ=7. This run
used 2000 zones to resolve a length scale of L=1095 fm, and was
allowed to run for slightly longer than two sound-crossing times so that
the shock leading the deflagration front propagates across the grid twice.
These tiled displays show a leading shock front
(in the temperature graphs) originating
from the right and traveling to the left with velocity
slightly higher than sound speed vs=1.03cs,
heating up the quarks just ahead of the cool hadron
phase (first row). The deflagration front (separating the two regions
with different scalar field values) travels to the left at a
speed slower than the leading shock front, vf/cs≈0.34.
Due to the reflection boundary conditions imposed on all our
calculations, a mirrored shock collision occurs when
the leading shock hits the left boundary (second row).
The reflected shock then travels
to the right, eventually colliding with the phase front
moving in the opposite direction (third row).
The shock heats up the fluid
and reduces the scalar field as it passes through
the hadron phase. Upon impact with the phase front, the
shock/front interaction also generates a rarefaction wave
traveling to the left which cools off the quarks in its path.
Finally the hadron matter at the right end is heated further as the
reflected shock hits the right boundary and collides with another incoming
shock (fourth row).

Figure 3: Scalar field (left column)
and temperature (right column) as a function of position
for a 1D deflagration
case with T=0.9Tc and κ=7. Top row: Early undisturbed
deflagration and shock fronts at t=3.2×10−21 s. Second row:
Shock collision in the quark phase at the left boundary at t=6.6×10−21 s.
Third row: Reflected shock interacting with deflagration front at
t=1.1×10−20 s. Fourth row: Interaction of two reflected
shocks in the hadron phase at the right boundary at t=1.4×10−20 s.

Figure 4 shows a similar series of displays
for a weak (supersonic post-shock velocity relative to the phase front)
detonation case with T=0.9Tc and κ=0.5,
run for almost two sound-crossing times. This run
used 2000 zones to resolve a length scale of L=1095 fm.
The outputs begin in the first row
with an early undisturbed state
showing the leading shock/detonation front originating at the
right end of the grid and traveling to the left
at speed vf/cs≈1.6. A rarefaction wave follows
the shock and separates
two hadron regions at two slightly different temperatures.
The second row shows the shock-shock collision which
occurs when the leading
shock front reaches and reflects off the left boundary. This interaction
heats up the fluid to temperatures higher than Tc, and
generates a quark region at the left boundary.
As the reflected shock travels to the right, the quark region grows.
The third row shows the reflected shock passing
through the rarefaction wave separating the two hadron
states. This rarefaction wave cools the newly formed
quark region, converting quarks back into hadrons.
The fourth row shows
the interaction of reflected rarefaction waves at the left boundary.
The complete reflected state
showing the thermal distortion and shock profiles after the
entire domain is converted to the hadron phase is displayed in the fifth row.
This compact structure is allowed to traverse the grid and interact
through each sequence again at
the right boundary (sixth row).
The complete reflected state from the right
boundary is shown in the final row.

Figure 4: Scalar field (left column) and temperature (right column)
as a function of position
for a 1D detonation
case with T=0.9Tc and κ=0.5. Top row: Early undisturbed
detonation front plus rarefaction at t=2.2×10−21 s.
Second row:
Shock-shock interaction at left boundary at t=4.8×10−21 s.
Third row: Reflected shock passing through the rarefaction wave at
t=5.7×10−21 s. Fourth row: Interaction of the rarefaction
waves at left boundary at t=6.6×10−21 s. Fifth row:
Reflected state from left boundary at t=9.7×10−21 s.
Sixth row:
Interaction of reflected states at right boundary at t=1.2×10−20 s.
Seventh row: Reflected state from right boundary at t=1.5×10−20 s.

Figure 5 shows outputs
from two phase front systems originating in separate
regions at different temperatures and are thus out of
initial thermal equilibrium: a deflagration front at T=0.9943Tc
at the left end, and a detonation front at T=0.9Tc at the right end.
This run
used 2000 zones to resolve a length scale of L=1806.1 fm.
The dissipation constant is κ=0.5 across both regions, and the
temperature discontinuity is initiated at the center of the grid.
The right end is initialized to the same state as the
previous detonation case in Figure 4.
The left end is similar to the weak deflagration case of
Figure 3, but
with a smaller dissipation constant and a higher initial temperature.
The first row illustrates an early stage where the deflagration
front has formed at the left end with vf/cs≈0.14,
the shock/detonation front and rarefaction wave have formed at
the right end with vs=vf≈1.6cs,
and a thermal shock front and rarefaction wave
are traveling in opposite directions from the
central thermal discontinuity at about the sound speed.
The second row shows the collision of
the detonation and thermal shocks.
The third row shows the passage of the thermal shock
through the detonation rarefaction, and the collision of
the deflagration shock with the thermal rarefaction wave.
The fourth row shows the interaction of the
deflagration front with the thermal rarefaction at the left
end, and the collision of the deflagration and detonation shocks
at the middle.
As the rarefaction wave passes through the hadron phase behind
the deflagration front at the left end, the phase transition is accelerated
to near sound speed.
The fifth row shows the accelerated shock and deflagration front features
traveling to the right at velocities vs≈0.67=1.16cs
and vf≈0.51=0.88cs.
Notice that the accelerated front velocity is close to, but slightly
faster than the deflagration velocity predicted by the
T=0.95Tc curve of Figure 2, which
yields vf=0.74cs=0.43.
The sixth row shows the interaction of the newly accelerated front with
the detonation propagating from the right end. This interaction
heats up and decomposes the hadrons back into quarks over
a substantially large region at the center.
The newly formed quark region gradually cools with the
shock passage and rarefaction wave interactions.
Profiles and distributions of quark and hadron regions
become increasingly complex as indicated in the final
row due to the various shock,
phase, and sound wave interactions as they reflect off the
boundaries and propagate through the grid.

Figure 5: Scalar field (left column) and temperature (right column)
as a function of position
for a 1D interaction of shocks and phase fronts
originating in two separate regions at different temperatures:
a deflagration system at T=0.9943Tc at the left end,
and a detonation configuration at T=0.9Tc at the right end.
The dissipation constant is
κ=0.5 everywhere, and the temperature discontinuity
is initiated at the grid center at t=0. Top row: Early undisturbed
deflagration and shock fronts (left), detonation front plus rarefaction
(right), and
thermal shock front plus rarefaction (center)
at t=1.8×10−21 s. Second row:
Interaction of detonation front and thermal shock at t=2.2×10−21 s.
Third row: Interaction of the detonation rarefaction and thermal
shock, and thermal rarefaction and deflagration shock at
t=3.1×10−21 s. Fourth row: Interaction of deflagration
front and thermal rarefaction, and deflagration shock and
detonation shock at t=4.8×10−21 s. Fifth row:
Acceleration of deflagration front by thermal rarefaction
at t=5.3×10−21 s.
Sixth row: Interaction of newly accelerated front and
detonation front at t=7.5×10−21 s. Seventh row: Interactions of multiple shocks resulting
in two separate quark regions at t=1.1×10−20 s.

We point out that, in general, Figure 2 cannot always
be used to predict reliably the accelerated velocity or
mode of propagation from rarefaction wave/front interactions.
For example, consider the same initial configuration used in producing
Figure 5, but with friction parameter κ=0.2.
In this case Figure 2 predicts, for a mean field temperature
T=0.95Tc, a supersonic front velocity of vf=1.56cs.
However, the actual (numerical) accelerated solution remains
a weak deflagration with subsonic velocity vf=0.54=0.94cs, only
slightly faster than for the κ=0.5 case.
This holds even in the limit κ→0:
the accelerated front velocity approaches but does not exceed
the sound speed as κ→0, and the propagation mode
remains a deflagration with clearly separated phase and shock front features.

iv.3 Two-Dimensional Results

We focus here on extending to two dimensions several of the
more interesting
parameter combinations presented in section §IV.2
that may potentially give rise to unstable behavior as
predicted by linear theory. In particular, we consider
six calculations: runs A and B are single deflagrating
fronts; run C is a single detonation transition;
runs D and E are interactions between a planar
and smaller spherical nucleating bubble of
deflagration (run D) and detonation (run E) types;
and run F is the interaction of a deflagration system with a detonation
front nucleating from two regions out of thermal equilibrium.
Table 4 summarizes
the parameters used in each of these runs.

Run

T/Tc

κ

λc

τ

Grid Dimensions

Stop Time

(fm)

(s)

(fm)

(zones)

(s)

A

0.9

7.0

1.2

1.3×10−21

278.03×6.14

1280×64

3.2×10−21

B

0.9943

1.5

31.

2.6×10−19

3059.83×152.99

1280×64

7.1×10−20

C

0.9

0.5

6.6

2.2×10−22

302.34×32.71

1024×256

3.5×10−21

D

0.9

7.0

1.2

1.3×10−21

1000×1000

512×512

1.2×10−20

E

0.9

0.5

6.6

2.2×10−22

1000×1000

1024×1024

1.2×10−20

F

0.9943/0.9

0.5

…

…

1806.1×451.53

2048×512

2.1×10−20

Table 4: Summary of parameters for the two-dimensional runs.

For all of these simulations we initialize the data with various
perturbations included. First, the planar fronts in each problem are
perturbed with a sinusoid of wavelength λ≈5λc,
where λc is the critical wavelength calculated from our 1D
studies. The amplitude of this perturbation is δx/x=0.05,
and the wavelength λ is used to set the
physical size of the grid parallel to the phase front.

Transverse and longitudinal
inhomogeneous fluctuations are also imposed on all the
initial data using the perturbation solutions discussed
in §III.
In particular, E=E0+δE, vx=δvx,
and vy=δvy, with
E0=3arT4+V(ϕ,T)−T∂TV(ϕ,T),
where

⎛⎜⎝(γ−1)δ(E/W)δvxδvy⎞⎟⎠=e−i(ωt+ky)∑jaj˜RjeIm(λ∗j)xe−iRe(λ∗j)x,

(7)

λ∗j are the eigenvalues,
and ˜Rj=Rj/|Rj|
are the normalized eigenvectors. The
expansion coefficients are defined as
aj=min(Avo,AE0) with
amplitude constant A=0.1 (if they are not set
to zero because of boundedness constraints as discussed
in §III).
v0 is the background average velocity in either
the quark or hadron regions as defined in
§III, but implied here to be
the velocity component orthogonal
to the interface boundary and measured in the rest frame of
the unperturbed surface. We take
v0=cs+0.5(1−c