Surface tension of dense matter at the chiral phase transition

If a first-order phase transition separates nuclear and quark matter at large baryon density, an interface between these two phases has a nonzero surface tension. We calculate this surface tension within a nucleon-meson model for domain walls and bubbles. Various methods and approximations are discussed and compared, including a numerical evaluation of the spatial profile of the interface. We also compute the surface tension at the other first-order phase transitions of the model: the nuclear liquid-gas transition and, in the parameter regime where it exists, the direct transition from the vacuum to the (approximately) chirally symmetric phase. Identifying the chirally symmetric phase with quark matter - our model does not contain explicit quark degrees of freedom - we find maximal surface tensions of the vacuum-quark transition $\Sigma_{\rm VQ}\sim 15 \, {\rm MeV}/{\rm fm}^2$, relevant for the surface of quark stars, and of the nuclear-quark transition $\Sigma_{\rm NQ}\sim 10 \, {\rm MeV}/{\rm fm}^2$, relevant for hybrid stars and for quark matter nucleation in supernovae and neutron star mergers.

To assess the possibility of the formation and evolution of quark matter in such systems one needs, besides the knowledge of the equation of state for nuclear matter under such extreme densities, information on the relevant time scales. At sufficiently high energy densities, one expects strongly interacting matter to become deconfined and essentially chiral [12][13][14], and thus chiral quark matter could provide the relevant degrees of freedom in the core of compact stars [15,16]. In fact, it was shown that deconfinement can happen at an early stage of a corecollapse supernova process, which could result not only in a delayed explosion but also in a neutrino signal of the presence of quark matter in compact stars [5]. However, the issue depends dramatically on the time scales of phase conversion as compared to the time during which the superdense region is probed [17,18].
Most model descriptions of strong interactions at high density and low temperature suggest a first-order phase transition for the chiral and the deconfinement transitions. If that is the case, a key ingredient is the surface tension, which sets the time scale -together with the * Electronic address: fraga@if.ufrj.br † Electronic address: hippert@if.usp.br ‡ Electronic address: a.schmitt@soton.ac.uk pressure difference at the phase transition and the temperature -for the phase conversion process [17]. Indeed, the surface tension is relevant for bubble nucleation of quark matter in supernovae [5,[17][18][19][20][21] and neutron star mergers [10,22]. It is also relevant for a possible quarkhadron mixed phase in the interior of neutron stars [23][24][25]. In this case, Coulomb effects together with the surface tension determine whether a mixed phase -being electrically neutral globally, but not locally -is preferred in a certain parameter regime over a globally and locally neutral homogeneous phase.
Ideally, the surface tension should be calculated from the underlying fundamental theory, Quantum Chromodynamics (QCD). However, current first-principle methods are not sufficient to determine the phase structure at large baryon chemical potentials, let alone to compute the surface tension at a possible first-order chiral or deconfinement phase transition. The only available rigorous methods for dense QCD are essentially effective theories for nuclear matter at relatively low baryon density and perturbative calculations for quark matter at ultra-high baryon density. The density regimes where the two approaches are valid are far apart, such that at a possible first-order transition at least one of them, very likely both, cannot be trusted [26]. Even if they happened to separately describe the low-and high-density phases reasonably well near the transition, this would be of little use for a rigorous calculation of the surface tension. The reason is that calculating the surface tension requires the knowledge of the entire potential, i.e., a single approach that contains nuclear and quark matter is required. At large baryon density, no such approach within QCD is currently available. Therefore, we currently rely on simple estimates or model calculations that contain both phases in a more or less realistic way [27][28][29][30][31][32][33][34][35][36][37][38][39] 1 .
Previous estimates of the surface tension were either performed in the framework of chiral models that lack the nuclear matter ingredient or employed two different models for nuclear and quark matter, which are glued together at the phase transition. To fill this gap, we employ a nucleon-meson model [45,46] that contains realistic nuclear matter, in the sense that its parameters are matched to the known properties of nuclear matter at saturation density. It does not contain quark degrees of freedom, thus the chirally restored phase can only be considered as a very rough approximation of dense quark matter.
Our model exhibits three, not just two, potentially stable phases: vacuum, nuclear matter, and the phase where the approximate chiral symmetry is restored ("quark matter"). As a consequence, there are three possibilities for first-order phase transitions. Besides transitions from the vacuum to nuclear matter and from nuclear to quark matter, a direct transition from the vacuum to quark matter is possible, depending on the parameter set. We identify the parameter regimes for these possibilities and compute the surface tension for all three transitions. The vacuum-quark surface tension is of phenomenological interest for the crust of strange stars [47,48] or so-called "strangelet dwarfs" which may exist if the surface tension is sufficiently small [49].
Finally, we also discuss various methods of calculating the surface tension. This is of particular interest because our model contains two mesonic condensates which have to be determined dynamically. This renders the calculation more complicated compared to the case of a single condensate. We compare numerical calculations using the domain wall and bubble profiles with approximations that merely require a numerical integration rather than solving a system of differential equations.
Our paper is organized as follows. In Sec. II we introduce the model and the approximations we use. We discuss the homogeneous phases in Sec. III, as a necessary preparation for calculating the surface tension. The computation of the surface tension is divided into two sections. In Sec. IV we explain the calculation for planar and spherical geometries, present several ways to approximate the full numerical result, and show some selected results, including temperature effects. In Sec. V we give a complete survey of the zero-temperature results for the various surface tensions in the parameter space of the model. Sec. VI presents our summary and outlook. 1 The surface tension can be calculated from first principles within lattice QCD if quarks are assumed to be unphysically heavy such that there is a first-order phase transition at zero chemical potential [40][41][42][43][44].

A. Model Lagrangian and approximations
To model nuclear matter and chiral symmetry at low temperatures and high densities, we adopt the following Lagrangian [45,46], where ψ is the nucleon spinor, and σ, π, ω µ are the mesonic fields. We have denoted ω µν ≡ ∂ µ ω ν − ∂ ν ω µ , τ = (τ 1 , τ 2 , τ 3 ) are the Pauli matrices, and µ is the baryon chemical potential. The nucleon spinor has four degrees of freedom in Dirac space and two degrees of freedom in isospin space, corresponding to neutrons and protons. We restrict ourselves to isospin-symmetric nuclear matter, such that these two degrees of freedom are degenerate. The pion decay constant and the omega mass are given by f π 93 MeV and m ω 782 MeV, respectively. The model contains a mesonic potential with 5 parameters a 1 , . . . , a 4 , , and Yukawa interactions between the nucleons and the mesons with the coupling constants g σ and g ω . The matching of these 7 parameters will be discussed below. For = 0, the Lagrangian is chirally symmetric. A (small) nonzero introduces a (small) explicit symmetry breaking, while condensation of the σ field breaks the (approximate) chiral symmetry spontaneously. Notice that the mass of the nucleon is generated entirely by chiral symmetry breaking, since there is no nucleon mass parameter in the Lagrangian. As mentioned in the introduction, the choice of the model is mainly motivated by the possibility to simultaneously account for realistic nuclear matter and for a chirally symmetric phase, and by its relative simplicity, allowing for a thorough calculation of the surface tension. Other similar choices, albeit leading to more complicated calculations, would be the extended linear sigma model from Refs. [50][51][52] or extensions of the present model to isospinasymmetric nuclear matter [53].
In what follows, we employ the mean-field approximation, i.e., we neglect mesonic fluctuations. These are expected to become particularly important at large temperatures, whereas we work only at zero or small temperatures. Moreover, we work within the "no-sea" approximation, i.e., we simply drop the vacuum term in the pressure, while a more elaborate evaluation would include renormalization, with the vacuum contribution depending on the renormalization scale. This would not render the calculation much more difficult and might be considered more rigorous from a theoretical point of view. However, the model is phenomenological to begin with, so it is unclear whether the result of a more complete evaluation is more realistic in any sense. Finally, we make use of the Thomas-Fermi approximation when computing the surface tension. This approximation has been employed frequently in the literature (e.g., Refs. [34,54]) and simplifies the calculation tremendously. In this approximation, the spatial dependence of the mesonic condensates is ignored when the fermions are integrated out, as opposed to a full solution of the Dirac equation for the fermions. Strictly speaking, this approximation is only valid for small gradients of the condensates. We shall nevertheless explore the full parameter space, including regions with large gradients, having in mind that our approximation has to be taken with care in these regions. These approximations lead to a relatively simple setup as follows.
The mesonic condensatesσ andω µ are introduced via the usual shift σ →σ + σ and ω µ →ω µ + ω µ , where now σ and ω µ are fluctuations. We do not include the possibility of omega condensation in the spatial components, which would break rotational symmetry,ω i = 0, and then omit the subscript 0 from the omega condensate for simplicity,ω ≡ω 0 . After separating the condensates, the Lagrangian becomes where we have assumed the condensates to be timeindependent but kept their spatial dependence. Notice that, being analogous to the Coulomb field in electrodynamics, ω 0 -and thusω -is not an independent dynamical field; hence the different sign in front of its gradient contribution. Its Euler-Lagrange equation, which we shall use below, has thus to be understood as a constraint rather than a minimization of the free energy. The fluctuation terms, containing σ, ω µ and π, will be neglected from now on. As customary, we have introduced the effective nucleon mass and the effective baryon chemical potential, and abbreviated the potential for the σ condensate, The sigma and pion masses are read off from the quadratic terms in σ and π, not shown explicitly in Eq. (2). In the vacuum, requiringσ = f π yields m 2 σ = a 1 + f 2 π a 2 and m 2 π = a 1 . With the pion mass m π = 139 MeV, this fixes the parameter a 1 . The constraintσ = f π can also be used to fix : f π is a minimum of U (σ) if = m 2 π f π . We fix the scalar coupling through g σ = m N /f π 10.097 with the vacuum mass of the nucleon m N = 939 MeV. We are thus left with the 4 parameters g ω , a 2 , a 3 , a 4 , whose matching to nuclear matter properties we discuss in Sec. II B. Now, integrating over the fermionic fields within the Thomas-Fermi approximation and dropping the vacuum contribution, we arrive at the free energy density where we have separated the gradient terms, such that Ω 0 depends only on the condensates themselves, not their derivatives, where T is the temperature and ε k = √ k 2 + M 2 is the single-nucleon energy. The factor 4 accounts for the two spin and two isospin degrees of freedom, which are all degenerate. We have kept the anti-nucleon contribution, e = −, although it is negligibly small for the temperature regime we are interested in. The Euler-Lagrange equations for the meson condensates are then where n s and n B are scalar and baryon densities, with the Fermi distribution f (x) = (e x/T − 1) −1 and the nucleon Fermi momentum k F = µ 2 * − M 2 .

B. Matching parameters to nuclear matter at saturation
To fix the remaining parameters g ω , a 2 , a 3 , a 4 , we make use of the following properties of infinite, zero-temperature, symmetric nuclear matter at saturation, where the saturation density n 0 and the binding energy E bind are known to good accuracy, while for the effective mass M 0 and the incompressibility K at saturation only a certain range is known (see, e.g., Refs. [1,[55][56][57][58][59] and [60,61], respectively). Of course, the values we use here for the bounds of these ranges contain some arbitrariness and should not be taken too literally. We first notice that the calculation of g ω decouples: from the definition of µ * , we know that g ωω = µ 0 − k 2 F + M 2 0 at saturation, where the chemical potential is given by µ 0 = m N + E bind = 922.7 MeV and the Fermi momentum k F is expressed in terms of n 0 with the help of Eq. (8b). From Eq. (7b), evaluated in the homogeneous case, ∇ 2ω = 0, we obtain m 2 ωω = g ω n 0 . Putting this together yields Hence for a given M 0 (and given n 0 and E bind ) we can directly compute g ω without having to choose a value for K. Next, a 2 , a 3 , a 4 are computed from the following three coupled equations: Eq. (7a) with ∇ 2σ = 0; the condition that at the nuclear matter onset the free energy Ω of nuclear matter be identical to that of the vacuum, Ω = 0; and the condition that the incompressibility K be given by a chosen value, where the incompressibility can be written as [1] where ε is the energy density. Additionally, one may use the surface tension of nuclear matter at saturation, as a constraint for the parameter space [45,46,62]. Again, the range of this quantity is more or less known, but the numbers chosen here for the upper and lower boundaries are somewhat arbitrary. Even though this constraint introduces a 5th quantity for matching the 4 parameters g ω , a 2 , a 3 , a 4 , it is conceivable that all constraints are met because K, M 0 , and Σ 0 are only known within a nonzero range. However, in contrast to the other 4 saturation properties, the surface tension is prone to errors if computed from the given Thomas-Fermi approximation, and thus it has to be considered with some care if treated as a further constraint. In fact, we shall show later that, if taken literally, (9) and (12) cannot be fulfilled simultaneously. We might simply proceed by fixing M 0 and K each to a certain value consistent with the experimental data, and compute the surface tension for the resulting set of parameters. However, a single number for the surface tension in the given phenomenological model, and within the given approximations, would not be of particular significance since the relation to its actual value in QCD is unclear. Hence, we shall rather perform a more general study by investigating the phase structure and surface tension in the parameter space of the model. The idea is firstly to obtain a range of values for the surface tension, for instance asking for the maximal value that the model allows for, and secondly to search for qualitatively different scenarios. In doing so, we keep the range for M 0 and K from Eq. (9) in mind, but also include regions where M 0 and/or K are beyond their "allowed" regime. While a given scenario might be outside the physical region in the given model, it may well occur in other phenomenological models, in more elaborate versions or approximations of the present model, or in QCD.

III. HOMOGENEOUS PHASES AND PHASE TRANSITIONS
As a preparation for the calculation of the surface tension, let us start by discussing the homogeneous phases.
To that end, we solve the stationarity equations (7) for a given baryon chemical potential µ and at temperature T = 0 for constantσ andω, so that ∇ 2σ = ∇ 2ω = 0.
In Fig. 1 we show the effective nucleon mass as a function of µ, including all unstable and metastable branches, which arise due to the multivalued nature of this function 2 . Recall that the nucleon mass is, up to a constant factor g σ , the same as the sigma condensate. The omega condensate shows the same qualitative features that are relevant for the following discussion, and thus here we can restrict ourselves to showing the behavior of only one of the condensates.
We are mainly interested in the first-order phase transitions and the associated spinodal regions since this is where we shall compute the various surface tensions 3 . We  abbreviate the three stable branches by V (vacuum), N (nuclear matter), and Q (chirally restored phase, "quark matter"). The baryon onset -the transition between V and N -is a first-order phase transition by construction. The chiral phase transition between N and Q can be of first order (upper right and lower left panels) or a crossover (lower right panel). This crossover only occurs for extremely large (unphysical) values of the incompressibility. Since our model breaks chiral symmetry explicitly, there are no second-order phase transitions in the model.
For certain parameter choices there is a first-order transition between V and Q (upper left panel). Even if our model is, in principle, not suitable to account for realistic quark matter, since it contains no quark degrees of freedom and no strangeness, we can phenomenologically identify Q with the quark matter phase. Then, this scenario is a realization of the strange quark matter hypothesis, according to which QCD would feature a direct transition from the vacuum to strange quark matter, and nuclear matter would be metastable [63,64]. If we do not impose stability (as opposed to metastability) of nuclear matter as an additional constraint, the parameter regime where this hypothesis is realized has to be kept and included in the calculation of the surface tension. Surprisingly, as we shall see in Fig. 2, this regime overlaps significantly with the parameter regime allowed by the constraints (9).
As a result of this analysis, we can compute three different surface tensions: Σ NQ , Σ VQ , and Σ VN . We mainly Qualitatively different regimes at zero temperature in the plane of effective nucleon mass at saturation M0 and incompressibility at saturation K. The lower (black) solid line separates the parameter region where nuclear matter is stable from the region where it is unstable. The upper (blue) solid line separates the regions where the transition between nuclear and quark matter is of first order and where it is a crossover. The dashed lines are relevant for the calculation of the surface tension at the first order chiral transition: above the upper dashed line there is no phase coexistence between the vacuum and the chirally restored phase for any µ, and below the lower dashed line there is no phase coexistence between nuclear matter and the chirally restored phase for any µ. The shaded bands indicate the allowed physical values for M0 and K according to Eq. (9). Right panel: First order phase transition lines in the plane of incompressibility and chemical potential for a fixed M0. Since the transition between nuclear matter and the vacuum is held fixed, the VN phase transition line is vertical. The topology of this diagram is the same for all M0. In particular, for each M0 there is a K for which all three phase transitions occur at the same critical chemical potential, which defines the lower (black) solid line of the left panel.
do so at the phase transitions, but also discuss bubbles of the stable phase immersed in the metastable phase, for which the spinodal regions are relevant. Notice, in particular, that in the presence of two first-order phase transitions, VN and NQ, the spinodal regions can overlap (as in the upper right panel of Fig. 1). in which first-order NQ and VQ transitions exist, which is useful in the calculation of the surface tension. It is only below the upper dashed line, for instance, that there is phase coexistence between the V and Q phases. In other words, above this line there exists no chemical potential at which these phases have the same free energy. In the small regime between the upper dashed line and the lower (black) solid line, V and Q have the same free energy at a chemical potential at which nuclear matter N is the ground state, and where the free energy minimum corresponding to the ground state N is between the minima corresponding to the metastable states V and Q. Therefore, in this regime, there is no stable domain wall configuration for the VQ transition and we shall not compute Σ VQ there.
Analogously, it is only above the lower dashed line that there exists, for each pair (M 0 , K), a chemical potential at which N and Q coexist. Again, there is a small stripe in the M 0 -K plane where a third phase -here the vacuum V -is preferred at the point where N and Q have the same free energy. In this case, however, the ground state has an effective nucleon mass that is not in between the effective nucleon masses of the two metastable states, and thus we do find stable domain wall configurations.
Finally, we see that the surface tension of the NQ transition must approach zero for sufficiently large K and/or M 0 , since the first-order phase transition line turns into a crossover in that region. The right panel of Fig. 2 provides a complementary view of the different phases in the plane of incompressibility versus chemical potential.
In summary, the surface tensions exactly at the critical chemical potentials can be calculated in the following regimes: Σ NQ can be calculated anywhere between the lower dashed line and the upper (blue) solid line and Σ VQ can be calculated anywhere below the lower (black) solid line. Since there is always a first-order transition between V and N by construction, we can compute Σ VN throughout the parameter space, although this includes regions where V and N are metastable phases. As mentioned above, the "allowed" regime, given by the shaded rectangle where the two shaded bands intersect, contains a significant fraction (about half of its area), where nuclear matter does not exist as a stable phase.

IV. DOMAIN WALLS, BUBBLES AND THE SURFACE TENSION
The surface tension is obtained from classical configurations connecting different homogeneous phases, which satisfy Eqs. (7) [65][66][67][68]. It can be calculated in different geometries and using different methods. The geometries we consider here are domain walls and bubbles. For the former, the surface tension is unambiguously defined as the free energy difference per unit area between the domain wall configuration and the homogeneous configuration of either phase, i.e., the surface tension is the energy cost that arises from abandoning the first minimum and walking uphill and downhill through the potential to reach the second minimum. On the other hand, stable solutions for a bubble are obtained in the spinodal region, where both phases exist as local minima but have different free energies. In this case, the condensates interpolate between the "false vacuum" (far away from the bubble at r = ∞) and some value close to the "true vacuum" (in the center of the bubble at r = 0). As we approach the phase transition, the value inside the bubble approaches the true vacuum and the radius of the bubble, which has to be determined dynamically, approaches infinity, i.e., we approach the domain wall solution [69].
In the present section, we select specific parameter sets to present the domain wall and bubble profiles, discussing nonzero-temperature effects, approximations and their relation to the full numerical result. These approximations are addressed in some detail due to the existence of two condensates in the present model. If only a single condensate were present the calculation would be straightforward. Numerical calculations of the domain wall or bubble profiles in related models with more than one condensate can be found in Refs. [34,54].

A. Domain walls
In the domain wall geometry the Euler-Lagrange equations (7) become After having identified the first-order phase transitions in Sec. III, we can solve this system of differential equations at the phase transition with the boundary conditionsσ(x = ±∞) =σ ± ,ω(x = ±∞) =ω ± , where the pairs (σ − ,ω − ) and (σ + ,ω + ) are the solutions of the homogeneous equations, corresponding to the two phases that have the same free energy at the phase transition [69]. We solve the equations numerically via successive over-relaxation (a useful method employed in similar contexts, e.g., in the calculation of flux tube profiles [70,71]).
Once the numerical solution is obtained, the surface tension is computed from where Ω hom ≡ Ω 0 (σ − ,ω − ) = Ω 0 (σ + ,ω + ) is the free energy of either homogeneous phase far away from the domain wall. The form of the surface tension (14) follows immediately from the free energy density (5) because Σ is defined as the free energy difference per unit area between the domain wall configuration and the homogeneous phase. The brute force numerical calculation is the most direct way to compute the surface tension; within the given Thomas-Fermi, mean-field, and no-sea approximations (and up to negligibly small numerical inaccuracies), it corresponds to the exact result. We show the domain wall profiles for a certain choice of the parameters M 0 and K and various temperatures in Fig. 3. The parameter choice is such that nuclear matter is stable and at zero temperature the baryon onset is succeeded by a first order chiral transition.
In the upper left panel we extend these two first-order phase transitions into the T -µ plane to show the critical chemical potentials at which the surface tension is calculated. We plot the entire chiral phase transition line, up to its intersection with the temperature axis, where it remains of first order, in contradiction with lattice results for QCD. One can check that, within our approximation, this is the case for the entire parameter space (in contrast to the chemical potential axis, where the chiral transition can become a crossover, as discussed above). However, this apparent contradiction with QCD is irrelevant since we should not trust our mean-field calculation for large temperatures. In what follows, we restrict ourselves to temperatures below 100 MeV, having in mind that even at these temperatures fluctuations might yield significant corrections to our result.
In the upper right panel the surface tension is shown as a function of temperature. Not surprisingly, it decreases with increasing temperature. In particular, it vanishes in the case of the VN transition 4 at the critical point which, for the given parameters, occurs at around T 19.6 MeV, in good agreement with experiment and more sophisticated nuclear physics calculations (see, e.g., Ref. [72] and references therein). Notice that our zerotemperature surface tension Σ VN is smaller than that of real-world nuclear matter (12) for the given parameter set. We shall see below that one can only fulfill  Eq. (12) in a parameter regime where nuclear matter is metastable or by going beyond the regime given by Eq. (9). The dashed lines in the upper right panel, barely distinguishable from the solid lines, show the result of an approximation which we explain now.

Semi-analytical approximation
It is useful to implement a simpler approximation for the surface tension, which does not require a numerical solution of the coupled differential equations. As Fig. 3 suggests, the results of the following approximation are very close to the full numerical results. And, in any case, we should keep in mind that we have already employed various approximations to set up the profile equations, which thus are not exact to begin with.
The approximation is derived with the help of the first integral of motion so that we can write the surface tension (14) as Assumingσ andω to be monotonic functions of x, we have introduced the functionω(σ), denoted its derivative with respect toσ by a prime, and, in the last step, rewritten the spatial derivative ofσ with the help of Eq. (15). So far, this is merely an alternative way of writing the surface tension, and not much is gained yet because the functionω(σ), which appears in Ω 0 and whose derivative appears explicitly, can only be obtained from the full solution of the differential equations. Notice that this is different in the case of a single condensate, where Eq. (16) only requires knowledge of the potential itself and reduces the calculation of the surface tension to a numerical integration.
For two condensates, we can simplify the calculation of the surface tension by approximatingω(σ): we drop the gradient ofω in Eq. (13b), such that Eq. (13b) becomes m 2 ωω = g ω n B , which can be solved to findω(σ). This has to be done numerically, even at zero temperature, where the equation becomes a sixth-order polynomial inω. By taking the derivative with respect toσ on both sides of this equation, we obtain a semi-analytical expression for the derivative [containing the numerical functionω(σ)], Here we have added the subscript "app" to emphasize that this expression is approximate because it is not identical to the derivative one obtains from the full numerically calculated domain wall profiles. We can thus write our approximate result as We refer to this approximation as the "semi-analytical approximation". It is used for the dashed lines in the upper right panel of Fig. 3, which are in excellent agreement with the full result, and shall also make use of it in Sec. V.

One-condensate approximation
Another useful approximation is to replace Eq. (13b) by m 2 ωω = g ω n B from the beginning and use the resulting functionω(σ) in the Euler-Lagrange equation forσ (13a). We can then compute the domain wall profile by solving a single differential equation numerically or, equivalently, by rederiving an expression for the surface tension that only requires a numerical integral. Were we only interested in the domain wall geometry, this would not be necessary since, as we have just demonstrated, we can solve the full system numerically with the relaxation method or can employ the semi-analytical approximation (18). However, we were not able to make the relaxation method work for the coupled differential equations in the case of the bubble geometry, where it turned out to be difficult to prevent the iterative procedure to relax to the trivial solution. Therefore, we shall make use of this "one-condensate approximation" in the calculation of the bubble profiles. Nevertheless, let us first derive the surface tension within this approximation in the domain wall geometry, as a foundation to define the surface tension of bubbles below. The difference with respect to the previous derivation is that now Eq. (15) becomes Inserting this into the surface tension (14) yields Notice that once we have replaced the Euler-Lagrange equation forω, Eqs. (19) and (20) follow without further approximations. In particular, now there is only one functionω =ω app since we do not solve the twocondensate differential equations.

C. Bubbles
Assuming spherical symmetry, Eqs. (7) become To compute the profiles of stable bubbles, one first needs to identify the spinodal region around the firstorder phase transition. Then, a stable bubble profile is found by solving Eqs. (21) with the boundary conditions σ(∞) =σ ∞ ,ω(∞) =ω ∞ , dσ dr r=0 = dω dr r=0 = 0, where (σ ∞ ,ω ∞ ) is the solution of the homogeneous equations corresponding to the phase with the larger free energy ("false vacuum"). The values of the condensates in the center of the bubble have to be determined dynamically [69].
In the calculation of the bubble profiles we employ the one-condensate approximation discussed above. We solve the differential equation forσ straightforwardly with a shooting algorithm. The surface tension of the bubble is not uniquely defined because there is no homogeneous reference state, unless one defines a sharp bubble radius. We require the surface tension of a bubble to approach the surface tension of the domain wall at the phase transition. In the given one-condensate approximation, this  Table I: Overview of methods and approximations used to calculate the surface tension Σ in the planar (domain wall) and spherical (bubble) geometry, together with the equations that show the respective expressions for Σ. The first two lines are used for our main results in Sec. V and shown to be in very good agreement with each other. The approximation of the third line is only used in Fig. 4, since the methods from first and second lines fail for bubbles. "Relax" and "shoot" refer to full numerical calculations of the profiles with relaxation and shooting algorithms. "Integrate" refers to a numerical integration, without having to solve any differential equation.
surface tension is given by Eq. (20), so that we may define the bubble surface tension as whereω(r) =ω(σ(r)) withω(σ) from m 2 ωω = g ω n B and σ(r) determined numerically from solving for the bubble profile. One might argue that, even though we have approximated the problem by a single differential equation, we could as well use the original expression (16) for the surface tension. However, we have found that rederiving the surface tension consistently within the onecondensate approximation gives a much better approximation to the full result. To keep track of the various methods and approximations used in this paper, we have collected them in Table I. In Fig. 4 we plot the surface tension for the spinodal regions of the VN (left panel) and NQ (right panel) transitions for M 0 = 0.75m N and K = 300 MeV. These parameters are the same as in the upper right panel of Fig.  1. In the left panel, for chemical potentials smaller than the onset chemical potential µ 0 , stable vacuum bubbles exist in metastable nuclear matter N, while for chemical potentials larger than µ 0 nuclear matter bubbles exist in the metastable vacuum V. In the right panel, analogously, we have N bubbles in Q below µ c and Q bubbles in N above µ c . The (black) dots away from the phase transition are the results from the numerical calculation of the bubble profiles.
As the phase transition is approached, the numerics get more challenging because the bubbles become larger and the change in the condensate occurs in a very small range of r compared to the range that needs to be considered in the calculation. Therefore, there are no results close to the phase transition, and we have employed an interpolation to cover this regime (thin solid line). Exactly at the phase transition we can calculate the surface tension from the domain wall profile, either by solving the single differential equation in the planar geometry or, equiva-lently, by using the second line of Eq. (20). This result is also shown as a (black) dot. Its agreement with the interpolated result is a check for the shooting algorithm. We have also indicated the result at the phase transition from the full two-condensate calculation [(red) diamond] and from the semi-analytical approximation (18) [(blue) square]. We see that in the case of the chiral phase transition all results are in very good agreement, while for the baryon onset there is a deviation of the order of 10% 5 .
There is an obvious qualitative difference between the left and right panels of Fig. 4 regarding the behavior at the boundaries of the spinodal region. The "standard" scenario is that the spinodal region is bounded from both sides by points at which the second, metastable solution ceases to exist. Therefore, as we move to the edges of the spinodal region, the local minimum that is assumed far away from the bubble becomes more and more shallow. As a consequence, the bubble profile flattens out and as we approach the spinodal boundary the surface tension approaches zero. This scenario is borne out in the left panel.
The different behavior exhibited in the right panel is best understood by consulting the upper right panel of Fig. 1. We see that, as we move from the NQ phase tran- 5 In Ref. [46], the surface tension Σ VN 1.1 MeV fm −2 is quoted for a very similar parameter set as we use for Fig. 4. Our result (say from the full domain wall calculation, but also from the approximations) is significantly smaller compared to that value, almost by a factor 1/2. We found that if we work with the approximation of calculatingω(σ) from the homogeneous equation m 2 ωω = gωn B , but then completely ignore theω term in the surface tension, i.e., if we use we do obtain Σ VN 1.07 MeV fm −2 , in agreement with Ref. [46]. It is clear from the comparison with the full result that this approximation is too simplistic. sition towards larger values of the chemical potential, the metastable branch of nuclear matter terminates. Therefore, this part of the spinodal region shows the usual behavior and the surface tension of the bubble goes to zero. As we move towards smaller values of the chemical potential, however, it is the stable N branch, not the metastable Q branch that terminates first. More precisely, as we go to smaller values of µ, the N branch first becomes metastable itself, indicated by the solid (blue) vertical line in the right panel of Fig. 4, before this branch disappears, indicated by the left vertical dashed line. Therefore, between these two vertical lines the bubble profile interpolates between two metastable phases, while the actual stable phase is the vacuum V (i.e., one could also compute the surface tension of V bubbles immersed in Q in that region). Since the minimum of the metastable phase Q never becomes shallow in the spinodal region (this minimum only ceases to exist at a much smaller µ), the surface tension remains large throughout this segment of the spinodal region.

V. SURFACE TENSION OF COLD AND DENSE MATTER
While in the previous section we have picked specific parameter sets to discuss temperature effects, surface tension of bubbles, and various methods to calculate the surface tension, here we show results for the surface tension in a large region of the parameter space, restricting ourselves to the surface tension of domain walls. As Fig.  4 suggests, the surface tension of finite bubbles is not expected to differ significantly from the one of infinite bubbles (domain walls), except for the fact that it can become arbitrarily small at the edges of the spinodal region.
We plot the surface tension for all three possible firstorder phase transitions in Fig. 5 for two values of the effective nucleon mass at saturation, M 0 , as a function of the incompressibility at saturation, K. The parameters of the model are then, at each point, adjusted such that saturation density n 0 and nuclear binding energy E bind are held fixed to their physical values, as explained in Sec. II B. In both panels the values for K ranges from zero to the point where the transition between nuclear matter N and the chirally restored phase Q turns into a crossover. Therefore, Σ NQ goes to zero at the upper end of the K scale. The surface tensions are computed everywhere where they are defined, i.e. when the curves end there is no longer a first-order phase transition between the respective phases (see discussion of Fig. 2 in Sec. III). We have included results from the full numerical calculation (solid lines) as well as from the semi-analytical approximation from Eq. (18) (dashed lines). As already suggested from the selected results in the previous section, the approximation (18) is in very good agreement with the full result.   In Fig. 6 we plot the surface tensions in a different way, scanning the M 0 -K parameter space more systematically, but keeping the parameters constrained to the physical regime given by Eq. (9). For simplicity, we have used the semi-analytical approximation (18) for this figure, because computing all curves with the full numerical procedure would have been very laborious. Since we have seen that this approximation is very good, the full result is not expected to differ much (in most parts of the plot it would be indistinguishable by naked eye from the approximation). In both panels of this figure, the "allowed" shaded rectangle from the left panel of Fig. 2 is visible in a distorted way. In addition, it has acquired a "scar". This scar appears by calculating the surface tensions of both the VQ and NQ transitions at the point where all three phases V, N, Q have the same free energy [lower (black) solid line in the left panel of Fig. 2 and tricritical point in the right panel of Fig. 2]. There is no reason why the VQ and NQ surface tensions should be identical at this point, thus the distorted rectangle is broken.
From Figs. 5 and 6 we draw the following conclusions: • The surface tension between vacuum and quark matter is larger than that between nuclear and quark matter. This is easy to understand already from Fig. 1: the jump in the effective nucleon mass M -in fact the jump in both condensatesσ,ω -is larger for the VQ transition than for the NQ transition. In a domain wall (or a large bubble) this jump has to be bridged and thus a large difference in condensates inevitably leads to large gradients and/or a wide wall and thus to a large surface ten- sion.
• The constraint given by Eq. (12) for Σ VN cannot be fulfilled simultaneously with the constraints in Eq. (9) for K and M 0 , as the left panel of Fig. 6 shows.
In the left panel of Fig. 5 we see that smaller values of M 0 -stretching the allowed region somewhat, but perhaps still consistent with experiment [58,59] -are needed to reach the correct value for the vacuum-nuclear surface tension. In this parameter regime, nuclear matter is metastable.
• At fixed M 0 , the surface tension of the chiral transitions, VQ and NQ, decreases if nuclear matter at saturation is made stiffer (increasing K). In contrast, the surface tension of the VN transition increases. The surface tension of the chiral transition tends to be smaller if the transition occurs at larger baryon densities.
• For the surface tensions of the chiral transitions within the physical values of K and M 0 we find Σ VQ 15 MeV/fm 2 ∼ (80 MeV) 3 and Σ NQ 10 MeV/fm 2 ∼ (70 MeV) 3 . These maximal values become larger if we allow for softer nuclear matter, but not by much: even in the extreme limit K → 0, and using M 0 = 0.75 m N , the maximum value is only slightly larger, Σ VQ 20 MeV/fm 2 ∼ (90 MeV) 3 . These values are in agreement with several existing calculations that use some model description of quark matter, but do not contain realistic nuclear matter (such as the Nambu-Jona-Lasinio model or quark-meson models) [27, 28, 31, 33-36, 38, 39]. They are therefore complementary to our approach, which does contain nuclear matter but only a toy version of quark matter. There are also studies that predict larger surface tensions, up to about 100 MeV/fm 2 or even larger, but they either do not use a single model for the two phases at the phase transition [32,37,73] (and, additionally, a simple approximation for the surface tension [74,75]), or they are based on simple estimates from dimensional analysis, not on a microscopic calculation [29,30].

VI. SUMMARY AND OUTLOOK
We have calculated the surface tension of dense matter within a nucleon-meson model, which accounts for realistic nuclear matter, but does not contain quark degrees of freedom. We have investigated the parameter space of the model systematically to locate the possible firstorder phase transitions. In doing so, we have matched the parameters to properties of nuclear matter at saturation, which defines the allowed window in the parameter space.
First-order transitions can occur between the vacuum, nuclear matter, and quark matter (more precisely, the approximately chirally symmetric phase). In particular, we have identified the parameter regime where there is a direct transition between the vacuum and quark matter, i.e. where nuclear matter is metastable. In this regime, we have computed the vacuum-quark surface tension, which assumes values of Σ VQ (8 − 15) MeV/fm 2 . In the remaining parameter space we have calculated the nuclearquark surface tension, resulting in somewhat smaller values, Σ NQ (2 − 10) MeV/fm 2 . These relatively small values seem to favor the formation of quark matter in the core of neutron stars via nucleation during their for-mation in a supernova explosion [17,19].
We have also discussed in detail various methods of calculating the surface tension. Besides a numerical evaluation of domain wall and bubble profiles and the resulting surface tension from the free energy, we have discussed two approximations. The first one, which we called "semi-analytical approximation", reduces the calculation of the surface tension to a numerical integral, such that no differential equation has to be solved. In the case of a single condensate, this reduction is exact. In the case of more than one condensate (in the present setup there are two) the reduction involves a small-gradient approximation. We have shown that this approximation is in very good agreement with the full result. The second approximation, which we called "one-condensate approximation", was employed for the surface tension of bubbles. In this case, we have approximated the problem by a single-condensate system, such that a numerical calculation of the bubble profiles can be done more easily. This approximation has turned out to be slightly worse: in some cases 10% off the exact value if extrapolated to infinitely large bubbles.
There are several extensions and applications of our study. We have pointed out that the spinodal regions of the first-order phase transitions can overlap, which gives rise to unconventional profiles of domain walls and bubbles. We have not investigated these possibilities systematically, and it would be interesting to see whether they may have observable consequences. Moreover, we have evaluated the model in the mean-field and no-sea approximations, and have employed the Thomas-Fermi approximation to calculate the surface tension. Obvious extensions would thus be to go beyond any of these approximations. One should keep in mind, however, that our model is of phenomenological nature, and thus can only give us a rough idea of the surface tension of dense QCD matter (provided there is a first-order chiral phase transition in dense QCD), even if improvements of the present approximations are performed. It would, of course, be interesting to calculate the surface tension in a model that shows a first-order chiral (or deconfinement) transition while describing both nucleons and quark degrees of freedom. However, it is very difficult to construct such a model (see for instance recent progress along these lines in holographic studies [76][77][78]).
Although we have restricted ourselves to a theoretical calculation of the surface tension, there are obvious applications in astrophysics, where cold and dense matter can be found in neutron stars. For applications to supernova explosions it is interesting to compute the temperature dependence of the surface tension more systematically and to estimate the associated nucleation times. It would also be interesting to apply our results to a firstorder transition during a neutron star merger, possibly affecting the gravitational wave signal.