**2. Equivparticle Model**

As an example, in this work we adopt the equivparticle model to investigate the properties of quark matter and their nuggets. The strong interactions are included with density-dependent quark masses in the equivparticle model, while quarks are considered as quasi-free particles [53,62–71]. It is thus straightforward to write out the Lagrangian density, i.e.,

$$\mathcal{L} = \sum\_{i=u,d,s} \Psi\_i \left[ i\gamma^\mu \partial\_\mu - m\_i(n\_\mathbf{b}) \right] \Psi\_{i\prime} \tag{2}$$

where Ψ*<sup>i</sup>* represents the Dirac spinor of quark flavor *i*, *mi*(*n*b) the equivalent mass with *n*<sup>b</sup> being the baryon number density. The Coulomb interactions can also be included by adding photon field in the Lagrangian density [41].

The strong interaction among quarks are then reproduced by equivalent quark masses, where many different mass scalings were proposed. For example, for density dependent masses *mi*(*n*b) = *mi*<sup>0</sup> + *m*I(*n*b) with *mu*<sup>0</sup> = 2.2 MeV and *md*<sup>0</sup> = 4.7 MeV being the current masses of *u* and *d* quarks [72], the inversely linear scaling *m*<sup>I</sup> = *B*/3*n*<sup>b</sup> was obtained by reproducing bag model results in the limit of vanishing densities [73]. Considering the contributions of linear confinement and adopting linearization for the in-medium chiral condensates, an inversely cubic scaling *m*<sup>I</sup> = *Dn*−1/3 <sup>b</sup> was derived [74], while the one-gluonexchange interaction was later included with *m*<sup>I</sup> = *Dn*−1/3 <sup>b</sup> <sup>−</sup> *Cn*1/3 b [66]. Meanwhile, in the limit of large densities, perturbation theory suggests repulsive interactions among quarks and the quark mass scaling becomes *m*<sup>I</sup> = *Dn*−1/3 <sup>b</sup> <sup>+</sup> *Cn*1/3 b [68]. An isospin dependent term was also introduced to examine the impacts of quark matter symmetry energy, which was given by *m*<sup>I</sup> = *Dn*−1/3 <sup>b</sup> − *τiδDIn α* b e <sup>−</sup>*βn*<sup>b</sup> with *τ<sup>i</sup>* being the third component of isospin for quark flavor *i* and *δ* = 3(*n<sup>d</sup>* − *nu*)/(*n<sup>d</sup>* + *nu*) the isospin asymmetry [53].

In this work we adopt the following quark mass scaling

$$m\_{\rm I}(n\_{\rm b}) = D n\_{\rm b}^{-1/3} + \mathbb{C} n\_{\rm b}^{1/3} + \mathbb{C}\_{I} \delta^2 n\_{\rm b}.\tag{3}$$

The first term corresponds to linear confinement with the confinement parameter *D* connected to the string tension *σ*0, the chiral restoration density *ρ* ∗ , and the sum of the vacuum chiral condensates ∑*<sup>q</sup>* h*qq*¯ i<sup>0</sup> [74]. The second term represents the contribution of one-gluonexchange interaction for *C* < 0 [66] and the leading-order perturbative interaction for *C* > 0 [68], where in both cases *C* is connected to the strong coupling constant *α*s. Finally, we have added the third term in Equation (3) to account for the quark matter symmetry energy, which is increasing with *C<sup>I</sup>* . In fact, in addition to the kinetic contribution, it was shown that the formation of *u*-*d* quark Cooper pairs (2SC phase) could effectively enhance the quark matter symmetry energy [54]. The impacts of quark matter symmetry energy on compact star properties were extensively investigated in the past few years, e.g., those in References [53–57]. Note that in contrast to the mass scaling proposed in Reference [53], here *m*<sup>I</sup> is identical for different quark flavor, i.e., neglecting the isovector-scalar channel. Meanwhile, as will be discussed later, the third term leads to the isovector contributions in the vector potentials of quarks.

Adopt mean-field and no-sea approximations, the energy density *E* of quark matter at zero temperature is obtained with

$$E = \sum\_{i} \int\_{0}^{\nu\_{i}} \frac{\oint p\_{i}^{2} p^{2}}{2\pi^{2}} \sqrt{p^{2} + m\_{i}^{2}} \mathrm{d}p = \sum\_{i} \frac{\oint p\_{i} m\_{i}^{4}}{16\pi^{2}} f\left(\frac{\nu\_{i}}{m\_{i}}\right),\tag{4}$$

where *<sup>f</sup>*(*x*) = <sup>h</sup> *x*(2*x* <sup>2</sup> + 1) √ *x* <sup>2</sup> + <sup>1</sup> − arcsh(*x*) i , *g<sup>i</sup>* (= 6 for quarks and 2 for leptons) the degeneracy factor for particle type *i*, *ν<sup>i</sup>* the Fermi momentum and is linked to the number density

$$
\eta\_i = \langle \overline{\psi}\_i \gamma^0 \psi\_i \rangle = \frac{\mathcal{G}\_i \nu\_i^3}{6\pi^2}. \tag{5}
$$

Note that the masses of leptons take constant values with *m<sup>e</sup>* = 0.511 MeV and *m<sup>µ</sup>* = 105.66 MeV. The chemical potentials *µ<sup>i</sup>* and pressure *P* can then be obtained according to the basic relations of standard thermodynamics.

For *ud*QM without leptons, to investigate its saturation properties, we then expand the energy per baryon to the second order [26], i.e.,

$$
\varepsilon\_{2f}(\eta\_{\mathbf{b}\prime}f\_Z) = \varepsilon\_0 + \frac{K\_0}{18} \left(\frac{\eta\_{\mathbf{b}}}{\eta\_0} - 1\right)^2 + \varepsilon\_s (2f\_Z - 1)^2,\tag{6}
$$

where *n*<sup>b</sup> = (*n<sup>u</sup>* + *n<sup>d</sup>* )/3 is the baryon number density, *f<sup>Z</sup>* = (2*n<sup>u</sup>* − *n<sup>d</sup>* )/(*n<sup>u</sup>* + *n<sup>d</sup>* ) the charge fraction, *ε*<sup>0</sup> the minimum energy per baryon at saturation density *n*<sup>0</sup> and charge fraction *f<sup>Z</sup>* = 1/2, *K*<sup>0</sup> the incompressibility parameter, and *ε*<sup>s</sup> the symmetry energy. The corresponding values obtained with equivparticle model using various combinations of parameters are then indicated in Table 1, where the symmetry energy is found to be increasing with *C<sup>I</sup>* .

**Table 1.** The adopted parameter sets for the quark mass scaling in Equation (3) and the corresponding properties of *ud*QM.


#### **3.** *ud***QM Nuggets**

At fixed surface tension values, the energy per baryon of *ud*QM nuggets can be obtained with a liquid-drop type formula, i.e.,

$$\frac{M}{A} = \varepsilon\_{2f}(n\_0, f\_Z) + \frac{1}{5} \sqrt[3]{36\pi n\_0} a f\_Z^2 A^{2/3} + \sigma \left(\frac{A n\_0^2}{36\pi}\right)^{-1/3} \tag{7}$$

where the bulk (first) term is fixed by Equation (6). The second term represents the Coulomb energy with *α* being the fine structure constant and the third term corresponds to the surface energy. By minimizing Equation (7) with respect to *fZ*, one obtains the charge fraction for *β*-stable *ud*QM nuggets, i.e.,

$$f\_{\mathbf{Z}} = \frac{10\epsilon\_{\mathbf{s}}}{\sqrt[3]{36\pi m\_0}aA^{2/3} + 20\epsilon\_{\mathbf{s}}}.\tag{8}$$

Then we can constrain the parameters by comparing the binding energy of *ud*QM nuggets with finite nuclei, where *ud*QM nuggets should always be unstable at *A* . 300 and stable as *<sup>A</sup>* <sup>→</sup> <sup>∞</sup>. At this moment, the heaviest *<sup>β</sup>*-stable nucleus is <sup>266</sup>Hs with its energy per baryon being 931.74 MeV [75–77]. We thus require the energy per baryon obtained with Equation (7) at *A* = 266 should always be larger than that of <sup>266</sup>Hs, where a lower limit for the surface tension value *σ*min is found for various combinations of parameters. The corresponding constraints are then indicated in Figure 1 for *C<sup>I</sup>* = 0. The red curve indicates the boundary with the minimum energy per baryon *ε*<sup>0</sup> = 922 MeV, so that *ud*QM can be more stable than nuclear matter at *δ* = 0 in the lower left region. If we demand *ud*QM star matter to be more stable than neutron star matter, the stable-unstable boundary is expected to shift slightly to the lower left region. Meanwhile, the surface tension value should be larger than *σ*min, which favors the upper right region at a fixed *σ*. The combination of both constraints indicate an stability window for *ud*QM at a given surface tension value *σ*.

Finally, it is worth mentioning that the constraints in Figure 1 are insensitive to *C<sup>I</sup>* , where we have only presented the results at *C<sup>I</sup>* = 0.

**Figure 1.** Various constrains for the parameters of the mass scaling in Equation (3) and minimum surface tension value *σ*min for *ud*QM nuggets to be unstable comparing with <sup>266</sup>Hs. The open circles correspond to the choices of parameters in Table 1.

The liquid-drop type formula in Equation (7) is valid for small *ud*QM nuggets. If we want to investigate the properties of *ud*QM nuggets at larger baryon numbers, the effects of charge screening and electrons should be considered [78]. In such cases, we adopt Thomas–Fermi approximation as was done in References [26,58–61]. The total energy of the system is obtained with

$$M = \int\_0^\infty \left[ 4\pi r^2 E(r) + \frac{r^2}{2a} \left(\frac{\mathrm{d}\varphi}{\mathrm{d}r}\right)^2 \right] \mathrm{d}r + 4\pi R^2 \sigma \,\tag{9}$$

where the energy density *E* is given by Equation (4) and the electric potential *ϕ* by solving

$$r^2 \frac{\mathbf{d}^2 \varrho}{\mathbf{d}r^2} + 2r \frac{\mathbf{d} \varrho}{\mathbf{d}r} + 4\pi\alpha r^2 \left(\frac{2}{3}n\_\iota - \frac{1}{3}n\_d - n\_\varepsilon\right) = 0. \tag{10}$$

The most stable structure of an *ud*QM nugget is fixed by minimizing Equation (9), which follows the constancy of chemical potentials, i.e.,

$$
\mu\_i(\vec{r}) = \sqrt{\nu\_i(\vec{r})^2 + m\_i(\vec{r})^2} + V\_i(\vec{r}) = \text{constant.}\tag{11}
$$

Here the vector potential is given by

$$V\_{\rm i} = \frac{\text{d}m\_{\rm I}}{\text{d}m\_{\rm i}} \sum\_{i=u,d} n\_{\rm i}^{\rm s} + q\_{\rm i} \rho\_{\rm \prime} \tag{12}$$

where the scalar density is

$$m\_i^s = \langle \psi\_i \psi\_i \rangle = \frac{g\_i m\_i^3}{4\pi^2} g\left(\frac{\nu\_i}{m\_i}\right) \tag{13}$$

with *g*(*x*) = *x* √ *x* <sup>2</sup> + <sup>1</sup> − arcsh(*x*). The quantities *<sup>n</sup>* s *i* , *n<sup>i</sup>* , *ν<sup>i</sup>* , *m<sup>i</sup>* and *E* represent the local properties of *ud*QM and vary with the space coordinates. Note that a density derivative term is added to the vector potential in Equation (12), which is essential in order to maintain the self-consistency of thermodynamics [62,68,79–82]. Meanwhile, since in Equation (3) we have added the third term to account for the symmetry energy of quark matter, in contrast to our previous findings, the density derivative term takes different forms for *u* and *d* quarks, i.e.,

$$\frac{\mathbf{d}m\_{\mathrm{I}}}{\mathbf{d}n\_{\mathrm{u}}} \quad = \quad \frac{1}{3}\frac{\mathbf{d}m\_{\mathrm{I}}}{\mathbf{d}n\_{\mathrm{b}}} + \frac{\mathbf{d}m\_{\mathrm{I}}}{\mathbf{d}\delta}\frac{3-\delta}{3n\_{\mathrm{b}}},\tag{14}$$

$$\frac{\mathrm{d}m\_{\mathrm{I}}}{\mathrm{d}n\_{\mathrm{d}}} \quad = \quad \frac{1}{3}\frac{\mathrm{d}m\_{\mathrm{I}}}{\mathrm{d}n\_{\mathrm{b}}} - \frac{\mathrm{d}m\_{\mathrm{I}}}{\mathrm{d}\delta}\frac{3+\delta}{3n\_{\mathrm{b}}}.\tag{15}$$

In such cases, even though we have assumed the same equivalent masses for *u* and *d* quarks, their vector potentials take different forms, i.e., contributions to the isovectorvector channel.

Beside the constancy of chemical potentials in Equation (11), the dynamic stability of the quark-vacuum interface should also be fulfilled, which determines the radius of the quark core *R* with

$$P\_{\rm QM}(R) = \frac{2\sigma}{R}.\tag{16}$$

The obtained energy per baryon for *ud*QM nuggets are then presented in Figure 2, where we have adopted the parameters indicated in Table 1 as well as the smallest possible surface tension value *σ* = *σ*min. It is found that the energy per baryon for various cases crosses at *A* = 266 and is generally decreasing with *A*. However, if large symmetry energies with *C<sup>I</sup>* = 40 and 80 MeV/fm<sup>3</sup> were considered, large *ud*QM nuggets are destabilized and there may exists *ud*QM nuggets at *A* ≈ 1000 that are more stable than others. Particularly, there are even cases where the minimum energy per baryon of *ud*QM star matter is larger than 930 MeV but that of √ *ud*QM nuggets smaller than 930 MeV, e.g., *C* = −0.5, *D* = 176 MeV, and *C<sup>I</sup>* = 40 MeV/fm<sup>3</sup> . Such kind of scenarios occur when the surface tension value is smaller than the critical one *σ*crit, which can be roughly estimated with Equation (1) with the corresponding values presented in Table 1. In such cases, the surface of an *ud*QM star will fragment into a lattice of *ud*QM nuggets immersed in a sea of electrons, similar to the cases of strange stars with crusts of strangelets [51,52].

The charge-to-mass ratios of *ud*QM nuggets are indicated in Figure 3, which are decreasing with baryon number due to Coulomb repulsion. Additionally, for the cases with a larger symmetry energy (larger *C<sup>I</sup>* ), *ud*QM nuggets are more positively charged, which gives larger Coulomb energy and leads to a barrier of energy per baryon as we increase *A* for the cases with *σ* < *σ*crit.

**Figure 2.** Energy perbaryon for *ud*QM nuggets as functions of the baryon number *A*. The experimental data for the heaviest *β*-stable nuclei <sup>266</sup>Hs is obtained from the 2016 Atomic Mass Evaluation [75–77].

**Figure 3.** Same as Figure 2 but for the charge-to-mass ratio of *ud*QM nuggets.
