# Numerical Relativity as a Tool for Studying the Early Universe

Authors: David Garrison
Journal: Journal of Gravity
 (2014)

---

## Abstract

Numerical simulations are becoming a more effective tool for conducting detailed investigations into the evolution of our universe. In this paper, we show how the framework of numerical relativity can be used for studying cosmological models. The author is working to develop a large-scale simulation of the dynamical processes in the early universe. These take into account interactions of dark matter, scalar perturbations, gravitational waves, magnetic fields, and  turbulent plasma. The code described in this report is a GRMHD code based on the Cactus framework and is structured to utilize one of several different differencing methods chosen at run-time. It is being developed and tested on the University of Houston’s Maxwell cluster.


---

## Body

## 1. Introduction

Our knowledge of how the universe evolved comes primarily from observations of large structures such as stars, galaxies, clusters, and super-clusters of galaxies as well as from observations of the cosmic microwave background (CMB) radiation. Based on these observations, the standard model of cosmology was developed during the mid to late twentieth century. Some elements of this model include the existence of primordial metric perturbations, magnetic fields, and an early universe filled with nearly homogenous and isotropic plasma [1]. The perturbed Friedmann-Robertson-Walker (FRW) metric, which describes the spacetime curvature of the early universe, takes the following form:
(1)
d
s
2
=
-
(
1
+
2
ϕ
)
d
t
2
+
ω
i
d
t
d
x
i
+
a
(
t
)
2
[
(
(
1
+
2
ψ
)
δ
i
j
+
h
i
j
)
d
x
i
d
x
j
]
.Herea
(
t
) is the scale factor and ϕ, ψ, ω
i, and h
i
j are the scalar, vector, and tensor perturbation terms. Many cosmological models relate density fluctuations and variations in the CMB to perturbations in the FRW metric at the time of recombination. These perturbations start off small and grow as a power-law with time as the competing forces of universal expansion and gravitational attraction affect their growth [1]. Work by Kodama and Sasaki [2], Sachs and Wolfe [3], and Mukhanov et al. [1] all showed analytically how metric perturbations could cause density perturbations in a hydrodynamic fluid.Recently, beyond the standard model cosmological theories, [4] has suggested that primordial fluids and fields are potential sources of observable gravitational waves. This realization opens up exciting new possibilities, making primordial gravitational radiation an important source of information about the early universe. Taken together, this leads to the idea that there was a dynamical interaction between matter, electromagnetic, and gravitational fields in the early universe that affected the evolution of our universe. Signatures of these interactions may still be observable today. The objective of this paper is to show how the tools of numerical relativity can be used to study such an interaction and to introduce a computer code written for this purpose.This project uses the framework of numerical relativity to develop a computational laboratory to study the evolution of the early universe. The initial focus of this work is on deriving the spectrum of gravitational waves produced by relativistic turbulence in the early universe. Future work may involve a more advanced study of how this gravitational wave spectrum is affected by the presence of dark matter or a preexisting primordial gravitational wave field. Numerical relativity has been used for years to study the collisions of compact objects such as black holes and neutron stars and to predict the spectrum of the gravitational waves produced by their interactions. We propose to use this tool to provide detailed studies of cosmological events that may also one day be observable using either gravitational or conventional astronomy. In cosmology, numerical simulations are capable of providing more detail than the analytic calculations that have been performed to date. For example, a modern general relativistic magnetohydrodynamic (GRMHD) code is capable of performing such simulations by evolving both the spacetime and plasma field dynamically and can therefore represent chaotic processes such as turbulence more accurately than analytic calculations. This paper is a summary of the techniques used to develop such a simulation.


## 2. Development of Initial Conditions

Every numerical simulation consists of three parts: initial data, numerical evolution, and data analysis. In this section, we will focus on the initial data needed for cosmological studies.


### 2.1. Background and Hypotheses

Several different mechanisms for producing primordial gravitational waves have been identified and studied by researchers. These include quantum fluctuations during inflation, bubble wall motion, and collisions during phase transitions, cosmological magnetic fields, oscillating classical fields during reheating, cosmological defects, and plasma turbulence [4–14]. We assume that the early universe was in a metastable state during a first order phase transition. The false vacuum was separated from the true vacuum by a potential barrier or a scalar field. Quantum tunneling occurred across the barrier in finite regions of space resulting in true vacuum bubbles inside the false vacuum phase. As the universe expanded and cooled, the energy difference between the false vacuum and true vacuum got larger, making the phase transition more probable. Eventually, the probability of nucleating one critical bubble per Hubble time became high enough to cause the phase transition to begin. This defined the transition temperature, which is believed to be about 1 TeV [4]. The nucleated bubbles expanded and collided, eventually filling the whole universe. The collision of two or more of these bubbles broke spherical symmetry and released some of their energy as gravitational waves. Since the expansion of the bubbles was accompanied by macroscopic motions in the cosmic matter field, the collision of these bubbles also resulted in the anisotropic stirring of the field. This caused turbulent motions which provided a primary source of gravitational waves for this research.


#### 2.1.1. Primordial Magnetic Fields

Magnetic fields are believed to have played a large part in the dynamics of the universe’s evolution. Little is known about the existence of magnetic fields in the early universe. There are no direct observations of primordial magnetic fields. Theories also disagree on the amplitude of primordial magnetic fields. There are currently several dozen theories about the origin of cosmic magnetic fields [15, 16]. The main reason that we believe that primordial magnetic fields existed is because they may have been needed to seed the large magnetic fields observed today. Most theories of cosmic magnetic field generation fall into one of three categories [15–17]: (1) magnetic fields generated by phase transitions; (2) electromagnetic perturbations expanded by inflation; and (3) turbulent magnetofluid resulting in charge and current asymmetries.Most models calculate the magnitude of primordial magnetic fields by starting with the observed strength of galactic or intergalactic magnetic fields and calculating how this field should have been amplified or diffused by external effects such as the galactic dynamo and expansion of the universe [15, 16]. A major problem is that there does not appear to be a universal agreement of how efficiently a galactic dynamo could have strengthened seed magnetic fields. Estimates of the strength of these seed fields can vary by tens of orders of magnitude. Seed magnetic fields produced during inflation are predicted to have a current strength somewhere between 1
0
-
11 G and 1
0
-
9 G on a scale of a few Mpc [15, 16, 18]. Magnetic seed fields generated by phase transitions are believed to be less than 1
0
-
23 G at galactic scales [15, 16]. Some turbulence theories imply that magnetic fields were not generated until after the first stars were formed therefore requiring no magnetic seed fields [15].Given how little is understood about primordial magnetic fields and the general lack of agreement among theoretical predictions, it seems clear that the existence of primordial magnetic fields can neither be confirmed or ruled out. It seems that the best we can do is set an upper limit on the strength of primordial magnetic fields and utilize this limit as a starting point in developing models of cosmic turbulence. Observations of the CMB limit the intensity of the magnetic seed fields to a current upper limit of1
0
-
9 G [15, 16, 18, 19].It is well known that gravitational waves can interact with a magnetofluid in the presence of a magnetic field. Work by Duez et al. [20] showed how gravitational waves can induce oscillatory modes in a plasma field if magnetic fields are present. Work by Kahniashvili and others [4, 9–13] has shown how a turbulent plasma can yield gravitational waves. The result may be a highly nonlinear interaction as energy is transferred from the fluid to the gravitational waves and back.


#### 2.1.2. Turbulence in the Early Universe

Turbulence provides a particularly interesting GW source because it is not well understood analytically. This turbulence is a natural result of dynamics of the early universe resulting from bubble wall collisions and other chaotic events during the first order phase transitions. Analytic work done to date [4, 5, 12–14] summarizes the dynamics of the phase transitions using two quantities, α and β. α is traditionally defined as the ratio of false vacuum energy and plasma thermal energy density. This provides a measure of the transition strength. If α is much less than one, the transition is very weak. If α is larger than unity, the transition is very strongly first order. β is the rate of variation of the nucleation rate at the transition time. It fixes the time scale of the phase transition once the transition has begun. After a time interval β
-
1, the whole universe is converted to a true vacuum phase. Therefore the turbulent stirring should only last β
-
1.The amount of gravitational waves emitted by bubble collisions and turbulence generated in the plasma are also determined from two quantities,κ and v
b. κ is the fraction of vacuum energy transferred into fluid kinetic energy and v
b is the velocity of bubble wall expansion. Bubble walls can propagate via two modes, detonation and deflagration [4, 5, 12–14]. For detonation, the bubble walls are thin compared to the radius and they propagate faster than the speed of sound. This results in
(2)
v
b
(
α
)
=
1
/
3
+
(
α
2
+
2
α
/
3
)
1
/
2
1
+
α
,
κ
(
α
)
=
1
1
+
0.715
α
[
0.715
α
+
4
27
3
α
2
]
.
If the bubbles propagate by deflagration, the walls are thick and have a lower energy density. It is currently believed that for a relativistic plasma, the deflagration expansion mode is unstable so only the detonation modes will result.The number density of turbulent eddies within a Hubble radius should depend onv
b and β. The characteristic velocity perturbation of the turbulent fluid for the largest eddies at the stirring scale is given by
(3)
v
0
=
3
κ
α
4
+
3
κ
α
.
For κ
α
≈
1, which corresponds to a strongly first order phase transition, v
0 is about 0.65 at the time of the electroweak phase transition. We later use v
0 as the maximum velocity of fluid elements in our studies. This velocity is randomized in amplitude and direction in order to simulate the initial conditions for turbulence.


### 2.2. Computational Model

In order to study the interaction of the plasma field and the background spacetime dynamically (or separately) our team has written and is testing/improving a GRMHD code implemented using the Cactus framework [21]. We now describe the basic variables and equations that constitute this model.


#### 2.2.1. The Spacetime Evolution Model

The spacetime metric can be written as(4)
d
s
2
=
-
N
2
d
t
2
+
γ
i
j
(
x
→
,
t
)
(
d
x
i
+
N
i
d
t
)
(
d
x
j
+
N
j
d
t
)
.
Here N is the lapse, N
i is the shift vector, and γ
i
j is the spatial 3-metric [22, 23]. For this work, 3-metric and its “time-derivative,” the extrinsic curvature, “K
i
j” will be evolved using a strongly hyperbolic version of the BSSN formulation of numerical relativity [24].


#### 2.2.2. The General Relativistic Magnetohydrodynamic Model

The fluid and electromagnetic fields of the GRMHD equations are developed from several well-known equations [25]. They include the conservation of particle number, the continuity equation, the conservation of energy-momentum, the magnetic constraint equation, and the magnetic induction equation. For a system consisting of a perfect fluid and an electromagnetic field, the ideal MHD stress-energy tensor is given by
(5)
T
μ
ν
=
(
ρ
0
h
+
b
2
)
u
μ
u
ν
+
(
P
+
b
2
2
)
g
μ
ν
-
b
μ
b
ν
,
h
=
1
+
ϵ
+
P
ρ
0
,
b
μ
=
1
4
π
B
(
u
)
μ
,
B
(
u
)
0
=
1
α
u
i
B
i
,
B
(
u
)
i
=
1
u
0
(
B
i
α
+
B
(
u
)
0
u
i
)
.
Here, P is the fluid pressure, ρ
0 is density, B
i is magnetic field, u
μ is four-velocity, h is the enthalpy, ϵ is specific internal energy, and b
2 is the magnitude of the magnetic vector field squared. The addition of viscosity modifies the MHD stress-energy tensor by incorporating the viscous stress tensor:
(6)
T
μ
ν
=
(
ρ
0
h
+
b
2
+
Q
)
u
μ
u
ν
+
(
P
+
b
2
2
)
g
μ
ν
-
b
μ
b
ν
+
Σ
μ
ν
.
Here Q is artificial bulk viscosity and Σ
μ
ν is the viscous stress tensor for artificial shear viscosity. Artificial viscosity is used here as a way of handling shocks although we are working on integrating more advanced HRSC techniques. Our viscosity terms [26] are described and defined below and are only meant to be used when the divergence of the fluid flow is negative:
(7)
Q
=
I
n
Δ
l
∂
k
V
k
(
k
q
Δ
l
∂
k
V
k
-
k
l
C
s
)
,
Σ
j
i
=
I
n
Δ
l
(
k
q
Δ
l
∂
k
V
k
-
k
l
C
s
)
Sym
(
δ
j
V
i
-
∂
k
V
k
3
δ
j
i
)
,
I
n
=
(
ρ
+
ρ
ϵ
+
(
P
+
Q
+
b
2
)
)
N
.
In these equations V
i is the fluid velocity, C
s is the local speed of sound, Δ
l is the minimum covariant zone length, and k
q and k
l are constants multiplying the quadratic and linear contributions, respectively. The Sym(⋯) function in the shear viscosity equation is a symmetry operation.Dark matter can be added to the system using a two-fluid approach where the stress energy tensor for dark matter is added directly to the stress energy tensor for the magnetofluid, therefore, completing the right-hand side of Einstein’s equation.


#### 2.2.3. Initial Conditions: Plasma Field

In order to get a more accurate picture of the primordial gravitational wave spectrum, we directly simulate the turbulent primordial universe with the most realistic initial conditions possible. These should include not only a plasma field with a realistic EOS, they should also include elements such as magnetic fields, dark matter, and possibly even gravitational waves produced by other sources.The study will begin att
>
1
0
-
6 seconds after the Big Bang near the beginning of the Hadron epoch, when the primordial plasma field began to look like relativistic plasma and the strong force could be safely ignored. At this point the Debye length is about 1
0
-
16 m so even a small computational domain should demonstrate the dynamics of the plasma. The plasma at this time was composed mainly of electrons, positrons, neutrinos, and photons. Many of the initial conditions at this epoch are well known or can be fairly easily calculated using the available literature [27]. In addition, given that most cosmological models agree that over 80% of the matter in the universe is composed of dark matter, our cosmological simulations can also include some nonmagnetized “dark” fluid. We can take this into account by adding a second pressureless nonmagnetized fluid to our initial plasma field.The initial matter field will be taken to be homogenous and turbulent. We introduce turbulence into the system by randomly varying the initial velocities of the fluid elements up to the magnitude ofv
0, (3). In addition, we will be working to better establish the initial conditions resulting from the first order EWPT. Based on arguments in previous work [15–17, 28, 29], the initial magnetic field during this epoch should be less than or equal to 1
0
17 G.


#### 2.2.4. Initial Conditions: Spacetime

For this computational study the initial spacetime is constructed in such a way as to mimic the conditions present during the Hadron epoch (t
>
1
0
-
6 seconds). I choose to begin the simulation during the Hadron epoch because at that time the primordial plasma field appeared to look like a relativistic plasma field that could be modeled using a GRMHD code as opposed to a quark gluon plasma field. At this time the strong force could be safely ignored as electromagnetic effects would dominate the plasma’s dynamic motions. Also, at this time any EWPT would be complete.Although initial calculations will use a fixed spacetime background, where the background metric is not evolved, we may find it necessary to dynamically model the turbulent/spacetime interactions in order to fully understand the physics of the early universe. To do this, we need good initial conditions for the curvature of space during this epoch.The Robertson-Walker (R-W) spacetime metric for a flat universe can be written as(8)
d
s
2
=
-
d
t
2
+
a
2
(
t
)
g
-
i
j
d
x
i
d
x
j
,
where t is the time-like coordinate, g
-
i
j is the maximally symmetric three-dimensional space metric, and a
(
t
) is the scale factor. For calculating the initial spacetime, conformal time, τ
=
-
(
1
/
a
(
t
)
H
(
t
)
) is often used instead of cosmic time, t
=
t
0
a
2, to simplify the equations. Therefore,
(9)
d
s
2
=
a
2
(
τ
)
g
μ
ν
d
x
μ
d
x
ν
.
Here, the scale factor and Hubble parameter can be calculated based on the temperature and mass-energy density of the early universe relative to today. In order to add gravitational waves, the metric g
μ
ν is broken into two components: g
~
μ
ν the background metric plus a perturbation h
μ
ν. The metric can be written as
(10)
g
μ
ν
=
g
~
μ
ν
+
h
μ
ν
.
The background metric may take any form. It is also assumed that any perturbations are linear. For this study, a R-W metric is presumed. Initial perturbations may or may not be used for numerical experiments within this study. This metric may involve scalar, (ϕ
,
ψ) vector (ω
i), and tensor (h
i
j) perturbations, (1). The tensor perturbations are symmetric and transverse-traceless so ∂
i
h
i
j
=
0 and δ
i
j
h
i
j
=
0. The initial amplitude and spectrum of any of these fluctuations depend on the theory used to explain the generation of perturbations from inflation. Note that by comparing (4) and (1), it can be shown that scalar (ϕ) and vector (ω
i) perturbations can be related to the chosen lapse and shift in the same way that scalar (ψ) and tensor (h
i
j) perturbations can relate to the three-metric and extrinsic curvature. Because the focus here is on tensor perturbations, a geodesic slicing is used so the lapse, N, is set to unity and the shift vector, N
i, is set to zero. Since there are no singularities in either the proposed study or the tests described in Section 3, geodesic slicing should be sufficient for this work. Later, if scalar perturbations are included, they can be added by modifying the lapse, as well as the three-metric and extrinsic curvature.The process of generating tensor perturbations or gravitational waves from quantum fluctuations during inflation is similar to the process of generating scalar or vector perturbations. For example, work by Grishchuk [30–32] gives a basis for calculating the spectrum and amplitude of these waves for different slow-roll parameters.The possibility of all polarizations is included usingh
i
j
+
,
h
i
j
×
,
h
i
j
L, and h
i
j
R. Here, h
i
j
+ and h
i
j
× are the plus and cross-polarizations, respectively, and h
i
j
L and h
i
j
R are the left and right rotating polarizations defined by
(11)
h
i
j
L
=
1
2
(
h
i
j
+
-
i
h
i
j
×
)
,
h
i
j
R
=
1
2
(
h
i
j
+
+
i
h
i
j
×
)
.
According to work by Alexander and Martin [33], rotating polarizations should dominate in the early universe and satisfy the equations:
(12)
□
h
i
j
L
=
-
2
i
θ
h
˙
i
j
′
L
,
□
h
i
j
R
=
+
2
i
θ
h
˙
i
j
′
R
.
Here, the dot denotes a time-like derivative with respect to τ and the prime denotes a spatial derivative along the gravitational wave’s direction of propagation. If Alexander’s θ value is set to zero, the unpolarized gravitational wave signature is recovered [34]. For this work, rotational polarizations may have the added benefit of introducing extra vorticity into the homogeneous plasma. This may result in an increased magnetic field due to the dynamo effect. During inflation, a nonzero θ term results in a decrease in the amplitude of h
L and an amplification of h
R and, therefore, cosmological birefringence.Initial gravitational wave spectra from sources such as bubble collisions or phase transitions can be added linearly. Therefore, the initial gravitational wave spectrum can correspond to those predicted by supersymmetry, loop quantum gravity, string theory, deformed special relativity, variable speed of light theories, or many other theoretical models in order to provide potential tests of theoretical physics once the system is evolved with the turbulent matter field.


#### 2.2.5. Numerical Evolutions

The development of a stable and accurate GRMHD code has been a work in progress for the past several years. We now have a new GRMHD code in the testing and improvement stage. This code is designed to perform spacetime evolutions using either 2nd order finite, 4th order finite, or Fourier pseudospectral differencing methods with magnetohydrodynamic and pressure-less matter fields as well as various boundary and gauge conditions. The code can evolve the matter field independently of the spacetime in cases where a fully dynamic spacetime is not needed so the spacetime metric does not have to be evolved.We utilize the Cactus framework to develop this code. We developed an arrangement for Cactus that contains the GRMHD initial data, analysis, and evolution thorns. This code handles the physics, while Cactus does the IO and parallelization. The code is structured so that all the differencing is done outside of the main loops. This allows us to choose between several differencing techniques such as finite differencing or Fourier spectral differencing at run-time. We also used the Cactus method of lines routines to supply the time integrators. The spacetime (LHS of Einstein’s equations) can be evolved using a strongly hyperbolic form of the BSSN equations as defined by Brown et al. [24]. The matter field (RHS of Einstein’s equations) is evolved by the form of the GRMHD equations as defined by Duez et al. [23] with divergence cleaning and artificial viscosity. Periodic boundary conditions are also used so that the simulation domain can accurately represent a homogenous slice of a much larger universe.In addition to developing an accurate GRMHD evolution code, it is important to extrapolate the data so that it can be compared to cosmological observations. Gravitational waves can be calculated directly from the stress-energy tensor using its quadrupole moments. By doing this, the spectrum and relative amplitude of primordial gravitational waves created as a result of the turbulence in the matter field can be determined. Eventually, these results can be compared to stochastic gravitational wave data from GW observatories and observations of the cosmic microwave background.There are many challenges to evolving the numerical code. First, there are the standard difficulties of dealing with a nonlinear code. Speed and accuracy are the most important issues. Also, there are additional challenges because the GRMHD code utilizes a nonlinear primitive variable solver to recover elements of the stress-energy tensor from the MHD evolution variables, shock capturing techniques, and a technique called divergence cleaning to maintain physical values for theB-field. Optimizing and improving these solvers are essential to developing a fast, accurate, and stable code.Before running the experiments we thoroughly tested the code. The Duez et al. paper [23] suggested four tests of a GRMHD code; however, because of the limited scope of this study I felt that only the following tests are necessary: gravitational wave-induced MHD waves, Minkowski spacetime MHD tests, such as shock tests and consistency with the standard model of cosmology. I will not include tests of unmagnetized relativistic stars or relativistic Bondi flow in this paper because the spacetime that they are simulating lacks stars and black holes.This code is being developed and run on a variety of computing resources including: UHCL’s Athena cluster, University of Houston’s Maxwell cluster, and University of Texas’ Ranger cluster via the XSEDE network.


#### 2.2.6. Data Analysis

The data analysis part of this project will focus around determining the gravitational wave spectrum from the simulation. A Fourier analysis of the quadrupole moments of the stress-energy tensor should yield the spectrum of gravitational waves produced by the turbulent matter field. Much of the data analysis work consists of the addition and fine tuning of new analysis routines in the code.Visualization of Cactus-generated data is done using a variety of open-source software such as VisIt, xgraph, ygraph, and gnuplot. Each requires implementation scripts to be written. These scripts will tell the visualization program how to read the Cactus-generated data files. Modifications of the data include fast Fourier transforms (FFT) for spectral analysis.As these numerical experiments are being performed, the output is analyzed. The effects of variations in the density, temperature, magnetic field, and initial turbulence will be studied in the output data. A Fourier analysis of the perturbed quadrupole moments will be performed in order to extract the spectrum of gravitational waves. This spectrum will then be extrapolated to give the current observable values. The result will be several templates of GW spectra resulting from different initial conditions.


## 2.1. Background and Hypotheses

Several different mechanisms for producing primordial gravitational waves have been identified and studied by researchers. These include quantum fluctuations during inflation, bubble wall motion, and collisions during phase transitions, cosmological magnetic fields, oscillating classical fields during reheating, cosmological defects, and plasma turbulence [4–14]. We assume that the early universe was in a metastable state during a first order phase transition. The false vacuum was separated from the true vacuum by a potential barrier or a scalar field. Quantum tunneling occurred across the barrier in finite regions of space resulting in true vacuum bubbles inside the false vacuum phase. As the universe expanded and cooled, the energy difference between the false vacuum and true vacuum got larger, making the phase transition more probable. Eventually, the probability of nucleating one critical bubble per Hubble time became high enough to cause the phase transition to begin. This defined the transition temperature, which is believed to be about 1 TeV [4]. The nucleated bubbles expanded and collided, eventually filling the whole universe. The collision of two or more of these bubbles broke spherical symmetry and released some of their energy as gravitational waves. Since the expansion of the bubbles was accompanied by macroscopic motions in the cosmic matter field, the collision of these bubbles also resulted in the anisotropic stirring of the field. This caused turbulent motions which provided a primary source of gravitational waves for this research.


### 2.1.1. Primordial Magnetic Fields

Magnetic fields are believed to have played a large part in the dynamics of the universe’s evolution. Little is known about the existence of magnetic fields in the early universe. There are no direct observations of primordial magnetic fields. Theories also disagree on the amplitude of primordial magnetic fields. There are currently several dozen theories about the origin of cosmic magnetic fields [15, 16]. The main reason that we believe that primordial magnetic fields existed is because they may have been needed to seed the large magnetic fields observed today. Most theories of cosmic magnetic field generation fall into one of three categories [15–17]: (1) magnetic fields generated by phase transitions; (2) electromagnetic perturbations expanded by inflation; and (3) turbulent magnetofluid resulting in charge and current asymmetries.Most models calculate the magnitude of primordial magnetic fields by starting with the observed strength of galactic or intergalactic magnetic fields and calculating how this field should have been amplified or diffused by external effects such as the galactic dynamo and expansion of the universe [15, 16]. A major problem is that there does not appear to be a universal agreement of how efficiently a galactic dynamo could have strengthened seed magnetic fields. Estimates of the strength of these seed fields can vary by tens of orders of magnitude. Seed magnetic fields produced during inflation are predicted to have a current strength somewhere between 1
0
-
11 G and 1
0
-
9 G on a scale of a few Mpc [15, 16, 18]. Magnetic seed fields generated by phase transitions are believed to be less than 1
0
-
23 G at galactic scales [15, 16]. Some turbulence theories imply that magnetic fields were not generated until after the first stars were formed therefore requiring no magnetic seed fields [15].Given how little is understood about primordial magnetic fields and the general lack of agreement among theoretical predictions, it seems clear that the existence of primordial magnetic fields can neither be confirmed or ruled out. It seems that the best we can do is set an upper limit on the strength of primordial magnetic fields and utilize this limit as a starting point in developing models of cosmic turbulence. Observations of the CMB limit the intensity of the magnetic seed fields to a current upper limit of1
0
-
9 G [15, 16, 18, 19].It is well known that gravitational waves can interact with a magnetofluid in the presence of a magnetic field. Work by Duez et al. [20] showed how gravitational waves can induce oscillatory modes in a plasma field if magnetic fields are present. Work by Kahniashvili and others [4, 9–13] has shown how a turbulent plasma can yield gravitational waves. The result may be a highly nonlinear interaction as energy is transferred from the fluid to the gravitational waves and back.


### 2.1.2. Turbulence in the Early Universe

Turbulence provides a particularly interesting GW source because it is not well understood analytically. This turbulence is a natural result of dynamics of the early universe resulting from bubble wall collisions and other chaotic events during the first order phase transitions. Analytic work done to date [4, 5, 12–14] summarizes the dynamics of the phase transitions using two quantities, α and β. α is traditionally defined as the ratio of false vacuum energy and plasma thermal energy density. This provides a measure of the transition strength. If α is much less than one, the transition is very weak. If α is larger than unity, the transition is very strongly first order. β is the rate of variation of the nucleation rate at the transition time. It fixes the time scale of the phase transition once the transition has begun. After a time interval β
-
1, the whole universe is converted to a true vacuum phase. Therefore the turbulent stirring should only last β
-
1.The amount of gravitational waves emitted by bubble collisions and turbulence generated in the plasma are also determined from two quantities,κ and v
b. κ is the fraction of vacuum energy transferred into fluid kinetic energy and v
b is the velocity of bubble wall expansion. Bubble walls can propagate via two modes, detonation and deflagration [4, 5, 12–14]. For detonation, the bubble walls are thin compared to the radius and they propagate faster than the speed of sound. This results in
(2)
v
b
(
α
)
=
1
/
3
+
(
α
2
+
2
α
/
3
)
1
/
2
1
+
α
,
κ
(
α
)
=
1
1
+
0.715
α
[
0.715
α
+
4
27
3
α
2
]
.
If the bubbles propagate by deflagration, the walls are thick and have a lower energy density. It is currently believed that for a relativistic plasma, the deflagration expansion mode is unstable so only the detonation modes will result.The number density of turbulent eddies within a Hubble radius should depend onv
b and β. The characteristic velocity perturbation of the turbulent fluid for the largest eddies at the stirring scale is given by
(3)
v
0
=
3
κ
α
4
+
3
κ
α
.
For κ
α
≈
1, which corresponds to a strongly first order phase transition, v
0 is about 0.65 at the time of the electroweak phase transition. We later use v
0 as the maximum velocity of fluid elements in our studies. This velocity is randomized in amplitude and direction in order to simulate the initial conditions for turbulence.


## 2.1.1. Primordial Magnetic Fields

Magnetic fields are believed to have played a large part in the dynamics of the universe’s evolution. Little is known about the existence of magnetic fields in the early universe. There are no direct observations of primordial magnetic fields. Theories also disagree on the amplitude of primordial magnetic fields. There are currently several dozen theories about the origin of cosmic magnetic fields [15, 16]. The main reason that we believe that primordial magnetic fields existed is because they may have been needed to seed the large magnetic fields observed today. Most theories of cosmic magnetic field generation fall into one of three categories [15–17]: (1) magnetic fields generated by phase transitions; (2) electromagnetic perturbations expanded by inflation; and (3) turbulent magnetofluid resulting in charge and current asymmetries.Most models calculate the magnitude of primordial magnetic fields by starting with the observed strength of galactic or intergalactic magnetic fields and calculating how this field should have been amplified or diffused by external effects such as the galactic dynamo and expansion of the universe [15, 16]. A major problem is that there does not appear to be a universal agreement of how efficiently a galactic dynamo could have strengthened seed magnetic fields. Estimates of the strength of these seed fields can vary by tens of orders of magnitude. Seed magnetic fields produced during inflation are predicted to have a current strength somewhere between 1
0
-
11 G and 1
0
-
9 G on a scale of a few Mpc [15, 16, 18]. Magnetic seed fields generated by phase transitions are believed to be less than 1
0
-
23 G at galactic scales [15, 16]. Some turbulence theories imply that magnetic fields were not generated until after the first stars were formed therefore requiring no magnetic seed fields [15].Given how little is understood about primordial magnetic fields and the general lack of agreement among theoretical predictions, it seems clear that the existence of primordial magnetic fields can neither be confirmed or ruled out. It seems that the best we can do is set an upper limit on the strength of primordial magnetic fields and utilize this limit as a starting point in developing models of cosmic turbulence. Observations of the CMB limit the intensity of the magnetic seed fields to a current upper limit of1
0
-
9 G [15, 16, 18, 19].It is well known that gravitational waves can interact with a magnetofluid in the presence of a magnetic field. Work by Duez et al. [20] showed how gravitational waves can induce oscillatory modes in a plasma field if magnetic fields are present. Work by Kahniashvili and others [4, 9–13] has shown how a turbulent plasma can yield gravitational waves. The result may be a highly nonlinear interaction as energy is transferred from the fluid to the gravitational waves and back.


## 2.1.2. Turbulence in the Early Universe

Turbulence provides a particularly interesting GW source because it is not well understood analytically. This turbulence is a natural result of dynamics of the early universe resulting from bubble wall collisions and other chaotic events during the first order phase transitions. Analytic work done to date [4, 5, 12–14] summarizes the dynamics of the phase transitions using two quantities, α and β. α is traditionally defined as the ratio of false vacuum energy and plasma thermal energy density. This provides a measure of the transition strength. If α is much less than one, the transition is very weak. If α is larger than unity, the transition is very strongly first order. β is the rate of variation of the nucleation rate at the transition time. It fixes the time scale of the phase transition once the transition has begun. After a time interval β
-
1, the whole universe is converted to a true vacuum phase. Therefore the turbulent stirring should only last β
-
1.The amount of gravitational waves emitted by bubble collisions and turbulence generated in the plasma are also determined from two quantities,κ and v
b. κ is the fraction of vacuum energy transferred into fluid kinetic energy and v
b is the velocity of bubble wall expansion. Bubble walls can propagate via two modes, detonation and deflagration [4, 5, 12–14]. For detonation, the bubble walls are thin compared to the radius and they propagate faster than the speed of sound. This results in
(2)
v
b
(
α
)
=
1
/
3
+
(
α
2
+
2
α
/
3
)
1
/
2
1
+
α
,
κ
(
α
)
=
1
1
+
0.715
α
[
0.715
α
+
4
27
3
α
2
]
.
If the bubbles propagate by deflagration, the walls are thick and have a lower energy density. It is currently believed that for a relativistic plasma, the deflagration expansion mode is unstable so only the detonation modes will result.The number density of turbulent eddies within a Hubble radius should depend onv
b and β. The characteristic velocity perturbation of the turbulent fluid for the largest eddies at the stirring scale is given by
(3)
v
0
=
3
κ
α
4
+
3
κ
α
.
For κ
α
≈
1, which corresponds to a strongly first order phase transition, v
0 is about 0.65 at the time of the electroweak phase transition. We later use v
0 as the maximum velocity of fluid elements in our studies. This velocity is randomized in amplitude and direction in order to simulate the initial conditions for turbulence.


## 2.2. Computational Model

In order to study the interaction of the plasma field and the background spacetime dynamically (or separately) our team has written and is testing/improving a GRMHD code implemented using the Cactus framework [21]. We now describe the basic variables and equations that constitute this model.


### 2.2.1. The Spacetime Evolution Model

The spacetime metric can be written as(4)
d
s
2
=
-
N
2
d
t
2
+
γ
i
j
(
x
→
,
t
)
(
d
x
i
+
N
i
d
t
)
(
d
x
j
+
N
j
d
t
)
.
Here N is the lapse, N
i is the shift vector, and γ
i
j is the spatial 3-metric [22, 23]. For this work, 3-metric and its “time-derivative,” the extrinsic curvature, “K
i
j” will be evolved using a strongly hyperbolic version of the BSSN formulation of numerical relativity [24].


### 2.2.2. The General Relativistic Magnetohydrodynamic Model

The fluid and electromagnetic fields of the GRMHD equations are developed from several well-known equations [25]. They include the conservation of particle number, the continuity equation, the conservation of energy-momentum, the magnetic constraint equation, and the magnetic induction equation. For a system consisting of a perfect fluid and an electromagnetic field, the ideal MHD stress-energy tensor is given by
(5)
T
μ
ν
=
(
ρ
0
h
+
b
2
)
u
μ
u
ν
+
(
P
+
b
2
2
)
g
μ
ν
-
b
μ
b
ν
,
h
=
1
+
ϵ
+
P
ρ
0
,
b
μ
=
1
4
π
B
(
u
)
μ
,
B
(
u
)
0
=
1
α
u
i
B
i
,
B
(
u
)
i
=
1
u
0
(
B
i
α
+
B
(
u
)
0
u
i
)
.
Here, P is the fluid pressure, ρ
0 is density, B
i is magnetic field, u
μ is four-velocity, h is the enthalpy, ϵ is specific internal energy, and b
2 is the magnitude of the magnetic vector field squared. The addition of viscosity modifies the MHD stress-energy tensor by incorporating the viscous stress tensor:
(6)
T
μ
ν
=
(
ρ
0
h
+
b
2
+
Q
)
u
μ
u
ν
+
(
P
+
b
2
2
)
g
μ
ν
-
b
μ
b
ν
+
Σ
μ
ν
.
Here Q is artificial bulk viscosity and Σ
μ
ν is the viscous stress tensor for artificial shear viscosity. Artificial viscosity is used here as a way of handling shocks although we are working on integrating more advanced HRSC techniques. Our viscosity terms [26] are described and defined below and are only meant to be used when the divergence of the fluid flow is negative:
(7)
Q
=
I
n
Δ
l
∂
k
V
k
(
k
q
Δ
l
∂
k
V
k
-
k
l
C
s
)
,
Σ
j
i
=
I
n
Δ
l
(
k
q
Δ
l
∂
k
V
k
-
k
l
C
s
)
Sym
(
δ
j
V
i
-
∂
k
V
k
3
δ
j
i
)
,
I
n
=
(
ρ
+
ρ
ϵ
+
(
P
+
Q
+
b
2
)
)
N
.
In these equations V
i is the fluid velocity, C
s is the local speed of sound, Δ
l is the minimum covariant zone length, and k
q and k
l are constants multiplying the quadratic and linear contributions, respectively. The Sym(⋯) function in the shear viscosity equation is a symmetry operation.Dark matter can be added to the system using a two-fluid approach where the stress energy tensor for dark matter is added directly to the stress energy tensor for the magnetofluid, therefore, completing the right-hand side of Einstein’s equation.


### 2.2.3. Initial Conditions: Plasma Field

In order to get a more accurate picture of the primordial gravitational wave spectrum, we directly simulate the turbulent primordial universe with the most realistic initial conditions possible. These should include not only a plasma field with a realistic EOS, they should also include elements such as magnetic fields, dark matter, and possibly even gravitational waves produced by other sources.The study will begin att
>
1
0
-
6 seconds after the Big Bang near the beginning of the Hadron epoch, when the primordial plasma field began to look like relativistic plasma and the strong force could be safely ignored. At this point the Debye length is about 1
0
-
16 m so even a small computational domain should demonstrate the dynamics of the plasma. The plasma at this time was composed mainly of electrons, positrons, neutrinos, and photons. Many of the initial conditions at this epoch are well known or can be fairly easily calculated using the available literature [27]. In addition, given that most cosmological models agree that over 80% of the matter in the universe is composed of dark matter, our cosmological simulations can also include some nonmagnetized “dark” fluid. We can take this into account by adding a second pressureless nonmagnetized fluid to our initial plasma field.The initial matter field will be taken to be homogenous and turbulent. We introduce turbulence into the system by randomly varying the initial velocities of the fluid elements up to the magnitude ofv
0, (3). In addition, we will be working to better establish the initial conditions resulting from the first order EWPT. Based on arguments in previous work [15–17, 28, 29], the initial magnetic field during this epoch should be less than or equal to 1
0
17 G.


### 2.2.4. Initial Conditions: Spacetime

For this computational study the initial spacetime is constructed in such a way as to mimic the conditions present during the Hadron epoch (t
>
1
0
-
6 seconds). I choose to begin the simulation during the Hadron epoch because at that time the primordial plasma field appeared to look like a relativistic plasma field that could be modeled using a GRMHD code as opposed to a quark gluon plasma field. At this time the strong force could be safely ignored as electromagnetic effects would dominate the plasma’s dynamic motions. Also, at this time any EWPT would be complete.Although initial calculations will use a fixed spacetime background, where the background metric is not evolved, we may find it necessary to dynamically model the turbulent/spacetime interactions in order to fully understand the physics of the early universe. To do this, we need good initial conditions for the curvature of space during this epoch.The Robertson-Walker (R-W) spacetime metric for a flat universe can be written as(8)
d
s
2
=
-
d
t
2
+
a
2
(
t
)
g
-
i
j
d
x
i
d
x
j
,
where t is the time-like coordinate, g
-
i
j is the maximally symmetric three-dimensional space metric, and a
(
t
) is the scale factor. For calculating the initial spacetime, conformal time, τ
=
-
(
1
/
a
(
t
)
H
(
t
)
) is often used instead of cosmic time, t
=
t
0
a
2, to simplify the equations. Therefore,
(9)
d
s
2
=
a
2
(
τ
)
g
μ
ν
d
x
μ
d
x
ν
.
Here, the scale factor and Hubble parameter can be calculated based on the temperature and mass-energy density of the early universe relative to today. In order to add gravitational waves, the metric g
μ
ν is broken into two components: g
~
μ
ν the background metric plus a perturbation h
μ
ν. The metric can be written as
(10)
g
μ
ν
=
g
~
μ
ν
+
h
μ
ν
.
The background metric may take any form. It is also assumed that any perturbations are linear. For this study, a R-W metric is presumed. Initial perturbations may or may not be used for numerical experiments within this study. This metric may involve scalar, (ϕ
,
ψ) vector (ω
i), and tensor (h
i
j) perturbations, (1). The tensor perturbations are symmetric and transverse-traceless so ∂
i
h
i
j
=
0 and δ
i
j
h
i
j
=
0. The initial amplitude and spectrum of any of these fluctuations depend on the theory used to explain the generation of perturbations from inflation. Note that by comparing (4) and (1), it can be shown that scalar (ϕ) and vector (ω
i) perturbations can be related to the chosen lapse and shift in the same way that scalar (ψ) and tensor (h
i
j) perturbations can relate to the three-metric and extrinsic curvature. Because the focus here is on tensor perturbations, a geodesic slicing is used so the lapse, N, is set to unity and the shift vector, N
i, is set to zero. Since there are no singularities in either the proposed study or the tests described in Section 3, geodesic slicing should be sufficient for this work. Later, if scalar perturbations are included, they can be added by modifying the lapse, as well as the three-metric and extrinsic curvature.The process of generating tensor perturbations or gravitational waves from quantum fluctuations during inflation is similar to the process of generating scalar or vector perturbations. For example, work by Grishchuk [30–32] gives a basis for calculating the spectrum and amplitude of these waves for different slow-roll parameters.The possibility of all polarizations is included usingh
i
j
+
,
h
i
j
×
,
h
i
j
L, and h
i
j
R. Here, h
i
j
+ and h
i
j
× are the plus and cross-polarizations, respectively, and h
i
j
L and h
i
j
R are the left and right rotating polarizations defined by
(11)
h
i
j
L
=
1
2
(
h
i
j
+
-
i
h
i
j
×
)
,
h
i
j
R
=
1
2
(
h
i
j
+
+
i
h
i
j
×
)
.
According to work by Alexander and Martin [33], rotating polarizations should dominate in the early universe and satisfy the equations:
(12)
□
h
i
j
L
=
-
2
i
θ
h
˙
i
j
′
L
,
□
h
i
j
R
=
+
2
i
θ
h
˙
i
j
′
R
.
Here, the dot denotes a time-like derivative with respect to τ and the prime denotes a spatial derivative along the gravitational wave’s direction of propagation. If Alexander’s θ value is set to zero, the unpolarized gravitational wave signature is recovered [34]. For this work, rotational polarizations may have the added benefit of introducing extra vorticity into the homogeneous plasma. This may result in an increased magnetic field due to the dynamo effect. During inflation, a nonzero θ term results in a decrease in the amplitude of h
L and an amplification of h
R and, therefore, cosmological birefringence.Initial gravitational wave spectra from sources such as bubble collisions or phase transitions can be added linearly. Therefore, the initial gravitational wave spectrum can correspond to those predicted by supersymmetry, loop quantum gravity, string theory, deformed special relativity, variable speed of light theories, or many other theoretical models in order to provide potential tests of theoretical physics once the system is evolved with the turbulent matter field.


### 2.2.5. Numerical Evolutions

The development of a stable and accurate GRMHD code has been a work in progress for the past several years. We now have a new GRMHD code in the testing and improvement stage. This code is designed to perform spacetime evolutions using either 2nd order finite, 4th order finite, or Fourier pseudospectral differencing methods with magnetohydrodynamic and pressure-less matter fields as well as various boundary and gauge conditions. The code can evolve the matter field independently of the spacetime in cases where a fully dynamic spacetime is not needed so the spacetime metric does not have to be evolved.We utilize the Cactus framework to develop this code. We developed an arrangement for Cactus that contains the GRMHD initial data, analysis, and evolution thorns. This code handles the physics, while Cactus does the IO and parallelization. The code is structured so that all the differencing is done outside of the main loops. This allows us to choose between several differencing techniques such as finite differencing or Fourier spectral differencing at run-time. We also used the Cactus method of lines routines to supply the time integrators. The spacetime (LHS of Einstein’s equations) can be evolved using a strongly hyperbolic form of the BSSN equations as defined by Brown et al. [24]. The matter field (RHS of Einstein’s equations) is evolved by the form of the GRMHD equations as defined by Duez et al. [23] with divergence cleaning and artificial viscosity. Periodic boundary conditions are also used so that the simulation domain can accurately represent a homogenous slice of a much larger universe.In addition to developing an accurate GRMHD evolution code, it is important to extrapolate the data so that it can be compared to cosmological observations. Gravitational waves can be calculated directly from the stress-energy tensor using its quadrupole moments. By doing this, the spectrum and relative amplitude of primordial gravitational waves created as a result of the turbulence in the matter field can be determined. Eventually, these results can be compared to stochastic gravitational wave data from GW observatories and observations of the cosmic microwave background.There are many challenges to evolving the numerical code. First, there are the standard difficulties of dealing with a nonlinear code. Speed and accuracy are the most important issues. Also, there are additional challenges because the GRMHD code utilizes a nonlinear primitive variable solver to recover elements of the stress-energy tensor from the MHD evolution variables, shock capturing techniques, and a technique called divergence cleaning to maintain physical values for theB-field. Optimizing and improving these solvers are essential to developing a fast, accurate, and stable code.Before running the experiments we thoroughly tested the code. The Duez et al. paper [23] suggested four tests of a GRMHD code; however, because of the limited scope of this study I felt that only the following tests are necessary: gravitational wave-induced MHD waves, Minkowski spacetime MHD tests, such as shock tests and consistency with the standard model of cosmology. I will not include tests of unmagnetized relativistic stars or relativistic Bondi flow in this paper because the spacetime that they are simulating lacks stars and black holes.This code is being developed and run on a variety of computing resources including: UHCL’s Athena cluster, University of Houston’s Maxwell cluster, and University of Texas’ Ranger cluster via the XSEDE network.


### 2.2.6. Data Analysis

The data analysis part of this project will focus around determining the gravitational wave spectrum from the simulation. A Fourier analysis of the quadrupole moments of the stress-energy tensor should yield the spectrum of gravitational waves produced by the turbulent matter field. Much of the data analysis work consists of the addition and fine tuning of new analysis routines in the code.Visualization of Cactus-generated data is done using a variety of open-source software such as VisIt, xgraph, ygraph, and gnuplot. Each requires implementation scripts to be written. These scripts will tell the visualization program how to read the Cactus-generated data files. Modifications of the data include fast Fourier transforms (FFT) for spectral analysis.As these numerical experiments are being performed, the output is analyzed. The effects of variations in the density, temperature, magnetic field, and initial turbulence will be studied in the output data. A Fourier analysis of the perturbed quadrupole moments will be performed in order to extract the spectrum of gravitational waves. This spectrum will then be extrapolated to give the current observable values. The result will be several templates of GW spectra resulting from different initial conditions.


## 2.2.1. The Spacetime Evolution Model

The spacetime metric can be written as(4)
d
s
2
=
-
N
2
d
t
2
+
γ
i
j
(
x
→
,
t
)
(
d
x
i
+
N
i
d
t
)
(
d
x
j
+
N
j
d
t
)
.
Here N is the lapse, N
i is the shift vector, and γ
i
j is the spatial 3-metric [22, 23]. For this work, 3-metric and its “time-derivative,” the extrinsic curvature, “K
i
j” will be evolved using a strongly hyperbolic version of the BSSN formulation of numerical relativity [24].


## 2.2.2. The General Relativistic Magnetohydrodynamic Model

The fluid and electromagnetic fields of the GRMHD equations are developed from several well-known equations [25]. They include the conservation of particle number, the continuity equation, the conservation of energy-momentum, the magnetic constraint equation, and the magnetic induction equation. For a system consisting of a perfect fluid and an electromagnetic field, the ideal MHD stress-energy tensor is given by
(5)
T
μ
ν
=
(
ρ
0
h
+
b
2
)
u
μ
u
ν
+
(
P
+
b
2
2
)
g
μ
ν
-
b
μ
b
ν
,
h
=
1
+
ϵ
+
P
ρ
0
,
b
μ
=
1
4
π
B
(
u
)
μ
,
B
(
u
)
0
=
1
α
u
i
B
i
,
B
(
u
)
i
=
1
u
0
(
B
i
α
+
B
(
u
)
0
u
i
)
.
Here, P is the fluid pressure, ρ
0 is density, B
i is magnetic field, u
μ is four-velocity, h is the enthalpy, ϵ is specific internal energy, and b
2 is the magnitude of the magnetic vector field squared. The addition of viscosity modifies the MHD stress-energy tensor by incorporating the viscous stress tensor:
(6)
T
μ
ν
=
(
ρ
0
h
+
b
2
+
Q
)
u
μ
u
ν
+
(
P
+
b
2
2
)
g
μ
ν
-
b
μ
b
ν
+
Σ
μ
ν
.
Here Q is artificial bulk viscosity and Σ
μ
ν is the viscous stress tensor for artificial shear viscosity. Artificial viscosity is used here as a way of handling shocks although we are working on integrating more advanced HRSC techniques. Our viscosity terms [26] are described and defined below and are only meant to be used when the divergence of the fluid flow is negative:
(7)
Q
=
I
n
Δ
l
∂
k
V
k
(
k
q
Δ
l
∂
k
V
k
-
k
l
C
s
)
,
Σ
j
i
=
I
n
Δ
l
(
k
q
Δ
l
∂
k
V
k
-
k
l
C
s
)
Sym
(
δ
j
V
i
-
∂
k
V
k
3
δ
j
i
)
,
I
n
=
(
ρ
+
ρ
ϵ
+
(
P
+
Q
+
b
2
)
)
N
.
In these equations V
i is the fluid velocity, C
s is the local speed of sound, Δ
l is the minimum covariant zone length, and k
q and k
l are constants multiplying the quadratic and linear contributions, respectively. The Sym(⋯) function in the shear viscosity equation is a symmetry operation.Dark matter can be added to the system using a two-fluid approach where the stress energy tensor for dark matter is added directly to the stress energy tensor for the magnetofluid, therefore, completing the right-hand side of Einstein’s equation.


## 2.2.3. Initial Conditions: Plasma Field

In order to get a more accurate picture of the primordial gravitational wave spectrum, we directly simulate the turbulent primordial universe with the most realistic initial conditions possible. These should include not only a plasma field with a realistic EOS, they should also include elements such as magnetic fields, dark matter, and possibly even gravitational waves produced by other sources.The study will begin att
>
1
0
-
6 seconds after the Big Bang near the beginning of the Hadron epoch, when the primordial plasma field began to look like relativistic plasma and the strong force could be safely ignored. At this point the Debye length is about 1
0
-
16 m so even a small computational domain should demonstrate the dynamics of the plasma. The plasma at this time was composed mainly of electrons, positrons, neutrinos, and photons. Many of the initial conditions at this epoch are well known or can be fairly easily calculated using the available literature [27]. In addition, given that most cosmological models agree that over 80% of the matter in the universe is composed of dark matter, our cosmological simulations can also include some nonmagnetized “dark” fluid. We can take this into account by adding a second pressureless nonmagnetized fluid to our initial plasma field.The initial matter field will be taken to be homogenous and turbulent. We introduce turbulence into the system by randomly varying the initial velocities of the fluid elements up to the magnitude ofv
0, (3). In addition, we will be working to better establish the initial conditions resulting from the first order EWPT. Based on arguments in previous work [15–17, 28, 29], the initial magnetic field during this epoch should be less than or equal to 1
0
17 G.


## 2.2.4. Initial Conditions: Spacetime

For this computational study the initial spacetime is constructed in such a way as to mimic the conditions present during the Hadron epoch (t
>
1
0
-
6 seconds). I choose to begin the simulation during the Hadron epoch because at that time the primordial plasma field appeared to look like a relativistic plasma field that could be modeled using a GRMHD code as opposed to a quark gluon plasma field. At this time the strong force could be safely ignored as electromagnetic effects would dominate the plasma’s dynamic motions. Also, at this time any EWPT would be complete.Although initial calculations will use a fixed spacetime background, where the background metric is not evolved, we may find it necessary to dynamically model the turbulent/spacetime interactions in order to fully understand the physics of the early universe. To do this, we need good initial conditions for the curvature of space during this epoch.The Robertson-Walker (R-W) spacetime metric for a flat universe can be written as(8)
d
s
2
=
-
d
t
2
+
a
2
(
t
)
g
-
i
j
d
x
i
d
x
j
,
where t is the time-like coordinate, g
-
i
j is the maximally symmetric three-dimensional space metric, and a
(
t
) is the scale factor. For calculating the initial spacetime, conformal time, τ
=
-
(
1
/
a
(
t
)
H
(
t
)
) is often used instead of cosmic time, t
=
t
0
a
2, to simplify the equations. Therefore,
(9)
d
s
2
=
a
2
(
τ
)
g
μ
ν
d
x
μ
d
x
ν
.
Here, the scale factor and Hubble parameter can be calculated based on the temperature and mass-energy density of the early universe relative to today. In order to add gravitational waves, the metric g
μ
ν is broken into two components: g
~
μ
ν the background metric plus a perturbation h
μ
ν. The metric can be written as
(10)
g
μ
ν
=
g
~
μ
ν
+
h
μ
ν
.
The background metric may take any form. It is also assumed that any perturbations are linear. For this study, a R-W metric is presumed. Initial perturbations may or may not be used for numerical experiments within this study. This metric may involve scalar, (ϕ
,
ψ) vector (ω
i), and tensor (h
i
j) perturbations, (1). The tensor perturbations are symmetric and transverse-traceless so ∂
i
h
i
j
=
0 and δ
i
j
h
i
j
=
0. The initial amplitude and spectrum of any of these fluctuations depend on the theory used to explain the generation of perturbations from inflation. Note that by comparing (4) and (1), it can be shown that scalar (ϕ) and vector (ω
i) perturbations can be related to the chosen lapse and shift in the same way that scalar (ψ) and tensor (h
i
j) perturbations can relate to the three-metric and extrinsic curvature. Because the focus here is on tensor perturbations, a geodesic slicing is used so the lapse, N, is set to unity and the shift vector, N
i, is set to zero. Since there are no singularities in either the proposed study or the tests described in Section 3, geodesic slicing should be sufficient for this work. Later, if scalar perturbations are included, they can be added by modifying the lapse, as well as the three-metric and extrinsic curvature.The process of generating tensor perturbations or gravitational waves from quantum fluctuations during inflation is similar to the process of generating scalar or vector perturbations. For example, work by Grishchuk [30–32] gives a basis for calculating the spectrum and amplitude of these waves for different slow-roll parameters.The possibility of all polarizations is included usingh
i
j
+
,
h
i
j
×
,
h
i
j
L, and h
i
j
R. Here, h
i
j
+ and h
i
j
× are the plus and cross-polarizations, respectively, and h
i
j
L and h
i
j
R are the left and right rotating polarizations defined by
(11)
h
i
j
L
=
1
2
(
h
i
j
+
-
i
h
i
j
×
)
,
h
i
j
R
=
1
2
(
h
i
j
+
+
i
h
i
j
×
)
.
According to work by Alexander and Martin [33], rotating polarizations should dominate in the early universe and satisfy the equations:
(12)
□
h
i
j
L
=
-
2
i
θ
h
˙
i
j
′
L
,
□
h
i
j
R
=
+
2
i
θ
h
˙
i
j
′
R
.
Here, the dot denotes a time-like derivative with respect to τ and the prime denotes a spatial derivative along the gravitational wave’s direction of propagation. If Alexander’s θ value is set to zero, the unpolarized gravitational wave signature is recovered [34]. For this work, rotational polarizations may have the added benefit of introducing extra vorticity into the homogeneous plasma. This may result in an increased magnetic field due to the dynamo effect. During inflation, a nonzero θ term results in a decrease in the amplitude of h
L and an amplification of h
R and, therefore, cosmological birefringence.Initial gravitational wave spectra from sources such as bubble collisions or phase transitions can be added linearly. Therefore, the initial gravitational wave spectrum can correspond to those predicted by supersymmetry, loop quantum gravity, string theory, deformed special relativity, variable speed of light theories, or many other theoretical models in order to provide potential tests of theoretical physics once the system is evolved with the turbulent matter field.


## 2.2.5. Numerical Evolutions

The development of a stable and accurate GRMHD code has been a work in progress for the past several years. We now have a new GRMHD code in the testing and improvement stage. This code is designed to perform spacetime evolutions using either 2nd order finite, 4th order finite, or Fourier pseudospectral differencing methods with magnetohydrodynamic and pressure-less matter fields as well as various boundary and gauge conditions. The code can evolve the matter field independently of the spacetime in cases where a fully dynamic spacetime is not needed so the spacetime metric does not have to be evolved.We utilize the Cactus framework to develop this code. We developed an arrangement for Cactus that contains the GRMHD initial data, analysis, and evolution thorns. This code handles the physics, while Cactus does the IO and parallelization. The code is structured so that all the differencing is done outside of the main loops. This allows us to choose between several differencing techniques such as finite differencing or Fourier spectral differencing at run-time. We also used the Cactus method of lines routines to supply the time integrators. The spacetime (LHS of Einstein’s equations) can be evolved using a strongly hyperbolic form of the BSSN equations as defined by Brown et al. [24]. The matter field (RHS of Einstein’s equations) is evolved by the form of the GRMHD equations as defined by Duez et al. [23] with divergence cleaning and artificial viscosity. Periodic boundary conditions are also used so that the simulation domain can accurately represent a homogenous slice of a much larger universe.In addition to developing an accurate GRMHD evolution code, it is important to extrapolate the data so that it can be compared to cosmological observations. Gravitational waves can be calculated directly from the stress-energy tensor using its quadrupole moments. By doing this, the spectrum and relative amplitude of primordial gravitational waves created as a result of the turbulence in the matter field can be determined. Eventually, these results can be compared to stochastic gravitational wave data from GW observatories and observations of the cosmic microwave background.There are many challenges to evolving the numerical code. First, there are the standard difficulties of dealing with a nonlinear code. Speed and accuracy are the most important issues. Also, there are additional challenges because the GRMHD code utilizes a nonlinear primitive variable solver to recover elements of the stress-energy tensor from the MHD evolution variables, shock capturing techniques, and a technique called divergence cleaning to maintain physical values for theB-field. Optimizing and improving these solvers are essential to developing a fast, accurate, and stable code.Before running the experiments we thoroughly tested the code. The Duez et al. paper [23] suggested four tests of a GRMHD code; however, because of the limited scope of this study I felt that only the following tests are necessary: gravitational wave-induced MHD waves, Minkowski spacetime MHD tests, such as shock tests and consistency with the standard model of cosmology. I will not include tests of unmagnetized relativistic stars or relativistic Bondi flow in this paper because the spacetime that they are simulating lacks stars and black holes.This code is being developed and run on a variety of computing resources including: UHCL’s Athena cluster, University of Houston’s Maxwell cluster, and University of Texas’ Ranger cluster via the XSEDE network.


## 2.2.6. Data Analysis

The data analysis part of this project will focus around determining the gravitational wave spectrum from the simulation. A Fourier analysis of the quadrupole moments of the stress-energy tensor should yield the spectrum of gravitational waves produced by the turbulent matter field. Much of the data analysis work consists of the addition and fine tuning of new analysis routines in the code.Visualization of Cactus-generated data is done using a variety of open-source software such as VisIt, xgraph, ygraph, and gnuplot. Each requires implementation scripts to be written. These scripts will tell the visualization program how to read the Cactus-generated data files. Modifications of the data include fast Fourier transforms (FFT) for spectral analysis.As these numerical experiments are being performed, the output is analyzed. The effects of variations in the density, temperature, magnetic field, and initial turbulence will be studied in the output data. A Fourier analysis of the perturbed quadrupole moments will be performed in order to extract the spectrum of gravitational waves. This spectrum will then be extrapolated to give the current observable values. The result will be several templates of GW spectra resulting from different initial conditions.


## 3. Testing the Code

The first test that we performed involved generating Alfvén and magnetosonic modes by gravitational waves and comparing the results against the semianalytic solutions from the Duez et al. paper [20]. This semianalytic solution is only valid for a time much less than the dynamical collapse time of the unperturbed fluid. We began by using the same initial conditions as defined by Duez et al. ’s general example [23]:
(13)
h
+
(
t
,
z
)
=
h
+
0
sin
(
k
z
)
cos
⁡
(
k
t
)
,
h
×
(
t
,
z
)
=
h
×
0
sin
(
k
z
)
cos
⁡
(
k
t
)
,
P
(
0
,
z
)
=
1.29
×
1
0
-
9
,
ρ
0
(
0
,
z
)
=
2.78
×
1
0
-
9
,
v
i
(
0
,
z
)
=
0
,
B
i
(
0
,
z
)
=
(
1.09,8.26,14.4
)
×
1
0
-
5
,
h
+
0
=
h
×
0
=
1.18
×
1
0
-
4
,
where k is the wave number and we assume the fluid is unperturbed at t
=
0. Results are shown in Figure 1. The test was performed using second order, fourth order, and Fourier spectral differencing as a one-dimensional problem. The results shown used 200 grid points for the second order differencing, 50 grid points for the fourth order differencing, and 32 grid points for the spectral differencing. Additional tests were performed where the initial conditions were varied. For every variation the test proved successful, the analytic and numerical results proved almost identical to within a few percent. Runs were also conducted with 100, 25, and 16 grid points for the second order, fourth order, and spectral differencing, respectively, in order to test for convergence. As shown in Figure 1, the main source of errors in this test was phase errors between the analytic and numerical solutions. This made calculating convergence difficult. We were able to calculate convergence of the results at each time by dividing the L2 norm of the errors of the low resolution runs by that of the high resolution runs. The result was an oscillating pattern with a mean, after 5 crossing times, of around 2.4 for the 2nd order finite differencing, 4.67 for the 4th order finite differencing, and 2.7 for the Fourier Spectral differencing. A major factor affecting the convergence rate may be the larger time steps taken in the low resolution runs. Although we cannot show that the overall convergence rate matches that of the differencing method, we can show that the code does converge for each differencing method.Gravitational wave-induced MHD waves test results. The left column used 2nd order finite differencing, the middle column used 4th order finite differencing, and the right column used Fourier spectral differencing. Density is displayed here by calculatingδ
ρ
/
(
ρ
h
0
), velocity is displayed here by calculating δ
v
x
/
h
0, and B-field is displayed here by calculating δ
B
x
/
(
B
0
h
0
) as in the Duez et al. paper [23]. Time is shown on the x-axis. Because the physical size of the grid is 2 units, a time of 30 units represents 15 crossing times.
(a)
(b)
(c)For our second test we utilized a FRW spacetime that contained parameters (temperature, energy density, Hubble parameter, and scale factor) based with those accepted for the universe att
=
1
0
-
6 s. The goal of the test was to determine if it evolved consistently with the Friedmann equations. We set the initial scale factor a
=
6.6
×
1
0
-
14, the initial Hubble parameter H
=
1.4
×
1
0
26 km/s/Mpc, the initial temperature T
=
4.1
×
1
0
13 K, and the initial energy density ρ
0
ϵ
=
1.1
×
1
0
23 J/m
3. The system evolved until around t
=
86 s and the final value of these parameters was found at the end of the simulation. The code’s calculated change in each parameter was then compared to the values predicted from the Friedmann equations in order to determine how well the code’s results matched the analytic solution. The final time, 86 s, corresponded to physical age of the simulated universe after seven days of computing time and was not chosen to have any particular relationship to the physical size of the domain. There were no structures or inhomogeneities in the system so spacial resolution was not important in this test. This test proved successful with the errors of less than one percent as shown in Table 1.Table 1
Consistency of standard model of cosmology parameters test results.
t
T
ρ
ϵ
H
a
Initial
1
0
-
6
4.120
×
1
0
13
1.091
×
1
0
23
2.401
×
1
0
26
1.0
Final
8.629
×
1
0
1
4.446
×
1
0
9
1.480
×
1
0
7
2.792
×
1
0
18
9.289
×
1
0
3
Final/initial
8.629
×
1
0
7
1.079
×
1
0
-
4
1.357
×
1
0
-
16
1.163
×
1
0
-
8
9.289
×
1
0
3
t
a
8.629
×
1
0
7
8.586
×
1
0
7
8.586
×
1
0
7
8.601
×
1
0
7
8.629
×
1
0
7
T
a
1.077
×
1
0
-
4
1.079
×
1
0
-
4
1.079
×
1
0
-
4
1.078
×
1
0
-
4
1.077
×
1
0
-
4
ρ
ε
a
1.343
×
1
0
-
16
1.357
×
1
0
-
16
1.357
×
1
0
-
16
1.352
×
1
0
-
16
1.343
×
1
0
-
16
H
a
1.159
×
1
0
-
8
1.165
×
1
0
-
8
1.165
×
1
0
-
8
1.163
×
1
0
-
8
1.159
×
1
0
-
8
a
a
9.289
×
1
0
3
9.266
×
1
0
3
9.266
×
1
0
3
9.274
×
1
0
3
9.289
×
1
0
3
t
b
0.0
4.338
×
1
0
5
4.338
×
1
0
7
2.824
×
1
0
5
-
2.980
×
1
0
-
7
T
b
2.716
×
1
0
-
7
0.0
-
5.833
×
1
0
-
16
9.506
×
1
0
-
8
2.716
×
1
0
-
7
ρ
ε
b
1.361
×
1
0
-
18
2.933
×
1
0
-
27
0.0
4.773
×
1
0
-
19
1.361
×
1
0
-
18
H
b
3.804
×
1
0
-
11
-
2.051
×
1
0
-
11
-
2.051
×
1
0
-
11
0.0
3.804
×
1
0
-
11
a
b
1.637
×
1
0
-
11
2.338
×
1
0
1
2.338
×
1
0
1
1.521
×
1
0
1
0.0
t
c
3.351
×
1
0
-
3
T
c
1.678
×
1
0
-
3
ρ
ε
c
6.686
×
1
0
-
3
H
c
3.357
×
1
0
-
3
a
c
1.678
×
1
0
-
3
aFinal/initial for each column parameter calculated based on the Friedman equations.
bRatio predicted by Friedman equations (analytic) − ratio calculated by the GRMHD code (numerical).
cAverage error: (numerical − analytic)/numerical averaged among all nonzero columns.At this point we are prepared to add standing gravitational waves with a spectrum consistent to Grishchuk’s predictions [30–32]. The gravitational waves were given random phases in order to avoid large nodes and antinodes. We then ran simulations with spacetime perturbations and large (1
0
17 Gauss) magnetic fields. These simulations showed that spacetime perturbations and magnetic fields had no significant impact on the expansion rate of the spacetime and therefore adding them should still result in a simulation consistent with cosmological theory to within the same error as found in Table 1.Finally, we conducted shock tests using similar initial conditions to Duez et al. [23] and Komissarov [35] as shown in Table 2. The results shown are based on runs utilizing the second order finite differencing method. The fourth order finite differencing method gives similar results, as expected, although the fourth order tests were conducted using a lower resolution grid. These tests could not be completed with our Fourier spectral differencing method because the suggested shock tests are not periodic in nature. Future work may involve testing all three differencing methods using a periodic shock test, but the current results suggest that this may not be a worthwhile effort until HRSC techniques can be incorporated into the code.Table 2
Initial states for one-dimensional MHD testsa.
Test
Left state
Right state
Grid
t
final
Fast shock (μ
=
0
.
2
b)
u
i
=
(
25.0,0.0,0.0
)
u
i
=
(
1.091,0.3923,0.00
)
n
=
40
2.5
B
i
/
4
π
=
(
20.0,25.02,0.0
)
B
i
/
4
π
=
(
20.0,49.0,0.0
)
P
=
1.0, ρ
0
=
1.0
P
=
367.5, ρ
0
=
25.48
Slow shock (μ
=
0
.
5
b)
u
i
=
(
1.53,0.0,0.0
)
u
i
=
(
0.9571
,
-
0.6822,0.00
)
n
=
200
2.0
B
i
/
4
π
=
(
10.0,18.28,0.0
)
B
i
/
4
π
=
(
10.0,14.49,0.0
)
P
=
10.0, ρ
0
=
1.0
P
=
55.36, ρ
0
=
3.323
Switch-off fast rarefaction
u
i
=
(
-
2.0,0.0,0.0
)
u
i
=
(
-
0.212
,
-
0.590,0.0
)
n
=
150
1.0
B
i
/
4
π
=
(
2.0,0.0,0.0
)
B
i
/
4
π
=
(
2.0,4.71,0.0
)
P
=
1.0, ρ
0
=
0.1
P
=
10.0, ρ
0
=
0.562
Switch-on slow rarefaction
u
i
=
(
-
0.765
,
-
1.386,0.0
)
u
i
=
(
0.0,0.0,0.0
)
n
=
150
2.0
B
i
/
4
π
=
(
1.0,1.022,0.0
)
B
i
/
4
π
=
(
1.0,0.0,0.0
)
P
=
0.1, ρ
0
=
1.78
×
1
0
-
3
P
=
1.0, ρ
0
=
0.01
Alfvén wavec (μ
=
0.62
6
b)
u
i
=
(
0.0,0.0,0.0
)
u
i
=
(
3.70,5.76,0.00
)
n
=
200
2.0
B
i
/
4
π
=
(
3.0,3.0,0.0
)
B
i
/
4
π
=
(
3.0
,
-
6.857,0.0
)
P
=
1.0, ρ
0
=
1.0
P
=
1.0, ρ
0
=
1.0
Shock tube 1
u
i
=
(
0.0,0.0,0.0
)
u
i
=
(
0.0,0.0,0.0
)
n
=
200
1.0
B
i
/
4
π
=
(
1.0,0.0,0.0
)
B
i
/
4
π
=
(
1.0,0.0,0.0
)
P
=
1000.0, ρ
0
=
1.0
P
=
1.0, ρ
0
=
0.1
Shock tube 2
u
i
=
(
0.0,0.0,0.0
)
u
i
=
(
0.0,0.0,0.0
)
n
=
200
1.0
B
i
/
4
π
=
(
0.0,20.0,0.0
)
B
i
/
4
π
=
(
0.0,0.0,0.0
)
P
=
30.0, ρ
0
=
1.0
P
=
1.0, ρ
0
=
0.1
Collision
u
i
=
(
5.0,0.0,0.0
)
u
i
=
(
-
5.0,0.0,0.0
)
n
=
200
1.22
B
i
/
4
π
=
(
10.0,10.0,0.0
)
B
i
/
4
π
=
(
10.0
,
-
10.0,0.0
)
P
=
1.0, ρ
0
=
1.0
P
=
1.0, ρ
0
=
1.0
aIn all cases, the gas satisfies the Γ-law EOS with Γ
=
4
/
3. For the first 7 tests, the left state refers to x
<
0 and the right state x
>
0.
b
μ is the speed at which the wave travels.
cFor the nonlinear Alfvén wave, the left and right states are joined by a continuous function separated by 0.5 units.Overall, the artificial viscosity methods used in this code seemed to produce less accurate results than the high resolution shock capturing methods used by Duez et al. [23] and Komissarov [35].Fast and Slow Shocks: the shock fronts for both of these cases were more distorted than the results of Duez et al. [23] and Komissarov [35]. Also, the fast shock appeared to move slower than the slow shock which does not seem to agree with the standard results.Switch-on/off Rarefaction: while the switch-off (fast rarefaction) seemed to agree with the published results the switch-on (slow rarefaction) appeared distorted at the shock front.Alfvén Wave: our Alfvén wave results seemed to agree with the published results except for the extra dip for z
<
0.Shock Tubes 1 and 2: shock tube tests also produced distorted shock fronts, particularly for z
>
0.Collision: the collision test produced the best results when compared to the established results.


## 4. Preliminary Results

In order to present a relevant proof of concept on the use of GRMHD in cosmology, we evolve a turbulent plasma field, with conditions similar to the universe when it was1
0
-
6 s old. We use similar conditions as outlined in Section 2.2.3 and the consistency of standard model of cosmology parameters test with an initial uniform magnetic field of 1
0
15 G. No initial gravitational waves or dark matter was included in the system. We also introduced turbulence with a random initial velocity of 0.65.The data presented in Figure3 correspond to evolving the initial conditions stated above to around 3.5
×
1
0
-
5 s. Before calculating the PSDs of each of the quantities, they were normalized by dividing the perturbation amplitude by the mean value of the quantity. The results where plotted logarithmically for a frequency range from 0.0015 Hz to 0.12 Hz. Although the simulation was run to 3.5
×
1
0
-
5 s, the normalized PSDs did not seem to vary much after about a hundred iterations (see Figure 2). The normalized PSDs of the relevant parameters were plotted on a logarithmic scale for a frequency range relevant to a potential space based gravitational wave interferometer such as eLISA. The results show that for strong uniform initial magnetic fields, noticeable perturbations are generated in the spacetime metric, density, temperature, and magnetic field terms which are different from the perturbations in the velocity field. Perturbations in the metric are the most interesting of these results because they correspond to gravitational waves. These were vanishingly small for magnetic fields less than or equal to 1
0
12 G. Much more work is needed to more fully understand the dynamics of the interaction between GRMHD turbulence and gravitational waves, but the results so far clearly show that gravitational wave generation from primordial turbulence is possible.Shock test results shown using second order finite differencing and artificial viscosity att
final time. The pressure profiles are shown on the left and the z component of velocity profiles are shown on the right.
(a)
(b)Figure 3
Logarithms of the power spectral densities of several quantities after a turbulent relativistic plasma were evolved using the GRMHD code from1
0
-
6 s to around 3.5 ×
1
0
-
5 s. The perturbations were normalized using the mean amplitude of each quantity before the PSD was calculated.


## 5. Discussion

We now have several parameters to work with in developing numerical experiments. Assuming the expansion properties of the background spacetime and composition of the matter field are fixed, we can alter the metric perturbations (scalar and tensor), initial magnetic field strength, and the turbulence of the initial matter field. By altering these properties, we can perform a variety of numerical experiments to determine the effects of scalar perturbations and gravitational waves on structure formation, limits on primordial magnetic fields, properties of gravitational waves formed by turbulent plasma, the dynamics of a turbulent plasma in an expanding universe, and other interesting scientific properties.One of the most interesting of these experiments involves the interaction between gravitational waves and the primordial plasma. Work by Duez et al. [20] showed that gravitational waves can induce oscillatory modes in a plasma. According to Shebalin [36–38], large-scale coherent structures grow naturally out of MHD turbulence. Here, structure is defined as strengthening magnetic fields, permanent density and temperature variations, and secondary relic gravitational waves. One can assume that spacetime perturbations in the early universe (sometime after t
=
1
0
-
6 seconds) interacted with the primordial plasma and resulted in Alfvén and magnetosonic modes [20]. These modes then interacted dynamically, possibly resulting in turbulence and structure formation [36–38]. Using the techniques of numerical relativity, we can test this assumption.



---

*Source: 407197-2014-04-22.xml*