Surface and Curvature Tensions of Cold, Dense Quark Matter: A Term-by-Term Analysis Within the Nambu–Jona–Lasinio Model

Abstract

In this paper, we conduct a thorough investigation of the surface and curvature tensions, $\sigma$ and $\gamma$, of three-flavor cold quark matter using the Nambu–Jona–Lasinio (NJL) model with vector interactions. Our approach ensures both local and global electric charge neutrality, as well as chemical equilibrium under weak interactions. By employing the multiple reflection expansion formalism to account for finite size effects, we explore the impact of specific input parameters, particularly the vector coupling constant ratio ${\eta }_V$, the radius R of quark matter droplets, as well as the charge-per-baryon ratio $\xi$ of the finite size configurations. We focus on the role of the contributions of each term of the NJL Lagrangian to the surface and curvature tensions in the mean field approximation. We find that the total surface tension exhibits two different density regimes: it remains roughly constant at around $100\mathrm{MeV}{\mathrm{fm}}^{-2}$ up to approximately 2–4 times the nuclear saturation density, and beyond this point, it becomes a steeply increasing function of $n_B$. The total surface and curvature tensions are relatively insensitive to variations in R but are affected by changes in $\xi$ and ${\eta }_V$. We observe that the largest contribution to $\sigma$ and $\gamma$ comes from the regularized divergent term, making these quantities significantly higher than those obtained within the MIT bag model.

1. Introduction

One of the major unresolved issues in modern astrophysics is understanding the internal structure and composition of neutron stars. Current theoretical and observational constraints can be explained by significantly different models, suggesting that the core may be composed of nucleonic matter, the appearance of hyperons, quark matter, or even dark matter [1,2].

The possibility of quark matter existing in the core of neutron stars is quite compelling, as current understanding of the Quantum Chromodynamics (QCD) phase diagram suggests that, under high-density and high-temperature conditions, nuclear matter would deconfine, releasing its constituent quarks. Indeed, there is some evidence that a quark–gluon plasma has formed in heavy-ion collisions [3], supporting the idea that deconfinement could also occur under the high-density conditions inside neutron stars (see [4] and references therein). If this is the case, many neutron stars, especially those with higher masses and consequently higher central densities, could be hybrid stars, consisting of a core of quark matter surrounded by hadronic matter.

Although research on hybrid stars has advanced our understanding of their general properties over the past few decades, several critical unresolved issues persist. Among these are the structure and microphysical properties of the region separating the quark and hadron phases. Studies indicate that, if finite-size effects, such as surface tension, curvature tension, and Coulomb energy, are sufficiently small, the minimization of the system’s energy tends to favor global electrical neutrality over local neutrality [5,6,7,8,9,10]. This results in the formation of what is known as a mixed phase, consisting of electrically charged geometric structures within a background of the other phase with opposite electric charge [11,12,13,14]. This mixed phase, often referred to as the quark–hadron pasta phase, features structures reminiscent of various types of Italian pasta, such as gnocchi, spaghetti, and lasagna, embedded in a uniform “sauce” of the other phase that extends over a wide range of densities within the star. On the other hand, if finite-size effects have a high energy cost, a sharp boundary would separate the inner core of quark matter from the outer layers of hadronic matter.

Despite the critical role of surface and curvature tensions, their precise values in the dense environment of neutron stars remain uncertain, as it is not possible to determine them from first-principle calculations in the low-temperature, high-density regime. Consequently, while various calculations of these parameters exist in the literature, they currently only offer a general evaluation of their potential range, and many of them disagree with each other in their qualitative implications. For example, using the thin-wall formalism, it has been found that surface tension values are generally low ($\sigma <30\mathrm{MeV}/{\mathrm{fm}}^2$) across various equations of state (EOS) models, such as the NJL model [15,16,17], the linear sigma model [17,18,19], and the Polyakov–quark–meson model [20]. Similarly, the multiple reflection expansion (MRE) formalism also results in low surface tension values (see references in [21]) when quark matter is approximated as a free particle gas, which could be a reasonable approximation at extremely high densities due to the property of asymptotic freedom [22]. However, when various interaction channels are taken into account, which better represent the non-perturbative density regime of neutron star matter, the MRE formalism yields significantly higher $\sigma$ values, as seen in studies using the NJL model [23] and the vector MIT bag model [24]. In contrast, the study of curvature tension is less extensive, although it has been shown that it alters the density range over which each type of geometric structure prevails [25].

In this context, it is clear that understanding the role of strong interactions in surface and curvature tensions is crucial to fully characterize the behavior and properties of the mixed phase. Although previous works have partially addressed this issue, a comprehensive analysis of how each term in the effective NJL Lagrangian impacts $\sigma$ and $\gamma$, as undertaken in this study, has not yet been performed. Among the terms we will analyze, the one representing repulsive vector interactions is particularly relevant. Indeed, their impact on $\sigma$ and $\gamma$ has not yet been thoroughly evaluated in the context of the NJL model. These terms are important not only for their ability to stiffen the EOS, thereby explaining the existence of recently observed two-solar-mass pulsars, but also because they can significantly increase surface and curvature tensions, as recently demonstrated in the context of the MIT bag model [24]. To address finite size effects, we will employ the MRE formalism, focusing on matter that is in chemical equilibrium under weak interactions while taking into account both local and global electric charge neutrality. Given these conditions, we will investigate the influence of the vector coupling constant, the radius of quark matter droplets, and the charge-per-baryon ratio. The present manuscript builds upon an earlier study [21] but extends it in several key aspects, such as the analysis of the contribution of each term of the NJL Lagrangian to $\sigma$ and $\gamma$, together with an analysis of the role of non-vanishing electric charge of the geometric structures.

This paper is organized as follows: Section 2 discusses the NJL model in the bulk, including two types of vector interactions. Section 3 addresses finite size effects by adopting the MRE formalism, separating the contribution of each term in the NJL Lagrangian to the surface and curvature tensions. Section 4 presents our numerical results, examining the impact on the surface and curvature tensions of varying the vector coupling constant, the radius of quark matter droplets, and the charge-per-baryon ratio. Finally, Section 5 concludes the paper by summarizing our findings, discussing their implications, and comparing them with previous works.

2. The Model in the Bulk

We begin with an $SU{\left(3\right)}_{\mathrm{f}}$ NJL effective model at zero temperature, incorporating vector interactions. The specific form of the Lagrangian under study is the following:

$$\begin{array}{cc}\hfill \mathcal{L}=& \overline{q}(i\cancel{\partial }-m)q+\frac{1}{2}{G}_S\sum _{a=0}^8\left[{\left(\overline{q}{\lambda }^{a}q\right)}^2+{\left(\overline{q}i{\gamma }_5{\lambda }^{a}q\right)}^2\right]\hfill \\ \hfill & -G_D\left[det\overline{q}\left(1+{\gamma }_5\right)q+\mathrm{h}.\mathrm{c}.\right]\hfill \\ \hfill & -\left\{\begin{array}{c}{\displaystyle \frac{1}{2}}{g}_V{\left(\overline{q}{\gamma }^{\mu }q\right)}^2\hfill \\ {\displaystyle \frac{1}{2}}{G}_V\sum _{a=0}^8\left[{\left(\overline{q}{\gamma }^{\mu }{\lambda }^{a}q\right)}^2+{\left(\overline{q}i{\gamma }^{\mu }{\gamma }_5{\lambda }^{a}q\right)}^2\right].\hfill \end{array}\right.\hfill \end{array}$$

(1)

In this equation, $q_i(i=u,d,s)$ denotes the quark fields, which are characterized by three colors and flavors, each associated with a corresponding current quark mass $m_i$. We assume $m_u=m_d$. The term with coefficient $G_S$ represents a $U{\left(3\right)}_L\times U{\left(3\right)}_R$ symmetric four-Fermi interaction, where ${\lambda }^a$ are the Gell–Mann matrices, with ${\lambda }^0=\sqrt{2/3}I$. The interaction proportional to $G_D$, known as the Kobayashi–Maskawa–’t Hooft interaction, accounts for the breaking of $U{\left(1\right)}_A$ symmetry. Regarding the vector channel, we consider two types of contributions, termed Model 1 and Model 2. In Model 1, the term governed by $g_V(>0)$ provides a flavor-independent repulsive interaction, whereas in Model 2, the term associated with $G_V(>0)$ introduces a flavor-dependent repulsion, as discussed in Ref. [26].

Within the mean-field approximation, the thermodynamic potential per unit volume is expressed as the sum of the following components:

$${\mathsf{\Omega }}_{\mathrm{MFA}}={\mathsf{\Omega }}_{\mathrm{div}}+{\mathsf{\Omega }}_{\mathrm{free}}+{\mathsf{\Omega }}_{\mathrm{cond}}+{\mathsf{\Omega }}_{\det }+{\mathsf{\Omega }}_{\mathrm{vec}},$$

(2)

where ${\mathsf{\Omega }}_{\mathrm{div}}$ stands for the regularized divergent contribution, ${\mathsf{\Omega }}_{\mathrm{free}}$ has the same functional form as the $T\to 0$ limit of the free Fermi gas, ${\mathsf{\Omega }}_{\mathrm{cond}}$ represents the condensate term, ${\mathsf{\Omega }}_{\det }$ is associated with the Kobayashi–Maskawa–’t Hooft interaction, and ${\mathsf{\Omega }}_{\mathrm{vec}}$ represents the vector term. Although the thermodynamic potential per unit volume, expressed as a sum of various components, was explicitly presented in [21], we adopt here a different approach for the first and second terms. Specifically, we now consider the second term to correspond to that of a free Fermi gas. The explicit form of each contribution is

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_{\mathrm{div}}& =& -6\sum _{i=u,d,s}{\int }_0^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{E}_i,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_{\mathrm{free}}& =& -6\sum _{i=u,d,s}{\int }_0^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}({\mu }_i^{*}-E_i),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_{\mathrm{cond}}& =& G_S\left({\varphi }_u^2+{\varphi }_d^2+{\varphi }_s^2\right),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_{\det }& =& -4G_D{\varphi }_u{\varphi }_d{\varphi }_s,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_{\mathrm{vec},1}& =& -\frac{1}{2}{g}_V{\left(\sum _{i=u,d,s}{n}_i\right)}^2(\mathrm{Model}1),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_{\mathrm{vec},2}& =& -\frac{1}{2}{G}_V\sum _{i=u,d,s}{n}_i^2(\mathrm{Model}2).\hfill \end{array}$$

(8)

In Equation (3), $\mathsf{\Lambda }$ represents an ultraviolet cutoff introduced to regularize the divergent integral, and $E_i=\sqrt{k^2+M_i^2}$ is the on-shell energy of the quark, self-consistently evaluated for the constituent quark masses $M_i$, which are dynamically generated via the NJL interactions governed by the coupling constants $G_S$ and $G_D$ [27]:

$$M_i=m_i-2G_S{\varphi }_i+2G_D{\varphi }_j{\varphi }_k.$$

(9)

In this equation, ${\varphi }_i$ are the quark condensates associated with each flavor, where $i,j,$ and k are flavor indices with $i\ne j\ne k$. The value of ${\varphi }_i$ is ascertained through a self-consistent solution of the following equation:

$${\varphi }_i=-6\int {\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\displaystyle \frac{M_i}{E_i}}[1-\theta (E_i-{\mu }_i^{*})],$$

(10)

simultaneously with Equation (9).

In the preceding expressions, the effective quark chemical potentials ${\mu }_i^{*}$ account for the influence of vector interactions and are defined as follows:

$$\begin{array}{ccc}\hfill {\mu }_i^{*}=& {\mu }_i-g_V\sum _{j=u,d,s}{n}_j\hfill & \hfill (\mathrm{Model}1),\end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\mu }_i^{*}=& {\mu }_i-G_V{n}_i\hfill & \hfill (\mathrm{Model}2),\end{array}$$

(12)

where $n_i$ denotes the particle number densities, calculated as

$$n_i\left({\mu }_i^{*}\right)=6{\int }_0^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}.$$

(13)

The Fermi momentum ${\kappa }_i$, used in Equations (4) and (13), is defined by

$${\kappa }_i=\theta ({\mu }_i^{*}-M_i)\sqrt{{{\mu }_i^{*}}^2-M_i^2},$$

(14)

where $\theta$ is the Heaviside step function. To determine ${\mu }_i^{*}$, one must substitute Equation (13) into Equations (11) or (12) and then find a self-consistent solution for these equations, simultaneously solving Equation (9).

The portion of Equation (4) containing $E_i$ could be incorporated into ${\mathsf{\Omega }}_{\mathrm{div}}$ by altering the limits of integration. However, we opt to define ${\mathsf{\Omega }}_{\mathrm{free}}$ as specified in Equation (4) because, at finite temperature, this definition ensures that the term has the same functional form as a free Fermi gas.

The complete thermodynamic potential is obtained by including the contributions of electrons, denoted as ${\mathsf{\Omega }}_e$, and a vacuum constant, denoted as ${\mathsf{\Omega }}_{\mathrm{vac}}$. Electrons are modeled as a non-interacting gas of massless fermions, and their thermodynamic potential per unit volume is given by

$${\mathsf{\Omega }}_e(T,{\mu }_e)=-{\displaystyle \frac{{\mu }_e^4}{12{\pi }^2}}.$$

(15)

The vacuum constant ${\mathsf{\Omega }}_{\mathrm{vac}}$ is introduced to ensure that the pressure becomes zero at both zero temperature and zero chemical potentials. Therefore, the complete thermodynamic potential per unit volume is expressed as

$$\mathsf{\Omega }={\mathsf{\Omega }}_{\mathrm{MFA}}+{\mathsf{\Omega }}_e-{\mathsf{\Omega }}_{\mathrm{vac}}.$$

(16)

It is worth noting that ${\mathsf{\Omega }}_{\mathrm{vac}}$ does not exert a direct impact on the values of the surface and curvature tensions, as will be elucidated below.

3. Finite Size Effects

To incorporate finite size effects into the thermodynamic potential, we adopt the MRE formalism (as detailed in [21] and its references), which, for the case of a finite spherical droplet, modifies the density of states as follows:

$${\rho }_i\left(k\right)=1+{\displaystyle \frac{6{\pi }^2}{kR}}{f}_{S,i}\left(k\right)+{\displaystyle \frac{12{\pi }^2}{{\left(kR\right)}^2}}{f}_{C,i}\left(k\right).$$

(17)

Here, the surface contribution to the density of states is

$$f_{S,i}\left(k\right)=-{\displaystyle \frac{1}{8\pi }}\left(1-{\displaystyle \frac{2}{\pi }}arctan{\displaystyle \frac{k}{m_i}}\right),$$

(18)

and the curvature contribution is given by

$$f_{C,i}\left(k\right)={\displaystyle \frac{1}{12{\pi }^2}}\left[1-{\displaystyle \frac{3k}{2m_i}}\left({\displaystyle \frac{\pi }{2}}-arctan{\displaystyle \frac{k}{m_i}}\right)\right].$$

(19)

The MRE density of states for massive quarks exhibits a reduction compared with the bulk density, leading to negative values within a specific small momentum range. To address this non-physical behavior, an infrared (IR) cutoff in momentum space is introduced. Consequently, the following replacement is necessary to compute the relevant thermodynamic quantities:

$${\int }_0^{\mathsf{\Lambda }}\cdots {\displaystyle \frac{k^{2}dk}{2{\pi }^2}}⟶{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR}}}^{\mathsf{\Lambda }}\cdots {\displaystyle \frac{k^{2}dk}{2{\pi }^2}}\rho \left(k\right).$$

(20)

The infrared cutoff ${\mathsf{\Lambda }}_{\mathrm{IR}}$ is the largest solution of the equation ${\rho }_i\left(k\right)=0$ with respect to the momentum k, and it depends on the flavor and the drop’s radius as shown in

Table 1. Infrared cutoff ${\mathsf{\Lambda }}_{\mathrm{IR}}$ for different particle masses m and different values of the drop’s radius R. In the bulk case ($R=\infty$), ${\mathsf{\Lambda }}_{\mathrm{IR}}$ is zero.

Particles	$\mathit{m}\left[\mathbf{MeV}\right]$	$\mathit{R}\left[\mathbf{fm}\right]$	${\mathsf{\Lambda }}_{\mathbf{IR}}\left[\mathbf{MeV}\right]$
quarks u, d	5.5	3	51.77
	5.5	5	32.86
	5.5	10	18.34
quarks s	135.7	3	95.07
	135.7	5	62.15
	135.7	10	33.42

.

Introducing the MRE density of states in Equation (2), the complete expression for the thermodynamic potential of a finite spherical droplet can be expressed as follows (see [21] and references therein):

$$\begin{array}{c}\hfill {\mathsf{\Omega }}^{\mathrm{MRE}}V=-PV+\sigma S+\gamma C+\cdots ,\end{array}$$

(21)

where $V={\displaystyle \frac{4}{3}}\pi R^3$, $S=4\pi R^2$, $C=8\pi R$, and the dots represent higher order terms. The coefficients corresponding to the volume, surface, and curvature in Equation (21) are recognized as the pressure (P), surface tension ($\sigma$), and curvature tension ($\gamma$) of the droplet, respectively, and are defined as

$$\begin{array}{ccc}\hfill P& \equiv & -{\displaystyle \frac{\partial \left({\mathsf{\Omega }}^{\mathrm{MRE}}V\right)}{\partial V}}{|}_{\mu ,S,C},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \sigma & \equiv & {\displaystyle \frac{\partial \left({\mathsf{\Omega }}^{\mathrm{MRE}}V\right)}{\partial S}}{|}_{\mu ,V,C},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \gamma & \equiv & {\displaystyle \frac{\partial \left({\mathsf{\Omega }}^{\mathrm{MRE}}V\right)}{\partial C}}{|}_{\mu ,V,S}.\hfill \end{array}$$

(24)

Furthermore, considering the different terms that constitute the mean–field approximation thermodynamic potential given in Equation (2), it is possible to obtain the contribution of each one of them to the total P, $\sigma$, and $\gamma$:

$$\begin{array}{ccc}\hfill P_{\mathrm{tot}}& =& P_{\mathrm{vec}}+P_{\mathrm{div}}+P_{\mathrm{free}}+P_{\mathrm{cond}}+P_{\det },\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{tot}}& =& {\sigma }_{\mathrm{vec}}+{\sigma }_{\mathrm{div}}+{\sigma }_{\mathrm{free}}+{\sigma }_{\mathrm{cond}}+{\sigma }_{\det },\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{tot}}& =& {\gamma }_{\mathrm{vec}}+{\gamma }_{\mathrm{div}}+{\gamma }_{\mathrm{free}}+{\gamma }_{\mathrm{cond}}+{\gamma }_{\det }.\hfill \end{array}$$

(27)

In the following subsections, we will deduce each of the previously mentioned contributions in detail.

3.1. Divergent Term

By taking the divergent term from Equation (3) and adopting the procedure from Equation (20), we obtain

$${\mathsf{\Omega }}_{\mathrm{div},i}^{\mathrm{MRE}}=-6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{E}_i{\rho }_i.$$

(28)

Proceeding as in Equation (21), that is, by separating the terms proportional to V, S, and C, we identify the pressure, the surface tension, and the curvature tension as follows:

$$\begin{array}{ccc}\hfill P_{\mathrm{div},i}& =& 6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{E}_i,\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{div},i}& =& -6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{\mathsf{\Lambda }}kdk{E}_i{f}_{S,i},\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{div},i}& =& -6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{\mathsf{\Lambda }}dk{E}_i{f}_{C,i}.\hfill \end{array}$$

(31)

3.2. Free Term

Applying the procedure from Equation (20) to Equation (4), we obtain

$${\mathsf{\Omega }}_{\mathrm{free},i}^{\mathrm{MRE}}=-6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}({\mu }_i^{*}-E_i){\rho }_i.$$

(32)

Splitting the previous expression as in Equation (21), we find

$$\begin{array}{ccc}\hfill P_{\mathrm{free},i}& =& 6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}({\mu }_i^{*}-E_i),\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{free},i}& =& -6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}kdk({\mu }_i^{*}-E_i)f_{S,i},\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{free},i}& =& -6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}dk({\mu }_i^{*}-E_i)f_{C,i}.\hfill \end{array}$$

(35)

3.3. Condensate Term

The condensate term is given in Equation (5), where the condenstates ${\varphi }_i$ are defined in Equation (10). Using Equation (20), the condensate contribution reads

$${\mathsf{\Omega }}_{\mathrm{cond}}^{\mathrm{MRE}}=G_S\left[{\left({\varphi }_u^{\mathrm{MRE}}\right)}^2+{\left({\varphi }_d^{\mathrm{MRE}}\right)}^2+{\left({\varphi }_s^{\mathrm{MRE}}\right)}^2\right]$$

(36)

where

$$\begin{array}{ccc}\hfill {\varphi }_i^{\mathrm{MRE}}& =& -6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\displaystyle \frac{M_i}{E_i}}[1-\theta ({\mu }_i^{*}-E_i)]{\rho }_i\hfill \end{array}$$

(37)

$$\begin{array}{ccc}& =& -6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\displaystyle \frac{M_i}{E_i}}{\rho }_i+6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\displaystyle \frac{M_i}{E_i}}{\rho }_i\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& =& -6{\int }_{{\kappa }_i}^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\displaystyle \frac{M_i}{E_i}}{\rho }_i.\hfill \end{array}$$

(39)

Proceeding in the same way as in Equation (21), the condensates take on the following form:

$${\varphi }_i^{\mathrm{MRE}}={\varphi }_i^V+{\displaystyle \frac{S}{V}}{\varphi }_i^S+{\displaystyle \frac{C}{V}}{\varphi }_i^C,$$

(40)

where we define

$$\begin{array}{ccc}\hfill {\varphi }_i^V& \equiv & -6{\int }_{{\kappa }_i}^{\mathsf{\Lambda }}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\displaystyle \frac{M_i}{E_i}},\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill {\varphi }_i^S& \equiv & -6{\int }_{{\kappa }_i}^{\mathsf{\Lambda }}kdk{\displaystyle \frac{M_i}{E_i}}{f}_{S,i},\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill {\varphi }_i^C& \equiv & -6{\int }_{{\kappa }_i}^{\mathsf{\Lambda }}dk{\displaystyle \frac{M_i}{E_i}}{f}_{C,i}.\hfill \end{array}$$

(43)

Replacing Equation (40) into Equation (36), we obtain

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_{\mathrm{cond}}^{\mathrm{MRE}}=& G_S\sum _i\left[{\left({\varphi }_i^V\right)}^2+{\displaystyle \frac{S^2}{V^2}}{\left({\varphi }_i^S\right)}^2+{\displaystyle \frac{C^2}{V^2}}{\left({\varphi }_i^C\right)}^2\right.\hfill \\ \hfill & +\left.2{\displaystyle \frac{S}{V}}{\varphi }_i^V{\varphi }_i^S+2{\displaystyle \frac{C}{V}}{\varphi }_i^V{\varphi }_i^C+2{\displaystyle \frac{SC}{V^2}}{\varphi }_i^S{\varphi }_i^C\right].\hfill \end{array}$$

(44)

We now multiply the preceding expression by the volume V, retaining only the terms that are proportional to $R^3$, $R^2$, and R, and disregarding any other powers of the radius R. The outcome is as follows:

$$\begin{array}{ccc}\hfill V{\mathsf{\Omega }}_{\mathrm{cond}}^{\mathrm{MRE}}& =& G_S\sum _i\left[V{\left({\varphi }_i^V\right)}^2+S2{\varphi }_i^V{\varphi }_i^S\right.\hfill \\ & +& \left.C\left(\frac{3}{2}{\left({\varphi }_i^S\right)}^2+2{\varphi }_i^V{\varphi }_i^C\right)+\cdots \right].\hfill \end{array}$$

(45)

From the above equation, we can easily identify the pressure, surface tension, and curvature tension associated with the condensate term, which turn out to be

$$\begin{array}{ccc}\hfill P_{\mathrm{cond}}& =& -G_S\sum _i{\left({\varphi }_i^V\right)}^2,\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{cond}}& =& 2G_S\sum _i{\varphi }_i^V{\varphi }_i^S,\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{cond}}& =& G_S\sum _i\left(\frac{3}{2}{\left({\varphi }_i^S\right)}^2+2{\varphi }_i^V{\varphi }_i^C\right).\hfill \end{array}$$

(48)

3.4. Determinant Term

The contributions of the determinant term can be derived similarly to the preceding subsection. Starting with

$${\mathsf{\Omega }}_{\det }^{\mathrm{MRE}}=-4G_D{\varphi }_u^{\mathrm{MRE}}{\varphi }_d^{\mathrm{MRE}}{\varphi }_s^{\mathrm{MRE}},$$

(49)

replacing Equation (40), and multiplying by the volume, we arrive at

$$\begin{array}{cc}\hfill V{\mathsf{\Omega }}_{\det }^{\mathrm{MRE}}=& -4G_D\times \{V\left({\varphi }_u^V{\varphi }_d^V{\varphi }_s^V\right)\hfill \\ \hfill & +S\left[{\varphi }_u^V{\varphi }_d^V{\varphi }_s^S+{\varphi }_u^V{\varphi }_d^S{\varphi }_s^V+{\varphi }_u^S{\varphi }_d^V{\varphi }_s^V\right]\hfill \\ \hfill & +C[{\varphi }_u^V{\varphi }_d^V{\varphi }_s^C+{\varphi }_u^V{\varphi }_d^C{\varphi }_s^V+{\varphi }_u^C{\varphi }_d^V{\varphi }_s^V\hfill \\ \hfill & +\frac{3}{2}({\varphi }_u^V{\varphi }_d^S{\varphi }_s^S+{\varphi }_u^S{\varphi }_d^S{\varphi }_s^V+{\varphi }_u^S{\varphi }_d^V{\varphi }_s^S)]\hfill \\ \hfill & +\cdots \},\hfill \end{array}$$

(50)

where terms of order $R^0$ were disregarded. From the latter expression, we obtain

$$\begin{array}{ccc}\hfill P_{\det }& =& 4G_D{\varphi }_u^V{\varphi }_d^V{\varphi }_s^V,\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill {\sigma }_{\det }& =& -4G_D\left[{\varphi }_u^V{\varphi }_d^V{\varphi }_s^S+{\varphi }_u^V{\varphi }_d^S{\varphi }_s^V+{\varphi }_u^S{\varphi }_d^V{\varphi }_s^V\right],\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill {\gamma }_{\det }& =& -4G_D\left[{\varphi }_u^V{\varphi }_d^V{\varphi }_s^C+{\varphi }_u^V{\varphi }_d^C{\varphi }_s^V+{\varphi }_u^C{\varphi }_d^V{\varphi }_s^V\right.\hfill \\ & & +\left.\frac{3}{2}({\varphi }_u^V{\varphi }_d^S{\varphi }_s^S+{\varphi }_u^S{\varphi }_d^S{\varphi }_s^V+{\varphi }_u^S{\varphi }_d^V{\varphi }_s^S)\right].\hfill \end{array}$$

(53)

3.5. Vector Term in Model 1

To analyze the contributions of the vector term in Model 1, we begin by

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_{\mathrm{vec},1}^{\mathrm{MRE}}=& -\frac{1}{2}{g}_V{\left(n_u^{\mathrm{MRE}}+n_d^{\mathrm{MRE}}+n_s^{\mathrm{MRE}}\right)}^2\hfill \\ \hfill =& -\frac{1}{2}{g}_V[{\left(n_u^{\mathrm{MRE}}\right)}^2+{\left(n_d^{\mathrm{MRE}}\right)}^2+{\left(n_s^{\mathrm{MRE}}\right)}^2\hfill \\ \hfill & +2n_u^{\mathrm{MRE}}{n}_d^{\mathrm{MRE}}+2n_u^{\mathrm{MRE}}{n}_s^{\mathrm{MRE}}\hfill \\ \hfill & +2n_d^{\mathrm{MRE}}{n}_s^{\mathrm{MRE}}].\hfill \end{array}$$

(54)

In the above expression, the particle number density in the MRE formalism is given by

$$\begin{array}{c}\hfill n_i^{\mathrm{MRE}}=6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}}{\rho }_i.\end{array}$$

(55)

After substituting the MRE density of states ${\rho }_i$ and separating the contributions from volume, surface, and curvature, the particle number density reads:

$$\begin{array}{c}\hfill n_i^{\mathrm{MRE}}=n_i^V+{\displaystyle \frac{S}{V}}{n}_i^S+{\displaystyle \frac{C}{V}}{n}_i^C,\end{array}$$

(56)

where we define

$$\begin{array}{cc}\hfill n_i^V& \equiv 6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}{\displaystyle \frac{k^{2}dk}{2{\pi }^2}},\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill n_i^S& \equiv 6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}kdk{f}_{S,i},\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill n_i^C& \equiv 6{\int }_{{\mathsf{\Lambda }}_{\mathrm{IR},i}}^{{\kappa }_i}dk{f}_{C,i}.\hfill \end{array}$$

(59)

Replacing Equation (56) into Equation (54), we find

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{\mathrm{vec},1}^{\mathrm{MRE}}V=-P_{\mathrm{vec},1}V+{\sigma }_{\mathrm{vec},1}S+{\gamma }_{\mathrm{vec},1}C+\cdots \end{array}$$

(60)

being

$$\begin{array}{ccc}\hfill P_{\mathrm{vec},1}& =& \frac{1}{2}{g}_V[{\left(n_u^V\right)}^2+{\left(n_d^V\right)}^2+{\left(n_s^V\right)}^2\hfill \\ & & +2n_u^V{n}_d^V+2n_u^V{n}_s^V+2n_d^V{n}_s^V]\hfill \\ & =& \frac{1}{2}{g}_V{\left(\sum _{i=u,d,s}{n}_i^V\right)}^2,\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{vec},1}& =& -g_V\sum _{i,j=u,d,s}{n}_i^V{n}_j^S,\hfill \end{array}$$

(62)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{vec},1}& =& -g_V\sum _{u,d,s}{n}_i^V{n}_j^C-\frac{3}{4}{g}_V\sum _{u,d,s}{n}_i^S{n}_j^S.\hfill \end{array}$$

(63)

3.6. Vector Term in Model 2

The vector term in Model 2 is given by

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{\mathrm{vec},2}^{\mathrm{MRE}}=-\frac{1}{2}{G}_V\sum _{i=u,d,s}{\left(n_i^{\mathrm{MRE}}\right)}^2.\end{array}$$

(64)

Replacing the particle number density $n_i^{\mathrm{MRE}}$ from Equation (56), we obtain

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{\mathrm{vec},2}^{\mathrm{MRE}}V=-P_{\mathrm{vec},2}V+{\sigma }_{\mathrm{vec},2}S+{\gamma }_{\mathrm{vec},2}C+\cdots \end{array}$$

(65)

being

$$\begin{array}{ccc}\hfill P_{\mathrm{vec},2}& =& \frac{1}{2}{G}_V\sum _{i=u,d,s}{\left(n_i^V\right)}^2,\hfill \end{array}$$

(66)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{vec},2}& =& -G_V\sum _{i=u,d,s}{n}_i^V{n}_i^S,\hfill \end{array}$$

(67)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{vec},2}& =& -G_V\sum _{i=u,d,s}{n}_i^V{n}_i^C-\frac{3}{4}{G}_V\sum _{i=u,d,s}{\left(n_i^S\right)}^2.\hfill \end{array}$$

(68)

3.7. Global/Local Charge Neutrality and Chemical Equilibrium

In this work, we are interested in studying finite-size droplets in $\beta$-equilibrium that may form, for example, within the mixed phase of a hybrid star. In such cases, chemical equilibrium is maintained through weak interactions among quarks, such as $d\to u+e^{-}+{\overline{\nu }}_e$, $u+e^{-}\to d+{\nu }_e$, $s\to u+e^{-}+{\overline{\nu }}_e$, and $u+d\leftrightarrow u+s$. Since we are considering a system at zero temperature, neutrinos leave the system freely, resulting in a vanishing chemical potential for the neutrinos. Therefore, the chemical equilibrium conditions read

$$\begin{array}{c}\hfill {\mu }_d={\mu }_s={\mu }_u+{\mu }_e,\end{array}$$

(69)

which can be expressed as

$$\begin{array}{ccc}\hfill {\mu }_u& =& {\mu }_q-\frac{2}{3}{\mu }_e,\hfill \end{array}$$

(70)

$$\begin{array}{ccc}\hfill {\mu }_d& =& {\mu }_q+\frac{1}{3}{\mu }_e,\hfill \end{array}$$

(71)

$$\begin{array}{ccc}\hfill {\mu }_s& =& {\mu }_q+\frac{1}{3}{\mu }_e,\hfill \end{array}$$

(72)

where ${\mu }_q$ is the quark chemical potential. The electric charge per unit volume is

$$n_Q=\frac{2}{3}{n}_u^{\mathrm{MRE}}-\frac{1}{3}{n}_d^{\mathrm{MRE}}-\frac{1}{3}{n}_s^{\mathrm{MRE}}-n_e.$$

(73)

For convenience, we will write the charge density in terms of the charge-per-baryon ratio:

$$\xi \equiv {\displaystyle \frac{n_Q}{n_B}},$$

(74)

where

$$n_B=\frac{1}{3}(n_u+n_d+n_s)$$

(75)

is the baryon number density. Replacing $n_Q=\xi (n_u+n_d+n_s)/3$ in Equation (73) we obtain the electric charge conservation equation:

$$0=\left(\frac{2}{3}-\xi \right)n_u^{\mathrm{MRE}}-\left(\frac{1}{3}+\xi \right)n_d^{\mathrm{MRE}}-\left(\frac{1}{3}+\xi \right)n_s^{\mathrm{MRE}}-n_e.$$

(76)

The calculation of the surface and curvature tensions presumes that there is a boundary separating the quark from the hadronic phase. Across this boundary, the electron background is uniform, i.e., it is the same on the internal and the external side of the drop’s boundary, simply because electrons do not feel the strong interaction. Additionally, the hadronic phase is predominantly positive since it is constituted mainly of neutrons and protons. Therefore, due to global charge neutrality, the quark phase inside the drop has to be negative. Thus, $n_Q$ must be negative at the inner side of the drop, i.e., $\xi \le 0$ (cf. Figure 5 of Ref. [25]). From Ref. [25] we learn that typical values of $\xi$ are in the range $-0.75<\xi \le 0$ (see also Ref. [28]).

4. Numerical Results

In this study, we have adopted the HK parameter set as described in Ref. [26], with the following specific values: $m_{u,d}=5.5\mathrm{MeV}$, $m_s=135.7\mathrm{MeV}$, $\mathsf{\Lambda }=631.4\mathrm{MeV}$, $G_S{\mathsf{\Lambda }}^2$ = 3.67, and $G_D{\mathsf{\Lambda }}^5$ = 9.29. Using these parameters, we have computed the corresponding IR cutoff for each flavor and size, as detailed in Table 1. In the figures we will use the parameter ${\eta }_V$, which represents both ${\eta }_V=g_V/G_S$ (Model 1) and ${\eta }_V=G_V/G_S$ (Model 2).

The thermodynamic potential is determined for each ${\mu }_q$ as a function of the condensates and the chemical potentials. From this, the thermodynamic quantities of interest can be calculated. Therefore, we numerically solve Equation (9) self-consistently, along with Equations (11) or (12), supplemented with the global/local electric charge neutrality condition (Equation (76)) and the chemical equilibrium conditions (Equations (70)–(72)). Note that previous calculations were performed under conditions of local electric charge neutrality ($\xi =0$) as described in [21]. One of the main focuses of this section is to extend the analysis to cases with finite $\xi$, which are of much greater astrophysical interest due to their role in determining the geometric structures that comprise the quark–hadron mixed phase in the interior of hybrid stars.

It is important to remark that, in general, when we numerically solve the set of equations related to all the conditions discussed above, there might be regions for which there is more than one solution for each value of ${\mu }_q$. To choose the stable solution among all of them, we require it to be an overall minimum of the thermodynamic potential.

4.1. Masses and Densities

In Figure 1 we begin by presenting results that are already well known in the literature [27] for the bulk case ($R=\infty$), with local electrical charge neutrality ($\xi =0$) and without vector interactions (${\eta }_V=0$). These results will be subsequently used for comparison when we consider different inputs for the parameters R, $\xi$, and ${\eta }_V$. At ${\mu }_q={\mu }_c\sim 350\mathrm{MeV}$, we observe a first-order phase transition, where $M_u$ and $M_d$ decrease sharply from $335\mathrm{MeV}$ to a significantly lower value (see Figure 1a). At this transition, the particle number density of u and d quarks jumps from zero to a finite value of approximately $2-3n_0$, while the particle number density of s quarks remains zero (see Figure 1b). However, due to flavor mixing, $M_s$ does not remain constant but instead exhibits a small jump at ${\mu }_q={\mu }_c$. Beyond ${\mu }_{\mathrm{c}}$, the contributions of ${\varphi }_u$ and ${\varphi }_d$ to $M_s$ gradually shrink, and $M_s$ remains nearly constant until ${\mu }_q$ surpasses $M_s$, leading to a non-zero $n_s$ as well. At that point, we observe a smooth crossover above the strange quark threshold.

In Figure 2, we present the constituent quark masses $M_i$ as a function of the quark chemical potential ${\mu }_q$ for quark drops of finite radii ($R=3,5\mathrm{fm}$) and the bulk scenario ($R=\infty$). We examine conditions with local electric charge neutrality ($\xi =0$) and global charge neutrality (considering $\xi =-0.2,-0.4$). Additionally, the influence of vector interactions is explored with ${\eta }_V$ values of 0, $0.5$, and 1. These analyses are conducted for both Model 1 and Model 2, encompassing all three quark flavors. Up, down, and strange quarks are represented by blue, red, and green curves, respectively. Each panel illustrates the effect of varying a single parameter while holding the others constant. The curve characterized by $\xi =-0.2$, $R=5\mathrm{fm}$, and ${\eta }_V=0.5$ is depicted with a solid line and is consistently repeated across all the panels. As in the baseline scenario depicted in Figure 1a, all curves asymptotically converge toward the current mass of the quarks. The up and down quarks undergo a continuous partial restoration of chiral symmetry, in contrast to the baseline scenario of Figure 1a, which shows a first-order transition. On the other hand, the strange quark consistently restores chiral symmetry through a crossover mechanism, akin to the pattern observed in the baseline case.

In Figure 2a,d it is shown how the curves $M_i$ versus ${\mu }_q$ are affected by changes in the parameter $\xi$, keeping $R=5\mathrm{fm}$ and ${\eta }_V=0.5$. As $\xi$ is changed from zero to increasingly negative values, there is a shift in the behavior of the constituent masses across both models shown. In the chiral symmetry broken phase, at low chemical potentials, the constituent masses $M_i$ align with the baseline model across all $\xi$ values (up to ${\mu }_q\sim 300\mathrm{MeV}$). As ${\mu }_q$ increases, there is a common qualitative behavior for the d and s quark masses, distinct from that of the u quark mass. For the d and s flavors, as $\xi$ becomes more negative, the constituent masses decrease for a fixed ${\mu }_q$. Indeed, chiral symmetry restoration occurs earlier for more negative $\xi$ values. However, the behavior for the u flavor is the opposite. At a fixed chemical potential, the constituent mass of the u quark increases as $\xi$ becomes more negative. Notably, in the most extreme case of $\xi$, chiral symmetry partial restoration for the u quark occurs at very high chemical potentials, after the restoration of d and s quarks. This effect is pronounced in Model 1. As we will discuss later, this behavior is related to the electric charge of each flavor. By imposing a more negative $\xi$, we require a suppression of u quarks, which results in an increase in their mass.

In Figure 2b,e, we show the behavior of the constituent masses $M_i$ while varying R, keeping $\xi =-0.2$ and ${\eta }_V=0.5$ constant. We observe that as R decreases, the effective mass exhibits a reduction. This effect is more pronounced at lower chemical potentials and diminishes asymptotically. Thus, we see that chiral symmetry partial restoration for all flavors is slightly anticipated when the system size decreases.

In Figure 2c,f, we present the behavior of the constituent masses $M_i$ while varying ${\eta }_V$, with $\xi =-0.2$ and $R=5\mathrm{fm}$ held constant. We observe that as ${\eta }_V$ increases, the effective mass shows a slight increase. This effect is more pronounced at intermediate chemical potentials and diminishes asymptotically. As is well known, increasing the vector interaction delays the partial restoration of chiral symmetry.

In Figure 3, we show the quark number densities $n_i$ as a function of the quark chemical potential ${\mu }_q$ for the same parameter choices of Figure 2. As before, each panel illustrates the effect of varying a single parameter while holding the others constant. Similarly, the curve characterized by $\xi =-0.2$, $R=5\mathrm{fm}$, and ${\eta }_V=0.5$ is depicted with a solid line and is repeated across all the panels. As in the baseline scenario, the particle number densities are zero at sufficiently low ${\mu }_q$; however, they increase gradually without the discontinuity seen in Figure 1b.

Figure 3a,d illustrate how the $n_i$ versus ${\mu }_q$ curves are influenced by variations in the parameter $\xi$ while maintaining $R=5\mathrm{fm}$ and ${\eta }_V=0.5$. For $\xi =0,-0.2$, the onset of u and d quarks occurs at roughly the same value of ${\mu }_q$, while s quarks appear at larger ${\mu }_q$. For the more extreme case $\xi =-0.4$, u quarks appear at a much higher chemical potential, even larger than the s quark onset. For a fixed chemical potential, as the value of $\xi$ becomes more negative, the amount of d and s quarks increases. Comparatively, there are always more d quarks than s quarks simply because the former has a lower mass. The behavior of u quarks is the opposite. As $\xi$ becomes more negative, the amount of u quarks decreases. This behavior is related to the electric charge of each flavor. Imposing a more negative $\xi$ disfavors the u quarks, which are the only ones with a positive charge. In Figure 3b,e we show the role of the drop’s radius on the particle number densities. While changes in R affect the $n_i$, these effects are significant only around the onset for each flavor. As ${\mu }_q$ increases, all curves tend to the bulk case. Finally, Figure 3c,f illustrate how the curves change with variations in the parameter ${\eta }_V$. Changes in ${\eta }_V$ affect the $n_i$, but these effects are not significant around the onset for each flavor. However, as ${\mu }_q$ increases, the effect of vector repulsive interactions becomes substantial and, for a fixed ${\mu }_q$, reduces the particle number densities as ${\eta }_V$ increases. This effect is more pronounced for Model 1.

4.2. Surface and Curvature Tensions

In Figure 4 and Figure 5, we illustrate the behavior of the surface and curvature tensions as functions of the baryon number density $n_B$ for the same parametrizations of previous figures. In a recent article [21], we studied $\sigma$ and $\gamma$ for the electrically neutral case $(\xi =0)$, without considering the separate contributions of each term in Equations (26) and (27). Here, we present a detailed breakdown of the role of these contributions.

Within each panel of Figure 4, the total surface tension (black curves) is shown along with the contributions from the divergent, free, condensate, determinant, and vector terms, as specified in Equation (26). Note that ${\sigma }_{\mathrm{tot}}$ exhibits two density regimes. Up to approximately $2-4n_0$, it remains relatively constant at around $100\mathrm{MeV}{\mathrm{fm}}^{-2}$; beyond this density threshold, it becomes an increasing function of $n_B$. The total surface tension is roughly insensitive to variations in R but is affected by changes in $\xi$ and ${\eta }_V$ starting from $n_B/n_0\gtrsim 2-3$. As $\xi$ becomes more negative (see Figure 4a,d) and when ${\eta }_V$ increases (see Figure 4c,f), an increase in the total surface tension is observed. The effect of ${\eta }_V$ grows with density and is much more pronounced for Model 1.

In all six panels, we observe that the divergent term has the most significant impact on the total surface tension. The influence of ${\sigma }_{\mathrm{div}}$ decreases with increasing density, while the effects of all other terms grow with density. Notably, ${\sigma }_{\det }$ has the smallest impact and approaches zero beyond approximately $2n_0$. The role of ${\sigma }_{\mathrm{cond}}$ is quite significant because it has a large negative value, especially in the low-density region. As the density increases, this term asymptotically approaches zero and becomes practically negligible beyond ∼$7n_0$. The fact that ${\sigma }_{\mathrm{cond}}$ is always negative can be understood as follows: the quantity ${\varphi }_i^V$ is always negative, as seen in Equation (41). On the other hand, the quantity ${\varphi }_i^S$ is always positive because $f_{\mathcal{S},i}\left(k\right)$ is always negative (see Equation (42)). Therefore, ${\sigma }_{\mathrm{cond}}$ is always negative (cf. Equation (47)). As the ${\varphi }_i$ tend to zero with increasing density, this causes ${\sigma }_{\mathrm{cond}}$ to also asymptotically approach zero. The contribution from ${\sigma }_{\mathrm{free}}$ is always positive and increases approximately linearly with density. Notice that ${\sigma }_{\mathrm{free}}$ is practically insensitive to variations in R (Figure 4b,e) and ${\eta }_V$ (Figure 4c,f). However, as $\xi$ becomes more negative, ${\sigma }_{\mathrm{free}}$ increases significantly (Figure 4a,d). The contribution from ${\sigma }_{\mathrm{vec}}$ is always positive and increases approximately linearly with density, becoming negligible at densities lower than ∼$3-4n_0$. The ${\sigma }_{\mathrm{vec}}$ term is practically insensitive to variations in R (Figure 4b,e) and $\xi$ (Figure 4a,d). However, with increasing ${\eta }_V$, ${\sigma }_{\mathrm{vec}}$ increases substantially, especially in Model 1.

Finally, in Figure 5 we display the variation of the curvature tension as a function of $n_B/n_0$. The total curvature tension ${\gamma }_{\mathrm{tot}}$ is shown along with the contributions from ${\gamma }_{\mathrm{div}}$, ${\gamma }_{\mathrm{free}}$, ${\gamma }_{\mathrm{cond}}$, ${\gamma }_{\det }$, and ${\gamma }_{\mathrm{vec}}$, as specified in Equation (27). The total curvature tension remains relatively unaffected by changes in R (Figure 5b,e), but it is influenced by variations in $\xi$ (Figure 5a,d) and ${\eta }_V$ (Figure 5c,f). Differently from ${\sigma }_{\mathrm{tot}}$, as $\xi$ becomes more negative, ${\gamma }_{\mathrm{tot}}$ decreases. However, the dependence with ${\eta }_V$ follows the same trend as for the surface tension, i.e., there is an increase in ${\gamma }_{\mathrm{tot}}$ as ${\eta }_V$ increases. The effect of ${\eta }_V$ grows with density and is much more pronounced for Model 1.

Across all panels, it is evident that the divergent term dominates the total curvature tension. Similar to the behavior observed with surface tension, the influence of the divergent term diminishes as density increases, while the effects of other terms rise with $n_B$. Notably, the determinant term has the least impact and approaches zero beyond approximately $2n_0$. The condensate term is always negative, and its influence is more significant in the low-density region. The contributions from ${\gamma }_{\mathrm{free}}$ and ${\gamma }_{\mathrm{vec}}$ are always positive and increase with density. The influence of ${\gamma }_{\mathrm{vec}}$ is particularly relevant in Model 1, especially at high densities. All terms are quite insensitive to variations in R and ${\eta }_V$, except for ${\gamma }_{\mathrm{vec}}$, which strongly depends on ${\eta }_V$.

5. Conclusions

We have undertaken a comprehensive study of the surface and curvature tensions of three-flavor cold quark matter using the NJL model, incorporating the vector channel into the Lagrangian. Our study ensures both local and global electric charge neutrality, as well as chemical equilibrium under weak interactions. Utilizing the MRE framework to address finite size effects, we examined the influence of particular input parameters, specifically the vector coupling constant, the radius of quark matter droplets, and their charge-per-baryon ratio. Special attention was paid to analyzing how each term in the Lagrangian contributes to the surface and curvature tensions. To the best of our knowledge, this is the first time such an analysis has been performed.

We find that the total surface tension ${\sigma }_{\mathrm{tot}}$ has two density regimes: it is roughly constant with a value around $100\mathrm{MeV}{\mathrm{fm}}^{-2}$ up to ∼$2-4n_0$, and, above this density, it is a steeply increasing function of $n_B$. The low-density plateau extends while the s-quarks are absent in the system. As soon as the s-quark fraction becomes finite, the total surface tension starts to increase. This is a quite EOS-independent behavior associated with the MRE density of states, which is quite sensitive to the particle’s mass. As emphasized in Refs. [23,24], the contribution of s-quarks to the total surface tension is much larger than the one of u and d quarks. This behavior can be understood if one keeps in mind that for massless particles, $f_S=0$, while for $m\to \infty$, $f_S$ remains finite. This means that light quarks tend to have a smaller contribution to the surface tension.

Our results show that the total surface tension is roughly insensitive to variations in R but is affected by changes in $\xi$ and ${\eta }_V$ starting from $n_B\gtrsim 2-3n_0$. As $\xi$ becomes more negative and as ${\eta }_V$ increases, an increase in ${\sigma }_{\mathrm{tot}}$ is observed (see Figure 4). The effect of ${\eta }_V$ grows with density and is much more pronounced for Model 1. Regarding the role of each term in the Lagrangian in ${\sigma }_{\mathrm{tot}}$, we observe that the largest impact comes from the divergent term. This influence decreases with $n_B$, while the effects of all other terms increase with $n_B$. The determinant term, in particular, has a practically negligible effect. The influence of the condensate terms is quite significant due to their large negative value, especially in the low-density region. As the density increases, this term asymptotically approaches zero and becomes practically negligible beyond ∼$7n_0$. The contributions of ${\sigma }_{\mathrm{free}}$ and ${\sigma }_{\mathrm{vec}}$ are always positive and increase with density, becoming relevant above ∼$3-4n_0$.

In contrast to the total surface tension, the curvature tension does not exhibit different density regimes, behaving consistently as a monotonically increasing function of $n_B$. This feature can be understood by examining the $f_C$ contribution to the MRE density of states. As discussed in [23,24], for massless particles $f_C=-1/\left(24{\pi }^2\right)$ and for $m\to \infty$, $f_C$ it is also finite. Consequently, light quarks play a significant role in the curvature contribution, rendering ${\gamma }_{\mathrm{tot}}$ less dependent on the emergence of s-quarks.

The behavior of the total curvature tension is mostly unaffected by variations in R but is sensitive to changes in $\xi$ and ${\eta }_V$, as illustrated in Figure 5. Unlike ${\sigma }_{\mathrm{tot}}$, ${\gamma }_{\mathrm{tot}}$ diminishes as $\xi$ becomes more negative. The relationship with ${\eta }_V$ mirrors that seen with surface tension; specifically, ${\gamma }_{\mathrm{tot}}$ rises as ${\eta }_V$ increases. The influence of ${\eta }_V$ intensifies with density.

Regarding the influence of each term in the Lagrangian on the total curvature tension, the divergent term clearly dominates ${\gamma }_{\mathrm{tot}}$. Similar to the behavior observed with surface tension, the influence of ${\gamma }_{\mathrm{div}}$ diminishes as density increases, while the contributions from other terms rise with $n_B$. Notably, ${\gamma }_{\det }$ is negligible. The condensate term is always negative and becomes more significant in the low-density region. The contributions from ${\gamma }_{\mathrm{free}}$ and ${\gamma }_{\mathrm{vec}}$ are consistently positive and increase with density. The role of ${\gamma }_{\mathrm{vec}}$ is quite relevant in Model 1, especially at high densities. All contributions are relatively insensitive to variations in R and ${\eta }_V$, except for ${\gamma }_{\mathrm{vec}}$, which strongly depends on ${\eta }_V$.

When comparing Model 1 and Model 2, it is noteworthy that the forms of the vector terms differ significantly. In Model 1, the vector term incorporates the square of the sum of the $n_i$, resulting in six quadratic terms. Conversely, in Model 2, the vector term is directly the sum of the $n_i^2$, yielding only three quadratic terms. This distinction directly influences the growth of ${\sigma }_{\mathrm{vec}}$ and ${\gamma }_{\mathrm{vec}}$ with density. In fact, in Model 1, these quantities increase much more with density than in Model 2. Therefore, to achieve comparable increases in both models for a given density, the value of ${\eta }_V$ used in Model 2 should be significantly higher than in Model 1 to compensate for the missing terms.

Another interesting aspect that arises from our results relates to the behavior of ${\sigma }_{\mathrm{tot}}$ in the low-density regime. Indeed, in the density range up to approximately $3-4n_0$, ${\sigma }_{\mathrm{tot}}$ shows a plateau with a value of about $100\mathrm{MeV}{\mathrm{fm}}^{-2}$, which is quite insensitive to variations in the parameters $\xi$, R, and ${\eta }_V$. This result is quite interesting because it shows that the hypothesis of constant surface tension, used in several works in the literature to describe the quark–hadron mixed phase, could be a reasonable approximation in the low-density regime.

Finally, it is worth comparing the results from the NJL model with those obtained in previous works using the MIT bag model. Since the MIT bag model does not incorporate dynamic masses, this comparison is particularly meaningful at the highest possible densities, which is precisely where the dynamic masses in the NJL model approach their current masses. Comparing the NJL Lagrangian of this work with that of the vector MIT bag model from Ref. [24], we see that both the vector and the “free” terms are present in both Lagrangians. Indeed, the vector term used in Ref. [24] has the same functional form as in Model 1 of the present work. On the other hand, both the determinant term and the condensate term in the NJL model are practically negligible at high enough densities. Therefore, the main difference between both models at high densities arises from the divergent term of the NJL model. If we compare the value of $\sigma$ from the MIT bag model, as given in [24], with ${\sigma }_{\mathrm{free}}+{\sigma }_{\mathrm{vec}},1$ from the NJL model, we find that the values are quite similar, and the same applies to $\gamma$. In fact, at high enough densities, the major difference in the total surface and curvature tensions between these models stems from the dominant contribution of the divergent term in the NJL model.

Our results have implications for understanding the properties of the quark–hadron mixed phase that may exist at the core of hybrid neutron stars. Although calculations in [25] show that density-dependent surface and curvature tensions have a small impact on the star’s structure (mass, radius, and tidal deformabilities), the geometric structure of the quark–hadron mixed phase is highly sensitive to how $\sigma$ and $\gamma$ behave. Our calculations demonstrate that both $\sigma$ and $\gamma$ increase as the baryon number density rises. Consequently, based on [25], we expect that drops, rods, and slabs—which form at lower densities—would be less suppressed by surface tension compared with tubes and bubbles, which emerge at much higher densities. Therefore, in the absence of curvature effects, drops, rods, and slabs would dominate the mixed phase, while tubes and bubbles would be rarely observed [25]. We also expect the curvature effects within the mixed phase to mirror those described in [25]. Specifically, slabs are not influenced by curvature, whereas bubbles exhibit a reversal in curvature energy compared with drops. This implies that surface and curvature effects partially counterbalance each other for tubes and bubbles, facilitating their stability even at higher densities and, consequently, with increased values of $\sigma$ and $\gamma$. In contrast, for drops and rods, surface and curvature effects combine, resulting in these configurations being less prevalent in the mixed phase. The characteristics of the mixed phase have significant astrophysical ramifications. When a hybrid neutron star undergoes dynamic perturbations, the interfaces between pasta structures and the surrounding matter can trigger phase conversion reactions. The diversity and extent of these structures lead to considerable variations in both the effective surface area and the conditions that facilitate these reactions. These variations are crucial, as they can influence the damping of r-modes [29] and the overall dynamic stability of hybrid stars [30,31] as well as strange white dwarfs [32]. Understanding these processes is essential for accurately modeling the behavior and evolution of hybrid neutron stars. None of these aspects were analyzed using $\sigma$ and $\gamma$ values determined within the context of the NJL model, presenting a promising direction for future research.