{'fraction_non_alphanumeric': 0.0301936382779968, 'fraction_numerical': 0.010683354707959277, 'mean_word_length': 4.463248491663603, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Holography in Quark-Gluon Plasma and Neutron Stars

In this thesis, QCD is studied from three different directions, with one overarching theme: holography. The holographic duality allows certain strongly coupled QFTs to be described in terms of much simpler classical gravity in one dimension more. The first direction from which QCD is studied in this thesis is by examining the effects of an external magnetic field on a particular holographic model of QCD, yielding interesting qualitative insight. The second approach examines how, in the same model, one can describe dense baryonic configurations, providing a new way to study the matter composing neutron stars. Indeed, the equation of state produced in this way is subsequently used to compute several neutron star properties which are observable, or will be in the near future. The last direction contains no holographic computations per se, but does incorporate several qualitative insights from holography into a new heavy ion code called Trajectum. This will in the near future be used to perform a Bayesian analysis, whereby it is hoped that these qualitative insights from holography can be tested on experimental data, to see how well the ideas coming from holography match up with experiment.


QCD
QCD is the non-Abelian gauge theory of SU (3), which is minimally coupled to a number of quark flavors N f . In the standard model, there are of course 6 flavors. However, the three heaviest flavors are too heavy to be of importance for many observables, and can be safely neglected. 1 In this thesis, we need a slight generalization of QCD, namely to that of a gauge group SU (N c ), where now N c and N f can be freely chosen as theory parameters. The Lagrangian for this generalized QCD is with A µ ∈ su(N c ), and m i the mass of quark flavor i. This theory has the property that the beta function of the coupling α s = g 2 /4π is negative to first order in perturbation theory provided that N f /N c < 11/2 [1,2]. A consequence of this negative beta function is that the coupling constant decreases towards higher energies, a phenomenon known as asymptotic freedom, and increases towards lower energies. In figure 1.1, one can see that this is indeed also seen in experiments. It can be seen that around 1 GeV, the coupling constant becomes O(1), and the theory can no longer be accurately described by perturbation theory. This is a huge obstacle in the way of understanding QCD at low energy scales. One method by which one can still compute certain observables in the non-perturbative regime is lattice QCD. This method discretizes QCD on a Euclidean lattice, enabling the computation of non-dynamical observables by Monte Carlo integration of the euclidean path integral. Lattice QCD is a reliable method to compute not only thermal properties of QCD, but also hadron masses, which have been favorably compared to experimental values. It is not without its downsides though, as the Euclidean formalism makes the computation of dynamical processes extremely challenging. Also, for similar reasons, it turns out to be rather difficult to introduce a finite baryon chemical potential, an issue known as the sign problem [4]. A detailed discussion of lattice QCD is beyond the scope of this thesis, but excellent introductions can be found in [5,6]. Two features of QCD which appear in the low energy regime are confinement and chiral symmetry breaking. Confinement is the phenomenon that states have to be color neutral, implying that particles carrying color charge, such as quarks and gluons, can not occur as isolated particles. One example of this which can be computed in lattice QCD is the quark-antiquark potential, which is shown in the left panel of figure 1.2 [7]. This quark-antiquark potential quantifies the potential energy between a heavy quark paired with an antiquark of the same flavor, and can be seen to grow linearly at large separation. Assuming that all quarks in the theory are infinitely massive, this  [7]. Right: renormalized chiral condensate ψ ψ (labeled qq in the main text), renormalized such that ψ ψ (T = 0) = 0. Figure taken from [8].
means that one would have to spend an infinite amount of energy to separate the quark-antiquark pair. In the case of realistic quark masses, instead this implies that once the quarks are separated far enough, the potential energy stored in the gluon field will be large enough such that a new quark-antiquark pair can be created. The new quarks then each pair with one of the original quarks to create two color neutral mesons. Chiral symmetry breaking refers to the approximate global chiral symmetry which acts on (1.1) such that the SU (N f ) L × U (1) L generators act on the left-handed quark components by multiplication, while the SU (N f ) R ×U (1) R generators act on the right-handed components. This symmetry is only approximate in (1.1), but becomes exact in the limit where the quarks are massless. The QCD vacuum, however, breaks this approximate symmetry further spontaneously. The order parameter of this chiral symmetry breaking is the chiral condensate operator q i q i , which can be defined for each flavor i. In the right panel of figure 1.2, one can see the chiral condensate as a function of temperature, computed using lattice QCD [8]. Note that the 1.1. QCD renormalized chiral condensate is shown, which is defined by the subtraction of a constant such that the chiral condensate at zero temperature vanishes. One can clearly see that indeed the order parameter q i q i , which is small at large temperatures, grows for small temperatures. A subsequent question one can ask is whether the transition between a chirally symmetric phase without confinement, known as the quark-gluon plasma, at high temperatures and the chirally broken confining vacuum are separated by a cross-over or a first order phase transition. In the Columbia plot [9], shown in the left panel of figure 1.3, one can see that the answer to this question depends on the masses of the quarks. One can see that for  [9]. Right: Equation of state at vanishing baryon chemical potential as a function of temperature, showing the pressure p, the energy density e, and the entropy s. The transition can clearly be seen to be a cross-over. Figure taken from [10].
the physical values of the quark masses, the transition is a cross-over. In the right panel of figure 1.3, one can see the equation of state as a function of temperature for vanishing baryon chemical potential [10]. Here too, it is apparent that the transition is a cross-over.
In the presence of a finite baryon chemical potential, the situation may be different. In RHIC and LHC are shown, with center-of-mass energies indicated. The cross-over seen in figure 1.3 is seen in figure 1.4 along the vertical axis. As one moves to larger chemical potentials, the cross-over is expected to turn into a first order phase transition at a critical endpoint. The search for such an endpoint is the purpose of the Beam Energy Scan program at RHIC [12]. If we look at the low temperature, large chemical potential region of the phase diagram, we see nuclear matter indicated. At chemical potential values beyond this value, still at low temperatures, we enter the regime in which the matter making up neutron stars exists. It is not known whether the densities inside a neutron star are large enough to probe a potential phase transition as indicated in figure 1.4, but at some large density it is expected that yet a new state of matter forms, known as a color superconductor [13]. One of the main reasons why these features in the phase diagram are as of yet unknown is the aforementioned sign problem, which precludes a first principles calculation of these features. The two experimental probes into 1.2. RELATIVISTIC HYDRODYNAMICS the phase diagram, namely heavy ion collisions and neutron stars, are both described by relativistic hydrodynamics, which is what we will describe next.

Relativistic hydrodynamics
Relativistic hydrodynamics is an effective theory which describes the behavior of fluids in local thermal equilibrium. It can be described by the conservation of conserved quantities that the underlying theory has. This always includes the conservation of the stress-energy tensor, but can also include other conserved currents, such as baryon number density. 2 Let us take this theory as an example. We then have where J µ is the conserved current associated to baryon number density. This set of equations cannot be solved though, which can be seen by a simple counting argument. Indeed, we have 10 independent components in the stress-energy tensor, and 4 in the baryon current, whereas we have only 5 equations to constrain them. Further input is therefore needed. This comes as no surprise, as it would be rather strange if the behavior of the conserved quantities were completely determined by the conservation laws themselves, and had no dependence on the underlying microscopic theory. The extra input used to close the system of equations is called the constitutive relations, which determine T µν and J µ in terms of the temperature T , baryon chemical potential µ and fluid velocity u µ . The constitutive relations can framed in terms of an expansion in derivatives of T , µ and u µ . Below, we will discuss three examples of such relations, namely that of ideal hydrodynamics with a conserved baryon number density, first order Israel-Stewart theory without a conserved baryon number density, and second order hydrodynamics, also without any conserved quantities other than the stress-energy tensor. These three cases correspond exactly to the three cases which will be used in the remainder of this thesis. We will however restrict ourselves to just a description of these theories, giving just the basic idea of the ingredients used to derive them. There are many different detailed derivations available in the literature, see e.g. [14][15][16]. Before moving on to the examples of constitutive relations however, note that one can easily couple the which we can compare to the known result for a fluid at rest to deduce that we should interpret the arbitrary functions e, P and n as energy density, pressure and baryon number density, respectively. Relating these three quantities through the equation of state, we can close the system of equations, rendering it solvable. These constitutive relations are used for the neutron star merger simulations in chapter 3, where they are solved together with the Einstein equations.
For the remainder of this section, we will disregard J µ , and consider a theory with only a conserved stress-energy tensor. We will also add the first order in derivative corrections. To this end, let us first introduce the following derivatives: where ∇ ν is the gradient in the fluid rest frame, and D is the time derivative in the fluid rest frame. We now write the first equation of (1.2) as T µν = eu µ u ν − (P (e) + Π)∆ µν + π µν .

RELATIVISTIC HYDRODYNAMICS
Here we have removed the dependence on µ as there is no more baryon number density, and added the bulk pressure Π and the shear tensor π µν , where π µν is traceless (π µ µ = 0) and orthogonal (u µ π µν = 0). We have also rewritten e(T ) and P (T ) as a single function P (e). For the bulk pressure and shear tensor we have the following expression in terms of derivatives: Π = −ζ(e)∇ · u, π µν = 2η(e)σ µν , (1.4) which is the most general expression at first order in derivatives which satisfies the second law of thermodynamics [16]. Here ζ(e) and η(e) are the bulk viscosity and shear viscosity, respectively. Both η(e) and ζ(e) are required to be positive to respect the second law of thermodynamics [16]. The σ µν tensor is a symmetric tensor satisfying the same tracelessness and orthogonality conditions as π µν : (1.5) Here we define the angled brackets as symmetrizing a tensor, and at the same time removing the trace. There is one big problem with these constitutive relations though, namely that they allow for superluminal propagation, thereby violating causality. This can be solved in the following way, which may seem ad hoc, but can be derived in several different ways [14]. The solution is to replace the identifications in (1.4) by the following differential equations, called the Israel-Stewart equations [17]: [Π + ζ(e)∇ · u] , (1.6) where the projectors in front of Dπ αβ ensure that the differential equation preserves tracelessness and orthogonality, and the positive functions τ π (e) and τ Π (e) are called the shear relaxation time and bulk relaxation time, respectively. Summarizing, the Israel-Stewart equations give us four parameters, called transport coefficients: η(e), ζ(e), τ π (e), τ Π (e), where the dependence on the energy density, or equivalently the temperature, is indicated. These transport coefficients depend on the microscopic details of 15 CHAPTER 1. INTRODUCTION the theory, and hence encode information about the underlying theory. Since they also enter the equations governing the hydrodynamical evolution, they also have an influence on macroscopic observables, which can in principle be measured experimentally. Second order hydrodynamics is a generalization of the above discussion. The stress-energy tensor is still described by (1.3), but the relaxation equations for the bulk pressure (1.6) and shear stress (1.7) are expanded to the following form: [Π + ζ(e)∇ · u + δ ΠΠ (e)∇ · uΠ (1.8) − λ Ππ (e)π µν σ µν ] , [π µν − 2η(e)σ µν + δ ππ (e)π µν ∇ · u (1.9) − φ 7 (e)π µ α π ν α + τ ππ (e)π µ α σ ν α − λ πΠ (e)Πσ µν .
These transport coefficients can, just as the first order transport coefficients, in principle be derived from the microscopic theory, and they also can potentially be measured experimentally. Both the first and second order constitutive relations will be used in chapter 4, where they will be used to describe the quark-gluon plasma stage of simulations of heavy ion collisions.

Heavy Ion Collisions
At large temperatures, QCD matter undergoes a transition to a quark-gluon plasma (QGP) phase, as can be seen in figure 1.4. Such large temperatures can be achieved experimentally by depositing extreme amounts of energy inside a small volume, and letting the system equilibrate towards a thermal state. For this system to be able to reach a near-equilibrium state though, the spatial extent of the system must be sufficient such that the energy density dissipates away faster than the system takes to equilibrate. If this condition is met, a significant portion of the system's evolution will be described by a QGP, which can be described using hydrodynamics as discussed above. 16

HEAVY ION COLLISIONS
Systems in which this is possible are heavy ion collisions. In a heavy ion collision experiment such as those conducted at RHIC and the LHC, two atomic nuclei are accelerated in opposite directions, and collide inside a particle detector. The resulting matter produced in the collision 'hydrodynamizes' on a timescale of less than 1 fm/c, which is much smaller than the spatial extent of the resulting plasma if the colliding nuclei are large and collide 'head-on'. An interesting question is indeed how small a collision system can be for it to still form a QGP (See [18] for a recent review.). In the following paragraphs, we will describe in some detail the physical processes occuring during a heavy ion collision. We will subsequently end this section with a discussion of experimental observables. See also [19] for a recent review of this topic. For the discussion of the processes occuring during a heavy ion collision, let us focus on a specific example. At the LHC, lead-208 nuclei are collided at center-of-mass energy per nucleon of 2.76 TeV and 5.02 TeV. Different snapshots of an animation of such a collision can be seen in the left panel of figure 1.5. Here the numbers in the bottom of each snapshot indicate the Figure 1.5: Left: Snapshots from an animation of a heavy ion collision taken at 5 fm/c before the collision, at the moment of collision, and at 5 and 50 fm/c after the collision, where the collision is viewed from the side, and where quark-gluon plasma is indicated in red. In the final snapshot the approximate rapidity y (see (1.10)) of particles is shown for y = 0, 1, 6. Figure taken from [19], which was adapted from [20]. Right: Location of nucleons participating in the collision, viewed in the plane transverse to the beam. Nucleons from one nucleus are shown in blue, the ones from the other are shown in red. Nucleons participating in the collision are shown as solid circles, while dotted circles indicate 'spectator' nucleons. Figure taken from [21].

CHAPTER 1. INTRODUCTION
time in fm/c, with t = 0 being defined as the moment the collision occurs, and the collision is viewed from the side, i.e. the beam passes through the figure from left to right. Because the nuclei each have a very large energy, in the lab frame they appear extremely Lorentz contracted, as can be seen in the snapshot at t = −5 fm/c. The nuclei then pass through each other, interacting and leaving matter in their wake. This matter is what we will be discussing the evolution of below.
Before continuing the discussion on the different stages this matter goes through, one important point is that the two nuclei need not collide head-on, as can be seen in the right panel of figure 1.5. In that figure, we are looking in the direction of the beam, i.e. the two dimensions shown are transverse to the beam. Since the nuclei are very small compared to the size of the beam, the so-called 'impact parameter', or the transverse distance between the centers of the colliding nuclei, is essentially random. As a consequence of this, heavy ion collisions as measured in an experiment are not all of the same type, as collisions with a small impact parameter (called central events) are very different from those with a large impact parameter (called peripheral or off-central events). One difference is that the number of participants in the collision correlates with the number of particles measured in the final state, which causes central events to have more particles in their final states. Another difference is in the initial geometry. Lead-208 is spherical on average, and therefore central events are also to a good approximation spherical. Offcentral events like the one shown in the right panel of 1.5 instead are quite elongated.
Let us now consider with the discussion of what happens after the collision. As was mentioned, when the nuclei pass through each other, they leave matter in their wake. This matter is produced, to a good approximation, in a way which is invariant under boosts in the beam direction. At some time after the initial collision, the resulting matter can be described by hydrodynamics. This process is called 'hydrodynamization'. Note that this is different from thermalization, as the matter is at this stage not yet in equilibrium, which shows itself in the fact that the matter is not homogeneous, and large gradients of the stress-energy tensor exist. The process by which this hydrodynamization happens is poorly understood, and even how fast this happens is subject to debate. Kinetic theory suggests that hydrodynamization occurs after roughly 1 fm/c [22]. Holography, which will be discussed in section 1.5, predicts even earlier values of perhaps 0.35 fm/c [23]. Furthermore, different models for this pre-equilibrium stage describing this 1.3. HEAVY ION COLLISIONS hydrodynamization process give qualitatively different answers for the state of the hydrodynamic fluid immediately after hydrodynamization.
The next stage in the description of a heavy ion collision is the hydrodynamical evolution of the fluid created in the hydrodynamization process. For this stage, the physical description is well understood, namely viscous hydrodynamics. What is less well understood are the values of the various transport coefficients entering the evolution through (1.8-1.9). Of the transport coefficients listed, the ones with the most pronounced effect on the experimental observables are the shear and bulk viscosities. This makes sense, because hydrodynamics is a derivative expansion, where higher order derivatives are assumed to be less important for the evolution. The shear and bulk viscosities are the only first order coefficients in this expansion, expressing the fact that they have the most influence on the evolution of the fluid and hence on the final experimental observables. In figure 1.6, the results from a Bayesian analysis are shown, in which among other quantities both these viscosities were fitted to experimental data. Of theoretical interested is 150   In both figures the 90% credible interval is shown. Also shown is the holographic result η/s = 1/4π for the shear viscosity. Both figures taken from [24].
that these viscosities can be obtained by means of holography, which will be discussed in section 1.5. In particular, [25][26][27][28] obtained a surprisingly small 19 CHAPTER 1. INTRODUCTION value for the ratio of the shear viscosity to entropy density ratio: where one should note that the assumption of infinite coupling strength is an important ingredient in the holographic computation. As can be seen in figure 1.6, this value is compatible with the results from the Bayesian analysis for values near the QCD cross-over, where the effective coupling strength is expected to be large. This lends credibility to the idea that holography can be used to at least give qualitative insight into QCD. At some point after the collision, the fluid has cooled and diluted enough so that the interactions can no longer maintain local hydrodynamic equilibrium, and hydrodynamics no longer provides a good description of the fluid. Theoretical models reflect this change by switching to a particle description at a certain temperature called the freeze-out temperature T fr . Note here that this change in description depends not so much on the time, but instead on the temperature. This means that even though the language in this section conveys this process as occuring sequentially in time, the time at which freeze-out occurs is not the same for different regions of the plasma.
After the system has cooled enough so that it is no longer described by hydrodynamics, there are still interactions between the particles, which can be well described by solving a Boltzmann equation. However, as the system expands further, at some point the particles become far enough separated that they no longer interact. After this time, except for the decay of unstable particles, no further interactions occur, and the particles travel in straight trajectories until they are detected. In fact, these final particles are all that can be measured. None of the other processes mentioned can be directly observed, so all conclusions about the processes described above have to be inferred from the final state particles and their correlations. As one can imagine, this is an enormous obstacle towards understanding the processes involved, because the final state typically depends on all of the physical processes involved in the collision.
Before discussing the various observables one can define in terms of the final state particles, let us mention one more physical process, which will be neglected in the rest of this thesis, but should nevertheless be mentioned. During the initial collision, it is possible that two nucleons undergo a hard scattering, creating high-momentum particles. These particles form jets, which subsequently propagate through the medium. This process contains a 20 1. 3

. HEAVY ION COLLISIONS
wealth of information about the medium, but as we will not simulate jets in chapter 4, we will not discuss them in detail.
Let us now move on to the discussion of the observables which can be experimentally measured. This discussion will necessarily be limited to a small subset, but this should give a good impression of the main types of observables and the type of information they carry about the physical processes mentioned above. Here, note that this will be a general discussion, and the precise definitions of observables computed in this thesis and their comparison to experimental data will be done in chapter 4. To start, let us note that the spatial extent of the QGP is only a couple of femtometers, which is too small to be able to measure any spatial information. Hence all observables are defined in terms of the momenta of the final state particles, where subtle differences between observables can be made based on which particles to count, such as conditions on the momenta and particle species.
For these definitions, let us decompose the transverse momentum of each particle in the following way: where p T is the transverse momentum and φ is the azimuthal angle. In addition, we define the rapidity y and the pseudorapidity η: where E is the energy of the particle, and the z-component of the momentum points along the beam axis. Note that in the case of massless particles, we have y = η, and also note that while computing y requires knowledge of a particle's mass, η is a pure angle. Most observables are defined in terms of only particles satisfying certain constraints on their momenta. The main reason for this is that experimentally, detectors are not 100% efficient in detecting every single particle from an event, where efficiencies vary depending on particularly p T and η. One could try and correct for this, but it is easier to just exclude particles from the most inefficient regions from the analysis. Indeed, for theorists it is easy to simply apply the same cuts, and this allows for a cleaner comparison.
With the momentum decomposition in hand, let us now define centrality. As was mentioned in the beginning of this section, the amount of overlap of the initial nuclei is very important, where central events with a small 21 CHAPTER 1. INTRODUCTION impact parameter produce many particles in a roughly spherical manner, while more peripheral events produce fewer particle in a more anisotropic way. Unfortunately, there is no way to experimentally determine the impact parameter, and therefore the 'centrality' is defined in a different way. Since we know that central events produce more particles than peripheral ones, it makes sense to use the number of particles (most often the number of charged particles to be precise) produced by an event as a proxy for the impact parameter. In this way, we determine for each event how many particles it produced, and sort all of them from many particles to few. Then we define centrality by percentiles, i.e. the event with the most particles is by definition 0% central, while the event with the fewest is by definition 100% central, and the other events interpolate between these extremes.
Using the above discussion, we can already define a few observables, such as the number of particles produced per unit pseudorapidity dN/dη and the mean transverse momentum p T . As it turns out, the number of particles produced correlates well with the entropy produced in the initial stage of the collision, because viscous corrections in the hydrodynamical evolution are too small to generate appreciable amounts of entropy, and the final state entropy is proportional to the number of particles. The momenta of the particles produced in the final state are to some approximation those of a boosted thermal ensemble. Because of this, the mean transverse momentum is mostly sensitive to the freeze-out temperature and the velocity of the fluid at the freeze-out surface.
It was mentioned above that the initial geometry of the plasma is generically anisotropic. It turns out that this initial spatial anisotropy is translated by the hydrodynamic evolution into anisotropy in momentum space, specifically in the azimuthal distribution of the momenta of the final state particles. In particular, one can perform a Fourier decomposition of the particle distribution dN/dφ in an event: where v n are called the anisotropic flow coefficients, and Ψ n are the event plane angles. Averaged over a large number of events, the flow coefficients show a pronounced dependence on the centrality and hence on the impact parameter. The reason for this is that the v n , but also correlations between different Ψ n , inherit information about the initial geometry, and this depends 22

NEUTRON STARS
strongly on the impact parameter. This is however not the complete story.
The viscosities tend to smooth out spatial structure. As such, large values of the viscosities tend to lower the final anisotropy present, making v n a probe of especially the shear viscosity. These observables will be discussed in more detail in section 4.3, along with figures of these experimental results compared to theoretical predictions described in section 4.2. Next, we will examine a different corner of the QCD phase diagram, namely neutron stars.

Neutron stars
When a star burns through its supply of hydrogen, it reaches the end of its life. In a sequence in which it starts burning ever heavier elements, it eventually sheds its outer layers to leave behind a compact remnant. The nature of this remnant is determined mainly by the mass of the progenitor star. For stars like our sun, the remnant will be a so-called white dwarf: an object with a mass in the order of magnitude of one solar mass (1 M ) and a radius comparable with that of the earth. Unlike an ordinary main sequence star, a white dwarf is not held in static equilibrium by thermal pressure of gas resisting gravitational collapse. Instead, the degeneracy pressure due to the Pauli exclusion principle of the electrons is what resists further gravitational collapse. There is however a limit to how much mass such an object can have before electron degeneracy pressure becomes insufficient to maintain hydrostatic equilibrium. This is called the Chandrasekhar limit, and is equal to about 1.4 M .
Indeed, if the progenitor star is too massive, the resulting compact remnant is no longer a white dwarf, but a neutron star instead. For small pressures, neutrons are unstable, as they can undergo beta decay into a proton, an electron and an anti-electronneutrino. At extreme densities, however, it is thermodynamically favorable for the protons and electrons inside ordinary matter to merge and form neutrons and neutrinos. This explains the name neutron star, as a neutron star is extremely neutron-rich. Observationally, it is known that the masses of known neutron stars are typically about 1.4 M , but with masses of around 2 M also occuring. The radius depends on the mass, as we see below, and recent experimental constraints by NICER put the radius of a typical 1.4 M neutron star at around 13 km [29]. Finally, we note that neutron stars are cold as far as QCD is concerned. This may 23 CHAPTER 1. INTRODUCTION seem like an odd statement given the fact that they have temperatures on the order of 100 eV [30]. 3 However, since neutron stars consist of densely packed neutron-rich matter, for which the relevant physics is QCD, one should compare this temperature to the energy scale of QCD, which is around 1 GeV. In this unit, we can safely neglect the effects of temperature on the structure of the neutron star and assume T = 0. 4 This means that, given their enormous density and cold temperature, neutron stars occupy the low temperature, large chemical potential region of the QCD phase diagram.
Let us now show that indeed the mass and radius are related. Assuming a non-rotating neutron star, one can assume spherical symmetry. This, in combination with the assumption of hydrostatic equilibrium (T µν = diag(e, P, P, P )), allows us to write down the Tolman-Oppenheimer-Volkov (TOV) equations: where r is the radial coordinate, G the gravitational constant and where the mass m(r) enclosed within radius r satisfies dm dr = 4πr 2 e.
Supplying these two equations with the boundary conditions at the center of the star at r = 0 that m(0) = 0 and that the central density is some specified value ρ * , the differential equations can be integrated to yield P (r) and m(r).
From the solution, we can then identify the radius R of the star as the point where P (R) = 0, and subsequently we can also find the mass M = m(R). 5 In this way, we obtain R(ρ * ) and M (ρ * ) as parametric functions of the central density ρ * . Note though that the TOV equations can only be solved given an equation of state. In [31], a large number of equations of state was generated, which are shown in the left panel of figure 1.7. Here only equations of state were used which are causal, i.e. the speed of sound is less than the speed of light, and which simultaneously satisfy constraints from nuclear matter models at low densities and from perturbative QCD at high densities. In the right panel, the resulting mass to radius relations are shown for the same equations of state.

NEUTRON STARS
One can see several important features. First of all, each equation of state has a maximum allowed mass for neutron stars it supports, in much the same way as we saw above for white dwarfs. It is observationally not precisely known what this maximum mass is precisely, but it is known that a neutron star named J0348+0432 has a precisely measured mass of 2.01 ± 0.04 M [32]. 6 This means that all the equations of state colored blue in figure 1.7 are excluded by this observation, as these equations of state do not support a 2 M star. Similarly, the equations of state colored red are also excluded by observational constraints, this time by the tidal deformability Λ of the neutron stars involved in the binary neutron star merger GW170817 [34]. Let us next discuss neutron star mergers. In 2017, the first neutron star coordinates, which implies that the area of a star of radius R equals 4πR 2 . For the mass, we define that the gravitational mass, as measured by examining the Schwarzschild metric outside the star, is the same as that of a black hole of mass M . 6 An even more massive star, J0740+6620, was detected after the publication of [31], with a mass of 2.14 +0. 10 −0.09 M [33].

CHAPTER 1. INTRODUCTION
merger, GW170817, was discovered using gravitational waves [34], and was accompanied by an electromagnetic counterpart [35][36][37]. In the remainder of this section, we will discuss the various stages involved in such a merger, and how QCD enters the problem. For more detailed reviews, see [38,39]. In figure 1.8, a schematic overview is given of the different stages of a merger event. The first stage of the merger is by far the longest lasting, as it covers the inspiral. This inspiral can quite literally take billions of years, as gravitational wave emission circularizes the orbit and slowly but surely shrinks the size of the orbit. As the orbits shrink, the orbital period does too, and the amplitude of the emitted gravitational waves increases. Only in the last minute or so does this happen to a large enough extent such that observatories such as LIGO and VIRGO can detect the gravitational waves emanating from the source. During the merger, the two stars exert a tidal force on one another, which slightly deforms the stars. This produces an imprint in the gravitational wave emission, the size of which depends on the equation of state through the tidal deformability Λ [40]. When the stars touch, two things can happen. Either the stars are heavy enough that the densities immediately exceed what the equation of state can support, and they collapse to a black hole. In this case, very little material will be ejected, producing only a small electromagnetic counterpart. Also, gravitational wave emission dies down quickly, as the resulting black hole 26 1.4. NEUTRON STARS will ring down with its characteristic quasinormal mode frequencies. The other option is that the stars merge to form a highly deformed object, which loses energy by gravitational wave emission, and ejects a substantial amount of neutron-rich matter. This matter, no longer under enormous pressure, decays to form heavy elements, and emits electromagnetic radiation in the process. The gravitational waves emitted during this phase are characteristic of the equation of state.
As the deformed object circularizes over a timescale of around 10 ms, gravitational wave emission dies down, there is again the possibility of gravitational collapse to a black hole. If this does not happen, the merger remnant will slowly lose angular momentum due to various processes over the course of a few seconds. The angular momentum effectively contributes partly to the pressure preventing the star from collapsing, and therefore as the star spins down, there is again the possibility of collapse to a black hole. If the mass of the merger remnant is below the maximum allowed mass however, the remnant will be a heavier neutron star.
To end this section, let us discuss the methods used to theoretically compute a waveform. Neutron star mergers are described theoretically by relativistic hydrodynamics coupled to general relativity. In figure 1.9, one can see the amplitude of gravitational waves emitted as a function of frequency. The low frequencies are mostly produced during the inspiral phase, while the peaks at high frequencies are mostly produced during the post-merger phase. Also indicated are the methods used, namely analytical methods for most of the inspiral, and numerical relativity for the merger part, where there is an overlap interval in which both methods are applicable. Numerical relativity is immensely computationally expensive, so it would be impractical to have to compute the inspiral using numerical relativity. In this way, the two methods neatly complement each other. Explaining either of these methods in detail is beyond the scope of this thesis, however good introductions can be found in [41,42].

Holography
Holography is a duality between two at first sight very different classes of theories. To illustrate this, let us focus on the first constructed example, namely that of N = 4 super Yang-Mills (SYM) theory living in four spacetime dimensions. For the purpose of this section, this can be seen as a highly supersymmetric version of QCD, where we take the number of colors N c to be infinite. In [25], a convincing argument was made that this theory is the same as type-IIB string theory living in an anti-de Sitter (AdS) background with five spacetime dimensions. In other words, the string theory lives in one dimension more than the gauge theory that it is dual to.
Furthermore, what makes this duality particularly interesting, is that the string theory side of the duality simplifies if in addition to N c → ∞ we also take the 't Hooft coupling λ ≡ g 2 N c to be infinite, where g is the gauge theory coupling constant. When taking this limit, known as the 't Hooft limit, two things happen: The string coupling on the string theory side of the duality goes to zero, leaving us with a classical string theory. Additionally, also the string length vanishes, reducing the classical string theory further to a theory of point-like particles, which in this example is classical type IIB supergravity.
In this way, we obtain a duality between on one side a strongly coupled quantum field theory, which we say lives on the boundary, and on the other side a classical gravitational theory which lives in one dimension extra, which we say lives in the bulk. This is extremely useful, because this duality relates something difficult, namely strongly coupled QFT, to something relatively easy, namely classical general relativity. Note that in [25], the duality was only introduced for N = 4 SYM theory and related theories, and that there is no known general way to obtain a holographic dual for an arbitrary QFT. 28

HOLOGRAPHY
It is expected though, both from the large N c expansion in gauge theory [43] and from black hole thermodynamics [44][45][46], that the class of theories with holographic duals is larger. We will touch further upon this problem of constructing holographic duals in section 1.5.1. Before doing so however, let us discuss some more how the duality precisely works. Indeed, for the two theories on either side of the duality to be equal, one needs a precise dictionary for how to relate quantities and problems on one side to the corresponding quantities and problems on the other side. Such a dictionary has been developed over the years, and in the following paragraphs we will discuss a selection of this dictionary, introducing only the quantities that will be used in this thesis. This discussion will just describe the dictionary without going into the derivations. An excellent review which goes in more detail can be found in [28].
Let us start the discussion of the dictionary with a few thermodynamical quantities. The first of these is the free energy. The starting point for this is the observation that the partition functions of both sides of the duality are equal [47][48][49]. After performing a Wick rotation, and using the fact that the bulk theory is classical, one obtains that the free energy of the boundary theory F obeys where T is the temperature, and S is the on-shell action of the bulk theory. Note here that the on-shell bulk action is divergent towards the AdS boundary. This problem is similar in origin to the UV divergences originating in QFTs, and the solution is similar, namely holographic renormalization [50]. In holographic renormalization, one regularizes the divergence by introducing a cutoff in the bulk spacetime integral for the action. Subsequently, one then compares the desired action to that of a reference solution, after which one can take the limit → 0 for the difference of the two regularized actions. In this way one can compute free energies up to an overall constant, which is for most purposes enough.
Let us next discuss temperature and entropy. In a bulk geometry with a horizon, such as one with a planar horizon called a black brane, one can obtain the temperature as the Hawking temperature of the horizon, which can be expressed in terms of the local metric at the horizon by requiring that the Wick rotated geometry has no conical singularity at the horizon [51]. The entropy can also be obtained purely from horizon data, namely by use 29 CHAPTER 1. INTRODUCTION of the Bekenstein-Hawking formula [44]: where now S is the entropy, and A is the area of the black hole. Note that in the case of a black brane solution, the black hole is infinite in extent, in which case it makes more sense to divide out the volume on the boundary. Indeed, when one examines the entropy density, the result is still finite. The next important element in the holographic dictionary that we will need is the field-operator correspondence [47,48]. Imagine an operator O(x) in the boundary theory that we want to compute, and imagine introducing a source for that operator φ(x). Then the field-operator correspondence tells us that in the corresponding bulk theory there is a field Φ(x, r), with r the bulk coordinate, where r = 0 corresponds to the boundary of AdS. This bulk field Φ then has the following near-boundary expansion: where ∆ is the scaling dimension of the operator O. For most operators the first term will be non-normalizable, and the second will be normalizable. One can see that in this way, one obtains a way to evaluate expectation values of operators in the boundary theory by an equivalent computation in the bulk theory, namely by extracting the subleading behavior of the corresponding bulk field. This result also allows for the computation of Green's functions. For example, by considering the appropriate space-time dependent metric fluctuation δg µν as the source for the stress-energy tensor T µν , one can obtain the Green's function for the stress-energy tensor, leading to the famous result mentioned earlier, namely that the shear viscosity divided by the entropy density of a holographic fluid is equal to 1/4π [26][27][28]. Let us next move on to two non-local operators, namely the Polyakov loop correlator and the entanglement entropy. In QCD, the Wilson line operator where P denotes path ordering and C denotes a closed path, contains information on among other things confinement. The reason for this is that if one takes C to run in the time direction from −∞ to ∞, and if one then takes two 30 1.5. HOLOGRAPHY such loops at a constant distance L from each other, the expectation value of this Polyakov loop correlator is equal to the potential energy stored in the gluon field separating a heavy quark-antiquark pair. The holographic dual of this operator is the on-shell action of an open string in the bulk attached to the path C on the boundary [52,53]. As the action of a string is just the area measured in the string frame metric, the complicated non-perturbative problem of evaluating the expectation value of the Polyakov loop correlator itself is therefore replaced in the holographic dual by the much easier task of finding a minimal surface.
A computationally related quantity to the Polyakov loop correlator is the entanglement entropy. In a QFT, if we imagine dividing the spacetime into a region A and its complement A c , we can partition the Hilbert space as H = H A ⊗ H A c , and define the reduced density matrix for a pure state Ψ ∈ H by ρ A = tr A c (|Ψ Ψ|). We can subsequently define the entanglement entropy as For static spacetimes, [54] proposed that the holographic dual of entanglement entropy is, similarly to the Wilson loop, a minimal surface in the bulk with its ends attached to the boundary of the region A. Important differences with the Wilson loop are that in the case of entanglement entropy, the minimal surface is a codimension 2 surface in the bulk, whereas in the case of the Wilson loop, the minimal surface is a dimension 2 surface. Also, for the entanglement entropy we use the Einstein frame metric, whereas for the Wilson loop, one had to use the string frame metric. The proposition was later generalized to non-static spacetimes in [55], and both propositions were proven in [56] and [57], respectively.

Bottom-up holography: IHQCD and V-QCD
One issue that we have so far glossed over is the fact that even though holography allows for an enormous simplification of certain computations, the theories discussed so far are not QCD. For example, even though N = 4 SYM theory is in essence 'just' QCD with a large number of colors and a lot of supersymmetry, phenomenologically the two theories are quite different, most notably in the fact that N = 4 SYM theory is conformal, whereas QCD is not. In the construction of holographic models to describe QCD, there are two general classes of models. On the one hand, there are the 'top-down' approaches, which includes N = 4 SYM theory, but also the previously mentioned WSS model. In a top-down approach, one starts from a string theoretical construction, and in that way arrives at a precise holographic theory. This has the obvious advantage that in this approach, the holographic dictionary is precisely known, and the general amount of control over the computations is larger. The main disadvantage of such theories is that, like N = 4 SYM, the phenomenological resemblance to QCD is not very good.
An alternative approach is the so-called 'bottom-up' approach, where one tries to construct a holographic model without a derivation from string theory, where the aim is to make the model as phenomenologically accurate as possible. Early examples of this approach are the 'hard-wall' models [82,83], which were followed by the 'soft-wall' model introduced in [84]. In this subsection, we will focus on the IHQCD model, including its extension V-QCD, as this is the model we will be using in later chapters.
In Improved Holographic QCD (IHQCD), a holographic theory is constructed for the fields dual to the tr F 2 and T µν operators in QCD. These fields are the dilaton Φ and the metric g µν , respectively. Note though that in this thesis we will write Φ = log λ. The IHQCD action is the following [85,86]: with R the Ricci scalar, and V g (λ) a potential function. The choice of a non-trivial potential V g allows for breaking of conformality in IHQCD. One can see this as follows: The metric which solves the IHQCD action is, near the AdS boundary, of the form , where the scale factor exp A(r) = 1/r can be interpreted as the renormalization scale. On the other hand, λ can be interpreted as the QCD coupling 32 1.5. HOLOGRAPHY strength. With the appropriate choice for the small λ expansion of V g (λ), one can then make sure that dλ/dA is equal to the QCD β-function in the UV.
Another demand on the potential fixes the large λ behavior of the potential as well. By requiring that the theory is confining and simultaneously has a linear glueball spectrum, the large λ behavior of V g is restricted to be of the form The intermediate behavior of the potentials can still be freely chosen, but can in principle be fixed by computing observables also computed on the lattice. Doing this results in the conclusion that IHQCD can match very well lattice results for pure Yang-Mills [87,88]. IHQCD does not contain quarks. For this reason, it has been extended to include N f flavor D4 andD4 branes, yielding V-QCD [89][90][91][92][93][94]. The V in the name stands for Veneziano, as we take N f to be large, with x f ≡ N f /N c fixed, a limit known as the Veneziano limit [95]. We then obtain the following action in addition to the one for IHQCD (1.11): and where the covariant derivative for the tachyon T is given by Here M N are the corresponding field strength tensors. This action contains 3 new phenomenological potentials: V f , κ and w, which we will discuss shortly. While (1.12) is required in chapter 3, in chapter 2 we can make the simplifying assumption that the non-Abelian 33 CHAPTER 1. INTRODUCTION parts of the gauge fields are zero, and that T = τ (r)I N f , which simplifies S DBI to the following expression: which we will call the diagonalized action. Here V µν is the Abelian component of both the left and right gauge fields, which are also assumed to be equal.
As was done for IHQCD, the potentials are chosen to satisfy phenomenological properties of QCD. For the V f -potential, the UV (small λ) behavior is fixed by requiring that the beta function matches that of QCD for different values of x f . The UV behavior of the κ-potential is determined by the RG flow of the quark mass [93], as well as the behavior at large quark mass [96]. In the IR, the potentials are constrained to reproduce phenomenologically reasonable features in the phase diagram, as well as the properties of meson spectra [94,[97][98][99][100][101][102]. In [103], the potentials were fitted to lattice data, resulting in a holographic model of QCD which matches with known phenomenological constraints as much as possible.

Chapter 2 (Inverse) Magnetic Catalysis in Holographic QCD
Magnetic fields play an important role in two widely studied QCD systems, namely heavy ion collisions, and neutron stars. In peripheral heavy ion collisions, the spectator nucleons, which are charged and are moving close to the speed of light, induce a magnetic field of 10 15 T which, in the appropriate units of the pion mass squared, gives around eB/m 2 π ≈ 10 [104][105][106][107][108][109][110]. In the context of neutron stars, magnetars exhibit magnetic fields of potentially 10 11 T [111], which in units of the pion mass squared is about eB/m 2 π ≈ 10 −3 . Given that one of the salient features of QCD at vanishing magnetic field is its phase structure, it makes sense to study the effect of the magnetic field on the phase structure. In particular, one can study how the phase transition temperatures, as well as the associated order parameters, change as one applies a magnetic field. It was precisely in this context that a surprising effect was discovered.
As was discussed in section 1.1, at low temperatures QCD spontaneously breaks chiral symmetry. At low temperatures, one expects that the order parameter of this chiral symmetry breaking, the chiral condensate, increases as one applies a magnetic field [112][113][114][115]. This phenomenon is called 'magnetic catalysis', and the reason for this is that Landau quantization leads to an effective reduction from 3 + 1 to 1 + 1 dimensions. In lower dimensions, the gauge theory IR dynamics are stronger, leading to a strengthening of the chiral condensate.
However, when lattice studies were done, surprisingly the opposite effect was seen in around the crossover temperature, and this effect was named  [118]. Right: Chiral cross-over temperature as a function of the magnetic field B. Figure taken from [116]. function of B for fixed temperature. One can see that for small temperatures, one sees magnetic catalysis, whereas for larger temperatures the condensate instead decreases with B. In the right panel, a related quantity is shown, namely the cross-over temperature as a function of B. One can see that the cross-over temperature decreases with B, signalling inverse magnetic catalysis. Given that IMC seems to require strong coupling to exhibit itself, it is natural to study this effect in holography, where we can hope to obtain a qualitative understanding of the mechanism leading to IMC.
In this chapter, which is based on [120][121][122] as well as upcoming work with the authors of [122], we will address these questions. While the questions asked in each of these papers are different, the model and the methods used are quite similar. So similar in fact, that it is possible to write down a 'master' model, which contains each of the models used in the papers that this chapter is based on by taking appropriate limits. In section 2.1, we will discuss this master model, as well as how to obtain useful information from it. This should allow for a streamlined treatment of computations that would otherwise have to return in slightly different setups in sections 2.2 through 2.5. The strategy for solving the model parallels [99], where of course the 36

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
discussion is slightly modified because our master model is more general than the one considered there. Also, in a few places, it was necessary to make nontrivial adjustments to the analysis. Wherever this occurs this will be clearly stated.

Analysis of V-QCD model in the presence of a magnetic field and anisotropy
In this section, we will go through the computations necessary to obtain the relevant observables for sections 2.2 through 2.5. This section is written with the aim of providing the reader with a practical guide as to how to perform these computations. As such, it necessarily contains a lot of technical details, which are required for the computations. The rest of this chapter has been written to only use the results from this section, and not the computations themselves, so it should be possible to follow the rest of this chapter without having read this section.

Extending V-QCD to incorporate a magnetic field and anisotropy
To study magnetic fields in V-QCD, we consider the diagonalized action 1.13. A magnetic field in the x 3 -direction can then be introduced by changing the ansatz for the Abelian gauge field to [120,121] where we recall that Φ was dual to the baryon chemical potential. This ansatz makes one important assumption, namely that all quark flavors have are identical, and in particular, that they have the same electric charge. In nature, this is of course not the case, but this assumption allows us to consider only the Abelian part of the DBI action, greatly simplifying the analysis. In addition to a magnetic field, we will also be adding an axion field χ to the action [122,123], where we will use the following ansatz: In this way, we can introduce an anisotropy to the system in a way that is different from a magnetic field. With this ansatz, the axion is dual to a

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
space-dependent theta term. As the axion will only appear in the action through a derivative, the presence of the axion does not break translation symmetry. Instead, it breaks rotational symmetry. To determine the remaining symmetry left over from rotations, we can distinguish two cases: one in which the axion is parallel to the magnetic field (a = 0, a ⊥ = 0) or where there is no magnetic field, and all other configurations. In the first case, the remaining symmetry is given by axial symmetry around the x 3 axis, whereas in the latter case the rotational symmetry is completely broken. As will become clear below, it turns out that to be able to use a diagonal ansatz for the metric, we have to choose either a ⊥ = 0 or a = 0 if a non-zero magnetic field is present. Also, for notational convenience, whenever no magnetic field is present a ⊥ and a will both be denoted a, since in this case the orientation of the axion is irrelevant.
With the addition of the magnetic field and the axion, the V-QCD action becomes 3) where V µν is the electromagnetic field strength tensor for the gauge field given by (2.1), and the potentials V g , V f , κ, w and the axion potential Z are given in appendix A.1. We will keep using these potentials throughout this chapter. In the rest of this section, it will be explained how this model can be solved to obtain black hole solutions which are dual to a QGP-like phase.
It is sufficient to focus on solutions containing a black hole, as horizonless solutions must always be obtained from a black hole solution where a limit is taken that lets the horizon shrink to zero size. This requirement ensures that the IR singularity contained in such a horizonless geometry is of the 'good' type, as discussed in more detail in [124].

Equations of motion and boundary conditions
To obtain the equations of motion, we first choose the following ansatz for the metric: 1 which contains the anisotropy factors U and W . These are necessary in order for the Einstein equations to be consistent. 2 Note that this is also the reason that if a magnetic field is present, either a ⊥ or a must vanish, because otherwise one of the Einstein equations cannot be satisfied. This problem could be remedied by choosing instead a more general metric ansatz. However this would greatly complicate the analysis, while the additional physical insight from allowing for a general angle between a and B would probably be limited. Note further that the metric contains a blackening factor f . This allows for the existence of a black hole horizon at the location where f = 0.
Before stating the Einstein equations and the equations of motion for the dilaton, tachyon and Φ field, note that in principle B and a also have equations of motion, so we are not completely free to choose any ansatz for them that we want. It is important therefore that we check that our ansätze (2.1), (2.2) are consistent with these equations of motion, and it turns out that this is indeed the case. Using the metric ansatz (2.5) we can write the Einstein equations as follows: where we use a dot for derivatives with respect to r, a convention we will keep throughout this chapter. We also define wheren is an integration constant that arises from integrating the Φ equation of motion. Using the above definitions, the Φ equation of motion can now be written as what is essentially a simple integral.
Even though Φ can be integrated out in favor of the integration constantn, we still need to integrate its equation of motion, because as we will see below the value of the chemical potential equals the difference of Φ evaluated at the boundary and at the horizon, necessitating that we evaluate this integral. Lastly, to complete the system, we also have equations of motion for λ and

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
τ : Observe that we have 8 equations of motion for 7 degrees of freedom, so in principle this system could be overconstrained. One can check however that (2.7) is a constraint, by taking the derivative of the right hand side, and using the other equations of motion to show that the derivative of (2.7) is automatically zero. This implies that if (2.7) is satisfied for one particular r, it is satisfied for any r. Therefore this equation will be trivially solved, provided that we choose the proper boundary conditions. Before stating the boundary conditions, it is useful to write the equations of motion in another form. The reason why this is useful is that near the boundary, which in r-coondinates is located at r = 0, A(r) grows like log r. Numerically this poses a problem, as this behavior makes it difficult to satisfy the boundary conditions at the AdS boundary to a good accuracy, and as we shall see below, the observables that we are interested in require the boundary conditions to be precisely met. The solution for this is to use A

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
as the independent variable instead of r. We can do this as long as A is a monotonic function of r throughout the bulk. Interestingly, this is not always the case. For this reason, it is prudent to use r as the independent variable near the horizon, and to do a coordinate transformation at some point in the bulk in order to use A as the independent variable near the boundary.
Since A satisfies a second order differential equation in r-coordinates, to perform the transformation to A-coordinates one has to introduce q ≡ e A dr/dA. Using this definition, one obtains for the Einstein equations: where G is now given by and where we introduce the notational convention, which will be used for the rest of this chapter, that a prime denotes a derivative with respect to A. For the remaining equations of motion, one obtains:

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
In addition to satisfying the equations of motion, solutions must also obey boundary conditions. We will first discuss the boundary conditions that need to be imposed at the horizon, before moving on the the boundary conditions at the boundary. The first of these boundary conditions is of course that f h = 0, where a subscript h will from now on always denote a quantity evaluated at the horizon. This first boundary condition simply follows from the definition of a black hole horizon, as this makes an observer stationary at the horizon move on a lightlike trajectory. We will also assume A h = W h = U h = 0 andḟ = 1. 3 While it may seem strange to just assume this, it turns out one can do this without loss of generality. The reason for this is that the solutions are invariant under symmetries which can be used to rescale the solutions to satisfy the boundary conditions at the boundary. For this reason it is irrelevant which choice we make for these assumptions, as any different choice will later be absorbed by these symmetries. This will be discussed in more detail in the next subsection. The last remaining boundary conditions are consequences of the horizon being just a coordinate singularity. This can be imposed by requiring that all variables are smooth

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
at the horizon. Taking equation (2.6) as an example, this means that we must require in particular thatÄ is finite. Given that f h = 0 by definition, to make sure thatÄ is finite we must make sure that all terms inversely proportional to f cancel. This leads to the following condition: Similar arguments for the other equations of motion yielḋ (2.10) Note here that even though (2.7) does not contribute a non-trivial equation of motion in the bulk, it is important to take it into account in the boundary conditions, so in particular the constraint onλ h comes from (2.7). The last boundary condition at the horizon that is required, is that Φ h = 0. This is needed because Φ h is a component of the gauge field, and in the Euclidean geometry the gauge field would not be continuous unless Φ h = 0 [125]. Applying these boundary conditions one is left with the freedom to choose λ h and τ h . These, together with B,n, and either a ⊥ or a , determine the entire parameter space of allowed solutions. For the boundary conditions near the boundary, one needs to consider the asymptotic behavior of the variables near the boundary. This asymptotic behavior essentially boils down to that to leading order the geometry is AdS, and it turns out that one can analytically expand around this ansatz for small r. For the f , U and W , the result of this procedure is that these variables approach constant values. Subleading corrections come in at O(r 4 ), and for

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A
MAGNETIC FIELD AND ANISOTROPY all works described in this thesis approximating them as constants is good enough. In order to make sure that on the boundary the t, x i coordinates agree with the familiar coordinates for Minkowski space, we require that where a subscript b will from now on always denote a quantity evaluated at the boundary. The near-boundary expansion for A and λ are a bit more complicated, they are given by [94]: where L UV , V 1 and V 2 are determined by the potentials as discussed in appendix A, and where Λ is an overall energy scale. It turns out that all quantities one might want to compute scale with Λ to some power, so in practice we just put Λ = 1 and if desired we can rescale later. For the next subsections, it turns out that it is convenient to combine both of these equations, to write A as a function of λ. Doing this, one obtains [99]: where b 0 and b 1 are given by the potentials and can be found in appendix A. The last near-boundary expansion that we will need is that of the tachyon [94]: where m q is the quark mass and qq is the chiral condensate. γ 0 /b 0 is given by the potentials, and can be found in appendix A. In the following we will consider massless quarks, so we will be imposing m q = 0 as the UV boundary condition for the tachyon. How this can be achieved will be detailed in section 2.1.5. With this, we now have a complete list of all the equations of motion, as well as the boundary conditions that we will need. In the next section, we will describe a set of symmetries of these equations of motion that are necessary to compute solutions. After that all the ingredients are set to describe the algorithm for obtaining the geometries, and finally the observables that one is after.

Symmetries of the equations of motion
Numerically, one of the simplest and best known methods to solve ODEs such as the ones describing this model is to initialize a solver for a specific value of the independent variable (in this case r or A), and then integrate the equation from that value to the entire domain one is interested in. However, we have boundary conditions on both the horizon and the boundary, and initializing the system of equations at one of the two locations by no means guarantees that the boundary conditions at the other location will also be satisfied. While there are methods available for solving such problems, it turns out we can use symmetry properties of the equations of motion to mostly overcome this issue. This will allow us to initialize the system of equations at the horizon, integrate towards the boundary, and rescale the solution such that the boundary conditions at the boundary are also met.
The symmetries of the equations of motion are essentially the diffeomorphism invariance that is left over after choosing the metric ansatz (2.5). It can easily be verified that the following five transformations leave the equations of motion invariant: • Shift of r: r → r + δ r .
• Shift of A: • Shift of U : • Shift of W : • Scaling of f : Together, δ r , δ A , δ U , δ W and δ f will denoted 'symmetry parameters' for the remainder of this chapter. Note first that these symmetries indeed justify our assumptions that A h = W h = U h = 0 andḟ = 1, as these choices can just be absorbed into the various deltas defined above. Next, observe that after generating a solution that satisfies the horizon boundary conditions, one can choose δ f , δ U , δ W and δ r to ensure that f b = 1, U b = 0, W b = 0, and that r b = 0. In equation (2.12), only the left hand side transforms under these transformations, and it is easy to see that by using the appropriate δ A , one can make sure that (2.12) is satisfied. Summarizing, we can satisfy all boundary conditions except that of the tachyon, namely that the quark mass vanishes. This issue will be addressed in section 2.1.5. Lastly, note that there is a price to pay for using these symmetries. Quantities like B, a andn enter in the transformations. This implies that while we are guaranteed to get a solution that satisfies the correct boundary conditions on both sides, we have no direct control over for instance the value of the magnetic field we want to compute the solution of. This need not necessarily be a bad thing though, because we are guaranteed that together with λ h , all possible values of B,n and a before any rescaling happens span the space of all possible solutions. If it is feasible to produce solutions which explore this entire parameter space, then we will be guaranteed to find all possible solutions for all possible rescaled values of B,n and a as well. In some of the setups that will be described below this is indeed the case, but in other cases it is necessary to fine tune the unrescaled values to produce the desired rescaled values. The procedure for doing so will also be described in section 2.1.5.

Computing observables
Now that the equations of motion, their symmetries, and the boundary conditions are established, we can describe how to obtain the solutions, and how

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
to extract useful information from these solutions. As it turns out, many quantities of interest can be expressed purely in terms of quantities defined at the horizon, and the symmetry parameters described in the previous section. The reason for this is that these quantities do not require the solving of any additional differential equations or integrals on top of the ones already mentioned in section 2.1.2. 4 The quantities for which this is the case will from now on be called 'background' observables, and their computation will be discussed first. There are also observables which do require the solving of additional equations of motion. In principle, one could follow the same computation scheme as for the background observables, but as this would require solving the same background equations of motion multiple times with the same boundary conditions, it is more efficient to solve the background once, and then solve the additional equations afterward. The computations of these non-background observables will be discussed towards the end of this subsection.

Background observables
In this subsection, we will focus on the background observables, where for convenience we will introduce the notation that quantities with a tilde denote quantities before the symmetries are used to impose the boundary conditions at the AdS boundary, and quantities without a tilde do correspond to the quantities with the proper boundary conditions imposed. Before moving on to the computation of these obervables, observe that the property of only needing symmetries and horizon data is extremely useful. Because the required symmetry transformations only depend on the non-rescaled solution that one has computed at the boundary, there is actually no need to retain the information about the bulk geometry at all, and therefore it can be immediately discarded. This does away with the need for reading/writing a lot of data to memory, making the computation faster. Of course, if one is interested in one of the quantities that do not have this property, this shortcut cannot be taken, but since the quark-antiquark potential and the entanglement entropy will only be computed for zero temperature solutions, it turns out that for a majority of the computations in this chapter, the shortcut is possible.

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
From the above discussion, the strategy for computing observables becomes clear. Choose λ h , τ h ,ñ,B andã, and then initialize according to the horizon boundary conditions described above. Subsequently integrate the equations of motion, switching from r-to A-coordinates once close enough to the boundary, and then integrate further up to some large enough value of A. Then extract the symmetry parameters δ A , δ U , δ W and δ f , and use them to evaluate the desired observables. 5 The next few paragraphs will describe how to compute the following observables in this way: • Entropy density s, • Baryon chemical potential µ, • Magnetic field B, • Anisotropies a and a ⊥ , • Baryon number density n, • 'Anisotropization' M a . This is the analog of magnetization for the anisotropy; • Quark mass m q , • Chiral condensate qq .
The first 5 of these are straightforward, whereas the latter 5 require a bit more work. We will now go through each of these observables in order, after which there will be a discussion on how to accurately determine the required symmetry parameters. The temperature is given by the Hawking temperature associated to the horizon. Using the metric ansatz 2.5, this can be expressed as where the Λ-dependence enters because of δ A , as we will see below the discussion of the observables that we can only extract δ A − log Λ from the model. Moving on to the entropy density, note that the entropy is given by the area of the black brane. Once again using the metric ansatz 2.5 and dividing out the overall volume factor, one obtains the entropy density where one may note that this result is different from literature by a factor 4. This amounts to choosing M 3 N 2 c = 1/4π in the action (2.3, 2.4), and is done for notational convenience. If desired, desired values of M 3 N 2 c can be reinstated by multiplying s, n, M and M a by appropriate factors. The chemical potential is given by the value of Φ at the boundary. As Φ can be shifted by a constant due to gauge symmetry, naively one would say that this is ill-defined. However, as discussed before in section 2.1.2, Φ h = 0 to ensure continuity of the gauge field at the horizon. Therefore we obtain: The magnetic field and the anisotropy can be obtained as follows: Next, we move on to compute the number density, magnetization and anisotropization. These quantities have in common that they are computed by taking derivatives of the on-shell action. Recall that where Ω is the grand potential. By definition, the number density is given by

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
which means that we can express n in terms of the action. To make this more concrete, consider the variation of the action with respect to Φ: where the last equality holds because the action is on-shell, and where the V /T factor comes from integration over t and x i . 6 As it turns out, the integration constant that we defined asn obeyŝ allowing us to write where we keep δΦ h = 0, and recognize δΦ b as being an infinitesimal change in baryon chemical potential. Combining the last equation with (2.15) and (2.16), one obtains that n =n/4π. Applying the appropriate rescalings, one then obtains where in the last step we used (2.14).
For the magnetization, defined by M V = −dΩ/dB, one can perform a similar computation to obtain the following integral: In this case there is no way to analytically evaluate the integral, so it has to be integrated along with the other equations of motion. Note however that M factorizes into a factor containing all the rescalings, and an integral that depends only on unrescaled variables. This allows us to still perform the integral first, without the need to retain any intermediate quantities. One further point of interest is that M is a divergent quantity, as the integrand

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
diverges near the boundary. This means that in practice we can only compute differences in M between solutions which have the same B, by cutting off the integral at some non-zero value of A, and comparing the differences. 7 Also, defining the anisotropization M a V = −dΩ/da , we can obtain in a similar fashion that where we note that U = 0 in the second case, whereas in the first case both a non-trivial U and W can occur if both a ⊥ and B are non-vanishing. The last quantities we would like to extract are the quark mass and the chiral condensate. Both of these can be extracted from (2.13). Given that the term containing m q becomes dominant near the boundary, it can be extracted by simply ignoring the term containing qq , and rearranging: The procedure for extracting the chiral condensate is slightly more involved.
The problem is of course that the qq term gets smaller relative to the m q term as one gets closer to the boundary. In [120], we developed a relatively simple way to solve this problem. The trick is to divide (2.13) by r[− log(rΛ)] −γ 0 /b 0 , so that the m q -term is to leading order constant. Subsequently taking a derivative on both sides and rearranging, one obtains the following expression: This works well enough to extract the chiral condensate reliably, but one has to be careful not to do the extraction too close to the boundary. This

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
is because (2.19) essentially computes the difference of two numbers, where that difference tends to zero as one approaches the boundary. For this reason, one eventually runs into accuracy issues.
To conclude this subsection, we will discuss how to obtain the symmetry parameters required for the computation of the observables discussed above. As an example, consider δ f . The blackening factor f does not change under any of the symmetry transformations listed in section 2.1.3 except for the one parameterized by δ f . This implies that after having computed an unrescaled solution by shooting from the horizon, one can extract δ f solely by using information from the blackening factor. Since f approaches a constant value with subleading corrections only entering at O(r 4 ) as one approaches the boundary, one can simply divide the unrescaledf and rescale such that it becomes zero: As the anisotropy factors U and W approach constants in the same way f does, δ U and δ W can be extracted in an analogous way by demanding that U b and W b both vanish to obtain The remaining symmetry parameter that is important for computing observables, δ A , is slightly more complicated to obtain. The key is to use 2.12, which, by substituting A =Ã + δ A , can be written in the following form: Using this equation, we can find δ A − log Λ as long as we evaluate the right hand side close enough to the boundary. Note that this is the reason why in the observables computed above every power of exp δ A came with a factor of 1/Λ. We can now either choose to compute observables in units of Λ, which is something we will mostly do in this chapter, or we can try to use a physically reasonable value for Λ, which is something we will do in chapter 3.
We will conclude this section with an improvement to the computation of δ A as described above, developed in [120]. Note that in (2.20) the subleading corrections go like λ. Also recall that λ goes to zero on the boundary by Taking these two facts together, one can see that if we linearly extrapolate δ A as computed by (2.20) to λ = 0, then the subleading linear correction should cancel. We tested that this works in practice, against the code used in [99]. We compared specifically that with the extrapolation, where one extrapolates linearly using A = 18 and A = 20 to compute two values of δ A , one gets a 0.7% deviation from a solution computed up to A = 400 without extrapolation. Furthermore, the deviation is an overall deviation which is always the same regardless of input parameters like λ h . This means that such a deviation is acceptable, as it effectively amounts to a slight redefinition of Λ. The linear extrapolation procedure presents an immediate advantage in computation time because it is computationally a lot cheaper. However, as it turns out, for (2.19), the accuracy issues one always eventually runs into as one gets closer to the boundary imply that for this observable one has to do the extraction at a value of A for which subleading corrections are sizable. For this reason, the extrapolation procedure is also done for (2.18) and (2.19), as both of these have subleading behavior which is linear in λ.

Non-background observables
We will now continue with what we defined to be 'non-background' observables. These observables will all be computed on top of a background metric, dilaton field and tachyon, where we can assume that the backgrounds have been properly rescaled using the symmetries, so that all the boundary conditions are satisfied. We will discuss the following observables: • Helicity 2 glueballs, • Quark-antiquark potential, • Entanglement entropy.
Let us start by discussing the helicity 2 glueballs. To compute the spectral density associated to these glueballs, one needs to examine the behavior of the δg 12 metric perturbation: where by an appropriate coordinate transformation we may assume q µ = (ω, 0, 0, q), and where the sum in the plane wave term goes over the time and These equations can also be put in the Schrödinger form, by defining h(r) = e −3A(r)/2−W (r)/2 ψ(r): Here we can easily determine whether the spectrum is discrete by looking at the asymptotics of V s . In the UV, V S diverges as A ∼ − log r, and if in the IR V S diverges as well, the spectrum will be discrete, otherwise it will be continuous.
To extract the spectral density, we recall that it can be obtained from the correlator of the energy-momentum tensort T 12 , which we can obtain from h(r) by the following near-boundary expansion: Here h has to satisfy infalling boundary conditions in the IR ψ(r) ∝ e iωr . We then want to extract the imaginary part of G(ω, q), which is equal to the spectral density. In practice, numerically, it is easier to define h(r) = e k(r) , and solve the corresponding equations of motion for k. This simplifies the IR boundary conditions to k(r) = iωr, and more importantly it removes the oscillatory behavior that h will usually have. Furthermore, the O(r 2 ) terms in the near-boundary expansion drop out as they are real, and allow one to write: where we note that this expression has the added advantage of being insensitive to h having the correct normalization to 1 at the boundary. Note lastly that because the equations of motion for h are linear, we can easily demand that the constant term in the near-boundary expansion equals 1.
As mentioned in the introduction, the quark-antiquark potential indicates whether it is possible to pull a quark-antiquark pair apart. In holography, this quantity is computed by evaluating the on-shell Nambu-Goto action of

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
a static string hanging in the bulk from two points on the boundary [126]. In principle, this is a divergent quantity. However, the difference of such a solution with a solution of two strings extending infinitely deep into the IR with non-holographic coordinates fixed is finite. By varying the distance between the endpoints, one can then examine the free energy, and hence the potential, associated to each quark-antiquark distance. Throughout this discussion, it is important that we use the string frame metric for the computation, as this is what the string 'feels'. This amounts ro replacing the A metric factor by A S = A + 2 3 log λ. As it turns out, it is not necessary to solve equations of motion for the Nambu-Goto action, as it has been worked out that the on-shell action satisfies [86,127]: where T f is the string tension, and r F is the turning point of the string in the bulk. This allows us to compute V (L) as a parametric function by choosing a value for r F and integrating. Also note that while the integrals above are given for V , one can obtain the integral for V ⊥ by setting W = 0. The last observable we will be examining is the entanglement entropy of the following two regions: • A region, A, defined by 0 < x 3 < L, where we note that x 3 is in the parallel to both a and B. We denote the entanglement entropy of this region by S E, .
• A region, B, defined by 0 < x 1 < L, where we note that x 1 is perpendicular to any source of anisotropy, be it a or B. We denote the entanglement entropy of this region by S E,⊥ .
The entanglement entropy can be computed in holography by finding a minimal surface with its endpoints fixed to the boundary of region A or B, respectively. Such a surface is called a Ryu-Takayanagi (RT) surface [54]. Given the symmetries of both regions considered, this results in a similar

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
computation to that for the quark-antiquark potential, with the difference that for the entanglement entropy we need to minimize the surface in the Einstein frame instead of the string frame. Using the same techniques as for the quark-antiquark potential, one can find both S E, and L as functions of the turning point in the bulk r F : 8 with A V an infinite factor arising becaus the region is spatially infinite in two dimension. Analogously, one obtains for S E,⊥ and L:

Satisfying constraints
With the above discussion, it is now possible to obtain a solution and extract observables given the inputs λ h , τ h ,B,ñ and eitherã ⊥ orã , which together parameterize the available space of solutions. We will now move up a level of abstraction, treating the entire discussion above as a function which takes in these inputs and outputs the observables listed in section 2.1.4. Here we must also emphasize that for some of these inputs this function will fail to return an answer, because it is not guaranteed that every set of inputs correspond

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
to a valid solution. All of these inputs are given at the horizon, and there is no guarantee that the geometry near the horizon smoothly connects to a near-boundary geometry. For example, for some of these solutions A might grow to negative infinity instead of positive infinity, or some of the other variables might hit poles at some point in the bulk. One issue that arises now is that while we have an efficient way of evaluating this function, we would rather have control over some of the outputs. In particular, even one of the boundary conditions, namely m q = 0, is not satisfied automatically. Also, it is sometimes needed to for example look at solutions at fixed B or a, which requires to fine-tune the inputs to reproduce the desired outputs. In this subsection, we will explain the solutions to both of these problems, starting with fixing m q = 0.
Requiring that m q = 0 can be done in two ways, namely setting τ h = 0, and finding a non-zero τ h for which m q = 0. While this seems like a trivial distinction, making this distinction immediately leads to finding out whether the chiral condensate is zero or not, which is one of the order parameters we're interested in. In the case where the chiral condensate is non-zero, it can immediately be concluded that such a solution spontaneously breaks chiral symmetry since m q = 0, which is one of the main features of V-QCD. First, note that setting τ h = 0 immediately implies thatτ h = 0 by (2.10), which by the equation of motion for τ , (2.9), implies that τ = 0 everywhere. Then from the near-boundary expansion (2.13), it immediately follows that m q and qq both vanish.
Finding a non-zero τ h for which m q = 0, keeping the other inputs fixed, is a more complicated procedure. The algorithm for doing so must determine where a solution with m q = 0 exists at all. If it does, it must choose the largest such value of τ h , and then accurately determine the outputs for that value of τ h . The reason for choosing the largest such value is that one needs to find, given T , µ, B and a, the value with the smallest grand potential. It turns out that the value for τ h which corresponds to this most stable solution is the one with the largest τ h [93,94]. The procedure to find this value is done using the steps described below, where any failure in any of the steps will lead to the conclusion that no solution with m q = 0 exists for a non-zero τ h . 9

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
This algorithm is mostly the same as in [99], except for the last step, which is modified in a non-trivial way. Also, several minor things like the amount of iterations before an algorithm terminates have been changed. These changes turn out to have little effect on the accuracy of the algorithm, but have a large effect on execution time, as every time we generate a bulk geometry and extract the observables takes a few milliseconds. For this reason iterating over this process is very computationally expensive, and lowering the amount of iterations saves a lot of computation time.
1. The first step involves finding any τ h for which a bulk solution exists, without any requirements on m q . This is done by starting with τ h = 1, and doubling τ h until a solution is found, which we call τ h,exist . If no solution is found after 30 iterations, and the procedure terminates.
2. We next need to find a solution τ h,min for which m q < 0. Having found at least one solution τ h,exist , we check the sign of m q (τ h,exist ). If it is negative, this step is complete. If it is positive, we bisect the interval (0, τ h,exist ). If m q < 0 for this solution, we found τ h,min and continue to the next step. If m q > 0, we continue searching in the interval (0, τ h,exist /2). If for this value of τ h there exists no solution, we continue searching in the interval (τ h,exist /2, τ h,exist ). We keep bisecting in this way until we find τ h,min , or, if we haven't found τ h,min after 30 iterations, we terminate the procedure, and conclude that no solution with m q = 0 exists.
3. After finding a value of τ h for which m q < 0, we need to find a value τ h,max for which m q > 0. Doing so guarantees, by continuity, that we will at least find one solution with m q = 0, so after this step, the algorithm can no longer terminate with a failure to find a solution. If m q (τ h,exist ) > 0, this step is trivially completed, and otherwise we keep doubling τ h,exist until a solution with positive m q is found. As in the previous steps, if after 30 iterations no acceptable solution is found, the algorithm terminates with a failure.
4. The next step is meant to find values of τ h bracketing the largest value of τ h for which m q = 0. This means that we have to search for zeroes between the zero that we already found and some large value of τ h . It is known that for large values of τ h the quark mass increases asymptotically [94], so what we do is to increase τ h stepwise until we see this asymptotic behavior. These steps need to be small enough so as to not miss zeroes, but must also not be too small, as this would slow down the algorithm considerably. A good compromise is to take the current highest τ h value, named τ h,last , and compute m q at τ h,last , τ h,last + ∆τ h and τ h,last + 2∆τ h , for an appropriate ∆τ h . If m q changes sign between any of these three values, we found a zero, and the τ h -value larger than where the zero is located becomes the new τ h,last . If m q does not change sign, we take τ h,last + 2∆τ h to become the new τ h,last .
To completely describe this step, we need two more things, namely what to choose for ∆τ h , and when to terminate the search. At the start of this step, ∆τ h is initialized as Also, ∆τ h is updated every time we find a new τ h,last . If we found a new τ h,last because we found a new zero, then we set ∆τ h to the difference between the last found zero and τ h,min , where we divide by 5 so as not to miss any zeroes. If we found a new τ h,last without a new zero, i.e. because m q (τ h,last + n∆τ h ) has the same sign for n = 0, 1, 2, then we construct a parabola through these three points, and look at its zeroes. If it has zeroes, ∆τ h is set to be the difference between the two zeroes, if there are no zeroes ∆τ h is doubled from its last value.
Lastly, the search is terminated if the following three requirements are met at the same time: • m q has not decreased in the last iteration.
• m q is more than 100 times larger than the largest m q found in all previous iterations.
These three conditions correspond to the ones described in [99]. Sometimes, however, this never terminates, so for these cases I have added the condition that the search is also terminated if τ h > 10 5 , which is large enough that the search can be safely terminated.
5. Now that we have found values of τ h which bracket the largest zero of m q , it is possible to iteratively find the zero. For this, [99] uses Brent's

ANALYSIS OF V-QCD MODEL IN THE PRESENCE OF A MAGNETIC FIELD AND ANISOTROPY
method, which combines the best properties of Newton's method, which doesn't always converge, but converges fast if it does, and the bisection method, which is guaranteed to find a zero, but does so more slowly.
In [121], we improved this algorithm slightly. Instead of using Brent's method to reach some small enough value of m q , we use Brent's method to reach a much larger interval. We then compute at 10 intermediate values within this interval, and perform a least square fit to a linear function, of which we can then find the zero analytically. This enables us to not only obtain the desired zero of m q , but also obtain an error estimate for quantities like T and s. We shall see in the rest of this section that this is crucial to obtain quantities like the speed of sound numerically.
With the above discussion, one can obtain a solution satisfying the boundary condition m q = 0, either by setting τ h = 0 to find a chirally symmetric solution, or by following the algorithm described above to obtain a chirally broken solution. To conclude this section, let us go yet one level of abstraction higher, and view the result of the discussion above as two functions (one for chirally symmetric and another for chirally broken solutions) which take in λ h ,ñ,B and either a ⊥ or a , and output the desired observables for a solution which satisfies m q = 0. The discussion below will explain how we can iteratively call these functions to fix B or a to some specified value by fine-tuningB andã, respectively. As the algorithm works the same for both chirally symmetric and chirally broken solutions, I will not distinguish between the two.
The algorithm is similar in setup, albeit simpler, than the one described above. One starts by finding two values of the input parameters which bound the desired solution. This is simpler than the procedure for finding τ h such that m q = 0, as generically B(B) and a(ã) are single-valued functions. For B, this is done by starting fromB = 0, thenB = 1, and then doubling until two values bracketing the desired solution are found. For a, one starts from a = 0, and incrementsã by 0.1 until either one finds values bracketing the solution. For both these cases, aB orã for which no solution exists will be interpreted as potentially bracketing the true solution if its 'neighbor' in the algorithm does exist. Empirically, it turns out that this scheme usually is able to find values bracketing the true solution.
Once two values which potentially bracket the true solution are found,

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
one once again uses Brent's method to shrink the size of these brackets. 10 As soon as Brent's method has shrunk the brackets to a small enough size, we once again perform a fit over intermediate values in this interval, and this is where the value of being able to compute an error bar on observables in determining τ h shows its value. As it turns out, the accuracy of observables varies by as much as an order of magnitude for different values of τ h andB or a. This also means that the accuracy of observables obtained from the fit to m q = 0 can vary by an order of magnitude. Now, having error bars, one can perform a weighted fit, in which the 'bad' points will receive a lower weight, enabling us to obtain a more accurate result. This turns out to be crucial to determine quantities like the speed of sound, which ultimately depend on quantities like dT /dλ h . For such an observable to not be completely dominated by numerical noise, one needs to be able to determine T and other such observables to as high an accuracy as possible.

Obtaining a phase diagram
The above discussion enables us to construct both chirally symmetric solutions and chirally broken solutions, both satisfying all the required boundary conditions. We are also able to constrain either B or a to have a specified value, and we can extract all the required observables from the solutions. Only one thing remains now in order for us to determine what the observables are as functions of T , µ, B and a, which together parameterize the phase diagram of the theory. In principle, this is an easy task; one just has to evaluate the grand potential of every solution, and if there are multiple solutions with the same T , µ, B and a, choose the one with the smallest grand potential. In practice, however, this is more complicated. The grand potential diverges like A 4 near the boundary. In principle one could try to perform holographic renormalization like one does for the magnetic field, but this turns out to be very hard to do numerically. What works better numerically is to use the first law of thermodynamics and integrate along a family of solutions where one is able to continuously vary the input parameters. Of course, this brings with it an integration constant, which we choose such that the thermal gas geometry has grand potential equal to zero. This thermal gas geometry is obtained by taking λ h → ∞, which makes the horizon shrink to zero size. 11 This corresponds to a horizonless geometry with a 'good singularity', as discussed in [124].
By integrating the first law in this way, one can obtain the grand potential for every solution. Note that this is a non-trivial statement, as it requires the whole parameter space to be connected, and in particular it requires that the chirally broken solutions have a limit in which they approach a corresponding chirally symmetric solution. Both of these statements turn out to be correct [94].
After one has computed the grand potential throughout the parameter space in this way, one knows which geometry is the dominant one for every T , µ, B and a. To then label the resulting phase diagram appropriately, in this chapter we will use the following order parameters: • Chiral condensate qq : this distinguishes chirally symmetric from chirally broken solutions.
• Confinement: this is defined in this chapter as the quark-antiquark potential having a linear branch. Usually this is equivalent with the geometry being horizonless. In section 2.4, we will find that in the presence of anisotropy this is no longer always true. This has some interesting consequences, which will be explored there.
This concludes the introduction into the computations needed throughout the rest of this chapter. In the next sections, we will explore various limits of the general model described so far, starting with one where we have only a magnetic field, and hence set µ = 0 and a = 0.

Inverse magnetic catalysis due to a magnetic field
In the beginning of this chapter, the puzzle of inverse magnetic catalysis (IMC) was introduced. An interesting possible explanation was put forward in [128,129], namely that there are two competing effects at work, called 'valence' and 'sea' quark contributions. The valence quark contribution corresponds to theqq operators appearing in the path integral. This effect is found to increase the condensate as B increases. The sea quark contribution, instead, arises from the quark determinant. This contribution tends to decrease the condensate as B increases. The appearance of magnetic catalysis or its opposite are then explained as follows: in regions of the phase diagram where valence quark effects are dominant, one finds magnetic catalysis, whereas in regions where the sea quark effects dominate, one finds inverse magnetic catalysis. Note that the sea quark contribution is a backreaction effect. To see that this is true, consider the quenched approximation, or the approximation of infinitely massive quarks. In this approximation, the quark determinant becomes a constant, and hence the sea quark contribution vanishes. In other words, when one takes the quarks as non-dynamical probes, there is no sea quark contribution. In this section, which is based on [120], we study this problem in the holographic model introduced in section 2.1. The aim of this study is twofold: to reproduce IMC in holography, and to investigate whether in holography one can similarly isolate two competing effects in analogy to the findings in [128,129]. In earlier studies, holographic gauge theories in the presence of a magnetic field have been investigated either in the absence of flavors, or with N f N c [130][131][132][133][134][135][136][137], or with smeared backreacted flavor branes in the Veneziano limit, which leads to a different flavor global symmetry group [138]. In our model, we have the correct flavor symmetry group, as well as a fully backreacted flavor sector in the Veneziano limit. The model can be obtained from the 'master' model introduced in section 2.1 by setting a ⊥ = a = µ = 0. The baryon chemical potential can be set to zero by choosingñ = 0.

Varying the w potential
An important ingredient to reproduce IMC is the choice of potentials. As briefly mentioned in section 2.1, the potentials used in this chapter can be found in appendix A.1. These potentials are to a large extent the same as the ones used in [99], with one important difference. The coupling of the magnetic field to the flavor sector depends on the w-potential. In [99], the choice is made that w(λ) = κ(λ). This is a natural choice, as it is expected that κ and w have similar asymptotics in both the UV and the IR. Indeed, this assumption is consistent with the flavor vector current two point function

INVERSE MAGNETIC CATALYSIS DUE TO A MAGNETIC FIELD
and the asymptotics of the meson spectra [98,101]. Also, QGP conductivity and its diffusion constant [139], as well as thermal photon emission during the QGP phase of heavy ion collisions [140], are well described by assuming this assumption. A simple modification to w(λ) = κ(λ) is to insert a multiplicative constant b such that w(λ) = bκ(λ). Looking at the appearance of w in the action (2.4) however, such a modification would just lead to a trivial redefinition of B. Another simple modification which keeps the asymptotics of κ and w the same is the one we will use, namely By looking specifically at this quantity, it is possible to compare to lattice results, and indeed we find qualitative agreement [118]. At small temperatures, the condensate increases monotonically as a function of B, while at temperatures around the phase transitions the condensate first increases, then jumps, and then decreases. For even larger temperatures, the condensate decreases monotonically. This is in qualitative agreement with [118], with the notable exception that real QCD has no first order deconfinement transition, whereas the holographic model does. 14 In the right panel of figure 2.4, we show ∆Σ at vanishing temperature for different choices of c. It can be seen that the condensate behaves like where Dq q depends on c. It is clear that Dq q decreases as c increases. One final quantity we examine is the magnetic susceptibility at B = 0 as a function of T for different values of c, which is shown in figure 2.5. This is a useful quantity, as it can be used to determine the behavior of the deconfinement transition around B = 0. This can be seen as follows: Using the first law of thermodynamics, and using that the thermal gas has s = 0 where M B and s are now taken to be observables in the deconfined phase.
Together with the observation that Knowing that the w-potential does not influence the entropy density for vanishing B, we can determine the dependence on c of the behavior of the phase transition as a function of B purely from the dependence of the magnetic susceptibility on c. In this way, we can conclude that the deconfinement transition will decrease more sharply with B for small values of c.

Varying the number of flavors
Now that we have established that it is possible to obtain IMC in V-QCD by setting c < 1, we can further investigate the mechanism behind its appearance. To do this, we vary the number of flavors by tuning This has the effect of varying the amount of backreaction that the quark sector has on the gluon sector. Indeed, examining (2.4), one can see that x f multiplies the entire quark sector of the action, so that tuning x f 1 turns the quark sector into a probe, whereas x f ∼ 1 allows for significant backreaction. One of the first things one can investigate is how the phase structure changes for different values of diagram in figure 2.6 contains the following four phases: • A quark-gluon plasma phase at large temperatures. This phase exhibits chiral symmetry and is not confining.

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
• An intermediate temperature phase, exhibiting chiral symmetry breaking, but, like the QGP-like phase, not confining. In real QCD, where the phase transition is a crossover, both chiral symmetry breaking and confinement set in in a smooth way. In our holographic model, since the chiral transition and the deconfinement transition are second 15 and first order, respectively, these transitions can occur at distinct temperatures depending on the choice of potentials.
• A low temperature thermal gas phase, exhibiting confinement and chiral symmetry breaking.
• At x f = 0.1, another deconfined chirally symmetric phase exists at intermediate temperatures, separated from the one mentioned above by a first order phase transition. This is the same phase that was previously found at B = 0 in [94].
In addition to the phase structure, it can be seen that around B/Λ 2 ∼ 1, for all choices of x f we see inverse magnetic catalysis, chararcterized by a decrease in the chiral transition temperature. Also for all considered values of x f , we see that for still larger values of B, the chiral transition temperature increases, signaling magnetic catalysis. For small values of B/Λ 2 , the behavior of the chiral transition depends on x f , where we see that for x f = 1 there is inverse magnetic catalysis, whereas for x f = 2/3 there is only slight inverse magnetic catalysis, and for smaller values of x f the dip in transition temperature decreases to a point that it is no longer detectable. 16 Another interesting feature is the reappearance of the chirally broken deconfined plasma phase at large B/Λ 2 , which seems to appear for all considered x f , but for x f = 0.1 this could not be shown conclusively due to numerical inaccuracies. In figure 2.7, the magnetic susceptibility for c = 0.4 for the same values of x f considered in figure 2.6 is shown. Because in this case the entropy density at the deconfinement transition is not the same for all x f , as opposed to the situation for figure 2.5, we can not conclude anything about the behavior of the deconfinement transition from the magnetic susceptibility. Nevertheless, 15 Note that the chiral transition is first order whenever it concides with the deconfinement transition. It is only second order when it is a distinct phase transition. 16 The dip can never really go away, as χ/s > 0 for all deconfined solutions, which by (2.24) implies d 2 T d /dB 2 < 0. Since for x 2/3 the chiral transition coincides with the deconfinement transition, this implies that there sould always be slight inverse magnetic catalysis. interesting features can be discerned. Firstly, the susceptibility becomes larger as x f is increased. This reflects the observation that the coupling between the magnetic field and the gluon sector becomes stronger as x f is increased. Also, as one decreases x f , an inflection point can be seen to appear, which eventually forms a first order phase transition somewhere between x f = 0.1 and x f = 1/3. The last observable we will discuss is the chiral condensate. It turns out that by examining the x f -dependence of the chiral condensate, we can isolate two competing effects in analogy to [128,129]. The 'valence' quark effect will be identified with the direct coupling of the magnetic field to the tachyon field, which is dual to the chiral condensate. The 'sea' quark effect, in contrast, is identified with the indirect coupling of the magnetic field to the tachyon field through the metric. To investigate what effect these two contributions have on the chiral condensate, we need to find a way to change their relative contributions. Looking at the tachyon equation of motion (2.8), and different x f . Note that in both panels x f is denoted as x.

INVERSE MAGNETIC CATALYSIS DUE TO
condensate is larger than for B = 0, the difference ∆Σ is smaller for larger x f than for smaller x f . This can be seen more clearly in the right panel of figure 2.8, where ∆Σ is shown for different x f , T = 0, B = 3, with c = 0.4. Since at large B the 'valence' quark contribution is constant, we find that indeed the 'sea' quark effect, which we control through x f , tends to lower the condensate. This corroborates the findings of [128,129].
To conclude this section, we find that we can indeed reproduce inverse magnetic catalysis in a holographic model, provided that we take backreaction of the flavor sector onto the gluon sector into account. Furthermore, we find, similarly to the lattice QCD studies [128,129], that there are two competing effects that influence the chiral condensate. Of these, the backreaction effect associated with 'sea' quarks, tends to decrease the condensate, whereas the direct 'valence' quark, increases it. In the next section, we will discuss what happens to IMC in the presence of a finite baryon chemical

INVERSE MAGNETIC CATALYSIS IN THE PRESENCE OF A
NONZERO CHEMICAL POTENTIAL potential, which, due to the sign problem, is a hard problem to address with lattice QCD.

Inverse magnetic catalysis in the presence of a nonzero chemical potential
In this section, which is based on [121], we extend the analysis of the previous section to also include a finite baryon chemical potential. There are several reasons why such an extension is interesting. Of these, the most important one is that this region of the phase diagram is explored by the two systems mentioned at the beginning of this chapter. Indeed, low-energy heavy ion collisions such as those performed at RHIC occur in a region of both nonnegligible baryon chemical potential [141] and magnetic field [110], whereas neutron stars occupy the low temperature, high baryon density region of the phase diagram, and in particular magnetars have substantial magnetic fields. 17 Additionally, it is of theoretical interest what effect the chemical potential can have on inverse magnetic catalysis [142]. The latter is a question that cannot be addressed using lattice QCD, due to the sign problem [4]. Because of this, attempts to explore this question have been made both with effective models [112] and using holography [25,93,134,135,137,138,143,144]. Both the results of the previous section, as well as lattice studies [128,129], show that the effects of backreaction are important to capture the physics behind inverse magnetic catalysis. Therefore, we choose to study this problem in the same setup as in the previous section, with the addition of the baryon chemical potential. In terms of the 'master' model introduced in section 2.1, this amounts to setting a ⊥ = a = 0. As in the rest of this chapter, we will use the potentials from appendix A.1. Throughout this section, we will additionally fix c = 0.4, and we keep the number of flavors fixed to x f ≡ N f /N c = 1.

Phase diagram and thermodynamics
As before, we start the discussion of the results by examining the phase diagram, which is shown in figure 2.9. The phases featured in the phase diagram are the same as in section 2.2, with the exception that the extra chirally symmetric phase that appears for x f = 0.1 is not present, as we fix x f = 1 in this section. The phase structure shows an interesting dependence on the magnetic field. Firstly, the deconfined, chirally symmetric phase can be seen to move down for small chemical potential values, whereas it increases for larger µ. The behavior of the chiral transition will be examined in more detail in section 2.3.2. In contrast to the chiral transition, the deconfinement transition does not exhibit a large dependence on the magnetic field, with one notable exception. From section 2.2, we already know that in the absence of chemical potential around B/Λ 2 ≈ 1, the two phase transitions join. As the deconfined, chirally symmetric phase continues to exist for larger values such effects.

INVERSE MAGNETIC CATALYSIS IN THE PRESENCE OF A NONZERO CHEMICAL POTENTIAL
of B, there is a triple point. Interestingly, as this triple point forms around B/Λ 2 ≈ 1, the second order transition between both deconfined phases turns first order. After this, as can be seen in the inset in figure 2.9, another triple point forms, after which the first order transition becomes a transition between two deconfined, chirally broken phases, before ending in a critical point.
An interesting observable which encodes much of the thermodynamical information is the speed of sound. In figure 2.10, we show the speed of sound in the direction of the magnetic field. 18 The speed of sound c s is given by  Note that since c s depends on derivatives of thermodynamic variables, it is crucial to obtain accurate estimates for these quantities, as any numerical noise will cause the derivative to be unusable. For this reason, the fitting procedure described at the end of section 2.1.5 is essential to obtain the speed of sound, although even with the improved accuracy obtained from this method the speed of sound still has some numerical noise, which was removed artificially from the curves in figure 2.10. Rather surprisingly, the speed of sound squared is, for certain values of the chemical potential and magnetic field, larger than the conformal value of 1/3, which is surprising given [145,146], though it does not contradict their findings. Finally, note that at large T , outside the range of figure 2.10, the speed of sound approaches the conformal value from below, consistent with [145,146].

Chiral condensate and inverse magnetic catalysis
As the chiral transition is second order in a large part of the phase diagram, we can use the behavior of the chiral transition temperature as a function of chemical potential and magnetic field to study inverse magnetic catalysis.
In the left panel of figure 2.11, the chiral transition temperature is shown as a function of µ for different values of B. For small values of the chemical potential, the chiral transition temperature decreases as a function of B for the values of B considered in the left panel of figure 2.11. This signals inverse magnetic catalysis, and is in line with expectation from section 2.2. Interestingly, for larger values of the chemical potential, the exact opposite behavior is seen, signaling magnetic catalysis instead. Looking explicitly at the change in the condensate around B = 0, we can examine in more detail where we have (inverse) magnetic catalysis. We do this by looking at the sign of where a positive sign signals magnetic catalysis, a negative sign signals inverse magnetic catalysis, and zero signals no change. 19 Using this quantity, we obtain the right panel of figure 2.11. We can clearly see now that for temperatures well below the phase transitions, magnetic catalysis occurs, and that for a small region below the chiral transition inverse magnetic catalysis occurs, but only for small enough chemical potential.
For larger values of B, this behavior continues. This can be seen both from the left panel of figure 2.11 and from figure 2.12, where Σ is shown, as defined in (2.22). One can see that the chiral condensate always decreases with temperature, and for most values of T and µ, the chiral condensate increases with B. However, for certain values of T , µ and B, the condensate decreases. One can check that indeed that for small B this is consistent with the right panel of figure 2.11, but also for larger values B, one can find examples of inverse magnetic catalysis. Indeed, one can see that for example for T Λ = 0.12, µ/Λ = 0.05, the chiral condensate has a finite value for B/Λ 2 = 5, whereas the condensate vanishes for B/Λ 2 = 10.
In figure 2.13, we show the magnetization divided by the magnetic field as a function of temperature for different values of the chemical potential and magnetic field. In section 2.2, we saw that for the first order deconfinement transition, one can predict the slope of the transition temperature dT d /dB from the sign of the jump of the magnetization across the transition. This was studied in more detail in [143]. It turns out one can extend the argument to second order transitions as well. In this case one looks at the difference in entropy densities ∆s, which must be zero along the phase transition line. Then by using a Maxwell relation and using that ∂s/∂T is always larger for the higher temperature phase, one arrives at the following formula, where approaches zero from above. In other words, from the direction of the kinks of M (T ), we can infer whether the chiral transition moves up or down. Indeed, this formula correctly predicts the behavior of the chiral transition temperature, even though many of the kinks in figure 2.13 are too small to be visible.
In conclusion, we find that in our holographic model, at finite chemical potential inverse magnetic catalysis persists up to some value of the chemical potential, after which magnetic catalysis sets in. Also, the region of the phase diagram covered by deconfined, chirally broken matter, grows in size at finite B. In the next section, we will take a different look at inverse magnetic catalysis. Instead of applying an external magnetic field, we turn on a different source of anisotropy, and we examine how different this is from applying a magnetic field.

Inverse anisotropic catalysis due to an anisotropy
When we apply an external magnetic field to a quark-gluon plasma, this breaks rotational invariance of the system, as the magnetic field is a vector. It is then an interesting question to what extent phenomena associated to magnetic fields, like (inverse) magnetic catalysis, are in fact more general phenomena. In other words, are these phenomena caused by anisotropy being introduced into the system, or does that anisotropy need to specifically be a magnetic field in order for these phenomena to appear? In this section, which is based on [122], we will introduce anisotropy by means of the anisotropy parameter a introduced in section 2.1, which is dual to a space-dependent θ-term. Note that this way of introducing anisotropy is very different from a magnetic field, as a couples to the gluon sector directly, while the magnetic field couples indirectly through the quark sector.
In terms of the 'master' model from section 2.1, we set the magnetic field B to zero, while also settingñ = 0 to ensure a vanishing chemical potential. In this section, we will not distinguish between a and a ⊥ . This is possible because in this section we set the magnetic field to zero, which means that there is no other source of anisotropy, and we can choose the x 3 -axis to align with a. The choice of potentials is the same as in previous sections, namely those in appendix A.1. Most of these potentials have been used in previous sections, and their asymptotics have been motivated in section 1.5.1. This is not true for the Z potential though. We choose it to be of the same form as what was used in [144,147], where we have set the constant term [85] to 1, as an overall factor can always be absorbed into a. In the IR, we must have λ ∼ λ 4 to ensure a linear glueball spectrum for the 0 −+ glueballs [86,[147][148][149], and the choice is the simplest choice that satisfies these constraints. We subsequently choose α = 10 [150]. Before discussing physical observables, will next discuss the IR asymptotics of the geometry, as the geometry in the presence of a is drastically different from the geometry at vanishing a. Next, we will discuss the thermodynamics, and we will conclude by examining various observables, many of which turn out to have interesting properties as a consequence of the

INVERSE ANISOTROPIC CATALYSIS DUE TO AN ANISOTROPY
different geometry at non-zero a.

IR behavior
To examine the IR behavior of the geometry, we choose f = 1 so that there will not be a horizon capping the geometry off. We also set the tachyon τ to zero, as the tachyon decouples from the other variables in the IR [93,98]. Before continuing, it is also convenient to define λ = e φ , dr dA e A = q = −e p , and W = W + A. In terms of these quantities, the equations of motion become It is not difficult to see that these equations admit fixed point solutions provided that e 2p * V g (λ * ) = 9, 3a 2 Z(λ * ) = 2V g (λ * ), (2.25) with W a constant which we can set to zero, as it can be absorbed into a. First note that this solution describes an AdS 4 × R solution, as e A ∼ 1/r and e A+W is constant. Also, this solution turns out to be unstable. This can be seen from the second order dilaton EoM: The term in the second square brackets vanishes at the fixed point, but the term in the first square brackets does not. Because of this a small φ will cause φ to be non-zero as well, which implies an instability. The stability can be restored by requiring

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
in addition to (2.25), which can only be satisfied for specific a * . It turns out that this fixed point is realized by a chirally symmetric vacuum that exists for x f = 1/3, as we will see later. Note that in this case one has to replace V g in the analysis by To find the IR behavior satisfied by the chirally broken vacuum, we need to generalize the fixed point to a 'slow roll' solution around the fixed point. In other words, we assume that deviations from the fixed point solution are small and linear in A. To do this, we set all second derivatives to zero, which, from the argument above, requires that Z(φ) ∝ V g (φ) 3 , which, for our choice potentials, holds to good enough precision. Under these assumptions, one can obtain that where Furthermore, by substituting the asymptotic behavior of the potentials V g ∝ V IR e 4φ/3 √ φ, Z ∝ Z IR e 4φ , we can obtain analytically that This has been investigavted for more generic potentials in [151]. Indeed, the above asymptotics describe an approximately AdS 4 × R metric.
To check that these asymptotics indeed describe the IR behavior displayed by the model, we compare the solutions of (2.26-2.  well with the exact numerical solution of the equations of motion given in section 2.1.2. In the right panel of figure 2.14, the solution of (2.26-2.30) is compared to the analytical asymptotics obtained in (2.31-2.32). These can be seen to eventually agree, but only for extremely large values of −A.

Thermodynamics
We will continue the discussion by examining various thermodynamical properties of the system at different values of a. As was mentioned before, some of these properties are substantially altered by the modifications to the IR different values of x f show the same qualitative behavior, we will show only the result for x f = 0. One can see that generally there are three branches of black hole solutions: one small black hole branch (labeled I in figure 2.15), one big black hole branch (labeled II in figure 2.15), and an unstable branch (unlabeled). If the anisotropy parameter a vanishes, the result qualitatively agrees with earlier results [87,88]. 20 In this case, the small black hole branch in fact has infinitesimal area, and corresponds to the horizonless thermal gas solution. As one turns on a non-trivial a, however, the behavior is qualitatively different, as a black hole solution now always dominates in region I. This is due to the AdS 4 × R geometry in the IR.  potential has a linear branch. As the system is anisotropic, this need not happen in every direction. This leads to a phase with anisotropic confinement, which is confining in the direction parallel to the magnetic field, and deconfined perpendicular to it. Also note that with this definition of confinement, the transition between a confining geometry and a deconfined geometry is a crossover. We will discuss the confinement properties of the system in more detail in section 2.4.3, but it can already be seen using this definition that for x f 1/3, the geometries become deconfined, i.e. the quark-antiquark potential has no linear branch.
Another interesting observation is that the first order phase transitions The numerics were not stable enough to determine whether this also happens for the other values of x f considered. Also, note that for x f = 1/3, the chiral transition decreases to zero temperature around a/Λ ∼ 1, and at some larger a the chiral transition reappears. In between these two values, the vacuum is the chirally symmetric vacuum mentioned above.
The last feature that can be seen in figure 2.16 is the behavior of the chiral transition. Comparing to figure 2.6, we can see that in both cases the chiral transition temperature first decreases as a function of a and B, respectively. After this decrease, an increase follows in both cases. This tells us that a magnetic field and an anisotropy have similar effects on the resulting chiral transition temperature. We will examine the chiral condensate explicitly in section 2.4.3.
Lastly, in figure 2.17 we show the 'anisotropic susceptibility' χ a , which we define by

INVERSE ANISOTROPIC CATALYSIS DUE TO AN ANISOTROPY
with M a as defined in section 2.1.4. This definition agrees with the standard definition of a susceptibility in the limit where a → 0. In the same way as what we did in the previous two sections, we can derive for the first order transition that dT c da = − ∆M a ∆s .
Noting that the higher temperature phase at the phase transition always has the largest entropy, we can read off from the jump of χ a whether T c is increasing or decreasing with a. Indeed, this agrees between figures 2.16 and 2.17.

Observables
Next, we examine in some more detail the chiral condensate, after which we will look at some more consequences of the IR geometry. The chiral condensate is shown in figure 2.18. Just as in the previous sections, the chiral condensate decreases with temperature. Also, at x f = 0 the condensate always decreases with a, while for larger x f , the condensate first decreases as a function of a, and then increases again. The latter behavior can be studied in more detail by defining For all values shown, one can see first a decrease in ∆Σ, which is then followed by an increase. In particular, the decrease in ∆Σ is greatest for temperatures just below the chiral transition. This gives evidence that indeed an anisotropy as introduced in this section leads to similar physical consequences for the chiral condensate, leading to the claim that indeed the cause of inverse magnetic catalysis may be the presence of anisotropy induced by the magnetic field. Next, we will investigate the spectrum of helicity 2 glueballs. As was mentioned in section 2.1.4, we can infer from the Schrödinger potential (2.21) whether the spectrum is discrete or continuous. For the geometry without the anisotropy, the Schrödinger potential diverges towards the IR, creating a discrete spectrum. In the presence of a though, we have V S (r) ∼ 2/r 2 due to the AdS 4 ×R geometry, leading to a continuous spectrum. Such a dissociation of mesons due to anistropy has been studied earlier in holographic models in [152,153]. In the left panel of figure 2.20, the resulting spectral density is shown for x f = 0, normalized by ω 4 . It can be seen that for small a, the glueballs have small but finite widths, whereas they would be described by delta functions at a = 0. As one increases a, the widths grow larger, and the peaks melt. Still, the mass gap which is exact at a = 0 is still approximately present even at a = 1. Note that the peak appearing around ω = 0 arises because due to the AdS 4 × R geometry changing the power law behavior around the origin to G ∝ ω 3 . In the right panel of figure 2.20, the widths of can be seen to be getting wider with a power-law dependence on a.
Since the helicity 2 glueballs are not stable for a = 0, it is interesting to examine another quantity related to the stability of mesons, namely the quark-antiquark potential V . This observable quantifies the potential energy of a quark-antiquark pair at some distance, where we assume the quarks to be infinitely heavy. If the quark-antiquark potential grows infinitely with distance, it is impossible to pull the quarks apart, and in this case they are confined into mesons. 21 However, if the potential is bounded from above, at some distance the force between the quarks will vanish, and in this case it is possible to pull the quarks apart.
In [86,127], a simple criterion was found to describe whether there is a linearly growing branch of the quark-antiquark potential which, in figure 2.16, is used as the criterion to label confinement. In our case, this implies that if A S ≡ A+ 2 3 log λ has a local minimum, then the quark-antiquark in the direction perpendicular to the anistropy V ⊥ will have a linear branch. Similarly, if A S + W/2 has a local minimum, then the quark-antiquark potential in the direction of the anisotropy V will have a linear branch. In figure 2.21, we show the entire quark-antiquark potential for different values of a and with x f = 1. One can see that in both the parallel and perpendicular cases, the solutions have multiple branches, with in particular the perpendicular case showing multiple swallowtail structures. Of these branches, we are supposed to pick the smallest one, which is stable. A consequence of this is that even if A S resp. A S + W/2 have local minima, the linear branch which is then guaranteed to exist may not be stable. As with the previous quantities we examined, the reason behind this is the AdS 4 × R geometry, which causes the A S minimum to be only a local one. However, even if the linear branch is not stable, for small enough a there is still a large potential barrier for the string to decay to the stable configuration. This may indicate that the decay is not very fast, in line with the observations from the glueball spectrum that the glueballs are stable but long-lived.
The last observable we will examine is the entanglement entropy of regions A and B as defined in section 2.1.4. The result is shown in figure 2.22. First, note that there are multiple branches, the lowest of which is the stable one. Up to about a/Λ ∼ 0.1, the results are similar to isotropic results [154,155]. entangling region even when a is small. Specifically, for S E,⊥ , the curve of the connected surfaces always goes to (L = 0, S E,⊥ = 0, where the branch of disconnected solutions with S E,⊥ = 0 connects. However, for S E, the point (L = 0, S E, = 0) is never reached. Instead, the swallowtail structure present without anisotropy quickly gets smaller, and then vanishes. As a/Λ 1, the result becomes very different for the different entangling regions. S E, crosses zero for a smaller value of L, whereas S E,⊥ moves in the opposite direction. Such behavior has also been observed for a magnetic field [156]. Overall, the dependence of the entanglement entropy on the anisotropy is less pronounced than that of the quark-antiquark potential. A possible explanation for this is that the characteristic scale of the entanglement entropy is LΛ ∼ 1, causing the Ryu-Takayanagi surface to remain relatively close to the boundary, whereas the characteristic scale of the quark-antiquark potential is roughly LΛ ∼ 10. As the modification to the geometry due to the anisotropy is most pronounced in the IR, this implies that the quark-antiquark is more drastically modified even for small a, while the entanglement entropy needs rather large values of a for a modification to become pronounced.

Interplay between magnetic field and anisotropy
In the final section of this chapter, we will investigate the interplay of both sources of anisotropy introduced before. From the discussion in section 2.1, this gives us two options, namely to have a parallel to B, or to have it perpendicular to B. Recall that we were unable to put the two anisotropies at an arbitrary angle, but this is likely not an issue, as one expects to be able to already infer a lot of information from these two cases.
In terms of the 'master' model introduced in section 2.1, we setñ = 0 to ensure µ = 0, and we allowB,ã andã ⊥ to be non-zero. However,ã andã ⊥ can not both be non-zero at the same time. We also set x f = 1 throughout this section. The potentials will be kept the same as in the previous sections, i.e. the ones from appendix A.1. We do re-examine the c-parameter from the is always the stable branch. One can see that in this case, at B = 0.1, the quark-antiquark potential gets steeper perpendicular to the magnetic field, while parallel to it the potential gets slightly shallower. For larger values of B the quark-antiquark potential gets steeper in both directions. As the quark-antiquark potential has a linear branch in all the cases considered here, one can also compute the slope of the asymptotic linear behavior. This slope is called the 'string tension', and is shown in the right panel of figure 2.23, again for T = 0, a = 0 and x f = 1, but now as a function of B and for different values of c. One can see that as c decreases,

CHAPTER 2. (INVERSE) MAGNETIC CATALYSIS IN HOLOGRAPHIC QCD
the perpendicular string tension increases faster as a function of B, whereas the parallel string tension increases more slowly, and even decreases as a function of B for a part of the range of magnetic field shown for c = 0.4, and for all of the range shown for c = 0.25. In [157], it was found on the lattice that the perpendicular string tension increases as a function of B, while the parallel string tension decreases. This is most in line with the result for the holographic model for c = 0.25, and therefore in the rest of the section, we will show results for both c = 0. 25  transition is the only feature in the phase diagram. Therefore we do not show a separate phase diagram. We can see several qualitative features from the chiral transitions. Firstly, we observe that around a ∼ B, there is equality between T (a ⊥ ) and T (a ). Also, a seems to effectively decrease the value of B for a > B, which seems to be most pronounced for a .
In conclusion, in this chapter we constructed a bottom-up holographic theory based on V-QCD which includes two different sources of anisotropy, namely the magnetic field B and the axion a. We were able to show that in this theory at zero chemical potential and zero a, it is possible to obtain inverse magnetic catalysis. Furthermore, here are two contributions to the chiral condensate which have competing effects, in support of the analysis done on the lattice in [128,129], namely that inverse magnetic catalysis is

INTERPLAY BETWEEN MAGNETIC FIELD AND ANISOTROPY
caused by 'sea' quark effects, which in our model correspond to the backreaction of the gluon sector onto the flavor sector. Next, we generalized the discussion to finite chemical potential, thereby for the first time giving insight into the fate of inverse magnetic catalysis at finite µ. Also, we found evidence to support the claim first made in [123] that IMC is a more general effect caused by the anisotropy induced by the magnetic field. Lastly, we investigated the interplay between both sources of anisotropy, yielding interesting dependence on the angle between both sources.

Chapter 3 Holographic Baryons and Neutron Stars
In the study of neutron stars, some of the most basic questions are the following: • Given the mass, what is the radius of the star?
• What is the maximum mass?
• How large is the tidal deformability?
• What is the post-merger gravitational wave spectrum?
• Do neutron stars contain deconfined quark matter cores?
As was discussed in the introduction, these quantities depend on two things. The first is general relativity, which we understand well, and of which we believe it holds in the regime of neutron stars to a high degree of confidence. The second is the neutron star equation of state, something we cannot compute from first principles in a controlled way. This provides us with both a challenge and an opportunity. On the one hand we are obviously challenged to find what the true EoS is. On the other hand, neutron stars provide us with a novel laboratory to learn something about QCD. As more data will be collected in the future, some candidate equations of state will be ruled out, while others will remain compatible with observations. Subsequently one can start looking into the assumptions that went into these equations of state, to hopefully learn something about which of these assumptions are viable, and which are not. For this to work however, it is important that a wide range of such models become available.
In this chapter, we will construct such an equation of state from holography, which hence takes non-perturbative effects into account. In [103,158,159], holography has been used to model the deconfined phase which one expects at extreme densities, while these models have continued to rely on more traditional approaches for the nuclear matter phase, which is expected to make up the bulk of the neutron star. In this chapter, we in addition also model most of the nuclear matter phase using holography. The holographic model we will use for this purpose is, as in the previous chapter, V-QCD. For the choice of potentials, however, we make a different choice from the previous chapter. In [103], observables at vanishing baryon chemical potential were compared to lattice data, and the potentials were tuned to reach a good agreement. Recall also from the introduction that various qualitative features of QCD are matched by a judicious choice of asymptotic behavior of the potentials, such as that chosen in V-QCD. This means that by choosing a set of potentials from [103], we can do the computations in this chapter with a holographic model of QCD which matches real QCD to the largest extent possible at the time of this writing.
In particular, we will use potentials 7a from [103], which for convenience can be found in appendix A.2. In upcoming work which is not part of this thesis, the analysis from this chapter will be repeated for several other possible choices to investigate to what extent our prediction from this chapter is a generic prediction from V-QCD as opposed to one for this specific set of potentials.
By the procedure outlined in this chapter it turns out to indeed be possible to construct an equation of state, which, in some sense, one can then see as a sort of 'extrapolation' of lattice data to the neutron star regime. Moreover, it turns out that the EoS one obtains in this way is compatible with all presently known constraints, and is hence a viable candidate. One thing that should also be emphasised though, is that constructing this EoS can not yet be done from the V-QCD action plus the chosen potentials without making approximations. Wherever these approximations occur they will be pointed out, but the task of improving the approximations further will be left for future work.
In the first section of this chapter, which is based on [160], we will obtain the holographic equation of state directly from the V-QCD action. This EoS will, for reasons we will see below, be assumed to approximate the true EoS

BARYONS IN V-QCD
best in the regime of high density, and is expected to not work very well in the low density regime. In the second section, based on [161], we will introduce a matching procedure between the EoS from holography, and a nuclear matter equation of state coming from an effective Skyrme nucleon-nucleon interaction which is expected to work better in the low density regime [162,163]. We will then conclude this chapter by exploring the observational consequences of the resulting hybrid equation of state. For completeness, note that section 3 of [160] will not be covered, as it does not lead to reasonable neutron star physics.

Baryons in V-QCD
As was discussed in the introduction, baryons arise in holography as solitonic objects living in the bulk. In V-QCD, the fields responsible for this are the gauge fields (A µ L ) ij and (A µ R ) ij , where the i and j are flavor indices. These fields are dual to the globally conserved currents In the previous chapter, the trace part of A L + A R has already been used as the source for the conserved total baryon number, i.e. the baryon chemical potential. In that analysis however, the other components of A L/R were set to zero. In this section, we do need to take those components into account.
Before moving on to a brief summary of the steps required for performing this analysis, let us first fix a few conventions that will be used throughout this section. This is to keep notational consistency with [160]. Firstly, we will denote spacetime indices with capital latin characters, and we keep the following convention for the field strength tensors: Also, note that when refering to an expression which applies to both left and right handed gauge fields, we will use L/R to denote this. Furthermore, since the L and R labels are indeed just labels, we will not be consistent in putting them up or down. Instead they will be put in whichever place is visually convenient for a given formula. Lastly note that since we use capital latin characters for spacetime indices, the labels for left and right handed gauge fields will be denoted with brackets as (L) and (R) whenever there is a possibility for confusion.

CHAPTER 3. HOLOGRAPHIC BARYONS AND NEUTRON STARS
In the following subsections, we will first expand the non-Abelian part of the DBI action to first order in the gauge fields, after which we will examine the Chern-Simons terms. After that, we will introduce a homogeneous approximation, which reduces the degrees of freedom to a single field h. This field is discontinuous at some point in the bulk, where the discontinuity can be related to the baryon number density. To treat the discontinuity properly, one it is therefore hard to perform the computation in the grand canonical ensemble where we derive the action. For this reason we proceed to take a Legendre transform, and we finally derive an action involving h which, when integrated on-shell, gives one the free energy density (in the canonical ensemble). We will briefly discuss how this is done numerically, and finally the result is Legendre transformed back into the grand canonical ensemble, hence giving us the phase diagram and the equation of state. The last subsection will then be devoted to a discussion of these results.

Expansion of the DBI action
The starting point for this subsection is the full flavored DBI action (1.12) introduced in section 1.5.1, which we repeat here for convenience: where we have and the covariant derivative for T is given by Note that the non-Abelian part of this action is in principle ambiguous because the order of the gauge fields is not specified, and its full form is not completely known. The first few terms as a series in F are known precisely

BARYONS IN V-QCD
though [164][165][166][167], and hence we will here also use a series expansion in F . In fact, we will consider only the first non-Abelian correction on top of an Abelian background. Using this restriction, the ambiguities are absent, and the trace in 3.1 does not require a special prescription, and can be taken to be an ordinary trace. Obviously, using an expansion in the non-Abelian field strengths assumes implicitly that these fields are small. This is one of the points where the present analysis could potentially be improved in the future, as it would be interesting to take the next order in the expansion into account to see by how much the result changes. An additional simplifying assumption we make is to take the tachyon to be flavor independent just as what was done in the diagonalized DBI action used in chapter 2. The first step towards expanding in the non-Abelian part of the gauge field is to split up A L/R into an Abelian part, which contains the Φ field familiar from chapter 2, and a non-Abelian part. To do this, we employ a slight abuse of notation and substitute where we treat A L/R as small perturbations, while Φ is kept to all orders. This split is artificial, however, and to make the following well-defined we have to impose the following condition on the non-Abelian part of the gauge field:

After performing the substitution (3.3) into (3.2) and expanding for small
where we define A = A L − A R . 1 One can obtain a similar expression for A (R) M N . One next has to substitute this last expression into (3.1) to obtain an 1 In this chapter, the prime always denotes a derivative with respect to r.

CHAPTER 3. HOLOGRAPHIC BARYONS AND NEUTRON STARS
action quadratic in the non-Abelian gauge fields. To do this, it is convenient to define the 'effective metric': , which leads to the following two identities: This subsequently leads to the following expression, where we keep all terms quadratic in A L/R and F L/R : where we define the diagonal part of the effective metric as with the indices ordered as (t, x 1 , x 2 , x 3 , r).
Putting this all together, including the analogous term for A (R) , one obtains

BARYONS IN V-QCD
where the terms involving F (L/R) rt cancel due to (3.4), and where we split up the zeroth and first order in the expansion. With the last expression we have found the first part of the action for the non-Abelian gauge fields. In the next subsection, we will discuss the terms coming from the Chern-Simons section of V-QCD.

Obtaining baryon number from the Chern-Simons action
Before stating the Chern-Simons (CS) action, let us first discuss why this part of the V-QCD action is important for describing baryons. In the approximations that follow, we will not attempt to backreact the baryons onto the background geometry, which is taken to be a thermal gas. In particular, this means that Φ will be taken to be zero throughout the bulk in the computation of the background geometry. This would mean that the gauge fields don't respond to a change in chemical potential, which is something they should do since they carry the correct charge. The Chern-Simons action provide such a coupling, and indeed coupling Φ to A L/R is entirely appropriate, as we will see that the baryon number one can derive from the Chern-Simons action is a total derivative, and hence counts solitons, which is exactly how baryons arise.
We will now state the Chern-Simons action itself, which depends on a CPodd potential V a (λ, τ ) [150]. This dependence is highly non-trivial though, as the explicit form of the Chern-Simons action changes with different choices of V a . The reason for this is that as one changes V a in a non-trivial way, gauge invariance and other essential properties are usually spoiled. In this chapter, we choose a string motivated ansatz, namely where the constant b is introduced with the motivation that if one is to have regular IR solutions at a finite QCD θ-angle, the contributions from the DBI action need to dominate over those coming from the CS action, which is achieved provided that b > 1 [150]. We will now set b = 1 for notational simplicity. Later we will restore b by rescaling τ . 2 The CS action is now given by [92]:

CHAPTER 3. HOLOGRAPHIC BARYONS AND NEUTRON STARS
where Now, to extract the coupling between Φ and the non-Abelian terms in Ω 5 , we make the same substitution (3.3) as in the previous subsection, and collect all terms involving Φ. To do this, it turns out to be useful to first modify Ω 5 by a total derivative, since Ω 5 is only defined up to total derivatives. The modification made is the following: 12 (3.9) Now one can check that

BARYONS IN V-QCD
Using a variation of the action with respect to Φ just as in section 2.1.4, one obtains the number density where we have used the Φ equation of motion in the last equality. The baryon number is therefore now given by: where N b is the total baryon number. To conclude this discussion of the CS action, note that indeed, as was mentioned at the start of this subsection, that H (Φ) 4 is indeed exact: In the next subsection, we will take the next step towards including a holographic treatment of baryons in V-QCD, namely the introduction of a homogeneous approximation. This will reduce the degrees of freedom to one single field h.

The homogeneous approximation
In principle, we could now attempt to construct a soliton dual to a baryon in V-QCD by solving the action defined by (3.6) and (3.7), while using (3.10) as a boundary condition to impose a fixed baryon number. One could subsequently put many such solitons together to achieve a finite baryon number density. While this has been successfully done using approximate methods in for example the Witten-Sakai-Sugimoto (WSS) model [60][61][62], this is very challenging. The reason for this is that the required solutions have non-trivial spatial profiles, and hence to obtain them one needs to solve PDEs, and not only that, the requirement of imposing a fixed number of baryons turns the problem into a particularly hard one.
To make progress, we once again look to the WSS model for inspiration, as approximations have been developed there in a controlled setup [63][64][65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81]. The particular approach we will follow is a homogeneous approach [78,80], and as the name suggests, this reduces the problem to a spatially homogeneous one, which simplifies the problem to the solving of ODEs, for which we can use standard techniques. This is not the only possible way forward though. A similar approach to the one we will follow was developed in the WSS model in [81], and this approach can also be adapted for V-QCD using similar techniques to the ones that will be described below, but doing so is beyond the scope of this thesis.
The idea behind our homogeneous approach is to assume a high density of baryons. These will energetically prefer to be located at some coordinate distance r c in the bulk. 4 At r = r c , there will therefore be a very non-trivial and highly inhomogeneous configuration of gauge fields. Far enough away from this region, however, one can approximate the gauge fields as being homogeneous, just as when describing an electrically charged plate one does not have to worry that the plate is made up out of atoms if one is solely interested in the electric field far enough away.
The assumption we will now make is to divide the bulk into three regions: two regions away from r c where the gauge field is homogeneous, and a region around r c which we ignore, and which we take to be infinitesimally thin. We replace instead this region around r c by a discontinuity in h at r c . If one

BARYONS IN V-QCD
once again looks at the electric plate analogy, this means that if asked what for example the total energy stored in the electric field is, we choose to only include the field outside the plate, which we assume to be homogeneous, and we ignore the fields inside the plate, which we cannot easily describe. Of course this approach is not perfect, but it is expected that one does obtain a qualitatively correct description. In the same way, we will likely not obtain the correct on-shell action in this manner, but we can still hope that the dependence of the on-shell action on the baryon chemical potential is captured. Of course we have no way of checking whether this actually happens without doing the full calculation, but if this happens we can already learn valuable information about the equation of state this way. In particular, while the equation of state would only be reliable up to an overall factor, the speed of sound would be unaffected by such a deviation. Indeed, in section 3.2, we will make precisely this assumption, and take the speed of sound as input from the holographic model.
Concretely, the discussion from the previous chapter means that we set N f = 2 and hence take with σ i the Pauli matrices. Here it is important to note that (3.12) respects chiral symmetry and parity [168,169]. We take all integrals to mean Specifically, this means that there are no delta function contributions at r c , which would otherwise come from various derivatives appearing in equations of motion. Now that the prescription is known, we can start substituting (3.12) into (3.7) and (3.11). Doing this yields and dt ∧ H 14) respectively. Before collecting these two results into the complete action for h, note that, were it not for the discontinuity at r c , (3.14) would always evaluate to zero due to the asymptotic behavior of the tachyon, which diverges towards the IR, and h, which should vanish on the boundary because we do not want to introduce a source for the baryon charge. This is not surprising given that H is exact, but it does show that the baryon number in the solutions is indeed sourced at the discontinuity, as is appropriate in this approximation.
We conclude this subsection by substituting (3.14) into (3.10), and then combining this result with (3.6) and (3.13). This results in the action for the homogeneous gauge field h: This puts the action for the non-Abelian gauge fields in a relatively simple form, containing just a single field h. In the next subsection, we will continue the analysis by examining the discontinuity, as we will need to decide by how much h should be discontinuous.

The Legendre transformed action
In this subsection, we will derive a condition that needs to be satisfied at the discontinuity. In fact, this condition will completely fix the available freedom at the discontinuity, except for its location, which we'll get back to in the next subsection. We will also see that this condition relates the jump in h at the discontinuity to the baryon number density. This makes it hard to perform the calculation in the grand canonical ensemble, as we would need to extract the baryon number density from the solution, and simultaneously also impose it as a condition at the discontinuity. The solution to this will be to perform a Legendre transform, after which we eliminate Φ in favor of the baryon number density. This gives us a different action, namely one in the canonical ensemble, which we will then solve in the next subsection.
To start, we derive the baryon number density from (3.15) in the same

BARYONS IN V-QCD
way as before: where we recall that Ξ depends on Φ . Here we also introduced the abbreviation From the Φ equation of motion, we can obtain the following: Since we can ignore any delta function contributions at r c , ρ is continuous, and we can then derive where we can immediately read off as the baryon number density on the boundary, which is the physical baryon number density. Since ρ is continuous, this implies that h is discontinuous, and from the above equation one can derive that = 2N c π 2 e −b τ (rc) 2 (1 − 2b τ (r c ) 2 ) Disc h 3 (r c ), (3.19) where we define the discontinuity as Disc g(r) ≡ lim Next, we perform the Legendre transformation to the action S h , whose on-shell value will correspond to the free energy, as follows: where we notice that S CS = d 5 x Φρ .

CHAPTER 3. HOLOGRAPHIC BARYONS AND NEUTRON STARS
Putting these two equations together, one obtains To arrive at the final Legendre transformed action, there is one last non-trivial step. We need to eliminate Φ from the action. The way to do this is to solve (3.16) for Φ . This cannot be done analytically, but one can obtain Φ as a series expansion, where higher order terms correspond to higher powers of F L/R , which can by our assumptions be ignored. Performing this calculation leads to as in chapter 2. Finally, we can substitute (3.21) into (3.20) to obtain the Legendre transformed action: We have now arrived at the action that we will solve numerically. In the next section, some important points regarding how this should be done will be given.

Numerical method
With (3.22), we now have an action for h, which, by substituting (3.18), we can solve. In this subsection, we will discuss step by step how this is done, and how one can finally obtain the grand potential density of the baryonic

BARYONS IN V-QCD
phase as a function of the baryon chemical potential. Using this, one can construct the phase diagram and compute the equation of state. We start the process by computing a background metric and profiles for λ(r), τ (r) and Φ(r) in the same way as done in chapter 2. For this background we will use the thermal gas background, as this is expected to be the only type of solution that will be dominant in the phase diagram. 5 Recall that we will not backreact the non-Abelian gauge fields onto these background fields, which means that we will use the same pre-computed background for all the solutions for h that we will generate in the rest of this subsection.
The general idea for computing the solutions for h is that we want to minimize the baryon action (3.22) while keeping fixed in (3.18). In the continuous sections of the bulk, this means that we simply have to satisfy the equations of motion, but this also means that we have to minimize other free parameters other than . In particular, h has boundary conditions at the boundary and in the deep IR. Another remaining free parameter appears to be the location of the discontinuity. The location, however, is completely fixed by the two boundary conditions for h and (3.19).
At the boundary, the equations of motion imply that However a non-vanishing C 1 would imply a source for the non-Abelian gauge field, and hence we set C 1 = 0. This leaves us with a boundary condition at the boundary parameterized by C 2 . In the deep IR, it turns out that where C τ is a constant that can be expressed in terms of the potentials. This in principle gives us a second boundary condition h 0 . However, such solutions will be left out of the remainder of the discussion, as we have numerically checked that the value of h 0 which minimizes the action is consistent with zero. Hence, for simplicity, we will just set h(r) = 0 for r > r c . If desired, the discussion below can be easily extended to reinstate h 0 . With h 0 set to zero, (3.19) reduces to

CHAPTER 3. HOLOGRAPHIC BARYONS AND NEUTRON STARS
This means that in order to find r c , we just monitor (3.23) as we solve for h by shooting from the boundary, and when we detect a sign change, we iterate to find the r c where (3.23) is satisfied. Summarizing the previous discussion, we now need to shoot from the boundary, solving the equations of motion for h until we find r c by means of (3.23). All of this we have to do for different C 2 , where in the end we must minimize the on-shell action with respect to C 2 . For this, there is one hurdle remaining, namely that the action diverges towards the boundary, implying that we need holographic renormalization to evaluate it. There is a natural choice to use for the subtraction, and that is to subtract the on-shell action for h = 0. This corresponds to the on-shell action of a solution without baryons, namely the thermal gas. In this way, by subtracting the action with h set to zero from the on-shell action with non-zero h, we find exactly the difference in free energy density between the baryonic solutions and the thermal gas solution. 6 Using all of the steps above, one is able to compute the free energy F ( , C 2 ) for any desired and C 2 . In the final step towards minimizing the free energy, one performs this computation for a wide range of values for C 2 given a fixed . The minimal value found in this way then yields the true free energy F ( ), where we note that this minimal value is a stable minimum. Finally, performing a Legendre transform on this last result enables one to obtain the grand potential, or equivalently, the pressure p(µ), as desired.

Results
We are now ready to discuss the results. For this, we will fill in the value of M corresponding to potentials 7a from [103], and we will also set N c = 3. Note though that our value for N f is not entirely consistent between the backgrounds, which are used also for the deconfined phase, and the baryonic solutions. For the baryonic solutions we use N f = 2, while the backgrounds use N f = 3, which is the most appropriate given that the potentials were fitted to lattice data. This is another source of uncertainty, that ultimately ends up modifying the overall normalization of the pressure. Lastly, note that for all the results in this section, we use b = 10 for the parameter introduced in (3.8).

BARYONS IN V-QCD
As a first step in examining the consequences of including baryons into V-QCD, let us examine the resulting phase diagram, which is shown in using the methods described in this chapter. The other phases were computed using the methods described in chapter 2. The phase diagram itself contains three phases: • A thermal gas phase at small temperatures and chemical potentials. This phase features confinement and chiral symmetry breaking, similar to the vacuum of real QCD. This phase does not feature any temperature dependence though, as this is suppressed by powers of N c [102], and hence doesn't appear in the large N c limit.
• A baryonic phase at small temperatures and intermediate chemical potentials. This phase shares the features of the thermal gas phase that were mentioned above, and in addition also has condensed baryons.
• A deconfined, chirally symmetric phase at large temperatures and/or large chemical potentials. This phase is akin to the quark-gluon plasma phase found in real QCD. Due to the construction of the potentials, this phase also has an equation of state at vanishing chemical potential which is fitted to that of real QCD.
The resulting appears to be qualitatively reasonable within the limitations of the approach, such as the large N c limit. The baryonic phase appears in the right place in the phase diagram, and importantly, the baryonic phase eventually yields to a deconfined phase at some large chemical potential. This is in contrast to a similar approach in the WSS model [80]. Before moving on to the equation of state, it is a good check to look at the location of the discontinuity. This is shown in figure 3  given by with p and the pressure and energy density, respectively. While this is equivalent to figure 3.3 up to an integration constant, this is a convenient additional way of visualizing the equation of state. In the thermal gas phase, the pressure and energy density are both zero, so the speed of sound vanishes.
In the baryonic and deconfined phases, it can be seen that the speed of sound exceeds the conformal value in two places. This indicates a very stiff equation of state. As it turns out, a stiff equation of state is likely needed to pass astrophysical constraints [170,171].
Summarizing, in this section we have constructed baryonic solutions in V-QCD using a homogeneous approximation. In the next section, we will use the resulting equation of state to compute neutron star properties, and thereby investigate the potentially observable consequences from such a strongly coupled approach.

Hybrid neutron star equations of state
In the previous section, we constructed an equation of state in V-QCD using a homogeneous approach. This equation of state is likely to be most accurate in the regime of large chemical potential, as the homogeneous ansatz assumes a large baryon density. Also, we reasoned that because the homogeneous ansatz only really looks at the tails of the soltions, the on-shell action, and hence the pressure, is probably not correct by an overall factor. Motivated by this, in this section, based on [161], we will construct a hybrid equation of state, which is equal to a nuclear matter model at low densities, and has speed of sound equal to the one constructed in the previous section at high densities. We will then examine the potentially observable consequences from this hybrid equation of state, namely the mass-radius relation, the tidal deformability, and post-merger gravitational waveforms resulting from a binary neutron star merger.

Matching procedure
As briefly mentioned above, there are two ingredients that go into the hybrid equation of state, namely a low density nuclear matter model, and V-QCD.

HYBRID NEUTRON STAR EQUATIONS OF STATE
For the nuclear matter model, we will use the SLy equation of state [162,163]. This EoS is based on effective Skyrme interactions between nucleons. For the V-QCD part of the EoS, we develop a family of solutions, parameterized by varying the parameter b as defined in (3.8).
To perform the matching, we take the two equations of state, with b fixed. We then multiply the pressure belonging to the baryonic phase 7 in V-QCD with a constant c b to reflect that we take only the speed of sound as input from V-QCD, and then we demand that by choosing appropriate µ * and c b . These two conditions together guarantee that the phase transition between the SLy and V-QCD parts of the equation of state is second order. The final equation of state is then given by Of course, these two parts should both just be different descriptions of the same phase, so ideally there would be no phase transition at all. However, a second order phase transition is the best one can do if one is only allowed to change the location of the matching µ * and the normalization c b . One could probably make the transition even smoother, but that would require more ad hoc modifications of both EoS ingredients, for which we have no justification.
In this way, for every chosen value of b, one obtains a hybrid EoS. This resulting hybrid EoS is shown for several choices of b in figure 3.5. At low densities, one can see that the hybrid EoS agrees with constraints coming from effective low-density methods, indicated by the blue band. This is of course true by construction, as this part of the EoS is equal to the SLy EoS, which is itself constrained to satisfy this constraint. At high densities, the hybrid EoS can be seen to agree with constraints from perturbative QCD at  The blue and orange bands show the constraints coming from effective low-density methods and perturbative QCD, respectively. It can be seen that the hybrid construction satisfies both of these constraints.
extremely large densities, indicated by the orange band. This is also true by construction, because the asymptotics of the V-QCD potentials are constrained such that the deconfined phase of V-QCD satisfies this constraint. Another interesting feature that one can see is that as b increases, the matching chemical potential µ * increases as well. This effectively means that as b increases, the proportion of the neutron stars described by SLy increases, and equations of state with smaller b describe neutron stars with a relatively larger part described by holography.

Mass-radius relation and tidal deformability
Now that we have a family of equations of state with input from holography, one can start computing its observable consequences to compare how well the equations of state compare to presently known constraints. One of these observables is the mass-radius relation. Given an equation of state, the

HYBRID NEUTRON STAR EQUATIONS OF STATE
Tolman-Oppenheimer-Volkov equations can be solved to yield a relation between the mass and radius of a non-rotating neutron star. This mass-radius relation is shown for the hybrid equations of state in figure figure 3.5. Also shown as a dashed line is the observational lower bound on the maximum neutron star mass given by [32].
central density, the matter inside the neutron star undergoes a transition to deconfined matter, corresponding to the horizontal lines in figure 3.5. The kink that is visible in figure 3.6 corresponds to the central density at which this happens, and the solutions to the left of this kink have a deconfined quark matter core. The negative slope of the M -R relation for these neutron stars with quark matter cores implies that these solutions are unstable though, so the family of hybrid equations of state does not allow for stable quark matter cores. One observable that can be read off from the mass-radius relation is the maximum mass of neutron stars, which for our family of equations of state decreases with increasing b. This is an important observable, as the mass of a neutron star is a quantity that can be experimentally measured with a reasonably good accuracy. Indeed, [32] measured the mass of the CHAPTER 3. HOLOGRAPHIC BARYONS AND NEUTRON STARS pulsar J0348+0432 to be 2.01 ± 0.04 M . This immediately implies that any equation of state that does not permit a stable neutron star of at least about 2 M is incompatible with this observation. Of course, due to the uncertainty in the measurement of [32] it is not precisely clear where one should draw the line, but for our equations of state this means that the EoS corresponding to b = 10.65 is near the edge of exclusion at M max = 2.00 M . In 2019, [33] measured a potentially even heavier neutron star, J0740+6620, at M = 2.14 +0.10 −0.09 M . Including this measurement, the hybrid EoS for b = 10.6 is about one standard deviation lighter than the one measured in [33], so we will take b = 10.6 as our upper bound for b.
As it turns out, we can also obtain a lower bound for b. The reason for this is that using the equation of state one can also compute the tidal deformability Λ. 8 The tidal deformability is a dimensionless number, which for our family of hybrid equations of state decreases as a function of b. The tidal deformability, as the name suggests, describes how easy it is for another object to gravitationally deform the star, and in fact this quantity has an influence on the gravitational wave signal of the inspiral phase of binary neutron star mergers. From GW170817 as measured by LIGO/VIRGO, it was not possible to measure the magnitude of Λ, but it was possible to extract an upper bound of Λ 580 at 90% confidence level for a neutron star of mass 1.4 M [172]. For the hybrid equations of state, it turns out the b = 10.45 corresponds to Λ 680 for such a neutron star, which would be ruled out, but b = 10.5 has Λ 550.
In this way we obtain a family of equations of state with input from holography with 10.5 b 10.6 which are compatible with current (2020) observational data. It is possible that future observations can constrain these values further though. In particular, future results from NICER, which recently published its first results [29], could potentially constrain the radius, thereby putting another constraint on the equation of state. Also, with more data, it is possible that the maximum mass constraint becomes more stringent, and also tighter constraints on Λ could potentially rule out more, or even perhaps all, values of b.

Holographic neutron star mergers
We saw in the previous subsection that the hybrid EoS does not exhibit quark matter cores. 9 An interesting question is whether there are any circumstances in which the phase transition would lead to observable consequences. As there are no static solutions with deconfined matter, such circumstances would have to be fleeting moments in a dynamical process. One such process is a binary neutron star merger. The gravitational waves from such a merger can be detected in LIGO/VIRGO, and as these detectors are upgraded, and as new ones are built, detection capabilities continue to grow. The gravitational waveform contains information about the equation of state, both in the inspiral, as in the post-merger signal. In the inspiral part of the waveform, the equation of state manifests itself mainly through the tidal deformability Λ, which in fact we have already used to constrain the parameter b. Right after the merger, the single object is in a highly excited state, and vibrates with characteristic frequencies, which depend on the equation of state. These frequencies can in principle be detected in the corresponding part of the waveform, though as of 2020 the post-merger signal has not been detected.
In this subsection, we will perform simulations of equal mass neutron star mergers, with the aim of extracting observable consequences of the hybrid EoS from the resulting waveforms. We will restrict ourselves to equal mass mergers, and the stars will be initialized on quasi-circular orbits with a diameter of 45 km. Here the reason why the orbits are only quasi-circular instead of circular is that the system continuously loses energy due to gravitational wave emission, and the reason to assume vanishing eccentricity is that the same gravitational wave emission tends to have circularized the orbits well before the merger takes place [39,173]. To generate these initial conditions, we will use the publicly available LORENE pseudospectral code [174]. As explained in the introduction, the subsequent merger is described by general relativity plus general relativistic hydrodynamics [175]. To solve these equations, we use the Einstein Toolkit [176][177][178], where we use the highorder, high-resolution shock-capturing code WhiskyTHC [179][180][181][182][183], which solves the ideal general relativistic hydrodynamics equations in conservative form [184]. For the solution of the Einstein equations themselves, we use the fourth-order finite differencing McLachlan code [185,186], which solves the Einstein equations in the CCZ4 formulation [187]. Here we use a "1 + log" slicing condition, and a "Gamma-driver" shift condition [188,189].
In figure 3.7, one can see the result of such a simulation for a binary neutron star merger where the two stars have mass equal to 1.3 M . The One can see that the signal starts with an inspiral phase where the stars get increasingly close together as they radiate away energy in the form of gravitational waves. Subsequently, when the stars touch, a highly excited object forms, which continues to emit gravitational waves until it settles down to a stable state. In this simulation, 40 ms after the merger no horizon has yet formed, but it is unclear whether this simulation would eventually collapse to a black hole. To investigate this, one would have to continue the simulation for much longer past 40 ms, which is computationally very expensive. For this reason, it is not clear whether the final state of this merger simulation would be a neutron star or a black hole.
If the initial neutron stars are taken to be much heavier, an event horizon will form immediately after the stars touch, and in this case the gravitational wave signal dies down quickly after that. highly excited object forms just as with the lighter stars, but now at 7.6 ms, some of the matter in the middle of the resulting object reaches densities large enough to cross the phase transition into quark matter, which immediately causes a collapse to a black hole. This is an example of phase transition triggered collapse [190].
Together, these three cases describe the possible outcomes from a neutron star merger for the hybrid equations of state. 10 A quantity which contains a lot of information is the power spectral density (PSD) [191]: where the integrations are performed from −7 to 24 ms around the merger, where t = 0 is defined as the maximum of the gravitational wave amplitude. In figure 3.9, the PSD is shown for the three possible cases discussed above, each with the hybrid equation of state with b = 10.5. Generically, the postmerger signal contains three characteristic frequencies, labeled f 1 , f 2 and f 3 [38]. Here f 1 is a universal value which is determined by the compactness M/R [191]. This implies that from f 1 no information from the equation of state can be obtained that couldn't already be obtained from the mass-radius relation. The other frequencies f 2 and f 3 depend on the equation of state  in a more non-trivial way, and therefore these frequencies would in principle allow one to obtain additional information regarding the EoS. Note also that the M = 1.5 M case does not contain the characteristic frequencies in its PSD. The reason for this is that the characteristic frequencies are caused by oscillations of the dense matter, which are absent in this case, because the configuration promptly collapses to a black hole.
A final interesting thing one can do, is to compare mergers of equal mass neutron stars where one keeps the mass fixed at 1.3 M , while changing the equation of state. In figure 3.10, the resulting PSD of such a comparison is shown. Here we show the resulting PSD for the hybrid equations of state with b = 10.5 and b = 10.6, as well as the SLy equation of state for comparison, so that we can clearly see what the impact is of including input from holography into the equation of state. The resulting f 1 frequencies lie on the universal curve proposed in [191]. The f 2 peak shifts to significantly lower frequencies for the hybrid equation of state. Note that this shift is larger for the b = 10.5 case, which has a smaller matching density, and therefore contains more input from holography. Also, the f 3 peak has a smaller amplitude for the hybrid equations of state as compared to the pure SLy case.
In conclusion, the hybrid equation of state makes testable predictions  for various neutron star observables. The interesting aspect of this is that holography gives us for the first time an approach to the EoS which is inherently strongly coupled. In the future, one could refine the analysis from both sections of this chapter to perhaps improve the reliability of the result.

Chapter 4 Simulation of Heavy Ion Collisions with Trajectum
In the last chapter, we studied the low-temperature, large chemical potential region of the phase diagram using holography. This gave us insights on how an explicitly strongly coupled framework such as holography can yield observationally testable predictions. In this chapter, we will focus on another region in the phase diagram, namely that of high temperature and close to vanishing chemical potential, where a quark-gluon plasma (QGP) exists (see [192] and references therein). To be precise, we will take the chemical potential to be zero, but a small chemical potential could in principle be added to the analysis presented here in the future. QCD matter in this region of the phase diagram can be created in heavy ion collisions, such as the ones performed at RHIC and LHC. At the LHC, the collisions are performed at 2.76 TeV and 5.02 TeV, using lead-208 and xenon-129 nuclei, as well as protons. In the future it is considered to collide oxygen-16 and argon-40 as well, but as of 2020 no decision has been made.
Heavy ion collisions are able to provide us a unique insight into QCD, and many quantities of theoretical interest can in principle be observed in this way. Examples of these are the following questions: • What does the initial state of a heavy ion collision look like?
• How can the dynamics of the matter created at the collision be described before it can be described by hydrodynamics? After which time does hydrodynamics apply?

CHAPTER 4. SIMULATION OF HEAVY ION COLLISIONS WITH TRAJECTUM
• What are the values of QGP transport coefficients like the shear viscosity η and the bulk viscosity ζ? Can higher order transport coefficients and relaxation times like τ π and τ Π be measured?
One first thing to note is that the question as to the nature of the equation of state at vanishing chemical potential is not included. The reason for this is that the finite temperature equation of state is well known from lattice QCD computations [10,193]. While it is true that the dependence of the EoS on the chemical potential is less well known, 1 this is beyond the scope of this chapter, since we will not consider finite chemical potential. A second thing to note is that out of equilibrium properties, like how the initial state reaches local equilibrium, and even near-equilibrium properties like transport coefficients, are very difficult to compute on the lattice. In holography, however, many of these quantities can be computed, and often give surprising results. A famous example of this is the prediction that the ratio of the shear viscosity over the entropy density of the QGP is given by [25][26][27][28] Another result from holography is that the matter created in the collision can be described using hydrodynamics a very short time after the collision [23]. However, all of these results are computed under the assumption of infinite coupling and infinite number of colors, so it is an interesting question how well these results match the properties of an actual QGP. Measuring these quantities of interest, however, is not as easy as it may seem. Most observables in the final state of a heavy ion collision depend on more than just one of the quantities of interest, making it difficult to measure an observable and subsequently inferring something about the quantity of interest. What is possible though is to create a detailed model which simulates heavy ion collisions for different choices of the quantities of interest, which we will from now on call 'input parameters'. Such a model can then produce predictions for many different observables for each particular set of input parameters. By using such a model in a Bayesian analysis, it is then possible to fit the input parameters to real experimental data, allowing one to learn something about which input parameters fit the data best [196][197][198][199][200][201][202][203][204].

OVERALL DESIGN OF TRAJECTUM
In this chapter, we will construct such a model, called Trajectum, which also includes code to analyze the results. 2 While it will not perform any computations using holography, many of its inputs will be inspired by holographic results of which we would like to test the validity in experiments. The actual Bayesian analysis will not be performed in this chapter, but the results will be compared to section 5.3 of [24] as a check that Trajectum is able to reproduce known results. In the next section, we will explain the overall design of Trajectum, after which we will explain the various components that work together to create the simulation of the heavy ion collisions. Finally, we will explain the observables which can be extracted from the resulting collisions, and compare them to previous work.

Overall design of Trajectum
Of course, Trajectum is not the first code to perform the tasks described above. Many codes exist which perform parts of the computation necessary to simulate a collision [24,. These codes then need to made to work together though, where care must be taken to hand over the result of each step to the next step in the computation in the right format. 3 The aim of Trajectum is to collect all of these steps and to reimplement them into a framework which provides a standard interface between the various required components. In this way, one ends up with two executables, collide and analyze. The first of these computes some desired number of heavy ion collisions with input parameters according to the user's preference. At the time of writing, collide does not provide a hadronic afterburner, but the output of collide is in the correct format for use in UrQMD [228,229]. The output of UrQMD can subsequently be used in analyze, which is able to compute a range of observables, which will be discussed below.
The advantage of enforcing a standard interface between the components has one very clear advantage. For each component required for the simulation, several choices are available. As an example, the simulation requires initial conditions, for which several models exist. Of these, the following are • Monte Carlo Glauber [209][210][211][212], • Ohio State University [213,214], • T R ENTo [215], • Gubser flow [233][234][235].
When using collide, the user can choose which of these to use, and as they adhere to the standard interface, they are guaranteed to work together correctly with the other components. Similarly, all the other components of collide can be interchanged for various different options. The standard interface is implemented in C++ using polymorphism.
As was mentioned above, many of the components implemented in Trajectum are reimplementations of existing earlier work. In the sections below, we will discuss each of the various components in detail. In certain places, the original works have been extended or modified in order to address some of the holography-inspired questions posed above. Wherever the implementation of an algorithm in Trajectum differs significantly from the one in the corresponding earlier work, this will be clearly stated.

collide executable
As was mentioned in the introduction, a heavy ion collision is simulated in several stages. This consists of a pre-equilibrium stage which provides initial conditions for the hydrodynamical simulation of the QGP, followed by a hydrodynamical evolution. From the hydrodynamical evolution a freezeout surface is computed, which is defined as an isotherm called T fr . At this freeze-out surface, the code translates from a continuous fluid description to a discrete particle description. The particles must subsequently be written to a file, so that they can be further processed by the hadronic afterburner UrQMD, which simulates interactions between the particles produced up to the point in time where they can be taken to be non-interacting. The afterburner also simulates the decay of unstable particles.
These stages each translate into one or more components in collide. An event is simulated in collide as follows: The component responsible for the pre-equilibrium stage (called 'initial conditions') generates initial conditions 4.2. COLLIDE EXECUTABLE for the hydrodynamical evolution. This is then passed to the component for the hydrodynamical evolution. The component for the hydrodynamical evolution (called 'hydrodynamics model') then evolves the evolution in time step by step. For this it depends on two auxiliary components: • A component (called 'transport coefficients') containing the equation of state, as well as first and (if needed) second order transport coefficients. 4 • A component (called 'PDE solver') which implements an algorithm for solving partial differential equations.
By splitting off these components from the hydrodynamical evolution itself, it becomes easier to alter the simulation, as one can easily implement a new solver or set of transport coefficients. In principle, to do this it is not even necessary to understand how each of the other components works. One only needs to know how for example the transport coefficients should interface with the framework provided by collide, which then guarantees that the new set of transport coefficients correctly works together with the other components.
After each time step is computed, the hydrodynamics model hands the new state of the fluid over to the last component (called 'hadronizer'), which is responsible for computing the freeze-out surface and generating particles from the fluid. When the hadronizer determines that after the last computed time step there is no new addition to the freeze-out surface, the hadronizer causes the entire computation to terminate, and collide will move on to the next event. An important note here is that even though the wording in the last sentence may suggest that collide waits for each event to be completed before moving on to the next event, collide is actually multithreaded, and will compute 20 simulations simultaneously. 5 Summarizing, we have 5 different components which together simulate a collision: • Initial conditions, • Hydrodynamics models,
In the following subsections, each of these components will be explained. In particular, all currently available options that the user can choose from will be covered. Also, an important point to mention at this point is that we will assume boost invariance, i.e. we solve hydrodynamics in 2 + 1D. This is done by taking the following metric (called the Milne metric): which is related to the Minkowski metric by the coordinate transformation with τ the proper time, and η the pseudorapidity. With this metric, boosts correspond to shifts in η, and we can implement boost invariance by assuming that all of the variables constituting the fluid have no η-dependence. Also, for the rest of this section, we define the canonical order of variables (τ, x, y, η) = (0, 1, 2, 3).

COLLIDE EXECUTABLE
In [233,234], an analytical solution of viscous hydrodynamics was derived under the assumption of a conformal equation of state, constant η/s and ζ/s = 0. Gubser flow initial conditions implement initial conditions as a time slice of this analytical solution. This provides a good non-trivial check on the implementation of the numerical code to solve the hydrodynamical evolution [235]. The other three initial conditions are phenomenological models attempting to describe the initial state of a heavy ion collision as well as possible. They also include a model of the evolution from the proper time of collision (τ = 0 + ) to the moment that the hydrodynamical evolution is initiated (denoted τ fs ). In the remaining paragraphs of this subsection, we will discuss the various steps involved in all of these models. As the models are rather similar in setup, this will be done for all three models simultaneously, where we point out the differences where they occur. The first step in all three remaining models is to determine the positions of nucleons in both nuclei. For protons, this is simple. As there is only one nucleon inside the nucleus, the single nucleon just sits in the center of the nucleus. For the other nuclei impemented in collide, we assume that the nucleons are distributed according to a Saxon-Woods distribution [236,237]: 6 ρ(r, θ, φ) = ρ 0 1 + w(r/R(θ)) 2 where ρ is the probability density to find a nucleon at radial distance r, polar angle θ and azimuthal angle φ, ρ 0 is a normalization factor, w and a are parameters, and Here R, β 2 and β 4 are more parameters, and Y nm are spherical harmonics. The values for these parameters for the nuclei implemented in collide are given in  also have to satisfy the property that the distance between two nucleons has to be larger than a minimal distance d min . This is enforced as follows: One samples coordinates r, θ and φ according to (4.2). Then, starting from the nucleons with the smallest r, φ is resampled in such a way that the nucleon satisfies the minimal distance requirement with respect to all nucleons with smaller r than itself. For all initial conditions, the sampled nucleons are subsequently converted into Cartesian coordinates, where the z-coordinate which points along the beam is discarded. The remaining coordinates are given a random offset, reflecting the fact that the nuclei collide with a random impact parameter. The next step is to determine the locations of the constituents of each nucleon. Here there is a difference between the initial conditions models we are considering, as the Monte Carlo glauber model does not have nuclear substructure, the Ohio State University model has substructure with 3 constituents per nucleon, and the T R ENTo model only optionally has substructure, but with an arbitrary number of nucleons N const . For the initial conditions without substructure, one can simply skip this step, and consider a single constituent of each nucleon, located at the location of the corresponding nucleon which was computed in the previous step. The procedure below describes how constituents are sampled in the T R ENTo model for initial conditions. The procedure for Ohio State University initial conditions is parameterized differently, but is physically equivalent. In the T R ENTo model without substructure, nucleons are modelled as Gaussian blobs of density with width w. 7 In the version with substructure, there are N const constituents, which are modelled as Gaussians blobs of density with width v, sampled from a Gaussian distribution around the center of the nucleon, 4.2. COLLIDE EXECUTABLE with width r. As it is convenient that the width parameter w from the model without substructure corresponds roughly to the one with substructure, the model with substructure is parameterized in terms of w and v, where v < w. The width of the distribution of constituents r is then chosen as For each constituent in each nucleon, the position is sampled from a gaussian with width r around the center of the corresponding nucleon. However, the average of these sampled constituents does not necessarily correspond to the center of the corresponding nucleon. Therefore, the constituents are moved by such an amount so that the center does coincide with the center of the corresponding nucleon. Then, because of the choice of r from (4.3), the nucleons will on average have width w. 8 Now that we have the positions of all the nucleons and (depending on the choice of initial conditions) their constituents, we have to determine which nucleons become 'wounded'. Wounded nucleons are precisely the nucleons that participate in the collision. To determine which nucleons participate, we examine all pairs of nucleons, where we take the first nucleon (nucleon A) from the first nucleus (which we will call nucleus A), and the second nucleon (nucleon B)) from the second nucleus (which we will call nucleus B). For each such pair, we then compute the overlap function where x A,i are the x-coordinates of the constituents of nucleon A, y A,i the y-coordinates of the constituents of nucleon A, and similarly for x B,j and y B,j . For the initial conditions without substructure, this overlap function is to be interpreted to have a single constituent for each nucleon, where we also have to replace v by w. Now that we have the overlap function, the probability that these two nucleons participate in the collision is

CHAPTER 4. SIMULATION OF HEAVY ION COLLISIONS WITH TRAJECTUM
where σ gg is a constant. In accordance with this probability, the nucleons are both marked as being 'wounded' with probability P coll . 9 If the two nucleons are marked as wounded, the corresponding nucleon pair is also marked as being a 'binary collision pair'. We have here introduced the parameter σ gg . In principle, we could leave this as a parameter to be determined by the user. Instead of this, however, we choose σ gg precisely such that the average cross-section of proton-proton collisions is the same as the parameter σ N N , which is available to the user to specify. Knowing which nucleons participate in the collision, we can start constructing a 'density' to initialize the plasma. Note that this density is just a function T (x, y). Depending on the model, it can be interpreted as an entropy density, energy density, or specific component of the stress-energy tensor. How exactly the different models interpret T will be discussed towards the end of this section. To compute T , let us first define the 'thickness functions' T A (x, y), T B (x, y) and T AB (x, y): where x ij , and y ij denote the x and y-coordinates of constituent j of nucleon i, and the outer sum goes over all wounded nucleons in nucleus A. What the weight γ ij is, depends on the initial conditions. For the Monte Carlo Glauber model, it is equal to 1, whereas for both the Ohio State University model and the T R ENTo model it is a random number, which for each constituent is drawn from a Gamma distribution with mean equal to 1, and standard deviation equal to a parameter called σ fluct . An analogous expression holds for T B (x, y) in terms of the wounded nucleons in nucleus B. The last function, T AB (x, y), is not determined in terms of the wounded nucleons, but in terms of the binary collision pairs. As T AB will only be used by the Monte Carlo Glauber, and since this model does not have constituents, the expression below is given in terms of the nucleon positions themselves: where (x i , y i ) are the coordinates of the nucleon from nucleus A participating in the binary collision, and (x j , y j ) are the coordinates of the nucleon from 4.2. COLLIDE EXECUTABLE nucleus B. In words, the wounded nucleons deposit density into the thickness functions T A and T B at the location of the constituents, and each binary collision pair deposits density into the thickness function T AB . How the three thickness functions are subsequently combined into one density T again depends on the specific model. In the Monte Carlo Glauber model, we have with α a parameter. In the Ohio State University model, the T AB -contribution is ignored, and we have In the T R ENTo model, we have with p a parameter. Note that in the limit p = 0, this reduces to the geometric mean T (x, y) = T A (x, y)T B (x, y).
At this point, for all models considered, we have a function T (x, y) with the dimension of inverse area. We now have several options for how to interpret this function. In both the Monte Carlo Glauber and Ohio State University models, T is interpreted as an entropy density as follows: where n is a dimensionless number called the norm, and τ fs is the proper time at which the hydrodynamical evolution is initialized. Using the equation of state, this can then be used to determine energy density and pressure. This does not completely fix the initial condition for the hydrodynamical evolution though, as the stress-energy tensor has 7 independent components, and fixing the energy density only fixes one. 10 In collide, there are two

CHAPTER 4. SIMULATION OF HEAVY ION COLLISIONS WITH TRAJECTUM
options to determine the other coefficients. The first sets the velocities to u µ = (1, 0, 0, 0), sets bulk pressure Π = 0, and subsequently computes the shear tensor as 11 π µν = 2ησ µν , (4.4) where the code gets the shear viscosity η from the transport coefficients model chosen by the user, and σ µν is defined as in equation 1.5. The second option sets π µν = 0, Π = 0, and [23]: with e the energy density, x 1 = x and x 2 = y.
In the T R ENTo model, a different approach is taken. Here, the stressenergy tensor is initialized as [24]: where n is now a parameter, also called the norm, with dimension of inverse length. Forp µpν , we havê where the rows and columns of the matrix correspond to µ, ν = 0, 1, 2. At this point, we deviate from the original implementation of T R ENTo. In the original implementation, the parameter we call v fs is set to 1, in which case (4.5) corresponds some density of massless particles determined by T free streaming from time τ = 0 + to τ fs . Indeed, this reflects in the name of the initialization time, where the fs stands for free streaming. The parameter v fs ∈ [0, 1], which we call the free streaming velocity, generalizes this ansatz by allowing free streaming to occur with a speed lower than that of the speed of light.
Another deviation from the original T R ENTo is in the initialization of the shear tensor and bulk pressure. In [23], it was observed that using holography to model the initial stage leads to a smooth transition to hydrodynamics. The energy density and velocity profile is identical to first order in proper 4.2. COLLIDE EXECUTABLE time to that obtained by free streaming, but the shear tensor obtained from holography is given by (4.4). 12 In order so that a future Bayesian analysis can potentially be used to decide whether free streaming or holography fits experimental data best, we introduce a parameter α ∈ [0, 1] with which we can interpolate between free streaming and holography.
The procedure to obtain the stress tensor then becomes the following. We first construct the stress tensor as given by free streaming (4.5), and decompose the stress tensor in terms of e, u µ , π µν fs and Π fs . We then compute π µν holographic and Π holographic as We subsequently compute and reconstruct T µν using e, u µ , π µν and Π. In this way, by choosing α, we can interpolate between weakly coupled initial conditions (α = 0) and strongly coupled initial conditions (α = 1). A final feature that is new in collide, which is available for all initial conditions models discussed above, is the ability to optionally bias the distribution of events generated by the initial conditions. To illustrate what this is and why one would want this, consider proton-lead collisions. Typically, for such collisions, we are interested in collisions which generate a lot of particles in the final state. These events correlate fairly well with the initial entropy deposited in the plasma by the initial conditions. 13 However, events which generate a lot of particles in the final state are rare. This leads to the inconvenient situation that to get a good amount of data to be able to achieve small statistical uncertainties for the observables one is interested in, one has to generate even more events, the vast majority of which will not be useful. The bias that we introduce now is based on the following idea. Of the entire collision simulation, the initial conditions are by far computationally the cheapest to compute. Since the initial entropy correlates well with the number of particles in the final state, we can compute the initial conditions TRAJECTUM for the hydrodynamical evolution, and then accept the event for further computation with a probability P (s), which depends on its initial entropy. By choosing the acceptance function P appropriately, we can selectively compute events which will have a large amount of particles in the final state, thereby improving statistics for observables relating to such events without needing to spend much more in computation time. Of course, when computing these observables, we need to be careful that our bias in the events that we generate does not translate into a bias in the final computed observable. This can easily be achieved by making sure that each event analyzed by analyze is weighted with weight 1/P (s). Also, care has to be taken when computing centrality classes, as these weights also have to be taken into account when determining them.

Hydrodynamics models
The next component we will describe is the hydrodynamics model itself. In collide, two models are available: • First order hydrodynamics, which solves the first order Israel-Stewart equations with only first order transport coefficients [17,[238][239][240], • Second order hydrodynamics, which also includes some second order coefficients.
Of course, one can just obtain the first order hydrodynamics equations by setting all the second order coefficients to zero. However, there is still a good reason to have two separate classes, namely that setting the second order coefficients to zero when running the code is a lot slower than explicitly leaving those coefficients out of the code to be evaluated, and solving the simplified equations. Since the equations which determine the hydrodynamical evolution are used by the solver every time step at every grid point, their speed of execution greatly impacts the overall execution time of the whole program. However, to explain how these two models work, we will just explain how the second order hydrodynamics model works, as the workings of the first order model can be easily determined by setting the second order coefficients to zero.
In hydrodynamics, one of the first things we have to decide when solving the equations numerically is which variables to choose. To illustrate this, one could imagine solving for each of the 7 components of the stress-energy 4.2. COLLIDE EXECUTABLE tensor. However, one could equivalently solve for the energy density e, two velocities v x and v y , the bulk pressure Π and the components π xx , π xy , π yy . 14 A convenient choice of variables turns out to be [240] T τ τ ,T τ x ,T τ y , Π,π µν , (4.6) where we defineT µν = τ T µν , andπ ηη = τ 2 π ηη , withπ µν = π µν for the other components. Just by counting, one can see that this number of variables is larger than the number of independent variables, as we are evolving every non-zero component of the shear tensor. It turns out however, that for the first step of solving the hydrodynamical equations, which will be discussed below, it is easier if these redundant components are evolved as well, so that they can be used in that step. Another question one might ask is why certain variables are defined by incorporating a factor τ into their definition. The reason for this is twofold. On the one hand, it reduces the number of Christoffel symbols entering the equations, and on the other hand it removes some general behavior of the variables. To see what is meant by this, consider T τ τ . By energy-momentum conservation, T τ τ will behave roughly like 1/τ . By instead evolvingT τ τ , this dependence is removed, making it easier to obtain a good accuracy in the numerics. Solving the hydrodynamical equations themselves involves two steps. The second step depends for its computation on quantities like the velocity and the energy density, which are quantities that are not listed in (4.6). This immediately explains the need for the first step, which is to reconstruct these quantities, which we will call 'auxiliary variables' from now on, from the variables in (4.6), which we will call the 'primary variables'. What the second step does depends on which solver is used. The finite difference solver requires for its algorithm a function which takes in all current values of the primary variables, as well as the auxiliary variables mentioned above, and outputs the proper time derivatives of the primary variables. The MUSCL solver instead requires the equations to be put in the form where U k is an array of all variables in (4.6), v x = u x /u τ and v y = u y /u τ are the fluid velocities in the x and y directions, and where we call S the source

COLLIDE EXECUTABLE
f (e) = 0 and the computation of u µ requires the pressure to be known as a function of the energy density, which is a task which is performed by the transport coefficients component. Next, we describe the computation of the source terms. ForT τ µ , we have [240] For the bulk pressure Π, we have used not all the second order expressions, but just the ones also kept by [24] (The full second order expressions can be found in [241][242][243]): whereσ ηη = τ 2 σ ηη , withσ µν = σ µν for the other components similar to the definition ofπ µν , and Lastly, for the shear tensor, we have Here the second order coefficients, which are set to zero for the first order model, are the coefficients δ ΠΠ , λ Ππ , δ ππ , φ 7 , τ ππ and λ πΠ .

Transport coefficients
The task of the transport coefficients model is twofold. It needs to compute the pressure by means of an equation of state, and it should compute the transport coefficients appearing at the end of the last subsection. In collide, there are three available transport coefficient models:

CHAPTER 4. SIMULATION OF HEAVY ION COLLISIONS WITH TRAJECTUM
• Ideal gas equation of state with 'constant' transport coefficients.
• Lattice QCD equation of state with 'constant' transport coefficients.
• Lattice QCD equation of state with temperature dependence in some of the transport coefficients.
Here 'constant' means that the transport coefficients are described by a single parameter. For example, for the shear viscosity η, 'constant' transport coefficients do not imply that η(T ) itself is constant as a function of temperature, but rather that η(T )/s(T ) is, where s(T ) is the entropy density specified by the equation of state. We will now describe each of the available models in some detail, starting with the ideal gas equations of state with 'constant' transport coefficients.
For the ideal gas equation of state, we assume an equation of state of the form with P the pressure, and α a constant which can be specified by the user. Additionally, the shear and bulk viscosities can be specified by the user by means the following constant combinations: The shear and bulk relaxation times τ π and τ Π , respectively, are specified as the following constant combinations: The second order coefficients are given by the following constant combinations: with P the pressure. Note that if one specifies ζ/s = 0, τ π sT /η small enough, and all other coefficients zero, one can use this model for transport coefficients in combination with the Gubser flow initial conditions to obtain the analytical 4.2. COLLIDE EXECUTABLE solution [233][234][235], with which one can check the accuracy of the numerical solution.
For both models which use a lattice QCD equations of state, we actually use a hybrid of a hadron resonance gas (HRG) for temperatures below 165 MeV, an analytical fit to a numerically constructed lattice QCD equation of state for temperatures above 200 MeV, and a polynomial interpolation of the trace anomaly in the intermediate temperature regime [24,244]. The polynomial interpolation is taken such that at the matching points, 165 MeV and 200 MeV, the trace anomaly, as well as its first 4 derivatives, is continuous. The hadron resonance gas equation of state can be computed from with the + for fermions, and the − for bosons. Also, E i (p) = m 2 i + p 2 , where m i and g i are the mass 16 and number of degrees of freedom of species i, respectively. The species used in the sum are precisely the particle content of UrQMD, which is necessary for consistency between the different stages of both collide and UrQMD itself.
The lattice QCD part of the equation of state is described by the following parameterization [10]: and the fit coefficients c t = 3.8706, a n = −8.7704, b n = 3.9200, d n = 0.3419,

PDE solvers
The PDE solver performs the task of solving the hydrodynamics equations.
To do this it interfaces with the hydrodynamics model, which provides the necessary functions, as discussed above. The PDE solver then updates the state of the fluid from proper time τ to τ + ∆τ , where the size of ∆τ can be chosen by the user. In collide, there are two main types of PDE solvers: • Finite difference solver, • MUSCL solver.
Both of these solvers have slight variations implemented, which allow the user even more flexibility. Each of the solvers has their advantages and disadvantages, so one should pick carefully when using collide on a problem. The finite difference solver uses a very simple algorithm, as we will see below. This has the advantage that, due to its low complexity, it is faster by about a factor 2 as compared to the MUSCL solver. The disadvantage of the finite difference solver is that it is not guaranteed to be stable. Under certain circumstances, like the presence of shocks or when using transport coefficients with extremely small viscosities, numerical instabilities may appear, which grow exponentially. The MUSCL solver, instead, is guaranteed to be stable. Indeed, the algorithm that it uses is stable by construction, as it was designed with stability in mind. For this reason, one can use MUSCL even for an ideal fluid, with all viscosities set to zero. The disadvantage that accompanies this stability, however, is increased execution time. Depending on the type of problem, one can choose which of these benifits outweigh the associated costs. In the following few paragraphs, we will discuss both of these solvers, pointing out the available variations along the way. Now let us first discuss the finite difference solver. One can write the hydrodynamics equations for the 'primary variables' defined in section 4.2.2 as follows: with f a function, and where U k ij is the collection of primary variables at grid site (i, j). Here the first entry corresponds to the x-coordinate, and the second corresponds to the y-coordinate. One can then perform the following

COLLIDE EXECUTABLE
where a x is equal to v x , and where F x (U k ) = v x U k is the flux of U k at the cell's interface. At this point, an optional shortcut is available in collide. To properly evaluate the flux for a specific set of primary variables, one needs to find v x and v y , which implies that one has to find the auxiliary variables, which is numerically expensive. The shortcut taken is that we compute v + x and v − x in an analogous way to the computation of U k,+ and U k,− , thereby limiting the amount of times one has to solve for auxiliary variables. 19 For U k,+ and U k,− , we have: with the quantityŨ k being roughly equal to the second order accurate derivative with respect to x. If one would take this rough equality to be exact, the algorithm would become identical to the finite difference method described above. This would then however also inherit the stability issues of that method. It turns out that the crucial element to achieve stability is the addition of a flux limiter. In the MUSCL implementation in collide, we use the minmod flux limiter: where the minmod flux limiter function is defined as minmod(x, y, z) = minmod(x, minmod(y, z)), minmod(x, y) = [sign(x) + sign(y)] · min(|x|, |y|) 2 .
Finally, note that in collide, one has the option to either integrate the time derivatives using the forward Euler algorithm, which is first order accurate, or with the midpoint method, which is second order accurate. In the results shown in section 4.3, we used the KT solver with the shortcut for the solving of auxiliary variables.

Hadronizers
After each update that the PDE solver makes to the fluid, the hydrodynamics model hands the new fluid state to the hadronizer. The hadronizer is TRAJECTUM then tasked with generating particles, which are output in a format compatible with UrQMD. In collide, at the time of writing, only one hadronizer is available, namely Cooper-Frye [223]. The hadronization procedure consists of several steps. First, the hadronizer computes a freeze-out surface. Subsequently, particles are produced at the freeze-out surface according to a modified thermal distribution, where the modifications encode the presence of shear stress and bulk pressure into the final state particles. In the paragraphs below, we will explain each of these steps in more detail. The first step to turn the fluid into particles is to compute the freezeout surface. This is an isotherm of a specific temperature T fr that the user can specify. At the freeze-out surface, a number of particles are produced according to a Poisson distribution, which has the property that the sum of Poisson distributed processes again follows a Poisson distribution. As a consequence of this, we can subdivide the freeze-out surface into triangles, and sample particles from each triangle individually. After this, we can even discard the freeze-out surface from which particles have already been sampled, as it is no longer necessary. This is exactly how collide tackles the problem. When the hydrodynamics model computes the state of the fluid at proper time τ + ∆τ , it gives the state of the fluid to the Cooper-Frye hadronizer, which still has in its memory the state of the fluid at proper time τ . Using the Cornelius algorithm, it then computes a triangulation of the freeze-out surface [222]. 20 For each triangle, it also computes a surface normal ∆σ µ [24]. The next step is to generate particles for each triangle. This is done using the Cooper-Frye formula [223]: where g i is the number of degrees of freedom for species i. The idea behind this formula is that the hadronization procedure should be not a change in physical process, but rather a change in description of that process. The particle species used in the formula are precisely the ones that are used in the HRG equation of state which describes the fluid, and because they are sampled from a thermal distribution, they indeed describe precisely the 4.3. ANALYZE EXECUTABLE HRG. 21 This ensures that, in the absence of shear stress or bulk pressure, the stress energy tensor of the fluid will be on average the same as that of the sampled particles. 22 In this sense, the sampled particles form a different description of the fluid, and when one then only keeps particles which move out of the fluid, one thereby obtains a consistent transition between the fluid and the HRG. Everything described so far assumes that π µν = 0 and Π = 0. If this is not the case, the sampling of the particles needs to be adjusted so that the stress-energy tensor of the fluid still matches the average stress-energy tensor of the sampled particles. The way this is done in collide is by rescaling the momentum of each sampled particle in the fluid rest frame [24]: where (λ shear ) ij ∝ π ij with a proportionality constant that can be computed by a hadron resonance gas computation, and λ bulk depends in a non-trivial way on Π. 23 For bulk corrections, also the particle density obtained from the Cooper-Frye formula (4.12) needs to be modified, as with the rescaling of the particle momenta described above the energy density is also modified. The way to correct this is to modify the number density [24]. We now have produced particles in a way such that the stress-energy tensor is on average continuous across the freeze-out surface. There is one final step, however, namely decaying the f 0 (500) particle. This particle is not treated by UrQMD, but given its light mass is produced in large quantities during an event [24,247]. In collide, the f 0 (500) can be optionally taken into account when compiling Trajectum. If the user includes the f 0 (500), it is immediately decayed into pions, which is justified given its short lifetime.

analyze executable
The analyze executable is tasked with computing the observables according to their proper definitions from the result of the hadronic afterburner. Additionally, it obtains from the output of collide each event's weight. The executable does this in two passes over the data file. The first pass reads every event, and computes quantities that will allow us to group the particles into centrality and N off trk bins during the second pass. During the second pass, each event is read again, and based on its N off trk value is placed in a bin together with similar events, after which each event is analyzed by the observables that the user requested. In the next few paragraphs, each of these two passes will be described in some more detail, after which a selection of observables will be shown.
The first pass over the data has the task to determine the centrality classes. This is done by counting N off trk for each event, which is a measure of the number of particles produced in the event. The definition of N off trk differs from experiment to experiment, and is usually defined in such a way that only particles are counted which are charged, and have kinematics such that the detector is able to detect them efficiently. Of course, in the output of collide and UrQMD, we are able to detect any particle whatsoever with 100% efficiency, but for comparison with the experiments it is still important to make the same cuts, so that we are comparing the exact same observable. This will be a recurring theme throughout this section. As briefly mentioned already, N off trk is defined by the number of charged particles in an event, and the kinematic cuts usually require that for a particle to be counted, the transverse momentum p T needs to be between a certain lower bound and a certain upper bound, and that the pseudorapidity η needs to have an absolute value below some upper bound. For example, the cuts used in this section are 0.4 GeV ≤ p T ≤ 10 GeV, |η| ≤ 2.4.
After N off trk has been counted for each event, the events are sorted from the largest N off trk to the smallest. This ordering then determines the centrality bins, and for each centrality bin, analyze then saves between which values N off trk should be such that the event is within that particular bin. Of course, when events are weighted according to their initial entropy, it is important to take these weights into account when computing the centrality bins, so as not to bias the binning. For example, 5% central would no longer mean that 5% of the events have a larger N off trk , but instead it means that 5% of the total weights belong to events with a larger N off trk . After the centrality bins 24 have been determined, a second pass over all 24 In analyze it is also possible to bin particles based on N off trk directly. For notational convenience, we will just refer to both of these binning options as 'centrality' bins, as for 4.3. ANALYZE EXECUTABLE events is made, in which each observable that the user requested can examine each event. Each observable again, typically, makes some cuts on which particles it counts. All of these cuts, namely cuts on transverse momentum, pseudorapidity, but also which type of particles to count, can be set by the user, where different settings can be chosen per observable. This may seem like a lot of flexibility, but this is necessary, as different experimental measurements often have different cuts, and as was mentioned before, to be able to compare to experimental values it is extremely important to use the exact same cuts as the experiment in question. In this way, the flexibility given to the user allows to compute a large number of observables at the same time, where each observable is precisely defined, by a single use of analyze,i.e. there is no need to run the executable multiple times. Another advantage of this flexibility is that the cuts on which types of particles to count makes it possible to compute observables defined in terms of 'identified particles', i.e. one can compute for example the mean transverse momentum of only charged pions, or of only protons.
In the framework provided by analyze during the second pass, each observable is allowed to output an array of what we will define 'intermediate quantities', along with a weight for each such quantity. The framework provided by analyze will then compute a weighted average for each element in the array, where not only the weights provided by the observable are taken into account, but also the weights for each event based on its initial entropy. The weighted averages are then given back to the observable, along with errors and correlations between the elements in the array, which provide the observable with all the information it needs to compute its final outputs. Also, the individual values for each event are given, which allows to, for example, compute event-by-event v 2 distributions [248]. To illustrate this, let us take the quantities v 2 {2} and v 2 {4} as examples [249]. These are both computed by the 'v2' component in analyze. For each event, 'v2' computes quantities known as 2 and 4 , which come with weights w 2 = M (M − 1) and , respectively, where M is the number of particles which satisfy the various cuts. 25 After each event has been computed, analyze computes the averages 2 = i∈eventsw i w 2,i 2 i i∈eventsw i w 2,i , 4 = i∈eventsw i w 4,i 4 i i∈eventsw i w 4,i , the remainder of the discussion the distinction is not important. 25 Precise definitions of 2 and 4 will be given in section 4.3.3.

CHAPTER 4. SIMULATION OF HEAVY ION COLLISIONS WITH TRAJECTUM
along with the standard deviations of the averages and correlations between the two measured averages. Herew i stands for the event weight of event i based on its initial entropy, i.e.w i = 1/P (s i ). The averages are then given back to the 'v2' component, which computes the final outputs where care is taken to estimate the final errors by standard error propagation methods.
In the remainder of this section, we will use analyze to compute a selection of available observables, which we will compare to available experimental data. 26 The aim of this is twofold. Firstly, as Trajectum is a new code, checking its predictions against known results is a good check that the code is relatively free of serious mistakes. For this purpose, we choose the maximum a posteriori (MAP) values obtained for √ s = 2.76 TeV in [24]. For certain choices of parameters of collide, the predictions of Trajectum should be the same as, or at least compatible with, to the result from [24]. As we will see below, this is the case. The MAP values we are using for the computation have been obtained by means of a Bayesian fit to experimental data, and as such represent some of the most state-of-the-art knowledge concerning which inputs to the simulation fit the data best. However, we would like to improve this analysis further by adding additional experimental data to the analysis. This brings us to the second aim of the remainder of this section.
Any additional experimental data should have the following properties to be useful in the Bayesian analysis, where by useful we mean that the additional observable provides an additional constraint to data. Firstly, the observable should have small statistical uncertainties for the number of events that we can feasibly generate for each set of input parameters to collide. The results shown below have been generated using 10000 minimal-bias events, and with limited computation time available it is not feasible to perform the Bayesian analysis with significantly more events. This means that the error bars shown in the plots below are what we can reasonably expect to achieve in the Bayesian analysis. In such an analysis, we need to be able to decide whether a particular observable is well-described by a particular set of 4.3. ANALYZE EXECUTABLE inputs to collide, or whether there is enough tension to exclude a particular parameter set. If the error bars are too large, this can not be achieved, and hence the particular observable will not give an additional constraint.
Another requirement for an observable to be useful in the Bayesian analysis is that the constraints it gives on the input parameters should not coincide with the constraints which are already given by other observables. To give an obvious example of this, imagine we would like to use the mean transverse momentum of positive pions in the analysis. This gives a constraint on the input parameters. In particular, it will constrain the bulk viscosity ζ. However, given that the mean transverse momentum of all pions is already being used in the fit, and since there is no significant difference in particles and their antiparticles in collisions at 2.76 TeV, the constraint given by the positive pions will completely coincide with the constraint from all pions. With regards to this requirement, the most interesting outcome of the analysis of the MAP values would be to find an observable which is completely incompatible with the prediction from [24]. Such an outcome could mean two things. Either there is a subset of input parameters to collide which is compatible with the data used in [24], or the model used to simulate the events is missing an ingredient which is essential to describe the new variable.

Charged particle multiplicities and spectra
Let us now start with the comparison of Trajectum to experimental data. In figure 4.1, the charged particle multiplicity is shown, where we only include particles which satisfy |η| ≤ 0.5. The resulting multiplicity is then divided by the size of this pseudorapidity range, so that we arrive at a quantity which is relatively independent of the particular pseudorapidity cut. In figure 4.1, also the experimental result from ALICE is shown [250]. One can see that Trajectum underestimates the experimental result by about 5-10%. This is no cause for alarm though, as [24] shows a compatible underestimation.
An observable which was not used in [24] is the transverse momentum spectrum, which is shown in figure 4.2 for events in the 0-5% centrality interval. Note that this figure shows the transverse momentum spectrum of identified particles, which reflects in the fact that not a cut in pseudorapidity, but in rapidity (see (1.10)) has been made. Also, in order to define the spectrum in such a way that it does not scale significantly with the size of the transverse momentum bins and the rapidity cut, the multiplicity in each bin has been divided by both the bin size and the rapidity range. Even though  [24], and experimentally measured data from ALICE [250], are shown.
we don't compare the transverse momentum spectrum with experimental data, one can see that for heavier masses the momenta shift to higher values, something we will see quantatively in the next subsection. Also, one can clearly see that even for rare particles like Ξ or maybe even Ω, it is possible to obtain reasonable statistics for the transverse momentum spectra of identified particles.

Mean transverse momentum
In figure 4.3, we show the mean transverse momentum of charged identified particles, where we have a rapidity cut of |y| ≤ 0.5. The mean transverse momentum is quite simply defined as the mean of the transverse momentum of all particles satisfying the cuts in all events. To mold this definition into a weighted average over events, we define for each event with w the weight associated to the event, and M the number of particles in the event which satisfy the cuts. This definition guarantees that if we then average p T event over all events with the given weights, the method agrees with the experimental one. Figure 4.3 shows excellent agreement with [251]. This is no surprise however, as the result from [24] agrees with the experimental data equally well, because p T was used in their Bayesian analysis.
Another observable we will examine is the mean transverse momentum fluctuations, which is shown in figure 4.4, where we include only charged particles with pseudorapidity |η| ≤ 0.8. The mean transverse momentum fluctuations are defined as [252]: where we average over all particle pairs (i, j) occuring in the same event.
To write this in the form required by analyze, we define the following two quantities for each event, in addition to the mean transverse momentum:   [24], and experimentally measured data from ALICE [251], are shown.
where each of these quantities have weight M (M − 1) and where M is the number of particles in the event which satisfy the cuts. 27 With these definitions, we now have that where the angle brackets are the properly weighted averages of C 1 and C 2 over all events. The result from Trajectum shown in figure 4.4 is in good agreement with [252], which is again unsurprising, because this observable was also used in the Bayesian analysis of [24], with good results.

Anisotropic flow
The anisotropic flow measures how anisotropic the azimuthal distribution of particles emitted from an event is. As this can vary substantially on an eventby-event basis, one usually defines observables as an average over a large number of events in a centrality class. Additional observables look at the amount of variation of anisotropic flow within a centrality class, which will be  [24], and experimentally measured data from ALICE [252], are shown.
covered in the next subsection. In this section, we define the flow coefficients v n {k}. Several of these are shown in figure 4.5 for charged particles with 0.2 GeV ≤ p T ≤ 5 GeV and |η| ≤ 0.8. To define the flow coefficients, let us first define for each event the cumulants [253]: where M is again the number of particles in the event which satisfy the cuts, and φ i is the azimuthal angle of particle i. We subsequently define for each n:    [24], and experimentally measured data from ALICE [249], are shown.
If we now average 2 and 4 with these weights, then for each n we obtain 2 and 4 , in terms of which are defined v n {2} = 2 , v n {4} = 4 2 · 2 2 − 4 .
In figure 4.5, one can see that, as with the previous observables, there is excellent agreement between Trajectum and experimental data [249]. Again, this is in agreement with [24] as well, who used exactly these experimental findings in their analysis.

Event-by-event anisotropic flow
As briefly mentioned above, one can also look at the amount by which the anisotropic flow, in this case specifically v 2 , fluctuates event by event. This event-by-event v 2 is shown in figure 4.6 for charged particles with |η| ≤ 2.5 for events in the 30-35% centrality interval. The method used in this figure for computing the v 2 of a single event is called the single particle method [248]. It can be defined in terms of the cumulants (4.13) as Event-by-event v 2 of charged particles with |η| ≤ 2.5 for events in the 30-35% centrality interval, using the single particle method. Both the result from Trajectum using the 2.76 TeV MAP values from [24], and experimentally measured data from ATLAS [248], are shown.
where M is again the number of particles in the event which satisfy the cuts. We can then store the v 2 values for each event, and then compute a probability distribution around the average value as in figure 4.6. As can be seen, even though this observable was not used in the Bayesian fit in [24], there is agreement with the experimental value [248], but the error bars are quite large. Also, the tail of the distribution is expected to be the most constraining [254], especially to the initial condition, but this is entirely missing from the Trajectum prediction, and would require more events to compute. Maybe one could still achieve a non-trivial constraint on the inputs by choosing larger centrality bins though. Another potential issue with this observable is that the experimental value [248] has been modified using an unfolding procedure, which attempts to remove the bias inherent in defining v 2 through (4.14). This implies that the two observables shown in figure 4.6 are not the same, and should not be compared. However, in the limit of the number of particles per event going to infinity, this mentioned bias in defining v 2 through (4.14) vanishes. Therefore, by studying the dependence of the Trajectum result shown in figure 4.6 on the acceptance cuts, particularly on the pseudorapidity, it may be possible to extrapolate the result to this desired limit. This would allow the comparison to the ATLAS data, but we leave this for future work.

Discussion and Outlook
In this thesis, we have studied QCD from three different angles, corresponding to the three main chapters. A recurring theme throughout this has been the idea to use holography to study QCD-like theories in a genuinely strongly coupled setup. We saw that this works rather well in many cases. In the case of observables at zero chemical potential, one can obtain qualitatively reasonable results which match results from lattice QCD, in which case holography can be used to try to understand the mechanisms behind those phenomena seen on the lattice. In the case of finite chemical potential, one can use holography in a regime where no results from other strongly coupled methods exist. This also gives qualitatively reasonable results, which do not directly contradict known experimental results. In this chapter, we will briefly summarize the main results, and provide an outlook for how the work in this thesis could be improved or extended in the future.

(Inverse) magnetic catalysis in holographic QCD
In chapter 2, we studied the V-QCD model in the presence of a magnetic field B and anisotropy a, which could be either parallel or perpendicular to B. This model is constructed such that it matches qualitative features of QCD, like the running of certain coupling constants in the UV, glueball spectra, and the phase structure at finite temperature. In particular, V-QCD displays confinement and spontaneous chiral symmetry breaking at low temperatures, while these features are lacking at high temperatures. This makes V-QCD a good tool to study the behavior of the chiral condensate associated to chiral symmetry breaking in the presence of a magnetic field. At zero temperature in a perturbative setting, one expects the chiral condensate to increase as a function of B. However, in lattice QCD studies, the opposite effect was seen, called inverse magnetic catalysis. In section 2.2, we saw that with a judicious choice of the parameter c occuring in the w-potential of V-QCD, it is possible to obtain the same behavior of the chiral condensate as seen on the lattice. Indeed, we saw that at small temperatures, the condensate increases as a function of B, whereas in a region between the chiral phase transition and the deconfinement transition, the condensate decreases as a function of B. We also saw that at larger values of B, eventually inverse magnetic catalysis disappears, and the condensate starts increasing again. Furthermore, we were able to identify two competing effects, namely a direct effect of B on the condensate which tends to increase the condensate, and an indirect effect which tends to decrease it. This was in agreement with earlier studies on the lattice.
We subsequently extended the analysis to include a finite baryon chemical potential µ in section 2.3. In this region of the phase diagram, lattice QCD is affected by the sign problem. In holography, however, there are no major technical issues preventing us from exploring this region. We found that the region of chirally broken deconfined plasma which exists at zero B extends in seze at finite B, particularly at large values of the chemical potential. In contrast, in the region between the chiral transition and the deconfinement transition, the condensate decreases as seen earlier in section 2.2, but now we can also conclude that this region of inverse magnetic catalysis only extends to a finite µ, after which the region of inverse magnetic catalysis disappears. Section 2.4 explores a different approach, by adding an anisotropy a, which is dual to a space-dependent theta term. This way we can examine the effects that anisotropy has on QCD, without that anisotropy being a magnetic field. In other words, we can investigate whether the effects that a magnetic field has on the plasma are due to the magnetic field specifically, or whether they apply more generally for other sources of anisotropy. We saw several interesting effects. The chiral condensate behaved in much the same way as it did in the presence of a magnetic field, giving evidence for the conjecture first proposed in [123] that inverse magnetic catalysis is caused by the anisotropy induced by the magnetic field. We named this conjectured effect 'inverse anisotropic catalysis'.
Finally in section 2.5, we explored the interplay between a magnetic field

HOLOGRAPHIC BARYONS AND NEUTRON STARS
B and anisotropy a, where we considered both the B a and B ⊥ a cases. We find that the presence of a seems to effectively decrease the value of B, especially for the parallel case. We also see that the chiral transition temperature for B a seems to cross that of B ⊥ a roughly where B ∼ a.
In the future it would be very interesting to determine whether it is possible to find an explanation for the observed behavior in section 2.5. Another interesting avenue of research would be to see if other sources of anisotropy can be added to lattice studies, thereby testing the conjecture of 'inverse anisotropic catalysis'. Additionally, the holographic model itself can be improved. For the studies presented in this thesis, the quark flavors are assumed to all be identical. In particular, the baryon number is assumed to be equal to the electric charge. In the future, one could attempt to couple the magnetic field to the quarks in a different way for different flavors, thereby getting closer to a model for QCD. This is however likely a very challenging extension. A future study which is simpler to achieve is to use the existing model to study transport coefficients in the presence of a magnetic field. In [255,256], the complete set of hydrodynamic transport coefficients for relativistic magnetohydrodynamics was derived to first order, but the values that these coefficients take are model-dependent. In principle, by studying perturbations around the background space-times constructed in section 2.1, one can derive the values for the transport coefficients in V-QCD.

Holographic baryons and neutron stars
In neutron star physics, one of the most basic questions is one of the biggest question marks. This question is what the neutron star equation of state is. The equation of state governs observables like the maximum possible mass of a neutron star, the mass to radius relation and tidal deformability, and the post-merger spectrum of a binary neutron star merger. On the other hand the answer to the question what the equation of state is would give us valuable information about QCD, including potentially the density at which a phase transition to deconfined matter occurs. Many models exist which derive an equation of state, but none so far have used explicitly strongly coupled methods. For this reason, in chapter 3, we used holography to study the equation of state and the resulting model for neutron stars.
In the Witten-Sakai-Sugimoto model, several approximations exist to incorporate baryons into the holographic model. In section 3.1, we adapted one such approximation into V-QCD. As it turns out, this gives a qualitatively reasonable phase diagram, with baryons condensing in the region in the phase diagram where one would expect this to happen. We also saw that the distance into the holographic direction where the baryons are located is O(1), giving confidence that indeed the approximation is reasonable. Looking at the zero temperature equation of state itself, we see that it exceeds the conformal value for the speed of sound in two places, a feature that is difficult to obtain in weakly coupled theories. This lends confidence to the idea that the assumption of strong coupling leads to genuinely different predictions for the equation of state as compared to weak coupling.
In section 3.2, we used the speed of sound of the resulting equation of state, and used it as the high-density part of a hybrid equation of state, the low density part coming from the SLy equation of state. This hybrid equation of state was shown to be compatible with currently existing constraints, including the maximum mass, mass to radius relation, and tidal deformability. Subsequently a merger simulation using the hybrid equation of state was shown, which included an example of phase transition induced collapse. We also examined the post-merger spectrum, finding a value for the frequency f 1 which is compatible with a universal relation, and a value for f 2 which shifts to lower frequencies as the hybrid equation of state is changed to incorporate the holographic part of the equation of state up to lower densities.
There are many opportunities in the future to improve the analysis from chapter 3. One example is ongoing work to repeat the analysis for different choices of potentials for V-QCD, where the potentials are taken from [103]. This analysis will also include different choices for the part of the equation of state used for the low-density part. Another extension already being done is to extend the analysis to the presence of a magnetic field, to see whether the equation of state changes as a function of B. Further in the future, it is important to investigate how good the various approximations made in section 3.1 are and, where possible, to improve them. This will likely be a large effort, but this is worth doing if it increases the reliability of the results. The merger simulation can also be improved, for example by including neutrinos and increasing the resolution. One particular effect which is very important for the outcome of a merger simulation is the temperature dependence of the equation of state. This is hard to do in bottom-up holographic models such as V-QCD since this requires stringy corrections, which are hard to obtain in a bottom-up model. It may be possible though to obtain a reasonable approximation, which can then be used in a merger simulation.

Simulation of heavy ion collisions with Trajectum
In chapter 4, the new Trajectum framework was introduce, which provides a consistent interface between the various components necessary to simulate heavy ion collisions. This framework consists of two executables, collide and analyze, to simulate collisions and to analyze the result, respectively. Also, several examples of components which fit into the Trajectum framework.
Most of these are reimplementations of existing models, but some have been modified in non-trivial ways. The aim of these extensions is to incorporate results from AdS/CFT, such as the way the fluid should be initialized at the start of the hydrodynamical evolution. Subsequently, a Bayesian analysis will be performed to attempt to infer whether experimental data shows a preference for AdS/CFT, or whether it points into another direction. We then tested the code using the maximum a posteriori (MAP) values obtained from a Bayesian analysis in [24]. This yielded excellent agreement. We also added new observables, which were not used in the analysis of [24]. This therefore potentially adds new constraints to the Bayesian analysis to be performed, which can then be used to learn something about the various parameters used in the simulation in collide.
In the future, there are several ways in which Trajectum can be improved. An interesting option is to include thermal photon emission by the fluid. This could give interesting new constraints, as it is directly sensitive to fluid quantities in the center of the collision, in contrast to the quantities we are using now, which are all produced at the edge of the quark-gluon plasma. Other obvious extensions are to extend the code to optionally work in 3+1D, which provides a better physical description of the collisions, as well as adding conserved quantities other than the stress-energy tensor. The latter include baryon number density and electric charge. With the inclusion of electric charge, it would also become possible to further extend the hydrodynamics into full magnetohydrodynamics, which could be used for various interesting problems [110,257,258]. A final interesting possibility is to include critical fluctuations which occur around a potential critical point [259,260]. This would allow us to study what the effect of a critical point on observables is, and these observables could then also be used in a future Bayesian analysis.
Throughout these chapters, the overarching theme has been that holography can be used as a tool to gain qualitative insight into strongly coupled theories such as QCD. This results both in new explanations of existing observations and in new qualitative predictions in areas where holography is so far the only way to gain insight in an explicitly strongly coupled setting. Furthermore, by simulating neutron star mergers and heavy ion collisions, the ideas coming from holography can be compared to experimental observations. All of these avenues taken together have taught us a lot about strongly coupled physics, and will surely teach us much more in the future.

Acknowledgements
First of all, I would like to thank my supervisors, Raimond and Umut, for hiring me and giving me this great opportunity to conduct research in so many different areas. You have both helped a lot whenever I was stuck at something, provided interesting projects, but you also allowed me to work in other collaborations, something I've learned a lot from. You've also always been encouraging to go to conferences and summer schools to present our work and learn new things.
My supervisors were not the only people I have had the pleasure to work with during my time as a PhD. Therefore I would also like to thank my other collaborators, Ioannis, Matti, Juan, Takaaki, Christian and Wilke. Without you, the work presented in this thesis could not have been done.
I would also like to thank my officemates, Bernardo and Sonja, for insightful discussions both within and outside physics. I've had a good time in our office the past four years, and I think that you've always made sure there is a good atmosphere inside the office.
Also the entire string group has been great to discuss both physics and non-physics related topics with, so I'd like to give thanks to Brice, Chongchuo, Chris, Damian, Domingo, Eric, Huibert, Kilian, Koen, Miguel, Natale, Nava, Phil, Pierre, Ronnie, Sebastian, Stefan and Thomas.
Moving outside of physics itself, I would like to thank my two paranymphs, Lotte and Stephanie, both for the support they've given me during the writing of this thesis, and for the invaluable support they will give me during the promotion ceremony itself. Of my two paranymphs, I would like to especially thank my girlfriend Lotte, who has supported me throughout my time as a PhD. She has encouraged me to never give up when things went less well, and I thank her enormously for her support.
Also my parents plus their spouses, Annelies, Arie, Jacobien and Eveline, I would like to thank for their support. In addition, I'd like to thank the rest of my family, my in-laws, and my friends in Ceros, USBC and SUF. Lastly, let me thank two very special friends that Lotte and I have made, namely Domingo and Danna. We have both very much enjoyed our dinner and board game evenings together, and we hope to come and visit you in your new home.
In dit proefschrift wordt QCD bestudeerd uit drie verschillende richtingen, metéén overkoepelend thema: holografie. De holografische dualiteit maakt dat sommige sterk gekoppelde kwantumveldentheoriën beschreven kunnen worden in termen van veel eenvoudigere klassieke zwaartekracht inéén dimensie extra. De eerste richting van waaruit QCD bestudeerd wordt in dit proefschrift is door de effecten van een extern magnetisch veld op een specifiek holografisch model van QCD te bestuderen, wat interessant kwalitatief inzicht geeft. De tweede richting bestudeert hoe, in hetzelfde model, het mogelijk is om baryonische configuraties met grote dichtheid te beschrijven, wat een nieuwe manier oplevert om de materie te bestuderen waar neutronensterren uit bestaan. De toestandsvergelijking die op deze manier verkregen wordt wordt vervolgens inderdaad gebruikt om verscheidene eigenschappen van neutronensterren te berekenen die geobserveerd kunnen worden, of geobserveerd zullen kunnen worden in de nabije toekomst. De laatste richting bevat op zichzelf geen holografische berekeningen, maar bevat wel verscheidene kwalitatieve inzichten vanuit holografie die worden toegepast in een nieuwe zware-ionen code genaamd Trajectum. Dit zal in de nabije toekomst gebruikt worden om een Bayesiaanse analyse te doen, waar gehoopt wordt dat deze kwalitatieve inzichten uit holografie getest kunnen worden op experimentele data, om te zien hoe goed de ideeën uit holografie met het experiment in overeenstemming zijn. APPENDIX A. V-QCD POTENTIALS from which we can obtain For both of the sets of potentials considered below, we will have which matches the perturbative QCD beta function, as well as the perturbative anomalous dimension of the quark mass in QCD. Since the potentials from section A.2 are only meant to be used for x f = 1, these potentials only satisfy these properties at this precise value of x f .

A.2 Baryons
In chapter 3, we use the potential which were fitted to lattice data in [103].