s v
verbatim,colback=white,colframe=codeborder,boxsep=0mm,left=2pt,right=2pt,top=2pt,bottom=2pt
#2

# Possible Value Analysis based on Symbolic Lattice

Qi Zhan<EMAIL_ADDRESS>1234-5678-9012 Zhejiang UniversityHangzhouChina

###### Abstract.

We propose a new static program analysis called program behavior analysis. The
analysis aims to calculate possible symbolic expressions for every variable at
each program point. We design a new lattice, transfer function, and widening
operator to accommodate the analysis. Furthermore, we extend the analysis to
interprocedural.

program analysis

††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

Static program analysis is a well-studied field (Cousot and Cousot, 1977),
which aims to reason about the behavior of computer programs without actually
running them. Static program analysis is useful not only in optimizing
compilers but also provide static bug or vulnerability detection tools
(Balakrishnan and Reps, 2004; Arzt et al., 2014).

In this paper, we propose a new static analysis technique, called program
behavior analysis. At a high level, we combine symbolic representation (King,
1976) and abstract interpretation (Cousot and Cousot, 1977) to calculate
possible symbolic values for every variable and memory cell at each program
point. With symbolic representation, the whole value flow is tracked by the
initial function parameters and global variables and their corresponding
memory structure; with abstract interpretation, possible values are
represented as an element of a lattice, where the top element $\top$
represents any value. The definition of the lattice is interleaved with the
symbolic expression, which is capable of representing any symbolic value.

As the lattice is infinite height, we propose a special widening technique to
approximate the transfer function. In addition, the values from different
control flows are joined together naturally based on the join of lattice. By
running the classical worklist algorithm in static analysis, we are able to
get precise function behavior from fixpoint results. The semantic level of the
code diff can be obtained by comparing the behavior of the function.
Furthermore, we believe that the semantic model can also be used to describe
the behavior of the whole program rather than limited to code diff, so we
extend the analysis to interprocedural. We discuss their potential usage in
future work.

The main contribution of this paper is a new type program behavior analysis
based on symbolic representation and we extend it to interprocedural.

##  2\. Possible Value Analysis

###  2.1. Preliminary

In this section, we present some basic techniques and terminology widely used
in academic research for program analysis that constitute the foundation of
our work.

Abstract Interpretation (Cousot and Cousot, 1977) is a well-established
foundation for static program analysis. In abstract interpretation, we
abstract the value as an element in abstract domain, and each program
statement is given an interpretation over abstract domain. This framework
guarantees that the results of analysis overapproximate the concrete behavior
of the program, which is called “soundness”.

Pointer Analysis is a particular program analysis technique, that aims to
compute the points-to information of the program, which is used in compiler
optimization, program verification, and program understanding (Jensen et al.,
2009; Arzt et al., 2014; Chandra et al., 2009). A typical flow-sensitive
Andersen’s algorithm is shown in Fig. 1. IN${}_{v}^{a}$ and OUT${}_{v}^{a}$
represent the input and output points-to set of the variable $a$ at the
program point $v$, respectively. For store statement, we separate the weak
update and strong update (Wilson and Lam, 1995). In the case that the input
points-to-set of the left-hand side variable is singleton, we can update the
memory value directly. When not, we must update the points-to relation by
unioning the original set.

Statement | Constrains |
---|---|---
$a=\&b$ | $\text{OUT}_{v}^{a}=\\{b\\}$ |
$a=b$ | $\text{OUT}_{v}^{a}=\text{IN}_{v}^{b}$ |
$a=\ast b$ | $\text{OUT}_{v}^{a}=\bigcup_{b^{\prime}\in\text{IN}_{v}^{b}}\text{IN}_{v}^{b^{\prime}}$ |
$\ast a=b$ | $\forall a^{\prime}\in\text{IN}_{v}^{a},\text{OUT}_{v}^{a^{\prime}}=\text{IN}_{v}^{b}$ | $|\text{IN}_{v}^{a}|=1$
| $\forall a^{\prime}\in\text{IN}_{v}^{a},\text{OUT}_{v}^{a^{\prime}}=\text{IN}_{v}^{a^{\prime}}\cup\text{IN}_{v}^{b}$ | $|\text{IN}_{v}^{a}|>1$
Figure 1. A typical flow sensitive Andersen’s algorithm

Symbolic Execution (King, 1976) is a way of executing a program abstractly.
The main idea is to treat concrete inputs as symbols and return symbolic
expressions expressed in input values. An abstract execution can cover
multiple possible inputs of the program that share a particular execution path
through the code. Our approach relies on the idea that using symbolic
execution to capture the possible value of the program, while we do not
consider the solution of symbolic constraints to gain concrete inputs.

###  2.2. Overview

Input: IN and OUT state of every basic block.

1 foreach  _basic block B_ do

2 OUT[B] = $\emptyset$;

3

4 end foreach

5W $\leftarrow$ all basic blocks;

6 while  _W is not empty_ do

7 Pick a basic block B from W;

8 old = OUT[B];

IN[B] = $\bigcup_{P~{}{\text{is a predecessor of B}}}$ OUT[P] ;

/* Section 2.3.5 */

OUT[B] = transfer(IN[B]) ;

/* Section 2.4 */

OUT[B] = IN[B] $\nabla$ OUT[B];

/* Section 2.3.3 */

9 if  _old $\neq$ OUT[B]_ then

10 Add all successors of B to W;

11

12 end if

13

14 end while

Algorithm 1 Worklist Algorithm

Possible value analysis computes the abstract value of each variable and
memory cell at each program point. The overall framework of possible value
analysis is illustrated in Algorithm 1. It is an abstract interpretation to
find a safe approximation for the set of symbolic expressions that each data
object holds at each program point. The main procedure is same as common data
flow analysis, i.e. we use worklist algorithm described in Alg. 1 to compute
the least fix point. We assume the control flow graph of the function is
given, and the analysis is performed on the control flow graph. Worklist
Algorithm is a typical technique used in the monotone framework of data flow
analysis (Kam and Ullman, 1977). It maintains IN and OUT sets for each basic
block, and a work list iteratively updates the data flow information of the
basic blocks. When OUT state of a basic block is updated, all successors of
the basic block are added to the worklist. The algorithm terminates until the
fixpoint is reached, i.e. worklist is empty. As the algorithm can be used to
solve data flow analysis constraints, and possible value analysis belongs to
the category, all we need to do is to instantiate abstract domain and transfer
function in the algorithm:

  1. (1)

Abstract domains. In Section 2.3, we define the lattice used in the analysis,
as well as the widening technique 2.3.3, and the abstract domains are defined
in Fig. 5.

  2. (2)

Transfer function. It interprets every program statement in the abstract
domain, we discuss it in Section 2.4.

###  2.3. Formalization

####  2.3.1. Language.

We formalize our approach using a simple SSA form (Rosen et al., 1988) call-
by-value language in Fig. 2(a), following the previous paper (Li et al., 2011,
2018) in the program analysis community. The language’s syntax and semantics
are straightforward and contain the basic features of a programming language.
Compared to (Li et al., 2011, 2018), we add a global variable declaration to
the language, which is necessary to depict the behavior related to the global
variable. In addition, we assume that the last statement of a function is
always a return statement to facilitate the discussion and inference rule.

ProgramP $\displaystyle\coloneq(F\mid G)^{+}$ FunctionF $\displaystyle\coloneq
f(v_{1},\dots,v_{n})\\{B*;\\}$ BlockB $\displaystyle\coloneq S$ StatementS
$\displaystyle\coloneq v_{1}=v_{2}\mid v_{1}=\&v_{2}$ $\displaystyle\mid
v_{1}=\ast v_{2}\mid\ast v_{1}=v_{2}$ $\displaystyle\mid
v_{1}=v_{2}~{}\text{binop}~{}v_{3}\mid v_{1}=\text{unop }v_{2}$
$\displaystyle\mid\text{if}(v)~{}\text{then}~{}\text{goto}~{}B$
$\displaystyle\mid v=\text{call}~{}f(v_{1},\dots,v_{n})\mid\text{return}~{}v$
$\displaystyle\mid v=\phi(v_{1},\dots,v_{n})$ $\displaystyle\mid s_{1};s_{2}$
binop $\displaystyle\coloneq+\mid-\mid\ast\mid/\mid\cdots$ unop
$\displaystyle\coloneq-\mid\neg\mid\cdots$ Global $\displaystyle g\in G$ (a)
Syntax of the language

Expressione $\displaystyle\coloneq\text{primitive}~{}p$
$\displaystyle\mid\text{function}(l_{1},\dots,l_{n})$
$\displaystyle\mid\text{binop}(l_{1},l_{2})$ $\displaystyle\mid\text{unop}(l)$
$\displaystyle\mid\phi(l_{1},\dots,l_{n})$ Primitivep
$\displaystyle\coloneq\text{const}(i)$ $\displaystyle\mid\text{arg}(i)$
$\displaystyle\mid\text{mem}(o)$ $\displaystyle\mid\text{global}(o)$ binop
$\displaystyle\coloneq+\mid-\mid\ast\mid/\mid\cdots$ unop
$\displaystyle\coloneq-\mid\neg\mid\cdots$ Lattice Elementl
$\displaystyle\coloneq\top\mid\bot\mid\text{Expression}(e)$ (b) Syntax of
symbolic expression and lattice

Figure 2. Syntax of language and lattice

####  2.3.2. Lattice Design.

$\top$binop$(\top,\top)$binop(binop($\top$, $\top$),
$\top$)binop(binop($\top$, $\top$), $\bot$)binop($\top$, binop($\top$,
$\top$))binop($\bot$, binop($\top$,
$\top$))binop($\bot,\bot)$primitive($p$)unop($\top$)unop($\bot$)$\bot$ Figure
3. Illustration of a subset of the lattice. For simplicity, we only consider
the primitive $p$ and binary expressions and unop expressions. Function and
phi expression can be added in a similar way. As the lattice is infinite
height, we have to omit some elements, which are represented as a dashed line.

As with every static analysis method, we design a lattice to abstract the
possible values of the program, as shown in Fig. 2(b). In our work, the
definition of lattice $\mathcal{L}$ is unusual, as it recursively interleaved
with the definition of the expression $e$.

Lattice $l$ is the abstract value of possible symbolic expressions. $\top$ and
$\bot$ are the maximal and minimal elements of the lattice, respectively. For
example, $f(\top,\text{binop}(\top,\text{const}(1)))$ is a valid lattice
element. It can be understood as a specific function $f$, where the first
argument is any possible value, and the second argument is a binary operator
of any possible value and const. The information is more precise than
$f(\top,\top)$ since we can restrict the possible constructions for the second
argument.

Expression $e$ can be a primitive value, a function call, a binary operation,
a phi node, or a unary operation, which corresponds to the concrete syntax of
the language in Fig. 2(a). Memory $o$ here is a memory location consisting of
a base address and an offset, where the

The partial order $\sqsubseteq$ and the join operator of set $\mathcal{L}$ are
defined in Fig. 4 and Table 1 recursively. For the same constructor, the
partial order and join are defined as the partial order and join of the
corresponding fields, respectively. It is easy to check that the order set
equipped with the above operators is indeed a complete lattice. Every subset
has a least upper bound and a greatest lower bound, and the join and meet
operators are idempotent, commutative, and associative. In the discussion
below, we use $\mathcal{L}=\langle L,\sqsubseteq,\sqcup,\bot,\top\rangle$ to
represent the lattice.

It seems that the lattice is too simple and not expressive enough, as it
consists of only one non-trivial element $e$. The magic lies in the fact that
the expression can contain lattice elements. An expression $e$ can contain
lattice elements, and a lattice element can contain expression $e$, the
lattice can be any symbolic expression. The strong expressiveness of the
lattice makes it possible to not miss any single value of the program. On the
other hand, the height of the lattice is infinite, which does not satisfy the
termination requirement of the worklist algorithm directly. Fig. 3 illustrates
a subset of the lattice, which only consists of primitive value and binary
operator.

$\displaystyle\bot\sqsubseteq l$ $\displaystyle l\sqsubseteq\top$
$\displaystyle p=p^{\prime}$
$\displaystyle\Rightarrow\text{primitive}(p)\sqsubseteq\text{primitive}(p^{\prime})$
$\displaystyle l_{i}\sqsubseteq l_{i}^{\prime}\text{ for }i=1,\dots,n$
$\displaystyle\Rightarrow\text{function}(l_{1},\dots,l_{n})\sqsubseteq\text{function}(l_{1}^{\prime},\dots,l_{n}^{\prime})$
$\displaystyle l_{i}\sqsubseteq l_{i}^{\prime}\text{ for }i=1,\dots,n$
$\displaystyle\Rightarrow\phi(l_{1},\dots,l_{n})\sqsubseteq\phi(l_{1}^{\prime},\dots,l_{n}^{\prime})$
$\displaystyle l_{1}\sqsubseteq l_{1}^{\prime}\text{ and }l_{2}\sqsubseteq
l_{2}^{\prime}$
$\displaystyle\Rightarrow\text{binop}(l_{1},l_{2})\sqsubseteq\text{binop}(l_{1}^{\prime},l_{2}^{\prime})$
$\displaystyle l\sqsubseteq l^{\prime}$
$\displaystyle\Rightarrow\text{unop}(l)\sqsubseteq\text{unop}(l^{\prime})$
Figure 4. Partial order of lattice Table 1. Join of lattice element

left | right | join $\sqcup$
---|---|---
$\bot$ | $l$ | $l$
$l$ | $\bot$ | $l$
primitive($p$) | primitive($p^{\prime}$) | primitive($p$), if $p=p^{\prime}$; $\top$, otherwise
function($l_{1},\dots,l_{n}$) | function($l_{1}^{\prime},\dots,l_{n}^{\prime}$) | function($l_{1}\sqcup l_{1}^{\prime},\dots,l_{n}\sqcup l_{n}^{\prime}$)
$\phi(l_{1},\dots,l_{n})$ | $\phi(l_{1}^{\prime},\dots,l_{n}^{\prime})$ | $\phi(l_{1}\sqcup l_{1}^{\prime},\dots,l_{n}\sqcup l_{n}^{\prime})$
binop($l_{1},l_{2}$) | binop($l_{1}^{\prime},l_{2}^{\prime}$) | binop($l_{1}\sqcup l_{1}^{\prime},l_{2}\sqcup l_{2}^{\prime}$)
unop($l$) | unop($l^{\prime}$) | unop($l\sqcup l^{\prime}$)
otherwise | $\top$

meet $\sqcap$ is defined similarly, we omit it here.

####  2.3.3. Widening

Since the height of the lattice is not limited, the worklist algorithm may not
terminate. For example, consider a very simple add statement $a=a+1$ in a
loop, a will be evaluated as (+ a 1), (+ (+ a 1) 1), and so on. Widening
technique is introduced by (Cousot and Cousot, 1977) to overcome this
limitation. In this section, we propose a series of simple widening strategies
for our lattice.

The intuition for widening comes from the recursive definition of lattice. As
shown in Fig. 3, the lattice can be seen as a layer hierarchy.
$\top,\bot,\text{primitive}(p)$ consist the first layer of the lattice.
binop$(\top,\top)$, binop$(\top,\text{primitive})$ consit the second layer and
so on. It is easy to define a depth function $\mathcal{D}$ of an element in
lattice:

$\displaystyle\mathcal{D}(\top)=0,\mathcal{D}(\bot)$
$\displaystyle=0,\mathcal{D}(\text{primitive}(p))=0$
$\displaystyle\mathcal{D}(\text{function}(l_{1},\dots,l_{n}))$
$\displaystyle=1+\max(\mathcal{D}(l_{1}),\dots,\mathcal{D}(l_{n}))$
$\displaystyle\mathcal{D}(\text{binop}(l_{1},l_{2}))$
$\displaystyle=1+\max(\mathcal{D}(l_{1}),\mathcal{D}(l_{2}))$
$\displaystyle\mathcal{D}(\text{unop}(l))$ $\displaystyle=1+\mathcal{D}(l)$
$\displaystyle\mathcal{D}(\phi(l_{1},\dots,l_{n}))$
$\displaystyle=1+\max(\mathcal{D}(l_{1}),\dots,\mathcal{D}(l_{n}))$

Based on the depth function, we can define a series of widening operators
$\nabla_{i}(x,f(x))$, where we only consider $f(x)$ in the simple widening
strategy. Only the element with a depth less than or equal to $i$ is kept
unchanged, and the element with a depth greater than $i$ is widened to the
corresponding element in $i$-th layer. For examle,
$\nabla_{0}(x,\text{binop}(\top,\top))=\top.$ The height of the sublattice
$L_{i}$ is limited by $i$, and the worklist algorithm will always terminate in
$L_{i}$.

####  2.3.4. Abstract Domains.

Based on the definition of lattice, abstract domains can be defined as usual.
The symbols and abstract domains are listed in Fig. 5, and introduced as
follows:

  * •

Expression $e$ is a value expression in the lattice $\mathcal{L}$, variables
$v$ are used to represent the local variables, and memory cell $o$ is used to
represent the memory locations.

  * •

The local variable context $\Sigma$ is a mapping from variable names to value
expressions.

  * •

The memory $m$ is a mapping from memory locations to value expressions. We
regard global variables as a part of memory locations.

  * •

The points-to relation $\Delta$ is a mapping of memory locations to the set of
memory locations, representing the possible locations to which a memory cell
can point.

$\displaystyle\begin{array}[]{ll}\text{Expression}~{}e\in
E&\text{Object}~{}o\in\mathcal{O}\quad\quad\text{Lattice}~{}l\in\mathcal{L}\\\
\text{Variables}~{}v\in\mathcal{V}&\text{Local
context}~{}\Sigma:\mathcal{V}\to\mathcal{L}\\\
\text{Memory}~{}m:\mathcal{O}\to\mathcal{L}&\text{Points-to
relation}~{}\Delta:\mathcal{O}\to\mathcal{P}(\mathcal{O})\\\ \end{array}$
Figure 5. Abstract domains

####  2.3.5. Join Operator

$\sqcup$ is originally used to join two elements in lattice, and is also
overloaded to merge context environments for convenience. Join of lattice
element is the definition of the lattice; Join of local variable context and
memory can be reduced to the join of lattice element; Join of points-to
relation is the union of the two sets. $\sqcup_{L}$ in left is the join
operator of the lattice $\mathcal{L}$. we assume that an element not in the
domain of the mapping is $\bot$ or an empty set.

$\displaystyle\centering\text{Expression}(e_{1})\sqcup\text{Expression}(e_{2})\@add@centering$
$\displaystyle=\text{Expression}(e_{1}\sqcup_{L}e_{2})$
$\displaystyle\text{Local}(\Sigma_{1})\sqcup\text{Local}(\Sigma_{2})$
$\displaystyle=\text{Local}(\Sigma),\text{where}~{}\Sigma(v)=\Sigma_{1}(v)\sqcup_{L}\Sigma_{2}(v)$
$\displaystyle\text{Memory}(m_{1})\sqcup\text{Memory}(m_{2})$
$\displaystyle=\text{Memory}(m),\text{where}~{}m(o)=m_{1}(o)\sqcup_{L}m_{2}(o)$
$\displaystyle\text{Points-to}(\Delta_{1})\sqcup\text{Points-to}(\Delta_{2})$
$\displaystyle=\text{Points-
to}(\Delta),\text{where}~{}\Delta(o)=\Delta_{1}(o)\cup\Delta_{2}(o)$

###  2.4. Transfer Function

[S-BinOp]Σ; m; Δ; v_1 ⇒e_1
Σ; m; Δ; v_2 ⇒e_2Σ; m; Δ; v = v_1 binop v_2 ⇓Σ[v ↦e_1 binop e_2]; m; Δ

[S-UnOp]Σ; m; Δ; v ⇒eΣ; m; Δ; v = unop v ⇓Σ[v ↦unop e]; m; Δ

[S-Seq]Σ; m; Δ; s_1 ⇓Σ’ ; m’; Δ’;
Σ’; m’; Δ’; s_2 ⇓Σ”; m”; Δ”Σ; m; Δ; s_1;s_2 ⇓Σ”; m”; Δ”

[S-Call]Σ; m; Δ; v_i ⇒e_i, ∀i ∈[1,n]Σ; m; Δ; v = f(v_1,…,v_n) ⇓Σ[v
↦f(e_1,…,e_n)]; m; Δ (intra)

[S-Return]Σ; m; Δ; v ⇒e Σ; m; Δ; return v ⇓Σ; m; Δ (intra)

[S-Phi]Σ; m; Δ; v_i ⇒e_i, ∀i ∈[1,n]Σ; m; Δ; v = ϕ(v_1, …, v_n) ⇓Σ[v ↦ϕ(e_1, …,
e_n)]; m; Δ

[S-Assign]Σ; m; Δ; v_2 ⇒e
Δ; v_1 = v_2 ⇓_pa Δ’Σ; m; Δ; v_1 = v_2 ⇓Σ[v_1 ↦e]; m; Δ’

[S-Ref]Δ;v_1=& v_2 ⇓_pa Δ’Σ; m; Δ; v_1 = &v_2 ⇓Σ; m; Δ’

[S-Load]e = ⨆m[o_i], ∀o_i ∈Δ(v_2)
Δ; v_1 = * v_2 ⇓_pa Δ’Σ; m; Δ; v_1 = *v_2 ⇓Σ[v_1 ↦e]; m; Δ’

[S-Store]Σ; m; Δ; v_2 ⇒e
∀o_i ∈Δ(v_1)
Δ; *v_1 = v_2 ⇓_pa Δ’Σ; m; Δ; *v_1 = v_2 ⇓Σ; m[o_i ↦m[o_i] ⊔e]; Δ’

(a) Intraprocedural transfer function rules

[S-Call]Σ; m; Δ; v_i ⇒e_i, ∀i ∈[1,n]
Σ[v_1↦e_1,…,v_n↦e_n]; m; Δ; s ⇓Σ’; m’; Δ’
Σ’; m’; Δ’; v ⇒e Σ; m; Δ; w = f(v_1,…,v_n){s; return v} ⇓Σ[w↦e]; m’; Δ’
(inter)

(b) Interprocedural transfer function rules

Figure 6. Transfer function rules

As we already defined the abstract domain and corresponding operations, the
last step is to define the transfer function. To this end, we define the
transfer function that soundly approximates the concrete semantics. The
transfer function of intraprocedural program behavior analysis serves as the
core of the analysis, and it is represented as a series of semantic rules.

$\begin{array}[]{c}\Sigma;m;\Delta;v\Rightarrow e\\\
\hline\cr\Sigma;m;\Delta,s\Downarrow\Sigma^{\prime};m^{\prime};\Delta^{\prime}\end{array}$

The premise of the judgment is read as, “with local variable context $\Sigma$,
memory $m$, and pointer mapping container $\Delta$, varible $v$ evaluates to
$e\in\mathcal{L}$.” The conclusion of the judgment is read as, “with local
variable context $\Sigma$, memory $m$, and pointer mapping container $\Delta$,
statement $s$ produces a new local variable context $\Sigma^{\prime}$, memory
$m^{\prime}$, and pointer mapping container $\Delta^{\prime}$.” The details of
the semantic is shown in Fig. 6(a). We divide the semantic rules into two
parts: memory-related rules and non-memory-related rules, and discuss them
separately.

❶ Rule 6(a) and 6(a) are standard symbolic execution rules, we evaluate the
expression on the current environment, and substitute the target variable with
a new element according to the operator. Rule 6(a) is similar to the previous
two rules, except that we create a new symbol $f(e_{1},\dots,e_{n})$ to
represent the return value of the function call. Rule 6(a) is used to handle
the phi node, which evaluates to a new phi node based on the values of the
incoming variables. Rule 6(a) is just a composition of two statements that
transfer the state.

❷ Rule 6(a), 6(a), 6(a) and 6(a) are used to handle the memory related
statements. To make the semantic rule clear, we separate the transfer function
into two parts: 1) points-to relation $\Delta$ transfer; 2) and local variable
context $\Sigma$ and memory $m$ transfer. Symbol $\Downarrow_{pa}$ denotes
arbitrary pointer analysis transfer function. Since there exist a lot of
pointer analysis algorithms, we can use any existing algorithms to implement
the memory-related part. Fig. 1 is an instance of the pointer analysis. So we
do not discuss the details of the pointer analysis in this paper. Besides, the
two parts can be implemented separately. On the other side, pointer analysis
is a time-consuming process, it is reasonable to use a flow-insensitive
pointer analysis to calculate the points-to relation, which may not need to be
calculated at every program point. In that case, the points-to relation is
calculated as a pre-analysis, and there is no need to update the points-to
relation in the semantic rule.

6(a) and 6(a) are straightforward since they simply update the pointer mapping
$\Delta$ based on the statement. For load statement $v_{1}=\ast v_{2}$, we use
the points-to relation to obtain the possible memory location of variable
$v_{2}$, then join the values of these memory locations to get the value of
$v_{1}$. For the store statement $\ast v_{1}=v_{2}$, we also use the points-to
relation to get the possible memory location of variable $v_{1}$, then store
the value in those memory locations. Similarly to strong and weak updates in
pointer analysis, we use the join operator to update the memory value when
there are multiple memory locations.

It is worth noting that the pointer analysis is not a trivial part of this
analysis. In the ideal case, we do not need an extra pointer analysis to
manage the memory information; the points-to relation can be directly
calculated from the possible values of the variables. The problem arises when
we encounter statements like $*a=b$ and $a\Rightarrow\top$ in practice. The
only sound way to interpret the statement is to assign to every possible
memory location the value evaluated from $b$. As a result, the entire analysis
will be useless, since for all memory cells $o$ in $m$, $m(o)=\top$.
Therefore, in store and load statement, we look up the points-to relation for
the address variable from pointer analysis instead of evaluating the address.

As an intra-procedural analysis, we do not need to consider the effects of the
callee function. To be specific, we make two simplifications in the analysis.
We assume:

  * •

The function call does not modify the memory pointed by the global variable
and its arguments.

  * •

The memory region pointed by parameters and global variables, directly or
indirectly is disjoint.

The reader may doubt the soundness of the analysis after simplification. In
general, a function call $f(v_{1},\dots,v_{n})$ can modify the memory pointed
by global variable and its arguments, which is ignored in our rule 6(a). We
argue that in intraprocedure, it is a reasonable simplification to ignore the
side effects of the callee function. It does not affect our analysis, since we
are analyzing the behavior of the given function itself. All other function
calls are treated as a black box. We do not need to know the exact behavior of
the callee function to analyze the function considered in the intraprocedural.
In the other hand, the whole analysis can be too imprecise to be useful if we
consider the callee function and set all related memory cells to $\top$,
thereby losing track of all information. When we expand the analysis to
interprocedural, 6(a) is replaced by its actual semantic of the callee
function 6(b). The same reason is also applied to the second simplification.

To ensure the effectiveness of the worklist algorithm, we also need to
guarantee the monotonicity of the proposed transfer function. For statements
that do not involve memory, the monotonicity is obvious based on the partial
order rule in Fig. 4. Taking binary operation $a=b+c$ as an example, suppose
$b\Rightarrow l_{1}\sqsubseteq l_{1}^{\prime}$ and $c\Rightarrow
l_{2}\sqsubseteq l_{2}^{\prime}$, then $a\Rightarrow$ binop($l_{1},l_{2}$)
$\sqsubseteq$ binop($l_{1}^{\prime},l_{2}^{\prime}$). Suppose the points-to
information $\Delta$ is monotonic, the transfer function of store statement is
also monotonic since we just store the value to every related memory location.
Load statement is also monotonic due to the monotonicity of the join operator.

Suppose that the points-to-information $\Delta$ is calculated by a sound
pointer analysis. Program behavior analysis performs essentially a flow-
sensitive analysis, the join of program state is defined from the join of the
lattice. Thus, program behavior analysis is sound if the pre-analysis is
sound. In summary, our analysis results are safe approximations of the actual
behavior of the program.

Equipped with above techniques and language, we can give the following
definition:

###### Definition 2.1 (Possible Value Analysis).

For a given function $f$, intraprocedural possible value analysis (pva) aims
to compute the abstract value $l\in\mathcal{L}$ for every variable and memory
cell at every program point, i.e.
$\text{pva}:(\mathcal{V}\bigcup\mathcal{O})\to\mathcal{L}$.

###### Theorem 2.2 (Soundness).

Suppose $\alpha$ is the function that represents the transfer function and
join, $\llbracket f\rrbracket$ is the concrete semantic of the single
function, and $\llbracket f\rrbracket^{*}$ is the abstract semantic of the
function, Then we have

$\alpha(\llbracket f\rrbracket)\sqsubseteq\llbracket f\rrbracket^{*}$

Currently, our analysis is path-insensitive. Our approach does not record the
control flow information during the analysis like symbolic execution does,
since we do not rely on the path constraint to solve or add certain assertions
to the program. It is also possible to add an assert statement to the program
and constrain the behavior following the previous work (Møller and
Schwartzbach, 2018). Regardless of the path constraint, control flow
information can be calculated since the condition itself is a lattice element.

##  3\. Discussion

###  3.1. Interprocedural Analysis

We have discussed the intraprocedural analysis in the previous section, and
the analysis can be expanded the analysis to interprocedural, i.e. from
function level to program level. Similarly to compare the behavior of two
functions, suppose that we have two programs $P$ and $P^{\prime}$, and we want
to compare the behavior of the two programs. From the perspective of black
box, we do not care about internal function calls or memory modification since
they are not observable objects. In the most extreme case, all function calls
are not recorded, and the whole program behavior can be reduced to calls to
IO-related operations, which corresponds to the common meaning of side effects
exactly.

It is worth noting that the implementation of interprocedural program behavior
analysis is simpler than the intraprocedural one due to less information
required. The semantic rule of the interprocedural analysis is shown in Fig.
6(b). We mark the different parts of the rule with numbers. Pointer analysis
$\Downarrow_{pa}$ can be replaced with any global pointer analysis algorithm.
In interprocedural analysis, the modification to function parameters,
variables, and every function return value can be directly returned to the
function call point, rather than being recorded as a feature. Therefore, in
rule 6(b), instead of recording the function call, we evaluate the callee
function in the current memory context. The context of the local variable
$\Sigma$ is represented by $\\{v_{i}\mapsto e_{i}\\}$, where $v_{i}$ is the
parameter of the callee function and $e_{i}$ is the value of the parameter.
After the evaluation, we set the return value and the new memory context back
to the caller function. As common interprocedural analysis, we use the call
graph to control the analysis process and the return value and context should
merge when the function is called in different places. The details of context
sensitivity are omitted in the paper.

In practice, program is not only composed of user-defined functions, but also
the standard library and third-party library (TPL) functions. Our framework
can easily handle library functions since we can just record them as done in
the intraprocedural analysis. The grain of the analysis can be controlled at
any level by interacting using the interprocedural or intraprocedural rule.

###  3.2. Limitation

Many code changes still can be evaluated to $\top$ and are useless in
practice. Our approach can benefit from the improvement of various program
analysis techniques; we believe that the precision of the analysis can be
further improved by more sophisticated pointer analysis algorithms or control
flow analysis algorithms. Benefiting from SSA form, our analysis is not
bothered by multiple assignments to the same variable. However, the SSA form
in many implementations is partial such as LLVM (Lattner and Adve, 2004),
which means that the memory-related instruction is actually not ssa. The
precision of our approach can be limited by $e=\bigsqcup{m[o_{i}]},\forall
o_{i}\in\Delta(v_{2})$ in 6(a).

##  4\. Related Work

Program behavior analysis work is the main part of our work. To the best of
our knowledge, we are the first to propose a static analysis to represent the
possible value of every variable and memory cell at each program point in a
symbolic representation, and we used the representation to capture the
semantic of code diff. Interval analysis (Cousot and Cousot, 1977) is a
classical static analysis technique to compute the possible value of a
variable at each program point. and value set analysis (VSA) (Balakrishnan and
Reps, 2004) is a static binary analysis technique, which uses abstract
interpretation to safely approximate the set of values of each data object at
each program point. The main difference between VSA and interval analysis and
our work is that their concern is numerical value, so the analysis is on the
numerical domain; while we consider the symbolic argument and memory, so the
analysis is on a symbolic domain.

##  5\. Conclusion

In this paper, we propose a new type possible value analysis and design a new
symbolic lattice. In this paper, we propose a new methodology using which we
can capture the semantic of code diff by observing the behavior of the
function. To this end, we propose a new type of program analysis, program
behavior analysis, to calculate possible values at each program point.
Furthermore, we extend the analysis to interprocedural and discuss possible
future work and downstream tasks.

## References

  * (1)
  * Arzt et al. (2014) Steven Arzt, Siegfried Rasthofer, Christian Fritz, Eric Bodden, Alexandre Bartel, Jacques Klein, Yves Le Traon, Damien Octeau, and Patrick McDaniel. 2014.  Flowdroid: Precise context, flow, field, object-sensitive and lifecycle-aware taint analysis for android apps.  _Acm Sigplan Notices_ 49, 6 (2014), 259–269.
  * Balakrishnan and Reps (2004) Gogul Balakrishnan and Thomas Reps. 2004.  Analyzing Memory Accesses in x86 Executables. In _Compiler Construction_ , Evelyn Duesterwald (Ed.). Springer Berlin Heidelberg, Berlin, Heidelberg, 5–23.
  * Chandra et al. (2009) Satish Chandra, Stephen J Fink, and Manu Sridharan. 2009.  Snugglebug: a powerful approach to weakest preconditions. In _Proceedings of the 30th ACM SIGPLAN Conference on Programming Language Design and Implementation_. 363–374.
  * Cousot and Cousot (1977) Patrick Cousot and Radhia Cousot. 1977.  Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints. In _Proceedings of the 4th ACM SIGACT-SIGPLAN Symposium on Principles of Programming Languages_ (Los Angeles, California) _(POPL ’77)_. Association for Computing Machinery, New York, NY, USA, 238–252.  https://doi.org/10.1145/512950.512973
  * Jensen et al. (2009) Simon Holm Jensen, Anders Møller, and Peter Thiemann. 2009.  Type analysis for JavaScript. In _Static Analysis: 16th International Symposium, SAS 2009, Los Angeles, CA, USA, August 9-11, 2009. Proceedings 16_. Springer, 238–255.
  * Kam and Ullman (1977) John B Kam and Jeffrey D Ullman. 1977.  Monotone data flow analysis frameworks.  _Acta informatica_ 7, 3 (1977), 305–317.
  * King (1976) James C. King. 1976.  Symbolic execution and program testing.  _Commun. ACM_ 19, 7 (jul 1976), 385–394.  https://doi.org/10.1145/360248.360252
  * Lattner and Adve (2004) Chris Lattner and Vikram Adve. 2004.  LLVM: A compilation framework for lifelong program analysis & transformation. In _International symposium on code generation and optimization, 2004. CGO 2004._ IEEE, 75–86.
  * Li et al. (2011) Lian Li, Cristina Cifuentes, and Nathan Keynes. 2011.  Boosting the performance of flow-sensitive points-to analysis using value flow. In _Proceedings of the 19th ACM SIGSOFT Symposium and the 13th European Conference on Foundations of Software Engineering_ (Szeged, Hungary) _(ESEC/FSE ’11)_. Association for Computing Machinery, New York, NY, USA, 343–353.  https://doi.org/10.1145/2025113.2025160
  * Li et al. (2018) Lian Li, Cristina Cifuentes, and Nathan Keynes. 2018.  Precise and scalable context-sensitive pointer analysis via value flow graph. In _Proceedings of the 2013 International Symposium on Memory Management_ (Seattle, Washington, USA) _(ISMM ’13)_. Association for Computing Machinery, New York, NY, USA, 85–96.  https://doi.org/10.1145/2464157.2466483
  * Møller and Schwartzbach (2018) Anders Møller and Michael I. Schwartzbach. 2018.  Static Program Analysis.  Department of Computer Science, Aarhus University, http://cs.au.dk/~amoeller/spa/.
  * Rosen et al. (1988) B. K. Rosen, M. N. Wegman, and F. K. Zadeck. 1988.  Global value numbers and redundant computations. In _Proceedings of the 15th ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages_ (San Diego, California, USA) _(POPL ’88)_. Association for Computing Machinery, New York, NY, USA, 12–27.  https://doi.org/10.1145/73560.73562
  * Wilson and Lam (1995) Robert P. Wilson and Monica S. Lam. 1995.  Efficient context-sensitive pointer analysis for C programs. In _Proceedings of the ACM SIGPLAN 1995 Conference on Programming Language Design and Implementation_ (La Jolla, California, USA) _(PLDI ’95)_. Association for Computing Machinery, New York, NY, USA, 1–12.  https://doi.org/10.1145/207110.207111


# Tidal Deformation and Radial Pulsations of Neutron Star with Holographic
Multiquark Core

Sitthichai Pinkanjanarod<EMAIL_ADDRESS><EMAIL_ADDRESS>High Energy Physics Theory Group, Department of Physics,Faculty of Science,
Chulalongkorn University, Bangkok 10330, Thailand  Department of Physics,
Faculty of Science, Kasetsart University, Bangkok 10900, Thailand  Piyabut
Burikham<EMAIL_ADDRESS>High Energy Physics Theory Group, Department of
Physics,Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand
Supakchai Ponglertsakul<EMAIL_ADDRESS>High Energy Physics Theory
Group, Department of Physics,Faculty of Science, Chulalongkorn University,
Bangkok 10330, Thailand  Department of Physics, Faculty of Science, Silpakorn
University,Nakhon Pathom 73000, Thailand

###### Abstract

Tidal deformation of neutron star with multiquark core is calculated using
nuclear and multiquark equations of state. The equation of state of the
multiquark phase from the holographic Sakai-Sugimoto (SS) model is relatively
stiff in the low density region and becomes softer at high densities. The
values of Love number and dimensionless deformation parameter, $k_{2}$ and
$\Lambda$, are found to be within the physically viable range under the
present constraints. Radial pulsation frequencies of the multiquark core for
$n=0-5$ modes are calculated for the entire mass range. For $M_{\rm
core}\simeq 2M_{\odot}$, the fundamental-mode frequency is approximately $2.5$
kHz for the energy density scale $\epsilon_{s}=23.2037$ GeV/fm3 of the
holographic SS model, this frequency is proportional to $\sqrt{\epsilon_{s}}$.

##  I Introduction

A number of detections of massive neutron stars (NS) with masses approximately
larger than 2 solar mass ($M_{\odot}$) Antoniadis:2013pzd ; Miller:2021qha ;
Clark:2002db ; Romani:2012rh ; Romani:2012jf ; vanKerkwijk:2010mt ;
Linares:2018ppq ; Bhalerao:2012xe ; Nice:2005fi ; Demorest:2010bx ;
Freire:2007sg ; Quaintrell:2003pn  indicates that the equation of state (EoS)
of the core of the massive NS could be softer than the usual stiff baryonic
nuclear matter in the crust. Substantial scan of various forms of the EoS
based on the sound speed and adiabatic index reveals the possibility that the
massive NS has a quark matter core Annala:2019puf .

The equation of state of nuclear matter at high density is still an open
problem of the strong interaction. For low density and extremely high
temperature (around trillion K), interaction between quarks and gluons
described by quantum chromodynamics (QCD) is weak and the quarks and gluons
could become deconfined in a phase of quark-gluon plasma (QGP). However, as
the NS cooled down, QCD becomes strongly coupled at low temperature while the
baryon chemical potential/density could still be large. In particular, at the
core of the NS, the radius of confinement could become larger due to extreme
density and quarks could be effectively deconfined from the localized hadrons.
However, the quarks are yet to be confined by gravity within the star. There
is a number of nonperturbative QCD approaches for the EoS of the interior of
the NS. For instance, the bag models of free quarks inside the star are simple
and effective but they are based on the assumption of weakly coupled quarks.
Lattice QCD, which is very useful for nonperturbative QCD, however becomes
unreliable at large baryon chemical potential Lombardo:2006yc . Colour
superconductivity (CSC) assumes that quarks form diquark condensate at high
density which can be described by the Nambu-Jona-Lasinio type theory
Alford:2007xm .

There is a complementary powerful tool for studying strongly coupled gauge
theories that emerged from the discovery of the gauge/gravity duality, i.e.,
the AdS/CFT correspondence Gubser:1998bc ; maldacena ; witb . One of the most
astonishing features in the gauge/string duality is the two theories defined
in different dimensions are conjectured to be equivalent under a strong-weak
duality that is useful to study the strongly coupled QCD. In the holographic
QCD the mesons and baryons are dually described by higher dimensional
configurations consisting of open strings and D-branes in the specific
background (see e.g. Ref. Witten:1998xy ; Callan:1998iq ; Brandhuber:1998xy
for earlier works and Ref. Nakas:2020hyo  for more recent development). In
Ref. bch , holographic multiquark (MQ) phase in the deconfined background with
broken chiral symmetry based on the Sakai-Sugimoto model ss1 ; ss2 ;
Aharony_chiral ; Bergman:2007wp  is demonstrated to be thermodynamically
prefered than the QGP phase at high chemical potential/density for low to
moderate temperatures up to the QGP phase transition around trillion Kelvins.
In Ref. bhp , the EoS of the multiquark phase is derived and applied to the
hypothetical multiquark star and NS with multiquark core. Preliminary
estimation allows the existence of NS with multiquark core with masses as
large as $3M_{\odot}$.

With the prospect of massive NS above $2M_{\odot}$ containing quark matter
core from observations, it is interesting to investigate how the holographic
multiquark EoS fit into the observational constraints. In addition to the mass
$M$ and radius $R$, the ratio $M/R$ which represents compactness of the star
also plays a crucial role in determination of the exact EoS of the matter
inside the NS and possibly other kinds of extreme compact objects such as
quark stars. A good probe for the compactness is deformability of the star
parametrized by Love number $k_{2}$ Mora:2003wt ; Berti:2007cd ;
Flanagan:2007ix ; Hinderer:2007mb  and the dimensionless parameter $\Lambda$
defined in Sect. IV. Tidal deformabilities and other properties such as moment
of inertia, sound speed and EoS of hybrid stars in holographic models of quark
and nuclear matter are also calculated in Ref. Annala:2017tqz ; Jokela:2018ers
; Zhang:2019tqd ; Jokela:2020piw . In this work we calculate the Love number
$k_{2}$ and $\Lambda$ of the NS with multiquark core by adapting Chiral
Effective Field Theory (CET) EoS Tews:2012fj  and FYSS EoS FYSS_2017 ;
FYSS_2013 ; FYSS_2011  for the nuclear crust. The deformation parameters of
the multiquark star (MQS) are also calculated for comparison. We found that
for massive NS with multiquark core, the values of deformation parameters are
within the allowed constraints from LIGO/Virgo. Radial pulsation frequencies
of the fundamental and excited modes of the MQ core are calculated and the
zeroth-mode instability is verified to occur at the maximum masses.

This work is organized as the following. Section II reviews the holographic
multiquark model and the corresponding EoS. The adiabatic index and sound
speed of multiquark core are explored in Section III. The tidal deformation of
the NS with multiquark core and MQS are calculated in Section IV. The radial
pulsations of the stars are discussed in Section V. Section VI concludes the
work.

##  II Holographic multiquark and the EoS

According to gauge-gravity duality from superstring theories, bound states of
multiquarks can be described holographically by strings and branes in
D4/D8/$\overline{\text{D8}}$ system called Sakai-Sugimoto (SS) model ss1 ; ss2
. The SS model can capture most features of the QCD e.g. low-energy spectrum
of mesons, confinement/deconfinement transition and chiral symmetry breaking.
In this model, hadrons naturally exist in the confined phase however, another
kind of bound states of quarks can also exist in the deconfined phase at low
and intermediate temperatures when the chemical potential/density is
sufficient large i.e., the multiquark states bch ; bhp  (see e.g. Ref.
Burikham:2011zz  for a concise review of holographic multiquarks).

###  II.1 Equation of state

As explored in details in Ref. bhp , the equations of state for the multiquark
involving pressure $P$ and mass density $\rho$ are given in terms of the
number density $n$ which could be summarized in the dimensionless form as

$\displaystyle P$ $\displaystyle=$ $\displaystyle an^{2}+bn^{4},$
$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\mu_{0}n+an^{2}+\frac{b}{3}n^{4},$ (1)

for small $n$ (“mql”) and

$\displaystyle P$ $\displaystyle=$ $\displaystyle kn^{7/5},$
$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\rho_{c}+\frac{5}{2}P+\mu_{c}\left(n-n_{c}\right)$ (2)
$\displaystyle+kn_{c}^{7/5}-\frac{7k}{2}n_{c}^{2/5}n,$

for large $n$ (“mqh”) respectively. Note that all quantities are dimensionless
and the conversion table to physical quantities is given in the Appendix. The
parameter $\mu_{0}$ is the onset chemical potential of the multiquark phase.
The critical mass density, number density and chemical potential where the EoS
changes from large to small $n$ are denoted by $\rho_{c},n_{c}$ and $\mu_{c}$
respectively. EoS of the deconfined multiquark matter also depends on the
colour charge of the multiquark which is parametrized by $n_{s}$, see Ref. bch
; bhp  for more details. For
$n_{s}=0,\rho_{c}=0.0841077,n_{c}=0.215443,\mu_{c}=0.564374~{}(\text{mqh}),a=1,b=0,\mu_{0}=0.17495~{}(\text{mql})$.
For
$n_{s}=0.3,\rho_{c}=0.0345996,n_{c}=0.086666,\mu_{c}=0.490069~{}(\text{mqh}),a=0.375,b=180.0,\mu_{0}=0.32767~{}(\text{mql})$.
In both cases of $n_{s}=0,0.3$, we found unique value $k=10^{-0.4}$ which
reflects universal behaviour of the MQ at high density regardless of the
colour charges. The SS model has only two free parameters, the number fraction
of hanging strings representing the colour charges $n_{s}$ and the energy
density scale $\epsilon_{s}$ defined in Ref. bch ; bhp . The pressure and mass
density scale with $\epsilon_{s}$ as $P,\rho\sim\epsilon_{s}$. Consequently,
the mass and radius of the pure MQS have the same scaling
$M,R\sim\epsilon_{s}^{-1/2}$. Remarkably, the compactness $M/R$ is thus
independent of $\epsilon_{s}$.

As the density decreases with increasing radial distance from the center of
the star, the confinement phase transition occurs. The gluons and multiquarks
undergo a transition into the colour-singlet baryonic matter. As in Ref.
Pinkanjanarod:2020mgi , we use the CET stiff EoS for the baryonic crust of the
NS. Moreover, the FYSS EoS FYSS_2017 ; FYSS_2013 ; FYSS_2011  at temperature
$T=0.1$ MeV is also implemented as an alternative crust EoS for comparison.
The exact values of $(P,\mu)$ at the confinement phase transition can be found
by studying the $P-\mu$ diagram of the two phases, the star profile as well as
the full mass-radius relations are obtained, see Ref. Pinkanjanarod:2020mgi
for details. In this work, we will focus on the deformability and dynamical
pulsations of such hybrid stars and the MQS.

##  III Adiabatic indices and sound speed

The adiabatic indices describe how pressure of the system changes with the
density. The content of a star could either have the stiff EoS where the
adiabatic index is large or soft EoS when it is small. Large adiabatic index
of matter in the interior could make the star more massive due to the larger
sustaining pressure per mass. Additionally, it represents compressibility of
the matter inside the star which determines the pulsation frequencies and
dynamical stability of the star under perturbation (see Eq. (12)).

Specifically, one kind of adiabatic index could be defined with respect to the
number density $n$ as

$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\frac{n}{P}\left(\frac{dP}{dn}\right).$ (3)

This form of adiabatic index (3) reflects the fact that in a pulsating star,
the number density $n$ varies with the pulsation. However, when considering
the stiffness of the EoS, we shall refer to the conventional relativistic
adiabatic (polytropic) index of the form,
$\gamma=\displaystyle{\frac{\rho}{P}(\frac{dP}{d\rho})}$. For massive NS with
multiquark core at the maximum mass, the range of values of
$(\Gamma,\gamma,c_{s}^{2})$ within the star are given in Table 1. The high-
density multiquark EoS is apparently softer than the low-density EoS with
smaller $\gamma$. The low-density “mql” EoS is considered stiff since
$\gamma>5/3$ while the “mqh” EoS in the core is softer with $\gamma\simeq
1.1$.

| $(\Gamma,\gamma,c_{s}^{2})$
---|---
$\epsilon_{s}$ |
(GeVfm-3) |  |
| mqh |  mql (small to large $r$)
23.2037 |  |
& FYSS, | (1.4, 1.1, 0.43) | $(3.6-2.3,2.6-2.2,0.97-0.05)$
$Y=0.10$ |  |
23.2037 |  |
& FYSS, | (1.4, 1.1, 0.43) | $(3.6-2.3,2.6-2.2,0.97-0.06)$
$Y=0.14$ |  |
26 & CET | (1.4, 1.1, 0.43) | $(3.6-3.1,2.6-2.8,0.97-0.33)$
28 & CET | (1.4, 1.1, 0.43) | $(3.6-3.0,2.6-2.7,0.97-0.26)$
Table 1: Values of $(\Gamma,\gamma,c_{s}^{2})$ for high density multiquark
“mqh” in the inner core and low density “mql” in the outer core with CET/FYSS
nuclear crust at the maximum mass of the massive NS with MQ core.

##  IV Tidal Love number and tidal deformation

###  IV.1 Definitions of tidal Love number and tidal deformation

Consider a static, spherically symmetric star of mass $M$ placed in an
external static quadrupole tidal field $\epsilon_{ij}$. The star will deform
in response and obtain a quadrupole moment $Q_{ij}$. In the star’s local
asymptotic rest frame, the metric component $g_{tt}$ at large $r$ is given by
Thorne:1997kt

$\displaystyle-\frac{(1+g_{tt})}{2}=$ $\displaystyle-$
$\displaystyle\frac{M}{r}-\frac{3Q_{ij}}{2r^{3}}\left(n^{i}n^{j}-\frac{1}{3}\delta^{ij}\right)+O\left(\frac{1}{r^{3}}\right)$
(4) $\displaystyle+$
$\displaystyle\frac{1}{2}\epsilon_{ij}x^{i}x^{j}+O(r^{3})$

where $n^{i}\equiv x^{i}/r$. We are interested in the strong field case. For
an isolated body in a static situation, these moments are uniquely defined:
$\epsilon_{ij}$ and $Q_{ij}$ are the coordinate-independent quantities.
Correspondingly, to the linear order in $\epsilon_{ij}$, the induced
quadrupole will be given in terms of the tidal field by

$\displaystyle Q_{ij}$ $\displaystyle=$ $\displaystyle-\lambda\epsilon_{ij},$
(5)

where $\lambda$ is a constant associated with the $\ell=2$ tidal Love number
$k_{2}$ by Flanagan:2007ix

$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle\frac{3}{2}\lambda
R^{-5}.$ (6)

The multipole moments $Q_{ij}$ and $\epsilon_{ij}$ can be decomposed as a
series of spherical harmonic basis functions

$\displaystyle\epsilon_{ij}$ $\displaystyle=$
$\displaystyle\sum_{\text{m}=-2}^{2}\epsilon_{\text{m}}\mathcal{Y}^{2\text{m}}_{ij},$
$\displaystyle Q_{ij}$ $\displaystyle=$
$\displaystyle\sum_{\text{m}=-2}^{2}Q_{\text{m}}\mathcal{Y}^{2\text{m}}_{ij}.$

The tensor spherical harmonic $\mathcal{Y}^{2\text{m}}_{ij}$, which are
symmetric and traceless, are defined by

$\displaystyle Y_{2m}(\theta,\phi)$ $\displaystyle=$
$\displaystyle\mathcal{Y}^{2\text{m}}_{ij}n^{i}n^{j},$

where $Y_{2m}(\theta,\phi)$ is the spherical harmonic with $\ell=2$ and
$n\equiv(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$. Thus, from Eqs.
(5), the tidal deformation $\lambda$ for a particular mode $m$ ($\ell=2$) can
be written as

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle-\frac{\epsilon_{\text{m}}}{Q_{\text{m}}},$

and in the dimensionless form as Raithel:2018ncd ,

$\displaystyle\Lambda\equiv\frac{\lambda}{M^{5}}\equiv\frac{2}{3}k_{2}\left(\frac{R}{M}\right)^{5}.$
(7)

Note that in this work, we will use the geometric units $G=c=1$ in our
calculation of $k_{2},\Lambda$.

###  IV.2 Calculations of tidal Love number and tidal deformation

In the moderately strong field, e.g. long before the merging of the binary, we
assume the perturbation is small comparing to the background metric
$g^{(0)}_{\alpha\beta}$ which takes the standard spherically symmetric form

$\displaystyle g_{\alpha\beta}$ $\displaystyle=$ $\displaystyle
g^{(0)}_{\alpha\beta}+h_{\alpha\beta},$ $\displaystyle g^{(0)}_{\alpha\beta}$
$\displaystyle=$
$\displaystyle\text{diag}[-e^{\nu(r)},e^{\lambda(r)},r^{2},r^{2}\sin^{2}\theta],$
(8)

where $h_{\alpha\beta}$ is a linearized metric perturbation. We only consider
the $\ell=2$, static, even-parity perturbations in the Regge-Wheeler gauge
Regge:1957td ; Hinderer:2009ca . The linearized metric perturbation
$h_{\alpha\beta}$ could then be expressed in the diagonal form as
Thorne_Campolattaro ; Hinderer:2009ca

$\displaystyle h_{\alpha\beta}$
$\displaystyle=-\text{diag}[e^{\nu(r)}H(r),e^{\lambda(r)}H(r),r^{2}K(r),$ (9)
$\displaystyle r^{2}K(r)\text{sin}^{2}\theta]Y_{2m}(\theta,\phi).$

Figure 1: Tidal Love number $k_{2}$ versus mass of MQS/NS with MQ core. Figure
2: Dimensionless tidal deformation $\Lambda$ versus mass of MQS/NS with MQ
core.

The nonzero components of the perturbations of the stress-energy tensor are
$\delta T^{0}_{0}=-\delta\rho=-(dP/d\rho)^{-1}\delta P$ and $\delta
T^{i}_{i}=\delta P$. Insert these relations and the metric perturbation Eq.
(9) into the linearized Einstein equation $\delta
G_{\alpha}^{\beta}=8\pi\delta T_{\alpha}^{\beta}$, and after elimination of
$K$, we obtain the following differential equation Hinderer:2009ca

$\displaystyle H^{\prime\prime}$ $\displaystyle+$
$\displaystyle\frac{2}{r}\left(1-2\frac{m}{r}\right)^{-1}H^{\prime}\left\\{1-\frac{m}{r}-2\pi
r^{2}(\rho-P)\right\\}$ (10) $\displaystyle+$ $\displaystyle
2\left(1-2\frac{m}{r}\right)^{-1}H\left\\{2\pi\left(5P+9\rho+\frac{\rho+P}{dP/d\rho}\right)\right.$
$\displaystyle-$
$\displaystyle\left.\frac{3}{r^{2}}-2\left(1-2\frac{m}{r}\right)^{-1}\left(\frac{m}{r^{2}}+4\pi
rP\right)^{2}\right\\}=0,$

where ′ denotes $d/dr$ and $m=m(r)$ is the accumulated mass of the star. For a
given EoS, the Tolman-Oppenheimer-Volkoff (TOV) equation can be solved as the
leading order of the Einstein’s field equations. The resulting solutions are
taken as the star profile to be used to subsequently solve for the
perturbation $H(r)$ at the next-to-leading order of the field equations.

The differential equation of $H(r)$ would be integrated radially outward
starting from the center with boundary conditions: $H(r)=a_{0}r^{2}$ and
$H^{\prime}(r)=2a_{0}r$ as $r\rightarrow 0$. The constant $a_{0}$ indicates
how much the star is deformed due to the given external field and could be
chosen arbitrarily as it cancels out in the calculation of tidal Love number.
Outside the star, where $T_{\mu\nu}=0$, the general solution for $H(r)$ at
large $r$ can be written in terms of associated Legendre function
$Q_{2}^{2}(r/M-1)\sim r^{-3}$ and $P_{2}^{2}(r/M-1)\sim r^{-2}$. By matching
the interior and exterior solutions of $H(r)$ and their derivatives at the
star surface $r=R$, this solution is uniquely determined. Equation (4) can
then be used to identify the coefficients of the external solution with the
axisymmetric tidal field and quadrupole moment via $\epsilon
Y_{20}(\theta,\phi)=\epsilon_{ij}n^{i}n^{j}$, and
$QY_{20}(\theta,\phi)=Q_{ij}n^{i}n^{j}=-\lambda\epsilon_{ij}n^{i}n^{j}$, where
$\epsilon$ and $Q$ are the coefficients of the $\ell=2,m=0$ spherical harmonic
of the tidal tensor and quadrupole moment respectively. Finally, the $\ell=2$
tidal Love number is given by Hinderer:2009ca

$\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle\frac{8C^{5}}{5}\left(1-2C\right)^{2}\left[2+2C(y-1)-y\right]$
$\displaystyle\times$ $\displaystyle\\{2C\left[6-3y+3C(5y-8)\right]$
$\displaystyle+$ $\displaystyle 4C^{3}\left[13-11y+C(3y-2)+2C^{2}(1+y)\right]$
$\displaystyle+$ $\displaystyle
3(1-2C)^{2}\left[2-y+2C(y-1)\right]\text{ln}(1-2C)\\}^{-1},$

where the star’s compactness parameter $C\equiv M/R$ and the quantity $y\equiv
RH^{\prime}(R)/H(R)$, which can be obtained by numerical integration of Eq.
(10) in the region $0<r\leq R$.

EoS | $M(M_{\odot})$ |  $R$(km) | $k_{2}$ | $\Lambda$
---|---|---|---|---
| 1.8 | 14.1 | 0.106 | 302
| 1.9 | 14.2 | 0.100 | 255
mq+CET, $\epsilon_{s}=26$ |  |  |  |
| 2.0 | 13.8 | 0.0752 | 113
| 2.2 | 12.2 | 0.0285 | 14.0
| 1.8 | 12.8 | 0.0619 | 108
mq+CET, $\epsilon_{s}=28$ | 1.9 | 12.3 | 0.0453 | 50.0
| 2.0 | 11.8 | 0.0328 | 22.7
| 1.5 | 12.0 | 0.100 | 313
mq+FYSS, $Y=0.10$ |  |  |  |
$\epsilon_{s}=23.2037$ | 1.9 | 11.9 | 0.0577 | 51.8
| 1.5 | 11.9 | 0.0956 | 290
mq+FYSS, $Y=0.14$ |  |  |  |
$\epsilon_{s}=23.2037$ | 1.9 | 11.8 | 0.0550 | 48.3
Table 2: Tidal love number $k_{2}$ and dimensionless tidal deformation $\Lambda$ of NS with multiquark core and nuclear CET/FYSS crust. EoS | $M_{1}(M_{\odot})$ | $M_{2}(M_{\odot})$ | $\Lambda_{1}$ | $\Lambda_{2}$ | $\tilde{\Lambda}$
---|---|---|---|---|---
| 1.7 | 1.7 | 418 | 418 | 418
mq+CET, $\epsilon_{s}=26$ | 2.0 | 1.4 | 113 | 1130 | 391
| 1.87 | 0.93 | 247 | 7980 | 1280
| 1.7 | 1.7 | 260 | 260 | 260
mq+CET, $\epsilon_{s}=28$ | 2.0 | 1.4 | 22.7 | 1120 | 313
| 1.87 | 0.93 | 62.3 | 8080 | 1080
| 1.7 | 1.7 | 128 | 128 | 128
mq+FYSS, $Y=0.10$, | 2.0 | 1.4 | 32.2 | 498 | 158
$\epsilon_{s}=23.2037$ |  |  |  |  |
| 1.87 | 0.93 | 59.5 | 5670 | 776
| 1.7 | 1.7 | 118 | 118 | 118
mq+FYSS, $Y=0.14$, | 2.0 | 1.4 | 30.2 | 459 | 146
$\epsilon_{s}=23.2037$ |  |  |  |  |
| 1.87 | 0.93 | 55.4 | 5310 | 727
Table 3: Dimensionless tidal deformation $\Lambda_{1}$, $\Lambda_{2}$ and
combined dimensionless tidal deformation $\tilde{\Lambda}$ of binary systems
of NS with multiquark core and nuclear CET/FYSS crust.

For the binary system, each mass is deformed by companion star’s gravitational
field, the combined deformation parameter of the binary $M_{1},M_{2}$ is
defined as

$\tilde{\Lambda}\equiv\frac{16}{13}\frac{(M_{1}+12M_{2})M_{1}^{4}\Lambda_{1}+(M_{2}+12M_{1})M_{2}^{4}\Lambda_{2}}{(M_{1}+M_{2})^{5}}.$

The value of $k_{2}$ for MQS and NS with multiquark core is plotted with
respect to mass in Fig. 1. The dimensionless $\Lambda$ is plotted in Fig. 2.
We consider the MQS with $\epsilon_{s}=23.2037,26,28$ GeV/fm3 and the NS with
multiquark core continuing to nuclear CET (FYSS) crust for
$\epsilon_{s}=26,28~{}(23.2037)$ GeV/fm3. In Fig. 2, the Gravitational Waves
events detected by LIGO/Virgo, GW170817 and GW190425, are also presented. The
lower and upper mass limits of both the binary masses are given by
($1.17M_{\odot},1.60M_{\odot}$) for GW170817 at 90% confidence interval. This
interval is converted to $1-\sigma$ or 68.27% confidence interval
($1.25M_{\odot},1.52M_{\odot}$) in Fig. 2. Similarly for GW190425, the 90%
confidence interval ($1.46M_{\odot},1.87M_{\odot}$) becomes
($1.54M_{\odot},1.79M_{\odot}$) for $1-\sigma$ 68.27% confidence interval.
Constraints on $\Lambda$ for each gravitational event Abbott:2020uma ;
TheLIGOScientific:2017qsa  are presented as upper limits in Fig. 2, i.e.,
$\Lambda<800$ for GW170817 and $\Lambda<600$ for GW190425. Table 4 shows
$k_{2},\Lambda$ for NS with CET/FYSS crust and multiquark core for certain
large masses relevant to detections of Gravitational-wave Observatories
Abbott:2020uma ; TheLIGOScientific:2017qsa ; GWcollab . Table 3 explores the
combined $\tilde{\Lambda}$ for binaries with masses possibly in the mass range
constrained by LIGO/Virgo Abbott:2020uma ; TheLIGOScientific:2017qsa . The
loose constraints are $\tilde{\Lambda}<600,700-800$ for the binary with masses
$(1.5-1.7,1.6-1.9)M_{\odot}$ (GW190425, low-spin prior) and
$(0.9-1.4,1.4-2.3)M_{\odot}$ (GW170817) respectively. We choose the mass range
of each companion in the binary to be larger than $1.7M_{\odot}$ so that it
has multiquark core.

Yagi and Yunes discovered separate universal classes of compact stars with
nuclear and quark matter EoS using the characteristic relationship between the
moment of inertia, (rotational) Love number and the quadrupole moment
Yagi:2013bca . It is interesting to see which class the NS with multiquark
core/MQS would fit into. Instead of using the rotational deformability
parameter, we plot relationship between $\Lambda$ and the compactness $C\equiv
M/R$ in Fig. 3.

Figure 3: Dimensionless tidal deformation $\Lambda$ versus compactness $C=M/R$
of MQS/NS with MQ core.

Similar to the rotational deformation case (between typical nuclear NS and
strange quark stars), two separate deformability classes of pure MQS and pure
nuclear NS curves exist as shown in Fig. 3. Pure MQS has lower $\Lambda$ at a
given $C$ than the pure nuclear NS and the hybrid star interpolates between
the two classes. The inverted triangle and diamond plot in Fig. 3 show the NS
with MQ core transiting from pure MQS curve to the CET curve as $C$ decreases.
NS with high-$Y(\equiv n_{p}/n_{b}=0.10,0.14$, where $n_{p}~{}(n_{b})$ is the
number density of proton (baryon) respectively) FYSS crust sides more with the
nuclear CET neutron star curve. Remarkably due to the independence of
$\epsilon_{s}$ of the compactness $C$, all MQSs with any $\epsilon_{s}$
degenerate on the same curve.

##  V Linearized adiabatic radial pulsations of the holographic multiquark
star

Infinitesimal radial pulsations of star were first discussed under general
relativity framework by Chandrasekhar Chandrasekhar:1964zz . Following the
rewritten form in Ref. Chanmugam ; Gondek:1997fd , there are 2 crucial
quantities describing the infinitesimal radial pulsations: the relative radial
displacement, $\xi=\Delta r/r$, where $\Delta r$ is the radial displacement of
the matter content of the star, and $\Delta P$, the corresponding Lagrangian
perturbation of the pressure. The two essential differential equations
governing the radial pulsation can be expressed as

$\displaystyle\frac{d\xi}{dr}$ $\displaystyle=$
$\displaystyle-\frac{1}{r}(3\xi+\frac{\Delta P}{\Gamma
P})-\frac{dP}{dr}\frac{\xi}{(P+\rho)},$ (12) $\displaystyle\frac{d\Delta
P}{dr}$ $\displaystyle=$
$\displaystyle\xi\left\\{\omega^{2}re^{\lambda-\nu}(P+\rho)-4\frac{dP}{dr}\right\\}$
(13) $\displaystyle+$
$\displaystyle\xi\left\\{\left(\frac{dP}{dr}\right)^{2}\frac{r}{(P+\rho)}-8\pi
re^{\lambda}(P+\rho)P\right\\}$ $\displaystyle+$ $\displaystyle\Delta
P\left\\{\frac{dP}{dr}\frac{1}{(P+\rho)}-4\pi re^{\lambda}(P+\rho)\right\\},$

where $\Gamma$ is the adiabatic index defined in Section III, $\omega$ is the
eigenfrequency of the radial oscillation and the $\xi$ and $\Delta P$ are
assumed to have a harmonic time dependence $\propto e^{-i\omega t}$. A key
feature of this system of radial pulsation equations, Eqs. (12) and (13), is
that they do not involve any derivatives of the adiabatic index, $\Gamma$.

There is subtlety in the calculation of second-order perturbation (13)
involving $\nu(r)$ from the zeroth-order star profile. We need to choose
$\nu(r=0)$ such that

$e^{\nu(R)}=1-\frac{2M}{R},$ (14)

at the star surface. Even though the Einstein’s field equation has shift
symmetry under $\nu\to\nu+\nu_{0}$ for arbitrary constant $\nu_{0}$ at the
zeroth order, the gauge fixing of $\nu$ according to (14) is required for the
calculation of first and second-order perturbations. This is the matching
condition to the Schwarzschild metric in the exterior region of the star. For
$\lambda(r)$, we have used the substitution

$e^{-\lambda(r)}=1-\frac{2m(r)}{r},$ (15)

in every perturbation equation so that the matching to the Schwarzschild
metric is thus guaranteed.

To solve the two coupled linear differential equations, Eqs. (12) and (13),
the boundary conditions at the center and surface of the star must be
specified. The first natural condition is the coefficient of the $1/r$-term in
Eq. (12) must vanish to avoid divergence as $r\to 0$,

$\displaystyle\Delta P(r=0)=-3(\xi\Gamma P)|_{r=0}.$ (16)

The relative radial displacement is normalized to unity at the center,
$\xi(0)=1$. The surface of the star $R$ is defined by $P(r=R)=0$. This also
leads to

$\displaystyle\Delta P(r=R)=0.$ (17)

After solving the TOV for the star profile for a given EoS, the two coupled
linear equations can be subsequently solved by the shooting method for the
eigenvalue $\omega^{2}$ and the perturbation profile $\xi(r)$ and $\Delta
P(r)$ satisfying the boundary conditions, Eqs. (16) and (17). For a given EoS
model of static star, a set of the eigenvalues;
$\omega_{0}^{2}<\omega_{1}^{2}<\cdots<\omega_{n}^{2}<\cdots$, with
corresponding eigenfunctions; $\xi_{0},\xi_{1},\cdots,\xi_{n},\cdots,$ could
be obtained. The characteristic standing-wave eigenfunction $\xi_{n}$ has $n$
nodes within the star. The eigenfunctions of pressure perturbation $\Delta
P_{n}(r)$ are simultaneously obtained.

For $\omega_{n}^{2}>0$ for all $n$, the static stellar model is stable with
respect to small radial adiabatic perturbations. If the lowest eigenfrequency
$\omega_{0}=0$, the configuration is on the verge of stability/instability
Cox1980 . Notably, the star configuration can be stable with respect to a set
of higher-modes pulsations and become unstable for lower modes. In this work,
only the zeroth-mode instability is studied.

###  V.1 Eigenfrequencies and instabilities

According to the radial pulsation Eqs. (12) and (13) together with the TOV
equation, the eigenfrequencies could be determined numerically by the shooting
process and they have been plotted in Fig. 4 and Fig. 5. In the SS model, the
frequencies of the multiquark-core oscillations scale with
$\sqrt{\epsilon_{s}}$. The zeroth-mode frequency is about $2.5$ kHz for
massive NS with masses around $2M_{\odot}$ for $\epsilon_{s}=23.2037$ GeV/fm3.
Instability with respect to the zeroth mode appears when the central density
exceeds the critical value at the maximum mass and enters into the unstable
branch of the $M-R$ diagram as expected. For NS with masses above
$2.15M_{\odot}$ (for $\epsilon_{s}=23.2037$ GeV/fm3), higher-mode oscillations
show distinct wiggly pattern of frequencies with respect to the mass due to
the interplay between perturbing waves in the “mqh” and “mql” structures.

Figure 4: The eigenfrequencies $f_{n}=\omega_{n}/2\pi$ of $n=0,1,2,3,4,5$
modes versus mass of multiquark core at energy density scale
$\epsilon_{s}=23.2037$ GeV/fm3. Solid (Dashed) line represents the frequency
of the MQ core in the stable (unstable) branch. Figure 5: The eigenfrequencies
$f_{n}=\omega_{n}/2\pi$ of $n=0,1,2,3,4,5$ modes versus compactness $C=M/R$ of
multiquark core at energy density scale $\epsilon_{s}=23.2037$ GeV/fm3. Solid
(Dashed) line represents the frequency of the MQ core in the stable (unstable)
branch.

##  VI Conclusions and Discussions

The deformation parameters of the NS with multiquark core (or hybrid star) and
MQS have been calculated in the holographic SS model. With the nuclear CET and
FYSS crust EoS, the massive NS with multiquark core is consistent with
observational constraints from the LIGO/Virgo on the mass, radius
Pinkanjanarod:2020mgi , and deformation parameter $k_{2},\Lambda$.
Deformability of MQ-matter star is clearly distinguishable from the baryonic-
matter star when plotted against the compactness $C$. Hybrid stars such as the
massive NS with MQ core has the deformability $\Lambda-C$ profile that
interpolates between the pure MQS and the pure baryonic NS.

Characteristic eigenfrequencies of the radial pulsation of the multiquark core
are calculated for the entire mass range for $n=0-5$ modes. The fundamental
mode frequency is approximately $2.5$ kHz near the maximum masses around
$2M_{\odot}$. A unique characteristic of the fundamental frequency is that it
is slowly increasing for a wide mass range up until approaching the maximum
mass where it drops abruptly to small value as shown in Fig. 4. For $M_{\rm
core}$ close to the maximum mass, the higher modes have more wiggly pattern
which is originated from the interplay between radial waves in the high
(“mqh”) and low (“mql”) density region of the MQ core.

The distinctive pattern in $\Lambda$ vs. mass (and compactness) where there is
kink at the transition between MQ core and nuclear crust could help
identifying the existence of MQ core in the massive NS. The sets of pulsation
frequencies of the MQ core can be compared with observations and EoS of other
models and reveal the properties of nuclear matter at extreme densities. For
example, whether the multiquarks could form and remain in the fluid/gas phase
or diquark could form superfluid below certain CSC temperature or both kinds
of multiquark phases, MQ and CSC, could be found under the appropriate
situations inside the NS (see e.g. Ref. BitaghsirFadafan:2018iqr ;
BitaghsirFadafan:2020otb  for holographic approach to CSC where the tidal
deformability is also computed). It is also possible that MQ is prefered in
the most dense and warm temperature region in the core, then in the less dense
and cooler layer CSC phase is prefered, and subsequently the confined nuclear
is prefered in the crust. Constraints from the observation of massive NS
around $2M_{\odot}$ suggest that in such case CSC EoS should be relatively
comparable to the MQ EoS. It would be interesting to explore such possibility
of NS with mixed phase in the future work.

###### Acknowledgements.

S.P. (first author) is supported in part by the Second Century Fund: C2F PhD
Scholarship, Chulalongkorn University. S.P. is supported by Rachadapisek
Sompote Fund for Postdoctoral Fellowship, Chulalongkorn University.

##  VII Appendix

quantity | dimensionless variable | physical variable
---|---|---
pressure | $P$ | $\epsilon_{s}P$
density | $\rho$ | $\displaystyle{\frac{\epsilon_{s}}{c^{2}}\rho}$
chemical potential | $\mu$ |  $\epsilon_{s}\mu$ fm3
number density | $n$ |  $n$ fm-3
mass | $M$ | $\displaystyle{\frac{c^{4}}{\sqrt{G^{3}\epsilon_{s}}}M}$
radius | $r$ | $\displaystyle{\frac{c^{2}}{\sqrt{G\epsilon_{s}}}r}$
angular frequency | $\omega$ | $\displaystyle{\frac{\sqrt{G\epsilon_{s}}}{c}\omega}$
Table 4: Conversion table of physical quantities for $\epsilon_{s}$ in GeVfm-3
units.

## References

  * (1) J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon and T. Driebe, et al. Science 340, 6131 (2013) doi:10.1126/science.1233232 [arXiv:1304.6875 [astro-ph.HE]].
  * (2) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, W. C. G. Ho, J. M. Lattimer and M. Loewenstein, et al. [arXiv:2105.06979 [astro-ph.HE]].
  * (3) J. S. Clark, S. P. Goodwin, P. A. Crowther, L. Kaper, M. Fairbairn, N. Langer and C. Brocksopp, “Physical parameters of the high-mass X-ray binary 4U1700-37,” Astron. Astrophys. 392, 909-920 (2002) doi:10.1051/0004-6361:20021184 [arXiv:astro-ph/0207334 [astro-ph]].
  * (4) R. W. Romani, A. V. Filippenko, J. M. Silverman, S. B. Cenko, J. Greiner, A. Rau, J. Elliott and H. J. Pletsch, “PSR J1311-3430: A Heavyweight Neutron Star with a Flyweight Helium Companion,” Astrophys. J. Lett. 760, L36 (2012) doi:10.1088/2041-8205/760/2/L36 [arXiv:1210.6884 [astro-ph.HE]].
  * (5) R. W. Romani, “2FGL J1311.7-3429 Joins the Black Widow Club,” Astrophys. J. Lett. 754, L25 (2012) doi:10.1088/2041-8205/754/2/L25 [arXiv:1207.1736 [astro-ph.HE]].
  * (6) M. H. van Kerkwijk, R. Breton and S. R. Kulkarni, “Evidence for a Massive Neutron Star from a Radial-Velocity Study of the Companion to the Black Widow Pulsar PSR B1957+20,” Astrophys. J. 728, 95 (2011) doi:10.1088/0004-637X/728/2/95 [arXiv:1009.5427 [astro-ph.HE]].
  * (7) M. Linares, T. Shahbaz and J. Casares, “Peering into the dark side: Magnesium lines establish a massive neutron star in PSR J2215+5135,” Astrophys. J. 859, no.1, 54 (2018) doi:10.3847/1538-4357/aabde6 [arXiv:1805.08799 [astro-ph.HE]].
  * (8) V. Bhalerao, M. H. van Kerkwijk and F. Harrison, “Constraints on the Compact Object Mass in the Eclipsing HMXB XMMU J013236.7+303228 in M33,” Astrophys. J. 757, 10 (2012) doi:10.1088/0004-637X/757/1/10 [arXiv:1207.0008 [astro-ph.SR]].
  * (9) D. J. Nice, E. M. Splaver, I. H. Stairs, O. Loehmer, A. Jessner, M. Kramer and J. M. Cordes, “A 2.1 solar mass pulsar measured by relativistic orbital decay,” Astrophys. J. 634, 1242-1249 (2005) doi:10.1086/497109 [arXiv:astro-ph/0508050 [astro-ph]].
  * (10) P. Demorest, T. Pennucci, S. Ransom, M. Roberts and J. Hessels, “Shapiro Delay Measurement of A Two Solar Mass Neutron Star,” Nature 467, 1081-1083 (2010) doi:10.1038/nature09466 [arXiv:1010.5788 [astro-ph.HE]].
  * (11) P. C. C. Freire, “Super-Massive Neutron Stars,” AIP Conf. Proc. 983, no.1, 459-463 (2008) doi:10.1063/1.2900274 [arXiv:0712.0024 [astro-ph]].
  * (12) H. Quaintrell, A. J. Norton, T. D. C. Ash, P. Roche, B. Willems, T. R. Bedding, I. K. Baldry and R. P. Fender, “The mass of the neutron star in Vela X-1 and tidally induced non-radial oscillations in GP Vel,” Astron. Astrophys. 401, 313-324 (2003) doi:10.1051/0004-6361:20030120 [arXiv:astro-ph/0301243 [astro-ph]].
  * (13) E. Annala, T. Gorda, A. Kurkela, J. Nättilä and A. Vuorinen, “Evidence for quark-matter cores in massive neutron stars,” Nature Phys. (2020) doi:10.1038/s41567-020-0914-9 [arXiv:1903.09121 [astro-ph.HE]].
  * (14) M. P. Lombardo, “Lattice QCD at finite density: Imaginary chemical potential,”PoS CPOD2006, 003 (2006)[arXiv:hep-lat/0612017 [hep-lat]]
  * (15) M. G. Alford, A. Schmitt, K. Rajagopal and T. Schäfer, Rev. Mod. Phys. 80 (2008), 1455-1515 doi:10.1103/RevModPhys.80.1455 [arXiv:0709.4635 [hep-ph]].
  * (16) S. S. Gubser, I. R. Klebanov and A. M. Polyakov, “Gauge theory correlators from noncritical string theory,” Phys. Lett. B 428 (1998), 105-114 doi:10.1016/S0370-2693(98)00377-3 [arXiv:hep-th/9802109 [hep-th]].
  * (17) Juan M. Maldacena, “Wilson loops in large N field theories,” Phys. Rev. Lett. 80 (1998) 4859-4862, [arXiv:hep-th/9803002].
  * (18) E. Witten, “Anti-de Sitter space, thermal phase transition, and conﬁnement in gauge theories,” Adv. Theor. Math. Phys. 2, 505 (1998) [arXiv:hep-th/9803131].
  * (19) E. Witten, JHEP 07 (1998), 006 doi:10.1088/1126-6708/1998/07/006 [arXiv:hep-th/9805112 [hep-th]].
  * (20) C. G. Callan, Jr., A. Guijosa and K. G. Savvidy, Nucl. Phys. B 547 (1999), 127-142 doi:10.1016/S0550-3213(99)00057-7 [arXiv:hep-th/9810092 [hep-th]].
  * (21) A. Brandhuber, N. Itzhaki, J. Sonnenschein and S. Yankielowicz, JHEP 07 (1998), 020 doi:10.1088/1126-6708/1998/07/020 [arXiv:hep-th/9806158 [hep-th]].
  * (22) T. Nakas and K. S. Rigatos, JHEP 12 (2020), 157 doi:10.1007/JHEP12(2020)157 [arXiv:2010.00025 [hep-th]].
  * (23) P. Burikham, A. Chatrabhuti and E. Hirunsirisawat, ”Exotic multi-quark states in the deconfined phase from gravity dual models,” JHEP 05 (2009) 006, [arXiv:0811.0243 [hep-ph]].
  * (24) T. Sakai, S. Sugimoto, “Low energy hadron physics in holographic QCD,” Prog. Theor. Phys. 113, 843-882 (2005) [hep-th/0412141].
  * (25) T. Sakai, S. Sugimoto, “More on a holographic dual of QCD,” Prog. Theor. Phys. 114, 1083-1118 (2005) [hep-th/0507073].
  * (26) O. Aharony, J. Sonnenschein, S. Yankielowicz, “A holographic model of deconfinement and chiral symmetry restoration”,Ann. Phys. 322 (2007) 1420. arXiv:hep-th/0604161.
  * (27) O. Bergman, G. Lifschytz and M. Lippert, “Holographic Nuclear Physics,” JHEP 11 (2007), 056 doi:10.1088/1126-6708/2007/11/056 [arXiv:0708.0326 [hep-th]].
  * (28) P. Burikham, E. Hirunsirisawat and S. Pinkanjanarod, ”Thermodynamic properties of holographic multiquark and the multiquark star,” JHEP 06 (2010) 040, [arXiv:1003.5470 [hep-ph]].
  * (29) T. Mora and C. M. Will, “A PostNewtonian diagnostic of quasiequilibrium binary configurations of compact objects,” Phys. Rev. D 69 (2004), 104021 [erratum: Phys. Rev. D 71 (2005), 129901] doi:10.1103/PhysRevD.71.129901 [arXiv:gr-qc/0312082 [gr-qc]].
  * (30) E. Berti, S. Iyer and C. M. Will, “A Post-Newtonian diagnosis of quasiequilibrium configurations of neutron star-neutron star and neutron star-black hole binaries,” Phys. Rev. D 77 (2008), 024019 doi:10.1103/PhysRevD.77.024019 [arXiv:0709.2589 [gr-qc]].
  * (31) E. E. Flanagan and T. Hinderer, “Constraining neutron star tidal Love numbers with gravitational wave detectors,” Phys. Rev. D 77 (2008), 021502 doi:10.1103/PhysRevD.77.021502 [arXiv:0709.1915 [astro-ph]].
  * (32) T. Hinderer, “Tidal Love numbers of neutron stars,” Astrophys. J. 677 (2008), 1216-1220 doi:10.1086/533487 [arXiv:0711.2420 [astro-ph]].
  * (33) E. Annala, C. Ecker, C. Hoyos, N. Jokela, D. Rodríguez Fernández and A. Vuorinen, JHEP 12 (2018), 078 doi:10.1007/JHEP12(2018)078 [arXiv:1711.06244 [astro-ph.HE]].
  * (34) N. Jokela, M. Järvinen and J. Remes, JHEP 03 (2019), 041 doi:10.1007/JHEP03(2019)041 [arXiv:1809.07770 [hep-ph]].
  * (35) K. Zhang, T. Hirayama, L. W. Luo and F. L. Lin, Phys. Lett. B 801 (2020), 135176 doi:10.1016/j.physletb.2019.135176 [arXiv:1902.08477 [hep-ph]].
  * (36) N. Jokela, M. Järvinen, G. Nijs and J. Remes, Phys. Rev. D 103 (2021) no.8, 086004 doi:10.1103/PhysRevD.103.086004 [arXiv:2006.01141 [hep-ph]].
  * (37) 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 (2013) no.3, 032504 doi:10.1103/PhysRevLett.110.032504 [arXiv:1206.0025 [nucl-th]].
  * (38) Shun Furusawa, Kohsuke Sumiyoshi, Shoichi Yamada, Hideyuki Suzuki., Nucl. Phys. A  957, 188. (2017).
  * (39) S. Furusawa, K. Sumiyoshi, S. Yamada, and H. Suzuki, Astrophys. J. 772,95 (2013).
  * (40) S. Furusawa, S. Yamada, K. Sumiyoshi, and H. Suzuki, Astrophys. J 738,178 (2011).
  * (41) P. Burikham and E. Hirunsirisawat, “Holographic multiquarks in the quark-gluon plasma: A review,” Adv. High Energy Phys. 2011 (2011), 123184 doi:10.1155/2011/123184
  * (42) S. Pinkanjanarod and P. Burikham, “Massive neutron stars with holographic multiquark cores,” [arXiv:2007.10615 [nucl-th]].
  * (43) K. S. Thorne, “Tidal stabilization of rigidly rotating, fully relativistic neutron stars,” Phys. Rev. D 58, 124031 (1998) doi:10.1103/PhysRevD.58.124031 [arXiv:gr-qc/9706057 [gr-qc]].
  * (44) C. Raithel, F. Özel and D. Psaltis, “Tidal deformability from GW170817 as a direct probe of the neutron star radius,” Astrophys. J. Lett. 857, no.2, L23 (2018) doi:10.3847/2041-8213/aabcbf [arXiv:1803.07687 [astro-ph.HE]].
  * (45) T. Regge and J. A. Wheeler, “Stability of a Schwarzschild singularity,” Phys. Rev. 108, 1063-1069 (1957) doi:10.1103/PhysRev.108.1063
  * (46) T. Hinderer, B. D. Lackey, R. N. Lang and J. S. Read, “Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral,” Phys. Rev. D 81, 123016 (2010) doi:10.1103/PhysRevD.81.123016 [arXiv:0911.3535 [astro-ph.HE]].
  * (47) K. Thorne and A. Campolattaro, “Non-Radial Pulsation of General-Relativistic Stellar Models. I. Analytic Analysis for L $\geq 2$ ” Astrophys. J. 149, 591
  * (48) B. P. Abbott et al. [LIGO Scientific and Virgo], “GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$,” Astrophys. J. Lett. 892 (2020) no.1, L3 doi:10.3847/2041-8213/ab75f5 [arXiv:2001.01761 [astro-ph.HE]].
  * (49) B. P. Abbott et al. [LIGO Scientific and Virgo],“GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,”Phys. Rev. Lett. 119, no.16, 161101 (2017) doi:10.1103/PhysRevLett.119.161101[arXiv:1710.05832 [gr-qc]].
  * (50) R. Abbott et al. [LIGO Scientific, Virgo and KAGRA], “Observation of gravitational waves from two neutron star-black hole coalescences,” Astrophys. J. Lett. 915, L5 (2021) doi:10.3847/2041-8213/ac082e
  * (51) K. Yagi and N. Yunes, “I-Love-Q,” Science 341 (2013), 365-368 doi:10.1126/science.1236462 [arXiv:1302.4499 [gr-qc]].
  * (52) S. Chandrasekhar, “The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Astrophys. J. 140, 417-433 (1964) [erratum: Astrophys. J. 140, 1342 (1964)] doi:10.1086/147938
  * (53) G. Chanmugam, “Radial oscillations of zero-temperature white dwarfs and neutron stars below nuclear densities,” Astrophys. J. 217(3),799-808
  * (54) D. Gondek, P. Haensel and J. L. Zdunik, “Radial pulsations and stability of protoneutron stars,” Astron. Astrophys. 325, 217-227 (1997) [arXiv:astro-ph/9705157 [astro-ph]].
  * (55) J. Cox, “Theory of stellar pulsations,” Princeton University Press, Princeton. (1980)
  * (56) K. Bitaghsir Fadafan, J. Cruz Rojas and N. Evans, Phys. Rev. D 98 (2018) no.6, 066010 doi:10.1103/PhysRevD.98.066010 [arXiv:1803.03107 [hep-ph]].
  * (57) K. Bitaghsir Fadafan, J. Cruz Rojas and N. Evans, Phys. Rev. D 103 (2021) no.2, 026012 doi:10.1103/PhysRevD.103.026012 [arXiv:2009.14079 [hep-ph]].