# Adaptive Optimizer for Automated Hyperparameter Optimization Problem

Huayuan Sun<EMAIL_ADDRESS>
Peking University
Beijing, China

###### Abstract

The choices of hyperparameters have critical effects on the performance of
machine learning models. In this paper, we present a general framework that is
able to construct an adaptive optimizer, which automatically adjust the
appropriate algorithm and parameters in the process of optimization. Examining
the method of adaptive optimizer, we product an example of using genetic
algorithm to construct an adaptive optimizer based on Bayesian Optimizer and
compared effectiveness with original optimizer. Especially, It has great
advantages in parallel optimization.

##  I Introduction

The effect of machine learning model largely depends on the selection of super
parameters, which has attracted more and more attention in recent years. The
problem can be formulated as solving the following minimization problem

$\mathop{\arg\min}\limits_{x\in\Omega}f(x),$

where $f$ is a specific metrics function, $\Omega$ is the hyperparameter space
and $x$ is the parameter in $\Omega$.

There are many hyperparameter optimization methods applied in practice. Random
Search and Grid Search are simple and straightforward methods to find optimal
parameters hyperparameter. The advantage is that the algorithm is clear and
suitable for parallel operation. However, they need to try enough parameter
combinations, which is very time-consuming. The other is based on Sequential
Model-Based Global Optimization [2-4], such as Bayesian optimization (BO) [1]
and Tree-structured Parzen Estimator (TPE), which can use relatively few steps
to achieve better results. While, they have poor support for parallel
operation. In addition, different optimization models perform great
differences in different optimization tasks, which may be related to the
machine learning algorithm used, the distribution of training data, the
selection of parameter space and evaluation function. In addition, different
optimization models are suitable in different optimization stages. For
example, we except more exploration in the early stage of optimization, and
later we need more local search, which is similar to the decreasing learning
rate in machine learning. The same is true for different optimizer.

In this paper, we propose a general framework, which can construct an adaptive
optimizer and automatically adjust the appropriate algorithms and parameters
in the optimization process. Our approach allow to automatically adjust the
appropriate algorithm and parameters in the optimization process, and it has
good support for parallel operation. In short, we introduce the variable
$\phi$ into the optimization objective, and change original optimization
problem into the form

$\mathop{\arg\min}F(\phi,N)=\mathop{\arg\min}\limits_{1\leq n\leq
N}f(\phi(\theta_{n},n)),$

where $\phi:(\theta_{n},n)\to x\in\Omega$ represents the optimization function
which suggest parameter $x$ in round $n$. $N$ is the maximum number of
iterations. Thus, the core problem is to determine the parameters
$\theta_{n}$. Among other algorithms, $\phi$ and $\theta$ is a definite
function and parameter.

To show the benefit of our adaptive optimizer compared to fixed optimizer in
the process of optimization, we consider a parameter optimization for lightgbm
by simple Bayesian optimizers and our approach which is adaptive Bayesian
optimizers respectively. Bayesian optimizationis a classical and effective
method to solve the above problems, for the derivative of neural network
hyperparameters cannot be obtained and Bayesian optimization is gradient-free
global optimization. It contains two core components, probabilistic surrogate
model and acquisition function. The probabilistic surrogate model is to
approximatly estimate $f(x)$ based on historical evaluation, commonly
including Gaussian Processes (GP), Eextra Trees (ET), Random Forests (RF),
Gradient Boosted Regression Trees (GBRT).

And acquisition function measures the utility by trading off exploration and
exploitation, such as
(1) Lower Confidence Bound (LCB):

$x_{t}=\mathop{\arg\min}\mu_{t-1}-\beta^{1/2}\sigma_{t-1}$

(2) Expected Improvement (EI):

$x_{t}=\mathop{\arg\min}E_{f(x)\sim
N(\mu_{t-1},\sigma_{t-1}^{2})}\left[\min(f(x)-f_{t-1},0)\right]$

(3) Probability of Improvement (PI):

$x_{t}=\mathop{\arg\min}P(f(x)\leq f_{t}(x))$

In the following sections, we first give a brief introduction to our algorithm
framework in section two. Afterwards, we state our adaptive approach for
hyperparameters optimization problem by genetic algorithm on two
examples111https://github.com/shyjin/Example-of-GA-for-Adaptive-Optimizer .

Algorithm 1 Adaptive Optimizer

1\. Generate initial points2. Calculate rewards $f_{j}^{*}$ of each optimizer
$j$ in the set of meta optimizers $\\{\phi_{j}(\theta_{n},n)\\}$3\. Choose
different $N_{s}$ optimizers based on $f^{*}$ 4\. Suggest next point $x_{i}$
according to the optimizer $\phi(\theta_{i},step_{i})$, and get the value
$f(x_{i})$5\. Repeat step 2 to 4 until termination criterion is reached6.
Output the best record

##  II Adaptive Optimizer

In this section, we introduce our approach of solving a hyperparameters
optimization problem. For convenience, we will not distinguish between
algorithms and parameters in the future. As long as the algorithm or
parameters are not exactly the same, we regard them as different optimizers.

Our method is basically the same as the usual algorithm, except that the step
2 and 3 of adjusting the optimizer is added. We take the adjusted $f_{n}$
returned by the optimizer in step $n$ as a reward to optimize the optimizer.
Then, there is a problem that the effectiveness of the reward function is
different at different iterations. For example, in the 15th round, we use
optimizer A to get a set of parameters with model accuracy of 0.9, which is
completely different from that in the 90th round with optimizer B. In this
example, optimizer A is obviously better than optimizer B for optimizer B gets
more information, but the accuracy is not improved. Then, if optimizer B
achieves an accuracy of 0.92 in the 90th round, how to make a consistent
comparison between the two optimizers is the focus of our discussion.

We defined adjusted reward function as

$f^{*}=\max\left(\epsilon,\frac{f(x)-\frac{1}{n}\sum
f_{k}}{std(f_{1},f_{2},...,f_{n})}-\alpha n\right),$

where $\epsilon>0$ is the minimum weight , $n$ is the number of current
iteration and $alpha=\frac{c}{N}$, $c$ is an non-negative constant like 1.96
or 2.23 and $N$ is the maximum number of iterations. In fact, any non-
decreasing function relative to $n$ is available. Considering that the
marginal benefit of previous iterations is great, non-convex functions are
more suitable than linear function.

Suppose the number of suggestion $N_{s}=1$, which means that we try once on
objective function and get rewards, then we update model. The simplest way to
select the next round of optimizer is to select it on a probabilistic basis.
For any fixed optimizer, its probability of selection is based on the largest
weight it has ever received. If it is not selected, we can set the average
weight of the population. Sum the selection weights of all the optimizers and
calculate the corresponding proportion, which is the probability of selection.
We give an example in Section 3.

For more complex problems, we ask more than one set of parameters and sent to
objective function. Then, running the $N_{s}>1$ set of parameters optimization
tasks in $N_{s}$ multithreads. Intelligent optimization algorithms such as
genetic algorithm can be used. The core approach is to use the adjusted reward
function $f^{*}$ above as fitness function, and cross-over and mutation are
added when selecting the optimizer for the next iteration.

##  III Examples

In order to test the effect of our method, we choose to build two practical
hyperparametric optimization problems for testing. General test and parallel
test are carried out respectively. We choose the most downloaded dataset
”credit-g” on OpenML. This dataset classifies people described by a set of
attributes as good or bad credit risks. And we divide the train dataset and
test dataset according to the proportion of 70% and 30%. The prediction model
is lightgbm, and AUC of the model in the test set is taken as the evaluation
function. AUC is a commonly used evaluation metrics in classification
problems. Five parameters and ranges commonly used in lightgbm algorithm are
selected as the hyperparameter space. The hyperparameters are

  1. 1.

num_leaves: Maximum tree leaves for base learners. The space is [4, 100], int
type.

  2. 2.

min_child_samples: Minimum number of data needed in a child. The space is [1,
100], int type.

  3. 3.

n_estimators: Number of boosted trees to fit. The space is [1, 100], int type.

  4. 4.

subsample: Subsample ratio of the training instance. The space is [0.1, 1],
float type.

  5. 5.

colsample_bytree: Subsample ratio of columns when constructing each tree. The
space is [0.1, 1], float type.

The maximum number of iterations for each optimization is $N$, and the number
of test parameters requested for each iteration is $N_{s}$. We take the best
score achieved in all iterations as the final score of the model. Without
losing generality, as we return the score to the optimizer, we will use
$-f^{*}$ and change the problem into solving the minimum value. In addition,
we repeated the experiment 10 times, and selected different random seeds each
time.

###  III-A Example: Simple Task, $N_{s}=1$

In this experiment, we set the maximum number of iterations to $N=100$ and
return $N_{s}=1$ set of parameter verification results each time. Our base
optimizers are Bayesian optimizers with GP, RF, ET and GBRT surrogate models
respectively. In each round of iteration, we will allocate the weight of the
next round of selection according to the historical scores of the four base
optimizers, and then install a certain proportion to select one of them. The
specific steps have been described in section two.

We use four base optimizers to perform optimization tasks separately, and
other settings are the same. The results are as follows: the average AUC of
the final 10 experiments of GP, RF, ET and GBRT are 0.7613, 0.7591,0.7638 and
0.7505 respectively; The average AUC of our algorithm is 0.7649. Our algorithm
is better than the four base classifiers. It should be emphasized that
although we have only increased by 0.0011 compared with the best performing
base optimizer ET, our advantage is that we don’t need to deliberately select
the optimizer. We only need to put a group of optimizers into the model to
achieve the effect no less than any optimizer. In reality, the optimization
problem is very complex. It is difficult for us to know which optimizers are
suitable for this problem before the actual test. Even at different stages of
optimization, the performance effects of different optimizers are different.
The proposed algorithm can automatically adjust the optimizer in the iterative
process

| GP | RF | ET | GBRT | MyAlgo
---|---|---|---|---|---
seed=0 | 0.7372 | 0.7777 | 0.7758 | 0.7887 | 0.7727
seed=1 | 0.7552 | 0.7679 | 0.7793 | 0.7887 | 0.7812
seed=2 | 0.757 | 0.7815 | 0.7900 | 0.7797 | 0.7614
seed=3 | 0.778 | 0.7782 | 0.7927 | 0.7674 | 0.7676
seed=4 | 0.7805 | 0.7687 | 0.7668 | 0.7865 | 0.7129
seed=5 | 0.7786 | 0.7394 | 0.7917 | 0.7339 | 0.781
seed=6 | 0.7671 | 0.6817 | 0.7168 | 0.6878 | 0.7891
seed=7 | 0.7404 | 0.7684 | 0.7266 | 0.7887 | 0.7371
seed=8 | 0.7715 | 0.7465 | 0.7611 | 0.6465 | 0.7554
seed=9 | 0.7473 | 0.7809 | 0.7366 | 0.7365 | 0.7907
auc_mean | 0.7613 | 0.7591 | 0.7638 | 0.7505 | 0.7649
TABLE I: In this table, the algorithms to be compared with are the base
Optimizer with GP, RF, ET and GBRT surrogate model.

###  III-B Example: Parallel Tasks, $N_{s}>1$

In this experiment, we set the maximum number of iterations to $N=50$ and take
$N_{s}=3$ set of parameter verification results each time. Our base optimizers
are also Bayesian optimizers with GP, RF, ET and GBRT surrogate models
respectively, and the difference is we use four different acquisition function
as well. They are LCB, EI, PI and GP Hedge (probabilistically choose one of
the above three acquisition functions at every iteration), and we have 16 base
optimizers for crossing surrogate models and acquisition function. Another
difference is the method to use to sample multiple points. the constant liar
strategy is used after suggesting one point. With this strategy a copy of
optimizer is created, which is then asked for a point, and the point is told
to the copy of optimizer with some fake objective, the next point is asked
from copy, it is also told to the copy with fake objective and so on. The type
of lie defines different flavours of strategies. We use minimum fake objective
here, which means set fake objective as the minimum value.

Here, we use genetic algorithm as the selection algorithm of the next round of
iterative optimizer. Other parts are the same as before and will not be
repeated. We calculate fitness according to $f^{*}$. Here for base optimizer
$j$, we make a slight adjustment and use

$f^{*}_{j}=\max\left(\epsilon,\frac{f_{j}-\frac{1}{n}\sum
f_{k}}{std(f_{1},f_{2},...,f_{n*N_{s}})}-\alpha n-b\right),$

as fitness, where $b$ is a positive constant and we prefer to focus on the
optimizer which have good fitness if $b$ is larger. The weight of each base
optimizer is the larger of its corresponding fitness ratio and a fixed value.
$k$ is the number of choosing optimizer $j$. In particular, we only select the
historical maximum as the result of one base optimizer, and then select four
optimizers as parents according to the corresponding weight. After retention,
crossover and mutation, take the first three as the next round optimizer.

Algorithm 2 Adaptive Optimizer with Genetic Algorithm

1\. Generate initial points2. Calculate rewards $f_{j}^{*}$ of each optimizer
$j$ in the set of meta optimizers $\\{\phi_{j}(\theta_{n},n)\\}$3\. Generate
$N_{s}$ Optimizers for next iteration 3.1 Choose different $NG$ optimizers
based on $f^{*}$ as parents 3.2 Retention with probability $\epsilon_{r}$ 3.3
Crossover surrogate model and acquisition function 3.4 Mutate with probability
$\epsilon_{m}$ 3.5 Check and replace the duplicate child 4. Suggest next point
$x_{i}$ according to the optimizer $\phi(\theta_{i},step_{i})$, and get the
value $f(x_{i})$5\. Repeat step 2 to 4 until termination criterion is
reached6. Output the best record

To compare our algorithm, we test the 16 base optimizers respectively as well,
combination of four surrogate models and four acquisition functions. We
calculate some indicators to measure the results, including the mean and
median of AUC for the effect of the optimizer, and the variance of AUC for the
stability of the optimizer. In Table 2, the results are shown in order of mean
value from large to small. The mean and median results of our algorithm are
ranked 3rd, and the variance is ranked 4th. The minimum value can be used to
measure the effect of the optimizer in the worst case, ranking 3rd place among
all optimizers.

Generally speaking, compared with all base optimizers, our algorithm is above
average in optimization effect and stability, and can easily run parallel
tasks, that is, $N_{s}$ tasks executed in each round of iteration. For the
base classifier we use, these $N_{s}$ prediction points must be executed in
turn during each round of optimization. Therefore, our algorithm has certain
advantages in parallel efficiency. The effect of the model is affected in the
following aspects:

Optimizer | auc_mean | auc_std | auc_max | auc_min | auc_median
---|---|---|---|---|---
RF_LCB | 0.7918 | 0.0076 | 0.8055 | 0.7789 | 0.7904
ET_PI | 0.7895 | 0.012 | 0.8014 | 0.7582 | 0.7921
myalgo_nan | 0.7891 | 0.0069 | 0.7973 | 0.7772 | 0.7911
ET_LCB | 0.7886 | 0.0137 | 0.8003 | 0.7508 | 0.7915
GP_LCB | 0.7862 | 0.0049 | 0.7919 | 0.7748 | 0.7863
RF_PI | 0.7855 | 0.0084 | 0.7968 | 0.7737 | 0.7859
GP_gp_hedge | 0.7852 | 0.0052 | 0.7944 | 0.7787 | 0.7841
RF_gp_hedge | 0.785 | 0.0128 | 0.8038 | 0.7638 | 0.7872
ET_gp_hedge | 0.7842 | 0.0199 | 0.8038 | 0.7429 | 0.7928
GP_EI | 0.7821 | 0.0074 | 0.7927 | 0.7695 | 0.7849
GP_PI | 0.7814 | 0.0045 | 0.7878 | 0.7752 | 0.7813
GBRT_LCB | 0.7787 | 0.0212 | 0.7978 | 0.7292 | 0.7862
GBRT_gp_hedge | 0.7743 | 0.0326 | 0.8012 | 0.7059 | 0.7934
RF_EI | 0.7648 | 0.0661 | 0.7924 | 0.5778 | 0.788
GBRT_PI | 0.7631 | 0.053 | 0.8055 | 0.6219 | 0.7817
ET_EI | 0.7572 | 0.0394 | 0.791 | 0.6735 | 0.7658
GBRT_EI | 0.7474 | 0.0338 | 0.8055 | 0.7 | 0.7503
TABLE II: In this table, the algorithms to be compared with combination of
four surrogate models and four acquisition functions

  1. 1.

Differences of base optimizers: similar base optimizers will generate similar
recommended parameters in each iteration, which will greatly affect the effect
of the model. Adding different kinds of optimizers will greatly improve.

  2. 2.

Maximum number of iterations N: The selection of optimizer needs the support
of historical data. If N is relatively small, due to the insufficient amount
of data, it will be greatly affected by random shake.

  3. 3.

$N_{s}$: As mentioned in the first point, when $N_{s}$ is too large, it is
easy to make the selected point close in the same round. Enough different base
optimizers need to be added.

##  IV Conclusion

With this paper, we have introduced a new methodology of automated
hyperparameter optimization. Machine learning algorithms are generally time-
consuming, so it is difficult to test which algorithms are more suitable for
hyperparametric optimization tasks. By using our method, we can automatically
adjust the optimizer to adapt to the task in the optimization process.
Moreover, a series of algorithms that do not support parallel tasks can be
combined to establish parallel tasks. Using a single optimizer is prone to
particularly poor performance, and multiple validation takes time much take.
The method is more suitable for difficult and time-consuming optimization
tasks.

## References

  * [1] Pelikan M. Bayesian Optimization Algorithm. Technical Report Illigal, 1999.
  * [2] Bergstra J, Bardenet R , B Kegl, et al. Algorithms for Hyper-Parameter Optimization Advances in Neural Information Processing Systems. 2011.
  * [3] Hutter F, Hoos H H, Leyton-Brown K . Sequential Model-Based Optimization for General Algorithm Configuration. Springer, Berlin, Heidelberg, 2012.
  * [4] Bergstra J , Yamins D L K , Cox D D . Making a Science of Model Search: Hyperparameter Optimization in Hundreds of Dimensions for Vision Architectures. jmlr, 2013.


# How perturbative QCD constrains the Equation of State at Neutron-Star
densities

Oleg Komoltsev  Faculty of Science and Technology, University of Stavanger,
4036 Stavanger, Norway  Aleksi Kurkela  Faculty of Science and Technology,
University of Stavanger, 4036 Stavanger, Norway  Univ Lyon, Univ Claude
Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France

###### Abstract

We demonstrate in a general and analytic way how high-density information
about the equation of state (EoS) of strongly interacting matter obtained
using perturbative Quantum Chromodynamics (pQCD) constrains the same EoS at
densities reachable in physical neutron stars. Our approach is based on
utilizing the full information of the thermodynamic potentials at the high-
density limit together with thermodynamic stability and causality. The results
can be used to propagate the pQCD calculations reliable around $40n_{s}$ to
lower densities in the most conservative way possible. We constrain the EoS
starting from only few times the nuclear saturation density $n\gtrsim
2.2n_{s}$ and at $n=5n_{s}$ we exclude at least 65% of otherwise allowed area
in the $\epsilon-p$ -plane. These purely theoretical results are independent
of astrophysical neutron-star input and hence they can also be used to test
theories of modified gravity and BSM physics in neutron stars.

##  I Introduction

The rapid evolution of neutron-star (NS) astronomy in recent years — in
particular the recent NS radius measurements [1, 2], the discovery of massive
NSs [3, 4, 5], and the advent of gravitational-wave and multi-messenger
astronomy [6, 7] — is for the first time giving us empirical access to the
physics of the cores of NSs. Within the Standard Model and assuming general
relativity, the internal structure of NSs is determined by the equation of
state (EoS) of strongly interacting matter [8, 9]. With these assumptions, NS
observations can be used to empirically determine the EoS [10, 11, 12, 13, 14,
15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30] (for reviews,
see [31, 32, 33, 34, 35, 36, 37]). And if the EoS can be determined
theoretically to a sufficient accuracy, comparison with NS observations allows
to use these extreme objects as laboratory for physics beyond standard model
(see, e.g., [38, 39, 40, 41, 42, 43, 44, 45, 46]) and/or general relativity
(e.g., [47, 48, 49, 50, 51]).

For both of these goals, it is crucial to make use of all possible controlled
theory calculations that inform us about the EoS at densities reached in NSs.
While in principle the EoS is determined by the underlying theory of strong
interactions, Quantum Chromodynamics (QCD), in practice we have access to the
EoS only in limiting cases. In the context of low temperatures relevant for
neutron stars, the EoS of QCD can be systematically approximated at low- and
at high-density limits. At low densities, the current state-of-the-art low-
energy effective theory calculations allow to describe matter to densities
around and slightly above nuclear saturation density $n\approx
n_{s}=0.16/\textrm{fm}^{3}$ [52, 53], but become unreliable at higher
densities $n\sim 5-10n_{s}$ reached in the cores of massive neutron stars. A
complementary description of NS matter comes from perturbative QCD (pQCD)
calculations which become reliable at sufficiently high densities $\sim
40n_{s}$, far exceeding those realised in NSs [54, 55].

Figure 1:  Theoretical constraints to $\epsilon$ and $p$ arising from low-
(CET) and high-density input (pQCD). The green lines show the envelope of
allowed values of $p(\epsilon)$; dashed green lines arise from trivially
imposing causality of the EoS while the thick solid lines correspond to newly
identified integral constraints arising from imposing simultaneously the high-
density limit for $p$, $n$, and ${\mu}$. The blue regions correspond to
allowed values $\epsilon$ and $p$ at fixed density $n$; the red regions would
be otherwise allowed but are excluded by the high-density input. At
$n=5n_{s}$, $75\%$ of otherwise allowed values of $\epsilon$ and $p$ are
excluded by the high-density input.

Several works have used large ensembles of parameterized EoSs to study the
possible behavior of the EoS in intermediate densities between the
theoretically known limits. Many of these works have anchored their EoS to the
low-density limit only, while others have interpolated between the two orders
of magnitude in density separating the low- and high-density limits [56, 10,
16, 27, 57]. While the input from pQCD clearly constrains the EoS at very high
densities, how this information affects the EoS around neutron-star densities
has so far been convoluted by the specific choices of interpolation functions.

The aim of the present work is to make the influence of the high-density
calculations to the EoS at intermediate densities explicit, and in particular,
to derive constraints to the EoS that are completely independent of any
specific interpolation function. By using the full information available in
the thermodynamic grand canonical potential as a function of baryon chemical
potential, $\Omega({\mu})=-p({\mu})$, we find stricter bounds than works using
only the reduced form of the EoS, i.e., the pressure as a function of energy
density $p(\epsilon)$ appearing in the hydrodynamic description of neutron
star matter [58, 59, 51, 60, 61]. The effect is demonstrated in fig. 1, which
shows the region in $\epsilon-p$ -plane that can be reached with a causal and
thermodynamically stable EoS with information of the high-density limit, and
in particular compares the allowed values at different fixed densities
$n=2,3,5,$ and $10n_{s}$ with and without the pQCD input.

##  II setup

In the following we consider all possible interpolations of the full
thermodynamic potential at zero temperature and in $\beta$-equilibrium,
$\Omega({\mu})=-p({\mu})$, between the low-density limit ${\mu}=\mu_{\rm L}$
and the high-density limit ${\mu}=\mu_{\rm H}$. We assume that at both of
these limits the pressure $p({\mu})$ and its first derivative, baryon number
density $n({\mu})=\partial_{{\mu}}p({\mu})$ are know

$\displaystyle p(\mu_{\rm L})$ $\displaystyle=p_{\rm L},\quad p(\mu_{\rm
H})=p_{\rm H},$ (1) $\displaystyle n(\mu_{\rm L})$ $\displaystyle=n_{\rm
L},\quad n(\mu_{\rm H})=n_{\rm H}.$ (2)

This information is readily available from the microscopic calculations which
are assumed to be reliable at these limits; representative values from state-
of-the-art chiral effective theory (CET) [62] and pQCD calculations [55] from
the literature are reproduced in table 1.

Thermodynamic consistency requires that pressure is a continuous function of
${\mu}$. Similarly, thermodynamic stability requires convexity of
thermodynamic potentials $\partial_{\mu}^{2}\Omega\leq 0$ so that $n({\mu})$
is a monotonically increasing function $\partial_{\mu}n({\mu})\geq 0$. Density
$n({\mu})$ does not need to be continuous and it can have discontinuities
(increasing the density at fixed ${\mu}$) in the case of 1st-order phase
transitions; we place no additional assumptions on the number or strength of
possible transitions. Causality requires that the speed of sound is less than
the speed of light, $c_{s}^{2}\leq 1$, and it imposes a condition on the first
derivative of the number density

$\displaystyle c_{s}^{-2}=\frac{\mu}{n}\frac{\partial n}{\partial{\mu}}\geq
1.$ (3)

| CET | pQCD
---|---|---
| soft | stiff | $X=1$ | $X=2$ | $X=4$
$\mu$ [GeV] | 0.966 | 0.978 | 2.6
$n$ [1/fm3] | 0.176 | 6.14 | 6.47 | 6.87
$p$ [MeV/fm3] | 2.163 | 3.542 | 2334. | 3823. | 4284.
Table 1: Collection of predictions for the thermodynamic quantities at the
low- (CET) and the high- (pQCD) density limits. CET limit corresponds to the
”stiff” and ”soft” EoS of [62] while the pQCD values are from a partial N3LO
calculation [55] at three different values of the renormalization scale
parameter $X$. The uncertainty of the pQCD limit can be assessed by varying
$X$ in the range $X=[1,4]$. All figures correspend to the ”stiff” CET and
$X=2$ pQCD.

We construct the allowed EoSs by considering all possible functions $n({\mu})$
allowed by the above assumptions connecting the low- and high-density limits.
Causality imposes a minimal allowed slope of the number density, i.e.,
$\partial_{\mu}n({\mu})\geq n/{\mu}$, that any causal EoS passing though a
given point in the ${\mu}-n$ plane can have. This can be visualized as a
vector field on the ${\mu}-n$ -plane as depicted in fig. 2, where the arrows
at each point corresponds to tangent lines with constant $c^{2}_{s}=1$. This
requirement imposes two fundamental constraints on the ${\mu}-n$ plane.
Starting from the point $\\{\mu_{\rm L},n_{\rm L}\\}$ we can follow the arrows
until $\mu_{\rm H}$ by solving eq. 3 with $c^{2}_{s}=1$, leading to
$n({\mu})=n_{\rm L}{\mu}/\mu_{\rm L}$. This produces a maximally stiff causal
EoS and the area under this line cannot be reached from the low-density limit
with a causal EoS. Correspondingly, upper limit for the $n({\mu})$ can be
obtained starting from $\\{\mu_{\rm H},n_{\rm H}\\}$ and following the arrows
backward to $\mu_{\rm L}$; the high-density limit $n_{\rm H}$ cannot be
reached from any point above this line by a causal EoS. These previously known
bounds (e.g. [58, 51]) are represented as orange lines in fig. 2.

Figure 2: Density as a function of chemical potential determines the EoS. At
any given point $\\{{\mu},n\\}$, causality imposes a minimal slope to any EoS
passing that point $\partial_{\mu}n(\mu)\geq\frac{n}{\mu}$, denoted by the
arrows. Thermodynamic consistency demands that $n(\mu)$ is a monotonic
function, and that the area under the curve in the plot is given by $\Delta
p$. These conditions cannot be simultaneously fulfilled by any EoS that passes
to the red area at any density and are therefore excluded. The lines $\Delta
p_{\rm min}$ and $\Delta p_{\rm max}$ display the constructions defined in eq.
7 and eq. 14.

The simultaneous requirement of reaching both $p_{\rm H}$ and $n_{\rm H}$
imposes further constraints and fixes the area under the curve $n({\mu})$.

$\displaystyle\int_{\mu_{\rm L}}^{\mu_{\rm H}}n(\mu)d\mu=p_{\rm H}-p_{\rm
L}=\Delta p.$ (4)

In the example shown in fig. 2, this requirement imposes that the area under
any allowed EoS is approximately one third of the area of the figure.

At each point we can evaluate the absolute minimum and maximum area under any
EoS ($\Delta p_{\rm min/max}$) that can be reached at $\mu_{\rm H}$ if the EoS
goes through that particular point. If $\Delta p_{\rm min}$ ($\Delta p_{\rm
max}$) is bigger (smaller) than $\Delta p$ then such point would be ruled out.

To obtain the minimum area at $\mu_{\rm H}$ for any EoS going through a
specific point $\\{\mu_{0},n_{0}\\}$ consider the following construction shown
in fig. 2 as a dashed blue line

$\displaystyle n({\mu})=\left\\{\begin{array}[]{rl}n_{\rm L}{\mu}/\mu_{\rm
L},&\mu_{\rm L}<\mu<\mu_{0}\\\ n_{0}{\mu}/\mu_{0},&\mu_{0}<\mu<\mu_{\rm
H}.\end{array}\right.$ (7)

For $\mu<\mu_{0}$, the smallest possible area is determined by the maximally
stiff causal line (i.e. $c_{s}^{2}=1$) starting from $\\{\mu_{\rm L},n_{\rm
L}\\}$ which we follow up to $\mu_{0}$. At $\mu_{0}$, we have a phase
transition where the density jumps to $n_{0}$. After that, the EoS follows the
maximally stiff causal line starting at $\\{n_{0},\mu_{0}\\}$ until $\mu_{\rm
H}$, where the EoS has another phase transition to reach $n(\mu_{\rm
H})=n_{\rm H}$. The solution to the equation $\Delta p_{\rm
min}(\mu_{0},n_{0})=\int_{\mu_{\rm L}}^{\mu_{\rm H}}n({\mu})=\Delta p$ is
shown in fig. 2 as the top red line denoted integral constraints; any EoS
crossing this line is inconsistent with simultaneous constraint on $p_{\rm H}$
and $n_{\rm H}$. This yields a maximum density for given ${\mu}$

$\displaystyle n_{\rm
max}({\mu})=\left\\{\begin{array}[]{ll}\frac{{\mu}^{3}n_{\rm L}-{\mu}\mu_{\rm
L}(\mu_{\rm L}n_{\rm L}+2\Delta p)}{({\mu}^{2}-\mu_{\rm H}^{2})\mu_{\rm
L}},&\mu_{\rm L}\leq{\mu}<{\mu_{\rm c}}\\\ n_{\rm H}{\mu}/\mu_{\rm
H},&{\mu_{\rm c}}\leq{\mu}\leq\mu_{\rm H},\end{array}\right.$ (10)

where ${\mu_{\rm c}}$ is given by the intercept of the causal line and the
integral constraint, i.e., the two cases in eq. 10

${\mu_{\rm c}}=\sqrt{\frac{\mu_{\rm L}\mu_{\rm H}(\mu_{\rm H}n_{\rm
H}-\mu_{\rm L}n_{\rm L}-2\Delta p)}{\mu_{\rm L}n_{\rm H}-\mu_{\rm H}n_{\rm
L}}}.$ (11)

Similarly, the procedure to maximize area under any EoS going through the
point $\\{n_{0},\mu_{0}\\}$ is shown as a dashed green line in fig. 2,

$\displaystyle
n({\mu})=\left\\{\begin{array}[]{rl}n_{0}{\mu}/\mu_{0},&\mu_{\rm
L}<{\mu}<\mu_{0}\\\ n_{\rm H}{\mu}/\mu_{\rm H},&\mu_{0}<{\mu}<\mu_{\rm
H}.\end{array}\right.$ (14)

Correspondingly, solving for $\Delta p_{\rm max}=\Delta p$ gives a constraint
to the minimal $n$ that can be obtained for a given ${\mu}$, depicted as the
bottom red line in fig. 2. Then, for a given chemical potential the minimal
allowed density is

$\displaystyle n_{\rm min}({\mu})=\left\\{\begin{array}[]{ll}n_{\rm
L}{\mu}/\mu_{\rm L},&\mu_{\rm L}\leq{\mu}\leq{\mu_{\rm c}}\\\
\frac{{\mu}^{3}n_{\rm H}-{\mu}\mu_{\rm H}(\mu_{\rm H}n_{\rm H}-2\Delta
p)}{({\mu}^{2}-\mu_{\rm L}^{2})\mu_{\rm H}},&{\mu_{\rm c}}<{\mu}\leq\mu_{\rm
H}.\end{array}\right.$ (17)

Note that the maximally stiff causal lines intercept with the integral
constraints at the same point ${\mu_{\rm c}}$ for lower and upper limits. This
happens because the EoS following $n_{\rm min}({\mu})$ up to ${\mu_{\rm c}}$
obtains the correct area $\Delta p$ only if the EoS jumps at ${\mu_{\rm c}}$
from $n_{\rm min}({\mu_{\rm c}})$ to $n_{\rm max}({\mu_{\rm c}})$. This is
also the EoS exhibiting a phase transition with the largest possible latent
heat

$\displaystyle Q=\mu_{c}(n_{\rm max}({\mu_{\rm c}})-n_{\rm min}({\mu_{\rm
c}}))=$ $\displaystyle-\mu_{\rm L}n_{\rm L}+\mu_{\rm H}n_{\rm H}+2p_{\rm
L}-2p_{\rm H}.$ (18)

We further note that there is one EoS that has a constant $c_{s}^{2}=1$
throughout the whole region, connecting the $n_{\rm max}(\mu_{\rm L})$ and
$n_{\rm min}(\mu_{\rm H})$. We denote this special EoS as

$n_{c}({\mu})=n_{\rm max}(\mu_{\rm L}){\mu}/\mu_{\rm L}=n_{\rm min}(\mu_{\rm
H}){\mu}/\mu_{\rm H}.$ (19)

Any point along this line maximises the area to the left of the point
(${\mu}<\mu_{0}$) and minimizes to the right of the point (${\mu}>\mu_{0}$),
so that this line corresponds to the maximal pressure at fixed ${\mu}$.

##  III Mapping to the $\epsilon-p$ plane

Figure 3: 3D rendering of the constraints in the ${\mu}-n-p$ phase space. Each
triangle is a slice of the $p-n$ constraints for fixed ${\mu}$. Note that
$n_{max},n_{min}$ and $n_{c}$, which are defined in eq. 10, eq. 17, and eq.
19, respectively, are projection on the $\mu-n$ -plane of the lines denoted
with the same names in the figure above.

For every allowed point on the $\mu-n$ -plane, we can find minimal and maximal
pressures, $p_{\rm min/max}(\mu_{0},n_{0})$, that can be obtained at that
point $\\{\mu_{0},n_{0}\\}$. Note that this is different from the minimal and
maximal pressures ($\Delta p_{\rm min/max}$) obtainable at $\mu_{\rm H}$ by an
EoS passing though $\\{\mu_{0},n_{0}\\}$ as discussed in the previous section.

The minimal pressure is given by the EoS that follows the maximally stiff
causal line terminating at the point $\\{\mu_{0},n_{\rm min}(\mu_{0})\\}$,
i.e., $n({\mu})=\frac{{\mu}}{\mu_{0}}n_{\rm min}(\mu_{0})$, and has a phase
transition from $n_{\min}$ to $n_{0}$ at $\mu_{0}$. This construction leads to
a lower bound of the pressure as function of ${\mu}$ and $n$

$\displaystyle p_{\rm min}(\mu_{0},n_{0})$ $\displaystyle=p_{\rm
L}+\frac{\mu_{0}^{2}-\mu_{\rm L}^{2}}{2\mu_{0}}n_{\rm min}(\mu_{0})$ (20)

In order to find the maximal pressure for a given point, the $\mu-n$ plane
needs to be divided in two different regions. For $n<n_{c}({\mu})$ (see eq.
19), the maximal pressure is obtained by following the maximally stiff causal
line terminating at the point itself $n({\mu})=n_{0}{\mu}/\mu_{0}$

$\displaystyle p_{\rm max}(\mu_{0},n_{0})$ $\displaystyle=p_{\rm
L}+\frac{\mu_{0}^{2}-\mu_{\rm L}^{2}}{2\mu_{0}}n_{0},\quad n<n_{c}({\mu}).$
(21)

For $n>n_{c}({\mu})$, the above construction would lead to a pressure that is
inconsistent with the integral constraint. Instead, a bound can be obtained by
noting that the EoS that maximizes the pressure at $\\{\mu_{0},n_{0}\\}$ is an
EoS that minimizes the pressure difference between $\mu_{0}$ and $\mu_{\rm
H}$. Thus the maximum pressure in that point is given by the difference
between $\Delta p$ and the pressure for maximally stiff EoS following causal
line starting at $\\{\mu_{0},n_{0}\\}$

$\displaystyle p_{\rm max}(\mu_{0},n_{0})$ $\displaystyle=p_{\rm
H}-\frac{\mu_{\rm H}^{2}-\mu_{0}^{2}}{2\mu_{0}}n_{0},\quad n>n_{c}({\mu}).$
(22)

The simultaneous bounds for ${\mu}$, $n$ and $p$ are visualized in fig. 3.
These constraints can be easily translated into bounds on the $\epsilon-p$
-plane using the Euler equation $\epsilon=-p+{\mu}n$ for a fixed density $n$.
To obtain the envelope of the allowed values of $p(\epsilon)$ irrespective of
$n$ (green lines in fig. 1 and fig. 4) we consider lines of constant enthalpy
$h=\mu n=\epsilon+p$. Finding the maximum and minimum pressure along these
lines maps the bounds from $\mu-n$ to $\epsilon-p$ -plane; for details, see
Appendix A.

Figure 1 shows the allowed range of $\epsilon$ and $p$ values with and without
imposing the high-density constraint for fixed densities of $n=2,3,5,$ and
$10n_{s}$ using the $X=2$ values for pQCD and the ”stiff” values for CET. At
the lowest density of $n=2n_{s}$ the high-density input does not offer
additional constraint. However, strikingly, already at $n=3n_{s}$ the largest
pressures allowed by causality alone are cut off by the integral constraint.
At this density only 68$\%$ of the $\epsilon-p$ values that would be allowed
without the high-density constraint remain (given by the ratio of the blue and
combined blue and red areas in fig. 1 when plotted in linear scale, also see
fig. A1). At $n=5n_{s}$ only $25\%$ of the allowed region remains after
imposing the integral constraint. At $n=10n_{s}$, the allowed range of values
is significantly reduced, now also featuring a cut of the lowest pressures
leaving around 6.5$\%$ of the total area.

We have checked the stability of these results against variation of different
pQCD and CET limits in table 1. While varying the CET parameters has a very
small effect on the excluded areas (of the order of line width in fig. 1),
varying the pQCD renormalization parameter increases the area in $\epsilon-p$
-plane. The union of allowed areas in range $X=[1,4]$ excludes
$[13\%,64\%,92.8\%]$ of the otherwise allowed area (without pQCD) for
$n=[3,5,10]\ n_{s}$. For comparison the same numbers for fixed $X=2$ are
$[32\%,75\%,93.5\%]$, see fig. A1 in the Appendix. The lowest density at which
$p$ and $\epsilon$ are limited by the high-density input is $n_{\rm
min}(\mu_{\rm L})=[2.2,2.5]\,n_{s}$ for $X=[1,4]$. The lowest density where
the pressure is limited from below is $n_{\rm max}(\mu_{\rm
L})=[4.8,9.0]\,n_{s}$ for the same values of $X$.

##  IV Speed-of-sound constraints

The EoSs that render the boundaries of the allowed regions in fig. 1 and fig.
2 are composed of the most extreme ones allowed by the above conditions. They
contain maximal allowed phase transitions and extended density ranges where
the speed of sound coincides with the speed of light. While it is clear that
these most extreme EoSs are unlikely to be the physical one, they cannot be
excluded with the same level of robustness as those which break the above
criteria.

A possible way to quantify how extreme the EoSs are is by imposing a maximal
speed of sound $c_{s}^{2}<{c_{s,{\rm lim}}^{2}}<1$ that is reached at any
density within the interpolation region [27]. In Appendix B we give the
generalizations for $n_{\rm min/max}$, $p_{\rm min/max}$, and ${\mu_{\rm c}}$
for arbitrary limiting ${c_{s,{\rm lim}}^{2}}$. The analysis is similar to
that of the above with the modification of maximally stiff causal lines to
lines of constant speed of sound $n({\mu})=n_{0}({\mu}/\mu_{0})^{1/{c_{s,{\rm
lim}}^{2}}}$.

Figure 4: Constraints to $p(\epsilon)$ for different limiting values of the
sound speed. The green envelope corresponds to $c_{s}^{2}=1$ while the blue
envelope corresponds to the conformal limit $c_{s}^{2}=1/3$. At the specific
value of the sound speed $(c^{2}_{s})^{\rm min}_{\rm lim}$ the allowed region
degenerates in to a single line. The red areas are not excluded if trivially
interpolating the reduced form of the EoS $p(\epsilon)$.

As the maximum speed of sound is reduced, the allowed range of possible
$n({\mu})$ diminishes rapidly leading to tighter bounds on the EoS. This is
demonstrated in fig. 4, showing the range of allowed values of $p(\epsilon)$
at different values of ${c_{s,{\rm lim}}^{2}}$. Interestingly, we note that
while the so-called conformal bound — that is, ${c_{s,{\rm lim}}^{2}}=1/3$ —
is consistent with the high-density constraint, imposing this condition to the
speed of sound gives an extremely strong constraint to the EoS and forces the
EoS to have a specific shape at all densities.

Decreasing ${c_{s,{\rm lim}}^{2}}$ even further, we eventually close the gap
between the lower and upper bounds completely, shown in fig. 4 as red line. At
this specific value of ${c_{s,{\rm lim}}^{2}}=(c^{2}_{s})^{\rm min}_{\rm lim}$
the allowed region degenerates into a single line. This is the minimum
limiting value for sound speed for which the low- and high- density limits can
be connected in a consistent way. Thus $c_{s}^{2}>(c^{2}_{s})^{\rm min}_{\rm
lim}$ has to be reached by any EoS at some density. For different values of
$X=[1,2,4]$, we find the minimal limiting values to be just below the
conformal limit [0.32,0.30, 0.32 ], see Appendix B. It is curious to note that
these values are very close to their pQCD values approaching the conformal
limit from below even though the pQCD values did not enter the calculation of
$(c^{2}_{s})^{\rm min}_{\rm lim}$.

##  V Discussions

Our results demonstrate the non-trivial and robust constraints on the EoS at
NS densities arising from pQCD calculations. These constraints are obtained
only when interpolating the full EoS $p({\mu})$ instead of its reduced form
$p(\epsilon)$.

The results highlight the complementarity of the high- and low-density
calculations. This is further demonstrated in fig. 5 depicting the expected
effect of future calculations on the allowed region of $\epsilon-p$ -values at
$5n_{s}$ and $10n_{s}$. We anticipate future results by extrapolating either
CET (”stiff”) up to $n=2n_{s}$ using values tabulated in [62] and/or
extrapolating the pQCD results ($X=2$) of [55] down to $n=20n_{s}$
(corresponding to ${\mu}=2.1$ GeV). Improvements of this magnitude have been
suggested, e.g., in [63, 64]. We see that both the low- and the high-density
calculations have a potential to significantly constrain the EoS in the future
and that the greatest benefit is obtained by combining the both approaches.
Lastly, we note that merely a $7\%$ reduction of the pQCD pressure for $X=1$
(or using $X=0.95$) would necessitate the violation of the conformal bound.
This implies that an improved pQCD calculations could be used in the future to
cleanly rule out all EoSs that remain subconformal at all densities. These
results highlight the importance of pursuing both the low- and high-density
calculations in the future.

Figure 5: The expected impact of future calculations on the theoretical
constraints to $\epsilon$ and $p$ at NS densities obtained by extrapolating
the CET and pQCD calculations to intermediate densities. Numbers in square
bracket show the end point of extrapolation of the CET and pQCD calculations.
Blue areas denoted as $[1.1n_{s},40n_{s}]$ correspond to blue area in fig. 1.
Green areas are extrapolation of CET limit up to $2n_{s}$ [62], while purple
regions are extrapolations of pQCD down to $20n_{s}$ [55]. The combined effect
of extrapolating both limits is given by the overlapping region.

##  VI Acknowledgments

We thank Eemeli Annala, Tyler Gorda, Tore Kleppe, Kai Hebeler, Achim Schwenk,
and Aleksi Vuorinen for useful discussions. We also thank Sanjay Reddy for
posing questions on multiple occasions about the role of high-density
constraints that in part motivated this paper.

##  Appendix A Boundaries on the $\epsilon-p$ -plane

Utilizing the Euler equation $\epsilon=-p+{\mu}n,$ one can map bounds from
${\mu}-n$ -plane to the corresponding limits on the $\epsilon-p$ -plane. In
order to find the extremal allowed values of $p(\epsilon)$, consider lines of
fixed enthalpy $h=\epsilon+p={\mu}n$. These lines correspond on one hand to
diagonal lines $p(\epsilon)=-\epsilon+h$ on the $\epsilon-p$ -plane, and on
the other hand to hyperbolas $n({\mu})=h/{\mu}$ on the ${\mu}-n$ -plane.
Therefore, minimising or maximizing the pressure for a constant $h$ on the
${\mu}-n$ -plane gives the minimal and maximal pressure on the corresponding
isenthalpic line on the $\epsilon-p$ -plane.

Substituting $n=h/{\mu}$ in eq. 20 we readily observe that the smallest
minimal pressure along isenthalpic line, $p_{\rm min}({\mu},h/{\mu})$, is
obtained at smallest value ${\mu}$ allowed by the constraints, that is, at the
crossing of the isenthalpic line and $n_{\rm max}({\mu})$.

Similarly, substituting $n=h/{\mu}$ in eq. 21 and eq. 22 we observe that the
maximal pressure $p_{\rm max}({\mu},h/{\mu})$ obtains its largest value for
the largest value ${\mu}$ allowed, at the crossing of the isenthalpic line and
$n_{\rm min}({\mu})$.

Therefore the allowed range of values in the $\epsilon-p$ -plane is bounded
from below by the line $\\{\epsilon_{\rm max}(\mu),p_{\rm min}({\mu},n_{\rm
max}({\mu}))\\}$ with $\mu_{\rm L}<{\mu}<\mu_{\rm H}$, where

$\displaystyle\epsilon_{\rm max}({\mu})=-p_{\rm min}({\mu},n_{\rm
max}({\mu}))+{\mu}n_{\rm max}({\mu})$ (23) $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{cc}\frac{\left({\mu}^{2}+\mu_{\rm
H}^{2}\right)\left({\mu}^{2}n_{\rm L}+\mu_{\rm L}(2p_{\rm L}-\mu_{\rm L}n_{\rm
L})\right)-4{\mu}^{2}\mu_{\rm L}p_{\rm H}}{2\mu_{\rm L}({\mu}-\mu_{\rm
H})({\mu}+\mu_{\rm H})},&{\mu}<{\mu_{\rm c}}\\\ \frac{1}{2}(({\mu}^{2}n_{\rm
H})/\mu_{\rm H}+\mu_{\rm H}n_{\rm H}-2p_{\rm H}),&{\mu}>{\mu_{\rm
c}}.\end{array}\right.$ (26)

This corresponds to the lower bound in fig. 1 and fig. 4. In fig. 3, the line
$\\{{\mu},n_{\rm max}({\mu}),p_{\rm min}({\mu},n_{\rm max}({\mu}))\\}$
corresponds to the dashed green line marked $n_{\rm max}$.

Similarly the allowed $\epsilon-p$ values are bounded by above by
$\\{\epsilon_{\rm min}({\mu}),p_{\rm max}({\mu},n_{\rm min}(\mu_{B}))\\}$,
with

$\displaystyle\epsilon_{\rm min}({\mu})=-p_{\rm max}({\mu},n_{\rm
min}({\mu}))+{\mu}n_{\rm min}({\mu})$ (27) $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{cc}\frac{1}{2}(({\mu}^{2}n_{\rm
L})/\mu_{\rm L}+\mu_{\rm L}n_{\rm L}-2p_{\rm L})&{\mu}<{\mu_{\rm c}}\\\
\frac{\frac{{\mu}^{4}}{\mu_{\rm L}^{2}}\frac{n_{\rm H}}{\mu_{\rm
H}}+(\frac{{\mu}}{\mu_{\rm L}})^{2}(\mu_{\rm L}^{2}\frac{n_{\rm H}}{\mu_{\rm
H}}-\mu_{\rm H}n_{\rm H}+4p_{\rm L}-2p_{\rm H})+2p_{\rm H}-\mu_{\rm H}n_{\rm
H}}{2((\frac{{\mu}}{\mu_{\rm L}})^{2}-1)},&{\mu}>{\mu_{\rm
c}}.\end{array}\right.$ (30)

This corresponds to the upper bound of pressure in fig. 1 and fig. 4. In fig.
3, the line $\\{{\mu},n_{\rm min}({\mu}),p_{\rm max}({\mu},n_{\rm
min}({\mu}))\\}$ corresponds to the dashed green line marked $n_{\rm min}$.

Additionally, in fig. A1 we present the impact of the pQCD renormalization
parameter variation on the allowed area for the fixed number density $5n_{s}$
and $10n_{s}$. It can be readily seen from the figure that for $5n_{s}$ pQCD
input excludes 65$\%$ (for $X=1$) of the otherwise allowed area, and for $X=2$
only 25$\%$ of the $\epsilon-p$ values is allowed. The expressions for the
areas without pQCD constraints are easily obtained by sending ${\mu_{\rm
c}}\rightarrow\infty$ and $n_{c}\rightarrow\infty$ and so taking the
expressions corresponding to $\mu<{\mu_{\rm c}}$ and $n<n_{c}$.

Figure A1: The impact of the pQCD renormalization parameter variation on the
theoretical constraints to $\epsilon$ and $p$ at NS densities. Blue area
denoted as $X=2$ corresponds to blue area in fig. 1. Purple and green areas
correspond to allowed area for renormalization parameter $X=1$ and $X=4$,
respectively. Thin black line corresponds to allowed area without high-density
input.

##  Appendix B Boundaries with arbitrary ${c_{s,{\rm lim}}^{2}}$

In this appendix we give generalizations of the bounds discussed in the main
text for limiting speed of sound ${c_{s,{\rm lim}}^{2}}=1$ to arbitrary
${c_{s,{\rm lim}}^{2}}$.

The line with the constant speed of sound is given by

$n({\mu})=n_{0}({\mu}/\mu_{0})^{1/c_{s}^{2}}.$ (31)

The assumption that no EoS is allowed to be stiffer than this with the
additionally imposed integral constraints leads to the minimal and maximal
densities at fixed $\mu$

$n_{\rm min}=\left\\{\begin{array}[]{ll}n_{\rm L}\left(\frac{{\mu}}{\mu_{\rm
L}}\right)^{1/c_{s}^{2}},\mu<{\mu_{\rm c}}\\\
\frac{\left(\frac{{\mu}}{\mu_{\rm
L}}\right)^{1/c_{s}^{2}}\left(c_{s}^{2}\left({\mu}n_{\rm
H}\left(\frac{{\mu}}{\mu_{\rm H}}\right)^{1/c_{s}^{2}}-\mu_{\rm H}n_{\rm
H}+\Delta p\right)+\Delta
p\right)}{c_{s}^{2}\left({\mu}\left(\frac{{\mu}}{\mu_{\rm
L}}\right)^{1/c_{s}^{2}}-\mu_{\rm L}\right)},\mu>{\mu_{\rm
c}}\end{array}\right.$ (32)

and

$n_{\rm max}=\left\\{\begin{array}[]{ll}\frac{\left(\frac{{\mu}}{\mu_{\rm
H}}\right)^{1/c_{s}^{2}}\left(c_{s}^{2}\left({\mu}n_{\rm
L}\left(\frac{{\mu}}{\mu_{\rm L}}\right)^{1/c_{s}^{2}}-\mu_{\rm L}n_{\rm
L}-\Delta p\right)-\Delta
p\right)}{c_{s}^{2}\left({\mu}\left(\frac{{\mu}}{\mu_{\rm
H}}\right)^{1/c_{s}^{2}}-\mu_{\rm H}\right)},\mu<{\mu_{\rm c}}\\\ n_{\rm
H}\left(\frac{{\mu}}{\mu_{\rm H}}\right)^{1/c_{s}^{2}},\mu>{\mu_{\rm
c}}\end{array}\right.$ (33)

where

${\mu_{\rm c}}=\left(\frac{\mu_{\rm H}^{1/c_{s}^{2}}(c_{s}^{2}(\mu_{\rm
L}n_{\rm L}-\mu_{\rm H}n_{\rm H}+\Delta p)+\Delta p)}{c_{s}^{2}\left(n_{\rm
L}\left(\frac{\mu_{\rm H}}{\mu_{\rm L}}\right)^{1/c_{s}^{2}}-n_{\rm
H}\right)}\right)^{\frac{c_{s}^{2}}{c_{s}^{2}+1}}.$ (34)

For any given $\mu$ the largest allowed latent heat associated to a phase
transition is given by

$\displaystyle Q({\mu})=\mu(n_{\rm max}({\mu})-n_{\rm min}({\mu}))$ (35)
$\displaystyle=\frac{{\mu}\left(\frac{{\mu}}{\mu_{\rm L}\mu_{\rm
H}}\right)^{1/c_{s}^{2}}\left(c_{s}^{2}n_{\rm L}\mu_{\rm
H}^{\frac{1}{c_{s}^{2}}+1}-\mu_{\rm L}^{1/c_{s}^{2}}(c_{s}^{2}(\mu_{\rm
L}n_{\rm L}+\Delta p)+\Delta
p)\right)}{c_{s}^{2}\left({\mu}\left(\frac{{\mu}}{\mu_{\rm
H}}\right)^{1/c_{s}^{2}}-\mu_{\rm H}\right)}$

for $\mu<{\mu_{\rm c}}$. For $\mu>{\mu_{\rm c}}$ the expression is obtained by
changing $H\leftrightarrow L$ in the above expression and multiplying the
whole expression by -1. The phase transition with largest allowed latent heat
takes place at $\mu={\mu_{\rm c}}$ and has the latent heat

$Q=\epsilon_{\rm H}-\epsilon_{\rm L}-\frac{p_{\rm H}-p_{\rm L}}{c_{s}^{2}}.$
(36)

The line with largest pressure at fixed ${\mu}$ is given by

$n_{c}({\mu})=n_{\rm max}(\mu_{\rm L})\left(\frac{{\mu}}{\mu_{\rm
L}}\right)^{1/c_{s}^{2}}.$ (37)

The minimal and maximal pressures at given point $\\{\mu_{0},n_{0}\\}$
generalize simply to

$p_{\rm min}(\mu_{0},n_{0})=p_{\rm
L}+\frac{c_{s}^{2}}{1+c_{s}^{2}}\left[\mu_{0}-\mu_{\rm L}\left(\frac{\mu_{\rm
L}}{\mu_{0}}\right)^{1/c_{s}^{2}}\right]n_{\rm min}(\mu_{0})$ (38)

and for $n_{0}<n_{c}(\mu_{0})$

$p_{\rm max}(\mu_{0},n_{0})=p_{\rm
L}+\frac{c_{s}^{2}}{1+c_{s}^{2}}\left[\mu_{0}-\mu_{\rm L}\left(\frac{\mu_{\rm
L}}{\mu_{0}}\right)^{1/c_{s}^{2}}\right]n_{0}$ (39)

and for $n_{0}>n_{c}(\mu_{0})$

$p_{\rm max}(\mu_{0},n_{0})=p_{\rm
H}+\frac{c_{s}^{2}}{1+c_{s}^{2}}\left[\mu_{0}-\mu_{\rm H}\left(\frac{\mu_{\rm
H}}{\mu_{0}}\right)^{1/c_{s}^{2}}\right]n_{0}.$ (40)

Reducing the ${c_{s,{\rm lim}}^{2}}$, the allowed values of $n(\mu)$
eventually degenerate into a line. This line is either a line of constant
speed of sound starting from $\\{\mu_{\rm L},n_{\rm L}\\}$ with a phase
transition at $\mu_{\rm H}$, or a line of constant speed of sound terminating
to $\\{\mu_{\rm H},n_{\rm H}\\}$ with a phase transition at $\mu_{\rm L}$.
Which of these cases is realised depends on the input parameter values, the
criteria for it is given by

$\displaystyle\frac{p_{\rm H}-p_{\rm L}}{\epsilon_{\rm h}-\epsilon_{\rm
l}}=\frac{\log(\mu_{\rm H}/\mu_{\rm L})}{\log(n_{\rm H}/n_{\rm L})}.$ (41)

If the left hand side of the eq. 41 is bigger than right hand side than the
value for $(c^{2}_{s})^{\rm min}_{\rm lim}$ is given by the solution of the
equation

$\displaystyle\int_{\mu_{\rm L}}^{\mu_{\rm H}}d\mu\,n_{\rm
H}\left(\frac{{\mu}}{\mu_{\rm H}}\right)^{1/(c^{2}_{s})^{\rm min}_{\rm lim}}$
$\displaystyle=\Delta p.$ (42)

Otherwise the equation reads

$\displaystyle\int_{\mu_{\rm L}}^{\mu_{\rm H}}d\mu\,n_{\rm
L}\left(\frac{{\mu}}{\mu_{\rm L}}\right)^{1/(c^{2}_{s})^{\rm min}_{\rm lim}}$
$\displaystyle=\Delta p.$ (43)

In general, these equations need to be solved numerically. The another way of
computing $(c^{2}_{s})^{\rm min}_{\rm lim}$ would be finding which of these
cases happen first. Instead of using criteria eq. 41 one can solve both eqs.
(42),(43) and take maximum of the solutions in order to find $(c^{2}_{s})^{\rm
min}_{\rm lim}$.

## References

  * Miller  _et al._ [2021] M. C. Miller _et al._ , The radius of PSR J0740+6620 from NICER and XMM-Newton data, Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _et al._ , A NICER view of the massive pulsar PSR J0740+6620 informed by radio timing and XMM-Newton spectroscopy, Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Demorest  _et al._ [2010] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. H. T. Hessels, A two-solar-mass neutron star measured using Shapiro delay, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Antoniadis  _et al._ [2013] J. Antoniadis _et al._ , A massive pulsar in a compact relativistic binary, Science 340, 1233232 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ [2021] E. Fonseca _et al._ , Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620, Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ , GW170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ , Multi-messenger observations of a binary neutron star merger, Astrophys. J. Lett. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Tolman [1939] R. C. Tolman, Static solutions of Einstein’s field equations for spheres of fluid, Phys. Rev. 55, 364 (1939).
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, On massive neutron cores, Phys. Rev. 55, 374 (1939).
  * Annala  _et al._ [2018] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Gravitational-wave constraints on the neutron-star-matter equation of state, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Margalit and Metzger [2017] B. Margalit and B. D. Metzger, Constraining the maximum mass of neutron stars from multi-messenger observations of GW170817, Astrophys. J. Lett. 850, L19 (2017), arXiv:1710.05938 [astro-ph.HE] .
  * Rezzolla  _et al._ [2018] L. Rezzolla, E. R. Most, and L. R. Weih, Using gravitational-wave observations and quasi-universal relations to constrain the maximum mass of neutron stars, Astrophys. J. Lett. 852, L25 (2018), arXiv:1711.00314 [astro-ph.HE] .
  * Ruiz  _et al._ [2018] M. Ruiz, S. L. Shapiro, and A. Tsokaros, GW170817, general relativistic magnetohydrodynamic simulations, and the neutron star maximum mass, Phys. Rev. D 97, 021501(R) (2018), arXiv:1711.00473 [astro-ph.HE] .
  * Bauswein  _et al._ [2017] A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas, Neutron-star radius constraints from GW170817 and future detections, Astrophys. J. Lett. 850, L34 (2017), arXiv:1710.06843 [astro-ph.HE] .
  * Radice  _et al._ [2018] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, GW170817: Joint constraint on the neutron star equation of state from multimessenger observations, Astrophys. J. Lett. 852, L29 (2018), arXiv:1711.03647 [astro-ph.HE] .
  * Most  _et al._ [2018] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, New constraints on radii and tidal deformabilities of neutron stars from GW170817, Phys. Rev. Lett. 120, 261103 (2018), arXiv:1803.00549 [gr-qc] .
  * Dietrich  _et al._ [2020] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Multimessenger constraints on the neutron-star equation of state and the Hubble constant, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Capano  _et al._ [2020] C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De, S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy, Stringent constraints on neutron-star radii from multimessenger observations and nuclear theory, Nat. Astron. 4, 625 (2020), arXiv:1908.10352 [astro-ph.HE] .
  * Landry and Essick [2019] P. Landry and R. Essick, Nonparametric inference of the neutron star equation of state from gravitational wave observations, Phys. Rev. D 99, 084049 (2019), arXiv:1811.12529 [gr-qc] .
  * Raithel  _et al._ [2018] C. A. Raithel, F. Özel, and D. Psaltis, Tidal deformability from GW170817 as a direct probe of the neutron star radius, Astrophys. J. Lett. 857, L23 (2018), arXiv:1803.07687 [astro-ph.HE] .
  * Raithel and Özel [2019] C. A. Raithel and F. Özel, Measurement of the nuclear symmetry energy parameters from gravitational-wave events, Astrophys. J. 885, 121 (2019), arXiv:1908.00018 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2020] G. Raaijmakers _et al._ , Constraining the dense matter equation of state with joint analysis of NICER and LIGO/Virgo measurements, Astrophys. J. Lett. 893, L21 (2020), arXiv:1912.11031 [astro-ph.HE] .
  * Essick  _et al._ [2020] R. Essick, P. Landry, and D. E. Holz, Nonparametric inference of neutron star composition, equation of state, and maximum mass with GW170817, Phys. Rev. D 101, 063007 (2020), arXiv:1910.09740 [astro-ph.HE] .
  * Al-Mamun  _et al._ [2021] M. Al-Mamun, A. W. Steiner, J. Nättilä, J. Lange, R. O’Shaughnessy, I. Tews, S. Gandolfi, C. Heinke, and S. Han, Combining electromagnetic and gravitational-wave constraints on neutron-star masses and radii, Phys. Rev. Lett. 126, 061101 (2021), arXiv:2008.12817 [astro-ph.HE] .
  * Essick  _et al._ [2021] R. Essick, I. Tews, P. Landry, and A. Schwenk, Astrophysical constraints on the symmetry energy and the neutron skin of 208Pb with minimal modeling assumptions, arXiv:2102.10074 [nucl-th] (2021).
  * Paschalidis  _et al._ [2018] V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke, and A. Sedrakian, Implications from GW170817 and I-Love-Q relations for relativistic hybrid stars, Phys. Rev. D 97, 084038 (2018), arXiv:1712.00451 [astro-ph.HE] .
  * Annala  _et al._ [2020] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Evidence for quark-matter cores in massive neutron stars, Nat. Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Ferreira  _et al._ [2020] M. Ferreira, R. Câmara Pereira, and C. Providência, Quark matter in light neutron stars, Phys. Rev. D 102, 083030 (2020), arXiv:2008.12563 [nucl-th] .
  * Minamikawa  _et al._ [2021] T. Minamikawa, T. Kojo, and M. Harada, Quark-hadron crossover equations of state for neutron stars: Constraining the chiral invariant mass in a parity doublet model, Phys. Rev. C 103, 045205 (2021), arXiv:2011.13684 [nucl-th] .
  * Blacker  _et al._ [2020] S. Blacker, N.-U. F. Bastian, A. Bauswein, D. B. Blaschke, T. Fischer, M. Oertel, T. Soultanis, and S. Typel, Constraining the onset density of the hadron-quark phase transition with gravitational-wave observations, Phys. Rev. D 102, 123023 (2020), arXiv:2006.03789 [astro-ph.HE] .
  * Baym  _et al._ [2018] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, From hadrons to quarks in neutron stars: a review, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Gandolfi  _et al._ [2019] S. Gandolfi, J. Lippuner, A. W. Steiner, I. Tews, X. Du, and M. Al-Mamun, From the microscopic to the macroscopic world: from nucleons to neutron stars, J. Phys. G 46, 103001 (2019), arXiv:1903.06730 [nucl-th] .
  * Raithel [2019] C. A. Raithel, Constraints on the neutron star equation of state from GW170817, Eur. Phys. J. A 55, 80 (2019), arXiv:1904.10002 [astro-ph.HE] .
  * Horowitz [2019] C. J. Horowitz, Neutron rich matter in the laboratory and in the heavens after GW170817, Annals Phys. 411, 167992 (2019), arXiv:1911.00411 [astro-ph.HE] .
  * Baiotti [2019] L. Baiotti, Gravitational waves from neutron star mergers and their relation to the nuclear equation of state, Prog. Part. Nucl. Phys. 109, 103714 (2019), arXiv:1907.08534 [astro-ph.HE] .
  * Chatziioannou [2020] K. Chatziioannou, Neutron star tidal deformability and equation-of-state constraints, Gen. Rel. Grav. 52, 109 (2020), arXiv:2006.03168 [gr-qc] .
  * Radice  _et al._ [2020] D. Radice, S. Bernuzzi, and A. Perego, The dynamics of binary neutron star mergers and GW170817, Ann. Rev. Nucl. Part. Sci. 70, 95 (2020), arXiv:2002.03863 [astro-ph.HE] .
  * Goldman and Nussinov [1989] I. Goldman and S. Nussinov, Weakly Interacting Massive Particles and Neutron Stars, Phys. Rev. D 40, 3221 (1989).
  * Giudice  _et al._ [2016] G. F. Giudice, M. McCullough, and A. Urbano, Hunting for Dark Particles with Gravitational Waves, JCAP 10, 001, arXiv:1605.01209 [hep-ph] .
  * Ciarcelluti and Sandin [2011] P. Ciarcelluti and F. Sandin, Have neutron stars a dark matter core?, Phys. Lett. B 695, 19 (2011), arXiv:1005.0857 [astro-ph.HE] .
  * Li  _et al._ [2012] A. Li, F. Huang, and R.-X. Xu, Too massive neutron stars: The role of dark matter?, Astropart. Phys. 37, 70 (2012), arXiv:1208.3722 [astro-ph.SR] .
  * Xiang  _et al._ [2014] Q.-F. Xiang, W.-Z. Jiang, D.-R. Zhang, and R.-Y. Yang, Effects of fermionic dark matter on properties of neutron stars, Phys. Rev. C 89, 025803 (2014), arXiv:1305.7354 [astro-ph.SR] .
  * Tolos and Schaffner-Bielich [2015] L. Tolos and J. Schaffner-Bielich, Dark Compact Planets, Phys. Rev. D 92, 123002 (2015), [Erratum: Phys.Rev.D 103, 109901 (2021)], arXiv:1507.08197 [astro-ph.HE] .
  * Ellis  _et al._ [2018] J. Ellis, G. Hütsi, K. Kannike, L. Marzola, M. Raidal, and V. Vaskonen, Dark Matter Effects On Neutron Star Properties, Phys. Rev. D 97, 123007 (2018), arXiv:1804.01418 [astro-ph.CO] .
  * Del Popolo  _et al._ [2020] A. Del Popolo, M. Le Delliou, and M. Deliyergiyev, Neutron Stars and Dark Matter, Universe 6, 222 (2020).
  * Jiménez and Fraga [2021] J. C. Jiménez and E. S. Fraga, Radial oscillations of quark stars admixed with dark matter, (2021), arXiv:2111.00091 [hep-ph] .
  * Damour and Esposito-Farese [1993] T. Damour and G. Esposito-Farese, Nonperturbative strong field effects in tensor - scalar theories of gravitation, Phys. Rev. Lett. 70, 2220 (1993).
  * He  _et al._ [2015] X.-T. He, F. J. Fattoyev, B.-A. Li, and W. G. Newton, Impact of the equation-of-state–gravity degeneracy on constraining the nuclear symmetry energy from astrophysical observables, Phys. Rev. C 91, 015810 (2015), arXiv:1408.0857 [nucl-th] .
  * Aparicio Resco  _et al._ [2016] M. Aparicio Resco, A. de la Cruz-Dombriz, F. J. Llanes Estrada, and V. Zapatero Castrillo, On neutron stars in $f(R)$ theories: Small radii, large masses and large energy emitted in a merger, Phys. Dark Univ. 13, 147 (2016), arXiv:1602.03880 [gr-qc] .
  * Doneva  _et al._ [2018] D. D. Doneva, S. S. Yazadjiev, N. Stergioulas, and K. D. Kokkotas, Differentially rotating neutron stars in scalar-tensor theories of gravity, Phys. Rev. D 98, 104039 (2018), arXiv:1807.05449 [gr-qc] .
  * Lope Oter  _et al._ [2019] E. Lope Oter, A. Windisch, F. J. Llanes-Estrada, and M. Alford, nEoS: Neutron Star Equation of State from hadron physics alone, J. Phys. G 46, 084001 (2019), arXiv:1901.05271 [gr-qc] .
  * Tews  _et al._ [2013] I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Neutron matter at next-to-next-to-next-to-leading order in chiral effective field theory, Phys. Rev. Lett. 110, 032504 (2013), arXiv:1206.0025 [nucl-th] .
  * Drischler  _et al._ [2019] C. Drischler, K. Hebeler, and A. Schwenk, Chiral interactions up to next-to-next-to-next-to-leading order and nuclear saturation, Phys. Rev. Lett. 122, 042501 (2019), arXiv:1710.08220 [nucl-th] .
  * Gorda  _et al._ [2018] T. Gorda, A. Kurkela, P. Romatschke, M. Säppi, and A. Vuorinen, Next-to-next-to-next-to-leading order pressure of cold quark matter: Leading logarithm, Phys. Rev. Lett. 121, 202701 (2018), arXiv:1807.04120 [hep-ph] .
  * Gorda  _et al._ [2021] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Soft Interactions in Cold Quark Matter, Phys. Rev. Lett. 127, 162003 (2021), arXiv:2103.05658 [hep-ph] .
  * Kurkela  _et al._ [2014] A. Kurkela, E. S. Fraga, J. Schaffner-Bielich, and A. Vuorinen, Constraining neutron star matter with Quantum Chromodynamics, Astrophys. J. 789, 127 (2014), arXiv:1402.6618 [astro-ph.HE] .
  * Annala  _et al._ [2021] E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Multimessenger constraints for ultra-dense matter, (2021), arXiv:2105.05132 [astro-ph.HE] .
  * Rhoades and Ruffini [1974] C. E. Rhoades, Jr. and R. Ruffini, Maximum mass of a neutron star, Phys. Rev. Lett. 32, 324 (1974).
  * Koranda  _et al._ [1997] S. Koranda, N. Stergioulas, and J. L. Friedman, Upper limit set by causality on the rotation and mass of uniformly rotating relativistic stars, Astrophys. J. 488, 799 (1997), arXiv:astro-ph/9608179 .
  * Tews  _et al._ [2019] I. Tews, J. Margueron, and S. Reddy, Confronting gravitational-wave observations with modern nuclear physics constraints, Eur. Phys. J. A 55, 97 (2019), arXiv:1901.09874 [nucl-th] .
  * Lope-Oter and Llanes-Estrada [2021] E. Lope-Oter and F. J. Llanes-Estrada, Maximum latent heat of neutron star matter independently of General Relativity, (2021), arXiv:2103.10799 [nucl-th] .
  * Hebeler  _et al._ [2013] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Equation of state and neutron star properties constrained by nuclear physics and observation, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Fujimoto and Fukushima [2020] Y. Fujimoto and K. Fukushima, Equation of state of cold and dense QCD matter in resummed perturbation theory, (2020), arXiv:2011.10891 [hep-ph] .
  * Fernandez and Kneur [2021] L. Fernandez and J.-L. Kneur, All order resummed leading and next-to-leading soft modes of dense QCD pressure, (2021), arXiv:2109.02410 [hep-ph] .