text:
\subsection{The spontaneous scalarization mechanism}
\label{sec:scalarization_mechanism}

\subsubsection{Tachyonic instability and nonlinear quenching}
\label{tachyonquenched}

Before we discuss spontaneous scalarization in the context of gravity, it is
instructive to review the dynamics of a real scalar field $\varphi$ with a
quartic self-interaction in Minkowski spacetime.
The Lagrangian for this field is,
\begin{equation}
\label{phi4L}
{\cal L}=\frac{1}{2} \eta^{\mu\nu} \partial_\mu \varphi\partial_\mu \varphi+V(\varphi),
\end{equation}
where $\eta_{\mu\nu}$ is the Minkowski metric,
\begin{equation}
\label{eq:phi4P}
V(\varphi)=\frac{1}{2}\mu^2\varphi^2+\frac{1}{4}\lambda \varphi^4,
\end{equation}
$\mu$ is the bare mass, and $\lambda$ is a coupling constant.
The scalar then satisfies the following field equation
\begin{equation}
\label{eq:phi4}
\Box_{\eta} \, \varphi + \mu^2\varphi + \lambda \varphi^3=0,
\end{equation}
where $\Box_{\eta} \equiv \eta^{\mu\nu} \partial_{\mu} \partial_{\nu}$ is the
flat-spacetime d'Alembertian.
Clearly, $\varphi=0$ is a solution of this equation.
Consider now small perturbations, $\delta \varphi$, around $\varphi=0$.
By linearizing~Eq.~\eqref{eq:phi4}, we find that $\delta \varphi$ obeys
\begin{equation}
\label{eq:phi4lin}
\Box_{\eta} \, \delta\varphi-\mu^2\delta\varphi=0.
\end{equation}
The corresponding dispersion relation is $\omega^2=k^2+\mu^2$, where $\omega$
is the frequency and $k$ is the wavenumber.
If $\mu^2>0$, the solutions to this equation are plane waves and the perturbations
decay. If instead $\mu^2<0$, one encounters a \emph{tachyonic} instability and
the perturbations with \emph{small wavenumber} exhibit exponential growth.

This exponential growth seems catastrophic at first sight, but it does not have to be.
As $\varphi$ grows, the linear approximation used above will quickly become
invalid, and the nonlinear self-interaction $\lambda \varphi^3$ will become
important. It will be this interaction that will determine the endpoint of the
instability.
Assume that $\lambda>0$, in which case the potential has the
well-known ``Mexican hat'' shape. Then, Eq.~\eqref{eq:phi4} will admit a second
solution with constant $\varphi$, which we denote as $\varphi_{\rm min}$ as it
will correspond to the minimum of the potential. Eq.~\eqref{eq:phi4} implies that
$\varphi^2_{\rm min}=-\mu^2/\lambda$.
Thus, the tachyonic instability simply drives the scalar field away from the
unstable local maximum of the potential and towards a stable minimum. This is
sometimes referred to as tachyon condensation and associated to a phase
transition of the system.

The key message from this simple example is that linearized perturbations
around the unstable maximum capture the onset of the tachyonic instability but
they are oblivious to the shape of the rest of the potential and hence they
cannot determine the endpoint.
Nonlinear interactions, represented in this specific case by the $\varphi^4$-term
in the potential~\eqref{eq:phi4P}, eventually quench the instability and
drive the field to a different, stable configuration.

\subsubsection{Tachyonic instability in curved spacetime}
\label{sec:tachcurved}

In the previous section we considered a scalar field that exhibited a tachyonic
instability in flat spacetime.
The generalization to curved spacetime is simple. If we promote the Minkowski
metric to some general curved background described by a metric $g_{\mu\nu}$,
Eq.~\eqref{eq:phi4lin} becomes
\begin{equation}
\label{phi4lincurv}
\Box \, \delta\varphi-\mu^2\delta\varphi=0,
\end{equation}
where $\Box = g^{\mu\nu}\nabla_\mu \nabla_\nu$, with $\nabla_{\mu}$ being the covariant derivative.
The key difference here is that in curved space $\mu^2<0$ is no longer
sufficient to have a tachyonic instability.

To see this, let us take $g_{\mu\nu}$ to be the Schwarzschild metric, which
can describe either a nonrotating \bh or the exterior spacetime of a nonrotating
\ns in \gr.
The line element is,
\begin{align}
 {\rm d} s^2=-\left(1-\frac{2M}{r}\right){\rm d} t^2+ \left(1-\frac{2M}{r}\right)^{-1}{\rm d} r^2+r^2 {\rm d} \Omega^2,
\end{align}
where ${\rm d} \Omega^2={\rm d} \theta^2 + \sin^2\theta {\rm d} \phi^2$
and $M$ is the \bh's mass.
Because the spacetime is static and spherically symmetry, we can decompose the
scalar perturbation $\delta \varphi$ into spherical harmonics $Y_{\ell m} (\theta, \phi)$ and assume a
harmonic time-dependence,
\begin{equation}
\delta \varphi = \sum_{\ell m}\frac{\psi_{\ell m}(r)}{r}
Y_{\ell m}(\theta,\phi) e^{- i \omega t},
\end{equation}
and by substitution into Eq.~\eqref{phi4lincurv} we obtain a Schr\"odinger-like
equation,
\begin{equation}
\frac{{\rm d}^2 \psi_{\ell m}}{{\rm d} r_{\ast}^2}
+ \left[ \omega^2 - V_{\rm eff}(r) \right] \psi_{\ell m} = 0,
\label{eq:eom_pert_schw}
\end{equation}
where we introduced the tortoise coordinate $r_\ast$ defined as ${\rm d} r/{\rm d}
r_{\ast} \equiv 1-2M/r$ and $V_{\rm eff}$ is an effective potential given by
\begin{equation}
V_{\rm eff} = \left( 1 - \frac{2M}{r} \right)
\left[
\frac{\ell(\ell+1)}{r^2} + \frac{2M}{r^3} + \mu^2
\right].
\label{eq:eff_potential}
\end{equation}
which encodes information about the background, curved spacetime.
To have an instability in a Schwarzschild \bh spacetime it is sufficient, but not
necessary that
\begin{equation}
    \int_{-\infty}^{\infty} {\rm d} r_{\ast} \, V_{\rm eff}(r)  \leqslant  0,
\label{eq:bound_state_criteria}
\end{equation}
where $r_\ast=-\infty$ corresponds to the horizon radius in tortoise coordinates. For $M=0$ (flat spacetime) this condition is always satisfied when $\mu^2<0$ and,
unsurprisingly, Eq.~\eqref{eq:eom_pert_schw} yields the same dispersion relation
we discussed in the previous section.
The situation is different for $M\neq 0$. Although not immediately obvious from inspecting Eq.~\eqref{eq:eff_potential}, it turns out that  $\mu^2$ would have to be sufficiently negative for  the tachyonic instability to occur.
The main lesson here is that in curved spacetime, the threshold for the
tachyonic instability to happen clearly depends on the spacetime and we will
return to how one can determine it later.

Note that, although we used Schwarzschild spacetime as an example above, one
can rederive Eq.~\eqref{eq:eom_pert_schw} for a general static, spherically
symmetric background, provided that $r_\ast$ is chosen appropriately.

\subsubsection{Scalarization and gravity}

We have so far considered a scalar field with a negative bare mass squared,
which is obviously not well-motivated.
However, fields can acquire an \emph{effective mass} squared, $\mu_{\eff}^2$,
in specific situations due to their coupling to other fields.
As an example, consider a scalar field that is nonminimally coupled to gravity
and that a term $\varphi^2 R$ is present in the action, where $R$ is the Ricci
scalar.
We then expect a contribution proportional to $\varphi R$ to the scalar's field
equation and, hence, the Ricci scalar contributes to the field's effective
mass, that is, $\mu_\eff^2 \propto \mu^2 + R$.
Assume that the scalar field has no bare mass (i.e.,~$\mu = 0$) and that the coupling
to $R$ is the only contribution to its effective mass $\mu_\eff$.
Then, in flat spacetime, scalar field perturbations would be massless, whereas
in curved spacetimes (with $R\neq 0$) they would be massive and, in general,
also be position dependent.
Moreover, the sign of $\mu^2_{\rm eff}$ would be
controlled by the sign of $R$ in this case. Hence, in some situations, it would
be possible for $\mu^2_{\rm eff}$ to become sufficiently negative in some
spacetime region and cause the scalar field to become tachyonically unstable
despite this being impossible in flat spacetime.

Just like in the flat spacetime example of the scalar with negative $\mu^2$ and
quartic interactions in Sec.~\ref{tachyonquenched}, this tachyonic instability
does not have to be catastrophic. It can simply signal that the scalar needs to
transition to a different configuration once curvature exceeds some threshold.
The instability implies that the scalar field will grow, nonlinearities will
become important and, if they can quench the instability, then one can end up
with a stable, different configuration, for \emph{both the scalar field and the
spacetime}.

This is precisely the idea behind spontaneous scalarization,\footnote{To the best of our knowledge, the
expression ``spontaneous scalarization'' first appeared in print in~\textcite{Damour:1996ke}.} first proposed in
\textcite{Damour:1993hw}: in a given (generalized) scalar-tensor theory, a
configuration with a constant scalar field and a metric that solves Einstein's
equations describes all gravitating systems except some that exhibit strong
gravity. For the latter case, curvature becomes significant enough to render
the constant scalar configurations tachyonically unstable.
The tachyonic instability is eventually quenched by nonlinear effects, and
there exists a stable configuration with a nontrivial scalar and a spacetime
that is no longer a solution to Einstein's equations.

\subsubsection{Strong field phase transitions and weak field screening}

We have not yet shown that the mechanism of spontaneous scalarization, as
described heuristically above, can be at play within a consistent gravity
theory --- we will do so in detail in the next section.
Nonetheless, assuming that the proposal can be successfully implemented in some
model, the following key observations can already be made:
\begin{itemize}
    \item Scalarization is a sharp transition to a new configuration that can
        differ significantly from the \gr configuration for the same object,
        even when one is very near the threshold of the tachyonic instability.
        This is intuitive when one thinks of scalarization as a linear
        instability quenched by nonlinearity: even for a mild instability
        (large timescale) the scalar has to grow significantly for nonlinear
        effects to become important and manage to stop further growth.
    \item The onset of the instability can be controlled by curvature
        couplings. Above we considered as an example a coupling between the
        scalar field and the Ricci scalar $R$, but one envisage couplings with other curvature
        invariants, as we will see in detail in Sec.~\ref{models}. Hence, there
        can be models in which scalarization will only occur in the strong
        field regime (where curvature can become large), while objects that exhibit weak gravity will show no
        deviation from \gr (because the curvature is small).
        %
        Combined with the previous  point, this suggests that spontaneous
        scalarization can be thought off as a strong-field phase transition,
        whereby a field that is dormant in the weak field transitioning to a
        nontrivial configuration in the strong field. Alternatively, one can
        think of scalarization in the reverse way: as a screening mechanism
        that forces a scalar field to transition to a trivial configuration in
        weak field and hence explain why this field has managed to remain
        undetected so far.
    \item The previous argument is based on the rather naive expectation that
        curvature invariants are a good measure of how strong the gravitational
        interaction is. As an example of the failure of this expectation, recall
        that for a Schwarzschild \bh, the Ricci scalar is zero, however other
        curvature scalars, such as the Kretschmann scalar $R_{\mu\nu\rho\sigma} \, R^{\mu\nu\rho\sigma}$
        are nonzero.
        %
        In general, curvature invariants  can have a
        complicated dependence on the properties of compact objects. Hence more
        work is needed to understand what controls the threshold of the
        tachyonic instability and the onset of scalarization. This will be
        addressed in Sec.~\ref{types}.
\end{itemize}

\subsection{Models of scalarization}
\label{models}

\subsubsection{Tachyonic instability and the minimal action}
\label{sec:minimal_action}

We now turn our attention to the gravity theories that can exhibit spontaneous
scalarization. As discussed above in detail, at the perturbative level, the hallmark
of spontaneous scalarization is a tachyonic instability.
This motivates the question: \emph{can we construct a minimal gravity theory which
can have scalar field perturbations that are tachyonically unstable?}
To do so, consider a gravity theory with a metric $g_{\mu\nu}$ and a scalar
field $\varphi$. Assume that the theory is such that
\begin{enumerate}[label=\textsc{A.\arabic*}]
    \item\label{itm:a1_grsol} Spacetimes that are solutions of Einstein's equations,
        potentially with a cosmological constant, and a constant scalar field
        are admissible solutions of this theory as well;
    \item\label{itm:a2_quad} All terms in the action are at least quadratic in $\varphi$.
    \item\label{itm:a3_2ndordereom} The field equations are second-order partial differential equations.
\end{enumerate}
Under these requirements, the equation governing the dynamics of scalar
perturbations $\delta \varphi$ on \gr spacetimes can be cast in the form
\begin{equation}
\label{genLinEq}
g_{\rm eff}^{\mu\nu}\nabla^{(0)}_\mu\nabla^{(0)}_\nu\delta\varphi - \mu_{\rm eff}^2 \, \delta\varphi+{\rm NLC}=0,
\end{equation}
where $g^{\mu\nu}_{\rm eff}$, $\nabla^{(0)}$, and $\mu_{\rm eff}^2$ are all computed in the
background spacetime $g^{(0)}_{\mu\nu}$, and NLC denotes nonlinear corrections.
Here, $g^{\rm eff}_{\mu\nu}$ is an effective metric which can differ from
$g^{(0)}_{\mu\nu}$ for certain types of nonminimal couplings between the metric and
the scalar field, and $\mu_{\rm eff}^2$ may contain not only a bare mass term, but also
other contributions.

Neglecting nonlinearities, and assuming that $g_{\rm eff}$ is nondegenerate and
has a Lorentzian signature, Eq.~\eqref{genLinEq} becomes a curved-spacetime
version of Eq.~\eqref{phi4lincurv}.
This means that one can identify all theories with a single scalar field that are
expected to lead to spontaneous scalarization by considering which couplings
between a scalar and the metric can contribute to $g_{\rm eff}$ and $\mu_{\rm eff}^2$
while still satisfying the assumptions above.
The benefit of taking into account all possible such terms is that it
would allow one to fully explore the onset of scalarization and identify a
class of gravity theories which result in a scalarized
spacetime.

Assumption~\ref{itm:a3_2ndordereom} ensures that there are no unwanted degrees of freedom, as
would be the case if the equations contained higher-order derivatives.
This assumption does limit the possibilities of the terms one can consider.
For example, a coupling term of the type $\varphi^2
\, R^{\mu\nu\lambda\sigma}R_{\mu\nu\lambda\sigma}$ in the action would
contribute to $\mu_{\rm eff}^2$ but leads to higher-order field equations (in
the absence of suitable counterterms).
To deal with this potential pitfall, one can follow the lines
of~\textcite{Andreou:2019ikc} and start from the Horndeski action
\cite{Horndeski:1974wa,Deffayet:2009mn}, also known as generalized scalar-tensor
theory\footnote{See also~\textcite{Motohashi:2018wdq} for a classification
of a broader class of scalar-tensor theories according to their \bh solutions,
including those of \gr.}:
\begin{equation}
\label{eq:Action_Horndeski}
S=\frac{1}{16 \pi G}\sum_{i=2}^{5}\int {\rm d}^{4}x \, \sqrt{-g} \, \mathcal{L}_i + S_{{\rm m}}[\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{equation}
and
\begin{subequations}
\begin{align}
\label{L2}
\mathcal{L}_2 &=  G_2(\varphi,X),\\
\label{L3}
\mathcal{L}_3 &= -G_3(\varphi,X) \, \Box \varphi,\\
\label{L4}
\mathcal{L}_4 &=  G_4(\varphi,X)R+G_{4X} \left[(\Box \varphi)^2-(\varphi_{\mu\nu})^2 \right],\\
\nonumber
\mathcal{L}_5 &= G_5(\varphi,X) G^{\mu\nu} \varphi_{\mu\nu} \nonumber  - \tfrac{1}{6} \, G_{5X} \left[ \left(\Box \varphi \right)^3 - 3 \, \Box \varphi \, (\varphi_{\mu\nu})^2 \right.\\
              &\quad \left. + 2 \, (\varphi_{\mu\nu})^3\right],
\label{L5}
\end{align}
\end{subequations}
where $G_{2}$, $G_{3}$, $G_{4}$, and $G_{5}$ are arbitrary functions of the scalar field $\varphi$ and
its kinetic term $X=-\nabla_\mu \varphi \nabla^\mu\varphi/2$.
Here, $G_{iX}=\partial G_{i}/\partial X$ ($i = 4, 5$), $G^{\mu\nu}$ is the Einstein tensor, and we also introduced the
notation $\varphi_{\mu\nu} = \nabla_{\mu} \nabla_{\nu} \varphi$ so, for example,
$(\varphi_{\mu\nu})^2 = \varphi_{\mu\nu} \varphi^{\mu\nu} = \nabla_\mu\nabla_\nu\varphi \, \nabla^\mu\nabla^\nu\varphi$.
Finally, $S_{{\rm m}}$ is the matter action, with matter fields collectively denoted by $\Psi_{{\rm m}}$.
This is the most general action for a metric and a scalar field that leads to
second-order field equations in four dimensions upon direct variation
(see~\textcite{Kobayashi:2019hrl} for a review).
We will assume for the moment that matter couples minimally to the metric only.
This means that the choice of fields $g_{\mu\nu}$ and $\varphi$
correspond to the so-called Jordan frame, but we will return to this issue later on.

Imposing assumptions~\ref{itm:a1_grsol} and~\ref{itm:a2_quad} on the action~\eqref{eq:Action_Horndeski} places
restrictions to the $G_{i}$ functions as we will see.
We refer to \textcite{Andreou:2019ikc} for a detailed discussion.
For our purposes, it is sufficient to say that by perturbing around an
arbitrary spacetime that is assumed to be a solution of Einstein's equations
with a constant scalar field, we can identify \emph{all} of the terms that
contribute, at the linear level, to $g^{\mu\nu}_{\rm eff}$ and $\mu^2_{\rm eff}$ as
defined in Eq.~\eqref{genLinEq}.
These terms amount to the following action
\begin{align}\label{eq:ActionCaseI}
    S_{\rm min} & = \frac{1}{16 \pi G} \int {\rm d}^{4}x \, \sqrt{-g} \, \left[R-\tfrac{1}{2}(\gamma_{1}+\gamma_{2}R)\nabla_\mu \varphi\nabla^\mu \varphi  \right.
                \nonumber \\
                &\left.\quad + \, \gamma_{2}R_{\mu\nu}\nabla^\mu\varphi\nabla^\nu\varphi - \tfrac{1}{2} \mu_\varphi^2\varphi^2 - \tfrac{1}{4}\beta\varphi^2 R + \tfrac{1}{2}\alpha\,\varphi^2\,\mathscr{G} \right.
                \nonumber \\
                &\left.\quad - \, 2 \Lambda \right] + S_{{\rm m}}[\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{align}
where $\mathscr{G}$ is the Gauss-Bonnet invariant, defined in terms of the Riemann
tensor and its familiar contractions as,
\begin{equation}
\mathscr{G} = R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}
- 4 R^{\mu\nu} R_{\mu\nu}
+ R^2,
\label{eq:def_gauss_bonnet}
\end{equation}
and where $\alpha$, $\beta$, $\gamma_i$, and $\mu^2_\varphi$ can be expressed
in terms of the $G_i$ functions and their derivatives evaluated in the
background configurations.\footnote{We are not following the notation of
\textcite{Andreou:2019ikc}, but have instead adapted it to match that of some
of the specific models we will later study. While expected, it is nontrivial
to show how the Gauss--Bonnet invariant emerges from Eq.~\eqref{eq:Action_Horndeski}.
This was first shown by~\textcite{Kobayashi:2011nu}.}
We will refer to this action, in a slight abuse of terminology, as the {\em
minimal action for scalarization}, in the sense that it contains all of the
terms that contribute to the \emph{onset of scalarization} manifesting as a tachyonic
instability. As such, it can be used to study and understand what triggers
scalarization and to determine the relevant instability thresholds.

Before we go further, it is worth examining what happens if we decide to drop
assumption~\ref{itm:a2_quad} altogether, but still impose assumptions~\ref{itm:a1_grsol} and~\ref{itm:a3_2ndordereom}.
Working along the same lines as before, one arrives at a different set of
terms, composing the action
\begin{align}\label{eq:ActionCaseII}
    S'_{\rm min} & = \frac{1}{16 \pi G} \int {\rm d}^{4}x \, \sqrt{-g} \, \left[R-\frac{1}{2}\frac{\gamma'_{1}+\gamma'_{2}R}{\varphi}\nabla_\mu \varphi\nabla^\mu \varphi \right. \nonumber \\
                 &\left. \quad + \, \frac{\gamma'_{2}}{\varphi}R_{\mu\nu}\nabla^\mu\varphi\nabla^\nu\varphi +\tau\varphi +\eta\,\varphi R +\lambda\,\varphi\,\mathscr{G} \right. \nonumber \\
                 &\left. \quad - \, 2 \Lambda \right] + S_{{\rm m}}[\Psi_{{\rm m}}; g_{\mu\nu}].
\end{align}
As before, $\tau$, $\eta$, $\lambda$, and $\gamma'_i$ can be expressed in terms
of the $G_i$ functions and their derivatives evaluated in the background
configurations. It might seem counter-intuitive that an action containing terms
linear in $\varphi$ leads to perturbation equations that are in the form of
Eq.~\eqref{genLinEq}, which has no source terms. This is due to the presence of
terms that are nonanalytic in $\varphi$.

The first impression may be that abandoning assumption~\ref{itm:a2_quad} has given rise to a second
minimal action for scalarization. However, action \eqref{eq:ActionCaseII}
is just a field redefinition away from action \eqref{eq:ActionCaseI}.
Indeed, one can start from Eq.~\eqref{eq:ActionCaseII}, introduce the redefinition\footnote{Here, and elsewhere in the text, when two actions are related by a field redefinition we will
not relabel the field in order to keep the notation lighter.} $\varphi \to \varphi^2$, and obtain action \eqref{eq:ActionCaseI}, with the following mapping of parameters: $\gamma_1=4 \gamma'_1$, $\mu^2_\varphi=-4\tau$, $\beta=-\eta$, and $\alpha=\lambda$ \cite{Andreou:2019ikc}.
This equivalence demonstrates that (i) assumption~\ref{itm:a2_quad} is redundant, and
(ii) that up to field redefinitons Eq.~\eqref{eq:ActionCaseI} is
sufficient to capture all terms that contribute to the onset of
scalarization and satisfy assumptions~\ref{itm:a1_grsol} and~\ref{itm:a3_2ndordereom}.

One can see by inspection that the $\gamma_i$ terms in
Eq.~\eqref{eq:ActionCaseI} will contribute to $g^{\mu\nu}_{\rm eff}$, while the
rest of the terms will contribute to $\mu^2_{\rm eff}$. Hence, if the latter
vanish, the former cannot trigger scalarization by themselves, as the effective mass would vanish.
Nevertheless, the $\gamma_i$ terms will affect the threshold of the tachyonic
instability we associate with scalarization (cf.~the discussion about the
tachyonic instability in curved spacetime in Sec.~\ref{sec:tachcurved}).
Additionally, $\mu^2_\varphi$ corresponds to the bare mass of the scalar field,
so it is expected to be positive. We then conclude that the terms which are expected to trigger
scalarization in the strong field regime are just the couplings of $\varphi$ to
$R$ and $\mathscr{G} $. In fact, we will see shortly that these are indeed the terms
present in the known models of scalarization.

To summarize, the minimal action \eqref{eq:ActionCaseI} can be used to study
the onset of spontaneous scalarization triggered by a nonminimal coupling to
gravity.
It contains all possible terms that contribute to the associated tachyonic
instability at the linearized level, so it could be used to study the threshold
and onset of this instability in full generality. As discussed earlier, as the
instability progresses, it is expected to be quenched nonlinearly, and the
endpoint will be a scalarized configuration.
The terms in Eq.~\eqref{eq:ActionCaseI} can contribute nonlinearly as well,
but one could add a plethora of other nonlinear interactions, ranging from
scalar self-interactions, e.g.,~$\varphi^4$, to nonminimal coupling terms that
do not contribute to linear perturbations around curved spacetime with constant
scalar, e.g.,~$\varphi^4 \, \mathscr{G} $.
That is, one can start from Eq.~\eqref{eq:ActionCaseI}, or even a subset of
terms therein, and construct different scalarization models. Models that differ
only by terms that are not in Eq.~\eqref{eq:ActionCaseI} will have the same
behavior in what regards the onset of scalarization, and hence the
configurations that one expects to \emph{not} scalarize, but they can differ in
the properties of scalarized solutions \cite{Andreou:2019ikc,Silva:2018qhn,Macedo:2019sem,Minamitsuji:2019iwp}.

Before proceeding on to discuss more specific known models, let us return
to the issue of the coupling to matter. We have so far assumed that matter couples
minimally to the metric only. This assumption is sufficient to ensure that the
theory is compatible with the \WEP~\cite{Will:2018bme}.
To satisfy the \WEP it is sufficient to have matter
couple minimally to some metric, but this does not need to be the same metric
(or choice of other fields) for which the theory has second-order field
equations.
However, it is known that a disformal transformation~\cite{Bekenstein:1992pj}
of the form
\begin{equation}\label{eq:disformal}
g_{\mu\nu} \rightarrow C(\varphi)\left[g_{\mu\nu}+D(\varphi)\nabla_\mu\varphi\nabla_\nu\varphi\right],
\end{equation}
leaves the Horndeski action~\eqref{eq:Action_Horndeski} formally
invariant~\cite{Bettoni:2013diz,Zumalacarregui:2013pma}. It was shown
in~\textcite{Andreou:2019ikc} that coupling matter minimally to a metric that
is related to $g_{\mu\nu}$ by such a disformal transformation as done, e.g., by
\textcite{Minamitsuji:2016hkk}, amounts to a redefinition of $\gamma_2$ in the linearized theory around spacetimes that are solutions of Einstein's equations.
Hence, such a coupling would be redundant when studying
the onset of scalarization using the minimal action~\eqref{eq:ActionCaseI}.
One could also entertain the idea of coupling matter to some composite metric
$\bar{g}_{\mu\nu}$ that depends on both $g_{\mu\nu}$ and $\varphi$ in a
different way compared to~Eq.~\eqref{eq:disformal}.
In such a case, it is likely that assumption~\ref{itm:a3_2ndordereom} would be violated and one
would have to start the analysis presented here with a generalization of the
action~\eqref{eq:ActionCaseI}.

\subsubsection{The Damour--Esposito-Far\`ese Model}
\label{sec:theory_def_model}

As we already mentioned, the concept of scalarization as we described it was
first discussed in \textcite{Damour:1993hw}. They considered the theory,
\begin{align}
S &= \frac{1}{16\pi G_{*}} \int{\rm d}^{4}x \, \sqrt{-g} \, \left[R -
2 \nabla_{\mu}\varphi \nabla^{\mu}\varphi
\right]
\nonumber \\
&\quad + S_{{\rm m}}[\Psi_{{\rm m}};\, {\cal A}^{2}(\varphi)g_{\mu\nu}],
\label{eq:action_st}
\end{align}
The action is said to be written in the Einstein frame, which means that,
contrary to our previous assumptions and conventions, the scalar field is coupled
\emph{minimally} to gravity and has canonical kinetic term.
The coupling with matter field $\Psi_{{\rm m}}$
is through the function ${\cal A}^{2}(\varphi)$.
Here, $G_{\ast}$ carries a subscript as it is not generally equal to $G$ used
so far. Variation of the action~\eqref{eq:action_st} with respect to $\varphi$
yields the field equation,
\begin{equation}
\label{DEFscalareq}
    \Box \varphi= -4\pi G_{*}\alpha(\varphi) T,
\end{equation}
where
\begin{equation}\label{eq:coupling_STT}
    \alpha(\varphi) = {{\rm d} \ln {\cal A}(\varphi)} / {{\rm d} \varphi},
\end{equation}
and $T = g_{\mu\nu} \, T^{\mu\nu}$ is the trace of the matter energy-momentum tensor in the Einstein frame
defined as $T^{\mu\nu} \equiv 2 (- g)^{-1/2} \delta S_{{\rm m}} / \delta g_{\mu\nu}$.
We see that $\alpha(\varphi)$ controls the coupling strength between the scalar field and matter.

If $\alpha(\varphi_0)=0$ for some
constant scalar-field value $\varphi_0$, the constant scalar configuration
with $T\neq 0$ will be an admissible solution of the theory.
It then follows from the generalized Einstein's equations,
\begin{equation}\label{eq:DEF_EinsteinEq}
    R_{\mu\nu} = 2 \nabla_{\mu} \varphi \nabla_{\nu} \varphi + 8 \pi G_{\ast} (T_{\mu\nu} - \tfrac{1}{2} g_{\mu\nu} \, T),
\end{equation}
that these will be solutions of \gr, since the first term in the right hand side vanishes.

At the same time, if we perturb Eq.~\eqref{DEFscalareq} linearly in $\varphi$ in a fixed
background metric which is a solution of \gr, and compare with Eq.~\eqref{genLinEq} we find that $\beta(\varphi_0) \equiv ({\rm d}
\alpha/{\rm d} \varphi)_{\varphi = \varphi_0} $ and $T$ determine the value and
sign of the effective mass square of the perturbations, namely,
\begin{equation}
    \mu_{\rm eff}^2 = -4\pi G_{*} \beta(\varphi_0) T.
    \label{eq:mueffsqr_def}
\end{equation}
For stars, one generally has $T < 0$ and hence, for
a negative sign of $\beta(\varphi_0)$ and the right magnitude of both quantities, the scalar can develop a tachyonic instability around a spacetime that describes stars in \gr, just as we have already discussed, which is studied in detail by~\textcite{Harada:1997mr}.
It was shown by~\textcite{Novak:1998rk} that this instability is quenched by
nonlinearities and the outcome is a \ns with a nontrivial scalar configuration.
These scalarized \ns{s} were shown by~\textcite{Damour:1993hw} to have
properties, such as their masses $M$ and radii $R$, which can be dramatically
different from their \gr counterparts.

In much of the literature considering scalarization, the function $\alpha(\varphi)$ is taken to have the
form\footnote{\textcite{Damour:1993hw} also studied the case in which ${\cal A} = \cos(\sqrt{6} \varphi)$
and hence $\alpha = - \sqrt{6} \, \tan(\sqrt{6} \varphi) \approx - 6 \varphi + 12 \varphi^3 + \dots$, that
includes higher powers in the scalar-field-matter interaction series~\eqref{eq:DEF_alpha} for
$\alpha_0 = 0$ and $\beta_0 = - 6$.}
\begin{equation}\label{eq:DEF_alpha}
  \alpha = \alpha_0 + \beta_0 \varphi \equiv \alpha_{{\mbox{\tiny DEF}}}\,,
\end{equation}
where $\alpha_0$ and $\beta_0$ are dimensionless constants. Sometimes this
choice, rather that the more general action of Eq.~\eqref{eq:action_st}, is
referred to as the \DEF model.
The constant $\alpha_0$ is then assumed to vanish to allow for constant
$\varphi$ solutions [cf.~Eq.~\eqref{DEFscalareq}], or it is assumed to be very small.
In the latter case, all stars will carry some nontrivial scalar field, but by
tuning down $\alpha_0$, any deviation from \gr would be undetectable until
scalarization kicks in.
In its original formulation, the \DEF model did not include a bare mass or
self-interactions for the scalar, but a potential $V(\varphi)$ can be added to
the action and this option has been considered in the literature, e.g., by
\textcite{Popchev2015,Chen:2015zmx,Ramazanoglu:2016kul}, as we will see in
Sec.~\ref{sec:ns_scalarization}.

So far it appears that the DEF model is not covered by our minimal action
\eqref{eq:ActionCaseI} because of our earlier assumption that the scalar does
not couple to the matter.
However, by defining ${\cal A}^2(\varphi) g_{\mu\nu}$ as a new metric and
rewriting the action \eqref{eq:action_st} in terms of that new metric, the
scalar field is no longer coupled to matter. This is referred to as the Jordan
frame. It was shown by~\textcite{Andreou:2019ikc} that, at linearized
level and after a suitable scalar field redefinition, the \DEF model is equivalent to
the action
\begin{align}\label{eq:ActionphisqR}
S &=\frac{1}{16 \pi G} \int {\rm d}^4x \sqrt{-g} \left[\left(1-\tfrac{1}{4}\beta_0 \varphi^2\right) R-\tfrac{1}{2}\nabla_\mu \varphi\nabla^\mu \varphi \right]
 \nonumber \\
&\quad + S_{{\rm m}}[\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{align}
This is indeed a particular case of the action \eqref{eq:ActionCaseI} in which $\gamma_1=1$, $\gamma_2=\alpha=\Lambda=0$, and $\beta=\beta_0$.
Hence, the \DEF model, in what regards the onset of the tachyonic instability
that leads to scalarization, is captured by the minimal action
\eqref{eq:ActionCaseI} and corresponds to one of the two couplings to curvature
that can trigger scalarization.

Before moving on, it is worth mentioning that the original formulation of the
\DEF model, which leads directly to Eq.~\eqref{DEFscalareq}, suggests that it
is the coupling to matter that controls and triggers scalarization. Indeed,
when $T=0$, as is the case for \bh{s}, Eq.~\eqref{DEFscalareq} becomes
$\Box \varphi= 0$, and admits only constant $\varphi$ solutions for stationary and asymptotically
flat configurations by virtue of a no-hair theorem by \textcite{Hawking:1972qk}
(this remains true when one includes a potential,
see~\textcite{Sotiriou:2011dz}).
However, our earlier analysis and the
correspondence between the \DEF model and action \eqref{eq:ActionphisqR} at the
linearized level, makes it clear that the \DEF model is part of a broader class
of theories in which scalarization is present and controlled by the couplings
to curvature, rather than matter, and this observation has been crucial for the
development of models that exhibit \bh scalarization.
It is the fact that $R$ and $T$ are related through the trace of the theory's
generalized Einstein's equation [cf.~Eq.~\eqref{eq:DEF_EinsteinEq}] which
allows for both interpretations in the \DEF model.

\subsubsection{Scalar-Gauss--Bonnet gravity}
\label{sec:theory_sgb}

It was first shown in \textcite{Doneva:2017bvd} and \textcite{Silva:2017uqg}
that theories described by the \sGB action
\begin{align}
\label{actionsgb}
    S &= \frac{1}{16 \pi G} \int {\rm d}^4x \sqrt{-g} \left[ R - \tfrac{1}{2} \nabla_{\mu} \varphi \nabla^{\mu} \varphi + f(\varphi) \, \mathscr{G} \right]
    \nonumber \\
    &\quad + S_{{\rm m}} [\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{align}
can exhibit \bh scalarization provided $({\rm d} f/ {\rm d}\varphi)_{\varphi = \varphi_0}=0$
for some constant $\varphi_0$. This is an existence condition for constant $\varphi$ configurations that are solutions of \gr.
As proven in \textcite{Silva:2017uqg}, \bh solutions of \gr are unique solution
to the theory~\eqref{actionsgb}, provided $({\rm d}^2 f / {\rm d} \varphi^2)_{\varphi = \varphi_0} \, \mathscr{G} < 0$.
To understand this, we can proceed as follows. By varying the action with respect to $\varphi$, we find
\begin{equation}
\label{eq:sgbscalar}
    \Box \varphi+ f_{,\varphi}(\varphi) \, \mathscr{G} = 0\,, \quad f_{,\varphi}(\varphi) = {{\rm d} f} / {{\rm d}\varphi}.
\end{equation}
Once again, we can consider linear perturbation of $\varphi$ on a fixed background and
compare with Eq.~\eqref{genLinEq}. We find that $({\rm d}^2
f/{\rm d}\varphi^2)_{\varphi = \varphi_0} \, \mathscr{G} $ plays the role of an effective mass square for
the scalar perturbations,
\begin{equation}
    \mu^{2}_{\eff} = ({\rm d}^2 f / {\rm d} \varphi^2)_{\varphi = \varphi_0} \, \mathscr{G}.
    \label{eq:mueffsqr_sgb}
\end{equation}
Hence, violating the condition $({\rm d}^2 f / {\rm d} \varphi^2)_{\varphi = \varphi_0} \, \mathscr{G} < 0$  is
necessary, but not sufficient, to develop a tachyonic instability that can lead
to scalarization.

In \textcite{Doneva:2017bvd} $f(\varphi)$ was chosen to be $f(\varphi) =
\lambda^{2}(1 - e^{-3/2 \varphi^2})/12$, whereas \textcite{Silva:2017uqg}
focused on  $f(\varphi)=\alpha\varphi^2/2$. Note that, in the linearized theory
around $\varphi=\varphi_0$, and provided the condition $({\rm d} f/{\rm d}\varphi)_{\varphi = \varphi_0}=0$ is satisfied, any choice of $f(\varphi)$ is equivalent to
$f(\varphi)=\alpha\varphi^2/2$.
Hence, all scalarization models described by action~\eqref{actionsgb} are
captured by the minimal action~\eqref{eq:ActionCaseI}, and in particular by the
coupling between the scalar and the Gauss--Bonnet invariant, on what regards
the onset of scalarization. However, different choices of $f(\varphi)$ will
exhibit different behavior in the nonlinear regime and hence scalarized \bh{s}
will generally have different properties.

Indeed, it was shown by \textcite{Blazquez-Salcedo:2018jnn} that the static,
spherically symmetric scalarized \bh{s} that were found by
\textcite{Silva:2017uqg} for $f(\varphi)=\alpha\phi^2/2$ are unstable against radial perturbations,
unlike their counterparts found by \textcite{Doneva:2017bvd} for $f(\varphi) =
\lambda^{2}(1 - e^{-3/2 \varphi^2})/12$. It was later demonstrated in
\textcite{Silva:2018qhn} by examining the case $f(\varphi)=\alpha\varphi^2/2 + \xi
\varphi^4$ that it is indeed the nonlinearity in $\varphi$ that controls the
stability of scalarized \bh{s}. It was further shown by
\textcite{Macedo:2019sem} and~\textcite{Minamitsuji:2019iwp} that a quartic scalar field
self-interaction would be sufficient to make scalarized \bh{s} stable
for the $f(\varphi)=\alpha\varphi^2/2$ case.
These results, which will be discussed in Sec.~\ref{sec:bh_scalarization}, are
a clear demonstration that although the onset of scalarization can be described
fully using action \eqref{eq:ActionCaseI}, the \emph{endpoint of the tachyonic
instability and the properties of the scalarized configurations will depend to
the nonlinear interaction between the scalar and curvature}, and are thus model
dependent.

We remark that models described by action \eqref{actionsgb} lead to
scalarization of compact stars as well for certain regions of their parameter
spaces, see~\textcite{Silva:2017uqg,Doneva:2017duq}.

\subsubsection{The Ricci-Gauss-Bonnet model}
\label{sec:mixed_model}
As we have seen, the \DEF model and the \sGB models of scalarization correspond
respectively to the $\varphi^2 \, R$ and $\varphi^2 \, \mathscr{G} $ terms (in
addition to $R$ and the canonical kinetic term for the scalar field) in the minimal
action~\eqref{eq:ActionCaseI}. We have also argued above that these two terms
are the only terms that can trigger scalarization as a tachyonic instability
around a spacetime that is a solution of \gr.
These facts together suggest considering the
following action \cite{Antoniou:2020nax,Antoniou:2021zoy}:
\begin{align}\label{action:RGmixed}
    S &= \frac{1}{16 \pi G} \int {\rm d}^{4}x \, \sqrt{-g} \, \left[R-\tfrac{1}{2}\nabla_\mu \varphi\nabla^\mu \varphi-\tfrac{1}{4}\beta\varphi^2 R
    \right. \nonumber \\
      &\quad \left. +\tfrac{1}{2}\alpha\,\varphi^2\,\mathscr{G}\right] + S_{{\rm m}}[\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{align}
This theory is interesting from the perspective of an \EFT.
The terms considered above can be seen as part of an \EFT in which the scalar
field enjoys reflection symmetry (i.e., invariance under $\varphi \to - \varphi$),
and shift symmetry (i.e., invariance under $\varphi \to \varphi + \rm{const.}$)
is broken only by the coupling to the curvature scalars.
This theory is not a complete \EFT, as there are other
operators compatible with these symmetries, such as $\varphi^4 \, R$ or
$G^{\mu\nu}\nabla_\mu \varphi\nabla_\nu \varphi$.
Nonetheless the theory is phenomenologically interesting for various reasons.

Firstly, it has \gr with a constant scalar as a late-time cosmic attractor for
$\beta>0$ \cite{Antoniou:2020nax}. To appreciate why this is important, recall
that one can think of scalization in terms of a tachyonic instability of compact object configurations that are solutions of \gr with
$\varphi=\varphi_0$.
Below the threshold of this instability, these configurations are expected to
be stable and exhibit no deviation from \gr. The attractive feature of
scalarization is that weakly gravitating systems will belong in this category
and hence scalarization can be a form of weak-field screening of the scalar
field. However, this argument assumes that $\varphi=\varphi_0$ everywhere in
the universe and only deviates from this value due to scalarization. If there
were another reason for $\varphi$ to evolve away from $\varphi_0$, this would
make even weakly gravitating objects develop a nontrivial scalar configuration.
Cosmic evolution can indeed cause such evolution as shown by~\textcite{Damour:1992kf}, to the extent that
weakly gravitating systems would become sufficiently scalarized to make the
\DEF model (see e.g.,~\textcite{Anderson:2016aoi}) and \sGB scalarization models
(see e.g.,~\textcite{Anson:2019uto,Franchini:2019npi}) fail weak-field and cosmological tests
of gravity without severely fine-tuning the initial conditions for cosmic
evolution.

The theory described by action \eqref{action:RGmixed} provides an elegant
solution to this problem \cite{Antoniou:2020nax}.
As pointed out in \textcite{Damour:1992kf}, the \DEF model has \gr as a cosmic
attractor for $\beta_0>0$, whereas scalarization requires $\beta_0$ to be
sufficiently negative.
It was already mentioned above that, at the linearized level around a \gr
background, the \DEF model is equivalent to action \eqref{eq:ActionphisqR}, so
it is reasonable to expect that the $\varphi^2 R$ term in action
\eqref{action:RGmixed} would tend to drive the scalar to a constant in late
time cosmology. Moreover, this term should be dominant over the $\varphi^2 \,
\mathscr{G}$ term at low curvatures, while the latter should dominate at high
curvatures and trigger scalarization.
It has indeed been shown in \textcite{Antoniou:2020nax} that cosmic evolution
in the model of action \eqref{action:RGmixed} tracks the cosmic evolution of
\gr from radiation domination onward, and $\varphi$ is driven to $\varphi_0$
very rapidly during matter domination.

It has also been shown in \textcite{Antoniou:2021zoy} and
\textcite{Ventagli:2021ubn} that for $\beta>0$ one can have a range of values
for $\alpha$ in which \bh{s} scalarize, but \ns{s} do not. This is interesting
because the strongest constraints on scalarization so far, and specifically the
\DEF model, are based on binary pulsars~\cite{Kramer:2021jcw,Zhao:2022vig}.
These constraints can be evaded if scalarization is limited to \bh{s}.
\textcite{Antoniou:2021zoy} has also provided strong indications that
scalarized \bh{s} should be radially stable for $\beta>0$, which is not the
case for $\beta=0$. Both of these issue will be discussed in
Sec.~\ref{sec:bh_scalarization}.

In summary, adding a $\beta \varphi^2 R$ term with $\beta>0$ to the simplest \sGB scalarization model addresses a series of concerns.
This term would anyway be there in an \EFT, as it has lower mass dimensions
than $\varphi^2 \, \mathscr{G}$. These considerations should act as a reminder that
scalarization theories are currently still toy models in need of a completion.

\subsubsection{Tensor-multi-scalar theories}

So far we have considered theories with a single field, but one could also study scalarization in models with multiple scalar fields.
These \TMST were studied by~\textcite{Damour:1992we} and
in their simplest form they are described by the action
\begin{align}\label{ActionTMST}
    S &= \frac{1}{16\pi G_{*}}\int {\rm d}^{4}x \, \sqrt{-g} \, [R - 2\gamma_{ab}(\boldsymbol{\varphi})\nabla_{\mu}\varphi^{a}\nabla^{\mu}\varphi^{b}   \nonumber \\
      &\quad - 4V(\boldsymbol{\varphi})] + S_{{\rm m}}[\Psi_{{\rm m}};\, {\cal A}^{2}(\boldsymbol{\varphi})g_{\mu\nu}],
\end{align}
where $\boldsymbol{\varphi}$ denotes a multiplet of $N$ scalar fields
$\{\varphi_1, \dots, \varphi_N\}$, and $\gamma_{ab}(\boldsymbol{\varphi})$ is a
target-space metric that mixes their kinetic terms.
This action can be seen as a generalization of action~\eqref{eq:action_st} to
multiple fields.
We will return to this theory in Sec.~\ref{sec:ns_scalarization}, as it has
been studied in the context of scalarized \ns{s}
by~\textcite{Horbatsch:2015bua} and~\textcite{Doneva:2020afj}.
However, we emphasize that one could consider further generalizations that
involve, e.g., couplings to the Gauss--Bonnet invariant or further derivative
interactions between the scalar fields.

\subsection{Types of scalarization}
\label{types}

One can classify types on scalarization based on which property of the compact
object triggers scalarization and controls the threshold of the tachyonic
instability.
As we have seen from our discussion of the model~\eqref{eq:ActionCaseI}, it is
the couplings between the scalar and curvature invariants that control the
onset of scalarization. So, if one thinks of the onset of scalarization as an
instability around a spacetime of \gr, what controls the onset reduces to how
curvature invariants depend on the properties of the object that curves
spacetime.

\subsubsection{Induced by compactness}

For static, spherically symmetric \bh{s}, there is a straightforward answer: in
\gr the exterior is described by the Schwarzschild spacetime, so $R=0$ and
$\mathscr{G}=48G^2M^2 /r^6$, where $M$ is the mass of the \bh and $r$ is the
areal radial coordinate. Hence, for a given $M$, $\mathscr{G}$ scales
monotonically with $r$ and whether the scalar will develop a negative enough
effective mass square outside the horizon is controlled by $M$, which also
determines the compactness of the \bh.\footnote{By compactness here, we are not referring to the definition as mass over horizon radius, used sometimes in the literature, as this is manifestly constant for Schwarzschild \bh{s}.}

For compact stars the issue is rather more complicated because typically the
tachyonic instability will be stronger in the interior of the star or only be
present there. In earlier papers focusing on the \DEF model, \ns scalarization
was said to be controlled by the compactness of the star.
This is because, in this model the trace of the energy-momentum tensor of the
star is what controls the onset, and for most \EOS, it shows a
steady increase as the radius decreases. Hence, the tachyonic instability is
stronger at the center of the star and how violent it is depends on how compact
the star is. However, exceptions do exist, as we will see in more detail later
on.

These spherically symmetric examples demonstrate that compactness plays a key
role in scalarization, and can be the trigger of a tachyonic instability.

\subsubsection{Induced by spin}
\label{sec:incded_by_spin}

Rapidly rotating \ns{s} in the context of scalarization were first studied in
\textcite{Doneva:2013qva}. They showed that for \DEF-like models, rapid
rotation can enhance scalarization by increasing the parameter space where
scalarization can occur.
Conversely, it is also known that larger spin leads to weaker scalarization for
\bh{s} in sGB gravity with ${\rm d}^2f/{\rm d}\phi^2 <0$ \cite{Cunha:2019dwb}.

In suitable theories, spin can in
itself induce a tachyonic instability that triggers scalarization~\cite{Dima:2020yac},
and spinning scalarized \bh{s} in these theories were
subsequently constructed explicitly by~\textcite{Herdeiro:2020wei} and~\textcite{Berti:2020kgk}.
The fact that spin can trigger \bh scalarization (rather than just control the
scalar charge), is important because it opens up the possibility that only
rapidly rotating \bh{s} carry scalar charge, irrespective of their mass.

These results, as well as further work on scalarized rotating compact objects
that will be discussed in Secs.~\ref{sec:ns_scalarization}
and~\ref{sec:bh_scalarization}, show vividly that spin, a ubiquitous
property of astrophysical objects, can play an
important role in the amount of scalar charge a scalarized object carries.

\subsubsection{Induced by matter or coupling to other fields}
\label{sec:matter_induced_other_field}

So far, we have linked scalarization to a tachyonic instability at the
perturbative level --- although it is fundamentally a non-perturbative effect
--- and we have focused on models in which that instability is controlled by
the nonminimal coupling to gravity.
However, as we saw above, in the \DEF model one can  think of $\mu_{\rm eff}^2$
as being controlled either by the trace of the stress-energy tensor of matter,
$T$, see Eq.~\eqref{eq:mueffsqr_def}, of by the Ricci scalar, $R$, see
action~\eqref{eq:ActionphisqR}.
The latter interpretation has the advantage of providing a unified framework of
scalarization of \bh{s} and stars as linked to a nonminimal coupling to
gravity, which is the perspective we followed. However, the former
interpretation highlights that $\mu_{\eff}^2$ in Eq.~\eqref{genLinEq} could
instead be attributed to any type of coupling between $\varphi$ and another
(matter) field.

For example,~\textcite{Stefanov:2007eq} considered a scalar field coupled to
nonlinear electrodynamics, while~\textcite{Herdeiro:2018wub} focused on
Einstein-Maxwell-scalar theory with the addition of the coupling $e^{-\lambda
\varphi^2}F_{\mu\nu}F^{\mu\nu}$, where $F_{\mu\nu}$ is the usual Faraday tensor.
In both cases it was shown that electrically charged \bh{s} can develop scalar
hair through scalarization. Further work in this direction will be summarized is
Sec.~\ref{sec:bh_scalarization_others}.

When $\mu_{\rm eff}^2$ in Eq.~\eqref{genLinEq} is thought of as introduced by a
coupling to matter, one is also led to consider whether surrounding matter,
such as an accretion disk, a companion, the galaxy, etc., could scalarize a \bh
even in models where \bh{s} cannot scalarize in vacuum. It was shown by~\textcite{Cardoso:2013fwa,Cardoso:2013opa}
that this can indeed occur in the \DEF model.

\subsubsection{Dynamical scalarization}
\label{sec:theory_dyn_sca}

So far we have discussed scalarization of isolated compact objects, but what
happens when they form a binary?
As we alluded to already, when we embed a \ns in an ambient scalar field
environment (in the \DEF model this is equivalent to introducing a small
nonvanishing $\alpha_0$), stars will always carry some small scalar charge
and they can still nonperturbatively develop large charge values when they
scalarize.
This is called \emph{induced scalarization}~\cite{Salgado:1998sg},
and binaries provide a natural scenario for it to occur.
Imagine two \ns{s}, each with its own compactness, such that one is scalarized
and the other is not.
As the system inspirals, at some point the nonscalarized \ns will start
experiencing the presence of the scalar field sourced by its companion,
and then induced scalarization will take place~\cite{Barausse:2012da}.

Another, perhaps more dramatic scenario, is that of \emph{dynamical scalarization}.
In this case, two nonscalarized \ns{s} can become scalarized once their orbital
separation becomes sufficiently small. Qualitatively, this can be quantified
by some measure of an ``effective compactness'' of the binary which, just for an isolated \ns,
can trigger scalarization once it reaches a certain threshold~\cite{Palenzuela:2013hsa,Shibata:2013pra,Taniguchi:2014fqa}.
In a quasicircular binary, this effective compactness only increases (it scales inversely with the
orbital separation), but that has not to be the case in an eccentric orbit.
In such cases, the effective compactness oscillates in time, being largest when the \ns{s} are closest.
This leads to a \emph{transient dynamical scalarization} of the system, whereby the two \ns{s}
continuously scalarize and \emph{descalarize} as the system inspirals.

What about \bh{s}? In \sGB models, the scalar field is sourced by the Gauss--Bonnet invariant,
and therefore a binary composed of two scalarized \bh{s} would in general result in an unscalarized
\bh remnant, since the latter has a larger mass and therefore a smaller spacetime curvature.
That is, the system descalarizes~\cite{Silva:2020omi}. On the other hand, depending on the initial
nonscalarized \bh{s}' spins and masses, one can have cases where the remnant scalarizes due to its large spin, i.e., there can be \emph{dynamical spin-induced scalarization}~\cite{Elley:2022ept}.

As we have seen, compact binaries lead to new manifestations of scalarization.
We will discuss these in more detail in Secs.~\ref{sec:ns_scalarization}
and~\ref{sec:bh_scalarization}.

\subsubsection{Beyond scalarization}
\label{sec:beyond_scalarization}

So far we have discussed scalarizations as a linear, tachyonic instability for a scalar field, that is then quenched nonlinearly and leads to a nontrivial scalar configuration. There are many ways to extend this paradigm and yet keep the key outcome: to have fields that undergo what resembles a phase transition --- a sharp change from a trivial to a nontrivial configuration --- in the strong field regime.

One direction is to generalize the mechanism to different fields, such as vectors, tensors or spinors, generating models of \emph{spontaneous vectorization},
\emph{tensorization} or \emph{spinorization}~\cite{Ramazanoglu:2017xbl,Ramazanoglu:2017yun,Ramazanoglu:2018tig,Ramazanoglu:2018hwk}. Another would be to construct models in which the transition is not triggered by a tachyonic instability, but some other linear instability~\cite{Ramazanoglu:2017yun}. Both of these directions will be discussed in Sec.~\ref{sec:other}.
A third direction comes from the possibility that the transition might not be triggered by a linear instability, but might instead be a fully nonlinear effect~\cite{Doneva:2021tvn}.
It was shown that if one chooses $f(\varphi)=\exp(\beta\varphi^4)$ in the action~\eqref{actionsgb}, one obtains a theory in which scalar perturbations are massless around Kerr \bh{s}, and hence there cannot be any linear tachyonic instability, and yet stable scalarized solutions still exist for certain masses and spins.

\subsection{Quantum aspects and classical analogues}

While most of this review will deal with scalarization from a \emph{classical}
field theory perspective, it is noteworthy to comment on the \emph{quantum
aspects} of this phenomenon.
In particular, \textcite{Lima:2010na} studied a quantum scalar field
nonminimally coupled to gravity living in the classical background of
a compact star spacatime (i.e., they worked in the semiclassical
approximation). They showed that some field modes can go through an exponential
growth causing the vacuum expectation value of the field operator
$\hat{\varphi}^2$ to grow, ultimately causing the vacuum
expectation value of the energy-momentum tensor $\hat{T}_{\mu\nu}$
of the field itself to grow in the same manner (see also~\textcite{Lima:2010xw}).
In this sense, one can say that classical curved spacetimes can ``awake the
vacuum'' state of quantum fields.
But what is the connection with spontaneous scalarization? The model studied
by~\textcite{Lima:2010na} is nothing but Eq.~\eqref{eq:ActionphisqR} for a
massive scalar field; the Jordan-frame equivalent of the \DEF model for a massive scalar.
In the parameter space spanned by the scalar-field-Ricci-scalar coupling
constant and stellar compactness $C \equiv M / R$ relevant for our discussion,
the regions where the instability occurs agree with the classical prediction to where
scalarization should happen as shown by~\textcite{Pani:2010vc}.
Thus, one can think of the classical ``scalar field perturbations'' which
we have so often spoken about so far to be seeded in Nature, more robustly, by
quantum field fluctuations~\cite{Landulfo:2014wra}.
We refer
to~\textcite{Landulfo:2012nz,Lima:2013uya,Mendes:2013ija,Mendes:2014vna,Santiago:2015bve} for
other works in the semiclassical approach.

Finally, while our review will focus on astrophysical implications of
scalarization, we remark that the realization of this effect in condensed
matter systems has also been studied.
More specifically, \textcite{Ribeiro:2019yzv} devised a classical analogue that
exhibit this phenomenon based on nonlinear optics of metamaterials and this
could, in principle, be observed experimentally.

\subsection{Black hole scalarization: vacuum spacetimes}
\label{sec:bhs_vac}

As we saw in Sec.~\ref{sec:theory_bg}, some gravity theories admit the
same vacuum \bh solutions as \gr, and yet can give rise to new branches of
solutions, with scalar hair, once certain conditions are met.
The prototypical example is described by scalar-Gauss-Bonnet theories, whose action, in the absence of matter, is given by Eq.~\eqref{actionsgb}, namely,
\begin{align*}
    S &= \frac{1}{16 \pi G} \int {\rm d}^{4}x \, \sqrt{-g} \, \left[ R - 2 g^{\mu\nu}\partial_{\nu} \varphi \partial_{\mu} \varphi + f(\varphi) \mathscr{G} \right].
\end{align*}
As we discussed in Sec.~\ref{sec:intro}, the first models exhibiting \bh
scalarization in vacuum were proposed in~\textcite{Doneva:2017bvd} and
\textcite{Silva:2017uqg}.
These models remain to date the most studied \bh scalarization models in
the literature and for this reason, they will be our main focus in this
section.

As explained in Sec.~\ref{sec:theory_sgb}, in theories described by the
action~\eqref{actionsgb}, a no-hair theorem guarantees that
the \bh solutions of \gr are unique to the theory, as long
as $({\rm d}^2 f / {\rm d} \varphi^2)_{\varphi = \varphi_0} \mathscr{G} < 0$, for some constant $\varphi_0$~\cite{Silva:2017uqg}.
If this condition is violated, scalar field perturbations can become
tachyonic unstable and the endpoint in expected to be a nonlinear, scalarized \bh~\cite{Ripley:2020vpk}.
Scalarized \bh{s} have been shown to form dynamically, as outcomes from
the core collapse of an initially unscalarized star~\cite{Kuan:2021lol}.

In what follows, we first review our understanding of isolated scalarized \bh{s}
in scalar-Gauss-Bonnet theories, focusing on their properties and their stability.
We then review what is understood about when they are found in binaries.
Finally, we give an overview of vacuum \bh solutions in models that
generalize the action~\eqref{actionsgb}.

\subsubsection{Spherically symmetric and rotating scalarized black holes}
\label{sec:bh_scalarization_sgb}

Spontaneously scalarized \bh solutions were first found in \textcite{Doneva:2017bvd,Silva:2017uqg}.
What distinguishes these works is the choice of the coupling function
$f(\varphi)$ that couples the scalar field to the Gauss-Bonnet invariant
[cf.~Eq.~\eqref{actionsgb}], respectively,
\begin{subequations}
\label{eq:gaussian_and_quadratic_couplings}
\begin{align}
    f &= (\lambda^2/12) \, [1 - \varepsilon \exp(-6\varphi^2) ],
    \label{eq:gaussian_coupling}
    \\
    f &= (\eta /8)  \, \varepsilon \, \varphi^2,
    \label{eq:quadratic_coupling}
\end{align}
\end{subequations}
where $\lambda^2$ and $\eta$ are coupling constants with dimensions of [length]$^2$ and $\varepsilon = \pm 1$,
which should not to be confused with the energy density of Sec.~\ref{sec:ns_scalarization}.
Both coupling functions agree in the small-$\varphi$ limit (i.e., when $\varphi \ll 1$)
and therefore make the same prediction for the onset of scalarization of \gr \bh{s}.
This threshold can be found by searching for bound state solutions, i.e., time-independent
solutions of the linearized field equation for $\varphi$ around a \gr \bh, which are
regular at the event horizon $r_{\rm h}$ and that vanish at spatial infinity.
Hence, the determination of the scalarization threshold reduces to a boundary value problem.
In the simplest case of a Schwarzschild \bh of mass $M$, the dimensionless quantity $\eta / M^2$ plays the role of the eigenvalue and can be determined numerically with standard techniques, such as the shooting method.
The smallest eigenvalue gives the threshold for the formation of the
``fundamental'' scalarized \bh solution (in the sense that the radial
profiles of the scalar field has $n=0$ nodes).
The other eigenvalues give the thresholds for the formation of ``excited states'',
that is, solutions with one or more nodes.
In Fig.~\ref{fig:decoupl_bh_scalarization} we show the results of such a calculation
in the spacetime of a Schwarzschild \bh and $\varepsilon = 1$,
as done in~\textcite{Silva:2017uqg}.
We remark that no bound states can form for $\varepsilon = -1$, because the effective
mass of scalar field perturbations is positive and therefore the effective potential
is positive definite. See the discussion in Sec.~\ref{sec:tachcurved}.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{figures/fig_decoupling.pdf}
\caption{Results of the numerical integration of the equation
of a monopole ($\ell = m = 0$) scalar perturbation $\delta \varphi$
on a Schwarzschild \bh background.
Top panel: the asymptotic value of $|{\rm d} \sigma / {\rm d} r|$ (with $\sigma \equiv r \, \delta \varphi$),
evaluated at a point $r_{\rm max} \gg r_{\rm h}$,
far away from the event horizon $r_{\rm h} = 2M$, as a function of $\eta / M^2$.
Cusps mark when bound state solutions form. They signal the transition
between stable Schwarzschild \bh{s} and scalarized solutions.
Bottom-left panel: effective potential $V_{\rm eff}$  (see Sec.~\ref{sec:tachcurved}) for
$(\eta / M)^2 = 0$ and $5$. In the latter case, $V_{\rm eff}$ develops
a negative region and can support bound states.
Bottom-right panel: radial profiles of $\delta \varphi$ for the first three
bound states, corresponding to $\eta /M^2 = 2.902$, $19.50$ and
$50.93$. These profiles have $0$, $1$ and $2$ nodes, respectively.
The properties of these solutions are shared between the models
of~\textcite{Doneva:2017bvd,Silva:2017uqg}.
Credit:~\textcite{Silva:2017uqg}.
}
\label{fig:decoupl_bh_scalarization}
\end{figure}

While both models~\eqref{eq:gaussian_and_quadratic_couplings} agree
in their prediction of the threshold for scalarization, they differ
significantly in their prediction of the properties of the nonlinear, scalarized solution, as already discussed in Sec.~\ref{sec:scalarization_mechanism}.
In Fig.~\ref{fig:scalarized_bh_D_M} we show the branches of scalarized \bh{s},
in a parameter space spanned by the dimensionless scalar charge and \bh mass.
The solutions were obtained by solving the full system of field equations (see e.g.~\textcite{Julie:2022huo} for detailed discussion).
The results are for the Gaussian coupling~\eqref{eq:gaussian_coupling} with $\varepsilon = 1$.
We show the fundamental, zero-node scalar field solution (dashed, red curves) and
also the excited scalarized \bh{s} (solid, green and blue curves).
We see that the branches of excited \bh{s} can terminate at finite masses $M$.
This is due to a violation of an inequality at the \bh horizon which must be satisfied for the scalar field to be real-valued and is typical of scalar-Gauss-Bonnet theories; see e.g.,~\textcite{Kanti:1995vq} and~\textcite{Antoniou:2017acq}.
The parameters in the coupling function~\eqref{eq:gaussian_coupling} are chosen in such a way that this disappearance of \bh solutions does not affect the fundamental branch.
For other couplings, such as Eq.~\eqref{eq:quadratic_coupling}, even this branch
can violate the regularity condition shortly after bifurcation [see e.g.,~\textcite{Doneva:2018rou,Silva:2018qhn}].
Hence, the domain of existence of scalarized solutions in the Gaussian model~\eqref{eq:gaussian_coupling} is larger
with respect to that of the ``quadratic'' model~\eqref{eq:quadratic_coupling}.

\begin{figure}[t]
\includegraphics[width=0.84\columnwidth]{figures/D_M_BH_merged.pdf}
\caption{The scalar charge $D$ normalized to the coupling parameter $\lambda$ as a function of the normalized \bh mass $M/\lambda$ for sequences of scalarized \bh solutions. The branches with different colors corresponds to solutions with different number of scalar field nodes, where only the nodeless scalar field solutions can be stable \cite{Blazquez-Salcedo:2018jnn}. The \gr solutions have zero scalar field and thus zero $D$, and that is why they lie on the $x$-axis, depicted in magenta.
The scalarized branches are symmetric with respect to the $y$-axis due to the reflection symmetry ($\varphi \to -\varphi$) in the theory.
Figure adopted from~\textcite{Doneva:2017bvd}.
}
\label{fig:scalarized_bh_D_M}
\end{figure}

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{figures/Herdeiro_2021_Qvsj_RBHs.png}
\caption{Rotating scalarized \bh{s} in the Gaussian model~\eqref{eq:gaussian_coupling}.
The shaded regions represent the domain of existence of scalarized rotating BHs in the
plane spanned by the dimensionless scalar charge ($Q_{\rm s} / \lambda$) and spin ($j = a / M$).
These regions are bound by the Kerr family of solutions $Q_{\rm s}/\lambda = 0$ (solid blue line) and, in the nonrotating limit, by the Schwarzschild solution (black dot-dashed line),
and also when the regularity condition for the scalar field is violated (dotted curve).
The inset shows the case for positive coupling constant studied in~\textcite{Cunha:2019dwb}, while
the main panel shows nonlinear spin-induced scalarized \bh{s}.
Credit:~\textcite{Herdeiro:2020wei}.
}
\label{fig:rotating_scalarized_bhs}
\end{figure}

Asymptotically flat, spinning scalarized \bh{s} in the Gaussian models (with $\varepsilon = 1$) and quadratic theories were explored respectively by~\textcite{Cunha:2019dwb} and~\textcite{Collodel:2019kkx}.
In discussing these solutions, it is convenient to consider a plane spanned by the dimensionless scalar charge $Q_s / \lambda$ and dimensionless spin parameter $j \equiv J / M^2$, where $J$ is the angular momentum of the \bh.
Figure~\ref{fig:rotating_scalarized_bhs} shows the domain of existence of scalarized rotating \bh
solutions. In particular, focus on the inset which corresponds to $\varepsilon = 1$ in Eq.~\eqref{eq:gaussian_coupling}.
We see that as the spin increases, so does the existence domain (shaded region) of the solutions.
This result can be understood in terms
of the Gauss-Bonnet invariant not being positive-definite for the Kerr metric~\cite{Cherubini:2002gen}.
More specifically, a Kerr \bh develops increasingly larger regions where the Gauss-Bonnet
invariant is negative as the spin increases and hence suppressing scalarization.
As a consequence, the deviations relative to \gr predictions on physical observables
(e.g., the location of the \isco of massive particles or the \bh shadow) are smaller for rapidly rotating \bh{s}. See Fig.~\ref{fig:scalarized_shadows} for an example.
\textcite{Collodel:2019kkx} studied the quadratic coupling \eqref{eq:quadratic_coupling} and reached similar conclusions. They also obtained radially excited solutions (i.e., scalar field profiles with nodes in the radial direction) and also the existence of ``angularly'' excited rotating \bh{s} (i.e., scalar field profiles with nodes in the angular directions).

\begin{figure}[t]
\includegraphics[width=0.48\columnwidth]{figures/j0_24-rH0_25.jpg}
\includegraphics[width=0.48\columnwidth]{figures/j0_24Kerr.jpg}
\caption{Gravitational lensing and shadow produced by rotating \bh{s}
under similar observation conditions. The left panel corresponds to a scalarized
\bh with mass $M / \lambda \approx 0.237$ and spin $j = 0.24$. The right panel corresponds
to a comparable Kerr \bh.
Credit:~\textcite{Cunha:2019dwb}.
}
\label{fig:scalarized_shadows}
\end{figure}

\subsubsection{The stability of scalarizated black holes
and implications to model building}
\label{sec:bhs_stability_nr}

Having established the existence of scalarized \bh{s}, the natural next task is to study their
stability.
A first indication comes from the study of the thermodynamical properties of such \bh{s}.
\textcite{Doneva:2017bvd} found that the scalarized solutions belonging to the fundamental branch in the Gaussian model are thermodynamically favored relative to the Kerr solution~\eqref{eq:gaussian_coupling},
in the sense that they have a smaller Wald entropy~\cite{Wald:1993nt}.
This suggest the stability of the former, i.e., in a dynamical process \bh{s} in this model
would evolve towards a state of smaller entropy.
In contrast, in the quadratic coupling~\eqref{eq:quadratic_coupling} model always results in scalarized \bh{s} with higher entropy relative to their \gr counterparts.
These observations remain true in the case of rotating \bh{s}~\cite{Herdeiro:2020wei,Collodel:2019kkx}.

\begin{figure}[t]
\includegraphics[width=0.82\columnwidth]{figures/scalbhs_eta_m01.pdf}
\caption{
Scalar-charge-mass diagram for scalarized \bh{s} with a quadratic scalar-Gauss-Bonnet model
with a scalar field potential $V(\varphi) = (\mu^2 / 2) \varphi^2 + (\lambda / 2) \varphi^4$.
The figure uses dimensionless quantities obtained by factors of $\eta^{1/2}$: \bh charge $\hat{Q} = Q / \eta^{1/2}$  and mass $\hat{M} / \eta^{1/2}$, and scalar field parameters $\hat{\mu} = \mu \eta^{1/2}$ and $\hat{\lambda} = \lambda \eta^{1/2}$.
The vertical line marks the threshold for scalarization $\hat{M}_{\rm th}$.
\bh{s} with $\hat{M} \geqslant \hat{M}_{\rm th}$ are radially unstable, while the Schwarzschild solution
is stable. For large $\hat{\lambda}$, we can obtain solutions with $\hat{M} < \hat{M}_{\rm th}$ and form two branches.
The ones which are radially unstable solutions are denoted with a dashed line (upper branch) and those stable to radial
perturbations are denoted with a solid line (lower branch).
The dot marks are the solutions marginally stable, which correspond to a minimum mass, but maximally scalar charged
\bh for fixed $\hat{\mu}$ and $\hat{\lambda}$.
Credit:~\textcite{Macedo:2019sem}.
}
\label{fig:eft_stabilization}
\end{figure}

\begin{figure}[t]
\includegraphics[width=\columnwidth]{figures/AntoniouEtAl2021_Plot1_1_Charge-Mass_beta_p.pdf}
\caption{
Similar to Fig.~\ref{fig:eft_stabilization}, but now  in ``mixed'' model, where a
quadratic coupling between scalar field and Ricci tensor is added to the
quadratic scalar-Gauss-Bonnet interaction.
In this figure, similar to Fig.~\ref{fig:eft_stabilization}, we show the charge
$\hat{Q} = Q / \alpha^{1/2}$ and mass $\hat{M} / \alpha^{1/2}$ made dimensionless
by division by $\alpha^{1/2}$, for different values of scalar-Ricci coupling
$\beta$. See Eq.~\eqref{action:RGmixed} for the theory's action.
The black lines correspond to values of $\beta$ for which all scalarized \bh{s} have masses below the \gr instability
mass threshold, while the red lines mark values of $\beta$ that lead to scalarized \bh masses that are larger than the \gr mass
threshold. Past the turning points (marked with blue dots), the \bh{s} are expected to be unstable to radial perturbations.
Credit:~\textcite{Antoniou:2021zoy}.
}
\label{fig:ricci_stabilization}
\end{figure}

A complementary question is whether scalarized \bh{s} are stable to linear perturbations.
The first step in this direction was carried by~\textcite{Blazquez-Salcedo:2018jnn} which considered radial scalar field and metric perturbations.
They showed that Schwarzschild \bh{s} become unstable at the scalarization threshold and that the stability of the fundamental (i.e., the zero-node) scalarized solution depends on which model in Eq.~\eqref{eq:gaussian_and_quadratic_couplings} is considered: in the quadratic case all solutions are unstable, whereas stable solutions exist in the Gaussian model.

This important difference between the two models can be traced back to the higher-than-quadratic terms in the scalar field in the theory's action.
More concretely,~\textcite{Silva:2018qhn} and~\textcite{Minamitsuji:2018xde} showed the inclusion of a quartic-interaction,
say $(\zeta/16) \, \varphi^4 \, {\mathscr G}$, is sufficient to yield
stable \bh solutions, as long as, $\zeta$ is negative and satisfies $\zeta / \eta < -0.8$.
These two conditions are satisfied by the Gaussian model. Indeed, expanding Eq.~\eqref{eq:gaussian_coupling} to order ${\cal O}(\varphi^4)$, we find $\zeta / \eta = - 3/2$,
helping explain the stability of \bh{s} in this model.

It was later shown by~\textcite{Macedo:2019sem} that the original quadratic model can be made stable under radial perturbations by making the scalar field massive and self-interacting $V(\varphi) \propto \mu^2 \varphi^2 + \lambda \varphi^4$,
as discussed in Fig.~\ref{fig:eft_stabilization}.
From an \EFT perspective, these terms are of lower order compared to the
quartic scalar-Gauss-Bonnet interaction and hence should be included.
A similar analysis was also performed for the Ricci-Gauss-Bonnet model by~\textcite{Antoniou:2021zoy} and is show in Fig.~\ref{fig:ricci_stabilization}.
See also~\textcite{Doneva:2019vuh} for a study of scalarization with massive scalar fields.

The nonlinear stability of \bh{s} in this model have also been studied in the time domain by~\textcite{Ripley:2020vpk}.
They found evidence for regions in mass-coupling parameter space in which the end state
of the radial instability of the Schwarzschild \bh is a stable scalarized \bh.
For larger couplings however, regions where the time-evolution equations change character from
hyperbolic to elliptic appear outside the \bh horizon. This signals a regime in which
the theory does not have a well-posed initial-value problem.~(See~\textcite{Hilditch:2013sba,Ripley:2022cdh} for discussions.)
These results are in agreement with earlier findings by~\textcite{Blazquez-Salcedo:2018jnn}, but in the linear regime.
The radial stability of scalarized \bh{s} in the Ricci-Gauss-Bonnet model introduced in Sec.~\ref{sec:mixed_model}, was analyzed by~\textcite{Antoniou:2022agj}.
They found that although the Ricci scalar does not affect the scalarization threshold, it can, at
the nonlinear level, stabilize the \bh{s} in the quadratic model.

The problem of nonradial perturbations of scalarized \bh solutions,
which is relevant in the context of testing such theories with ringdown \gw
signals, was addressed in a series of papers by Bl\'azquez-Salcedo and
collaborators~\cite{Blazquez-Salcedo:2020rhf,Blazquez-Salcedo:2020caw}.
On the other hand, the stability of rotating scalarized \bh{s} remain as
an open problem, even in the slow rotating limit.

\subsubsection{Spin-induced scalarization}

One property of the Gauss-Bonnet invariant of the Schwarzschild metric is that it
is positive definite everywhere and, as a consequence, spontaneous scalarization
can only occur for positive values of the coupling constant.
Does this always have to be the case? In the case of the Kerr metric, it is known
that regions where $\mathscr{G} < 0$ outside the outer horizon are possible if
the \bh spins with $a / M \geqslant 0.5$~\cite{Cherubini:2002gen}
suggesting the possibility of a
\emph{spin-induced scalarization} when the coupling constant is negative, i.e.,
$\varepsilon = -1$ in Eqs.~\eqref{eq:gaussian_and_quadratic_couplings}.

Although conceptually simple, the problem is not straightforward to analyse because
(i) the analytic form of $\mathscr{G}$ forbids the separation of variables as done
in Sec.~\ref{sec:intro} and thus the determination of the scalarization threshold,
in principle, will require more sophisticated numerical methods compared to the Schwarzschild case and
(ii) the Kerr metric (with $a / M \geqslant 0.5$) has regions where $\mathscr{G}$ can be either negative or positive, where in the latter regions, for negative coupling constants,
the effective mass may trigger a superradiant instability~\cite{Brito:2015oca}.
To overcome these difficulties,~\textcite{Dima:2020yac} reduced the $(2+1)$-dimensional
nonseparable scalar field equation into a $(1+1)$-dimensional system of equations which are
coupled through spherical harmonic multipole indices $l$ and $m$. They then performed time-domain
numerical integration of the scalar field borrowing methods previously developed to study superradiant instability~\cite{Dolan:2012yt}.
\textcite{Dima:2020yac} showed that the tachyonic instability is the dominant one and charted the
parameter space in which spin-induced scalarization would occur. This is shown in Fig.~\ref{fig:domain_spin_induced}.
\textcite{Hod:2020jjy} worked analytically in the infinitely-large coupling limit ($\eta / M^2 \to \infty$) and
confirmed the expectation that $a / M \gtrsim 0.5$ is the minimal necessary spin value for which spin-induced scalarization occurs.
The corrections for large, but finite, coupling were obtained by~\textcite{Hod:2022hfm}.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{figures/DimaEtAlSpinInduced.png}
\caption{Instability timescale $\tau$ for spin-induced scalarization
in the dimensionless spin $a/M$ and dimensionless coupling constant $-\eta / M^2$ plane.
The instability threshold for the total reconstructed field is
shown by the solid green line, while the threshold when the
$m = 0$ modes are excluded is shown by a blue dotted line.
The red dashed line corresponds to the instability threshold for
the $m = 0$ odd modes, while the dot-dashed cyan line marks
the instability threshold for the spherical mode $l = m = 0$.
Credit:~\textcite{Dima:2020yac}.
}
\label{fig:domain_spin_induced}
\end{figure}

Nonlinear spin-induced scalarized \bh solutions were obtained by~\textcite{Herdeiro:2020wei}
and by~\textcite{Berti:2020kgk}.
They confirmed the existence of scalarized \bh{s} solutions in the parameter space region
in which linear theory predict the spin-induced tachyonic instability.
In Fig.~\ref{fig:rotating_scalarized_bhs} we show
the domain of the spin-induced scalarized \bh{s} in the Gaussian model~\eqref{eq:gaussian_coupling}.
We see that scalarized solutions only exist for $j \gtrsim 0.5$, and that their
scalar charge increases with spin.
Moreover, there are solutions that can violate the Kerr bound $j \leq 1$ on \bh spins.

A complementary approach to that of~\textcite{Dima:2020yac}, was used by~\textcite{Doneva:2020nbb,Doneva:2021dqn}. In these works the scalar field equation is evolved
in (2+1)-dimensions, i.e., without doing a multipolar decomposition.
The influence of a non-vanishing scalar field mass explored in~\textcite{Doneva:2020kfv}.
A similar time-domain study, but using the hyperboloidal foliation method~\cite{Zenginoglu:2007jw},
was performed by~\textcite{Zhang:2020pko} and a model with both the Pontryagin and Gauss-Bonnet invariants was studied by~\textcite{Myung:2020etf}.
\textcite{Annulli:2022ivr} studied the effect of strong magnetic fields
on spin-induced scalarization. They found that the magnetic field suppresses this effect, i.e., it pushes the scalarization threshold to larger values of the dimensionless spin
$j$ of the \bh.

The case of spin-induced scalarization has been much less studied in terms on its stability and on how additional nonlinear interactions affect the scalar charge.

\subsubsection{Scalarized black holes in binary systems}

Most work in \bh scalarization has focused on isolated \bh{s},
but to be able to confront these models against \gw observations
one has to turn to the two-body problem.
What phenomenology would we expect in \bh binaries?
To answer this question we will first consider works that explored the
nonlinear regime of the late inspiral and merger of \bh binaries, a regime
that relies on \emph{numerical relativity}.
Next, we will consider works that focus on the inspiral alone, a regime which
can  be modeled with \pn theory.

\begin{figure*}[t]
\includegraphics[width=\textwidth]{figures/log10phi_snapshot_single.pdf}
\caption{Dynamical descalarization in binary \bh head on collisions.
Scalar field and Gauss–Bonnet dynamics on the $xy$–plane of the \bh{s} with initial mass-ratio $q = m_1 / m_2 = 1/2$.
We show the amplitude of $\log_{10}|\Phi|$ (color map) together with the isocurvature levels
of the Gauss–Bonnet invariant at the beginning of the evolution (top left), during the \bh{s}'
approach (top right), shortly before the collision (bottom left) and shortly after
the merger (bottom right).
These levels correspond to
$1 M^{-4}$ (solid line),
$10^{-1} M^{-4}$ (dashed line),
$10^{-2} M^{-4}$ (dot-dashed line) and
$10^{-3} M^{-4}$ (dotted line).
Credit:~\textcite{Silva:2020omi}.
}
\label{fig:descalarization_snapshot}
\end{figure*}

The first study in this direction was done by~\textcite{Silva:2020omi}, which
studied the scalar field dynamics (i.e., in the decoupling limit) in head-on
collisions and nonspinning, quasicircular inspirals of binary \bh spacetimes obtained
from numerical relativity simulations. This work adapted the methods developed in \textcite{Witek:2018dmd}
for shift-symmetric scalar-Gauss-Bonnet theory to theories which allow for spontaneous scalarization.
In particular, because they were interested at phenomenology near the scalarization
threshold,~\textcite{Silva:2020omi} adopted the quadratic
model~\eqref{eq:quadratic_coupling}, with a coupling strength $\eta$ such that
both (or one, depending on the mass ratio) binary components or the remnant \bh
can carry a scalar field bound state through the simulation.
They showed that \bh binaries can either form a scalarized remnant or
dynamically descalarize by shedding off the initial scalar hair (i.e., the scalar
bound state configuration) depending on the values of the coupling constants and
the mass ratio between the two \bh{s}.
An example of dynamical descalarization is shown in Fig.~\ref{fig:descalarization_snapshot}.
Dynamical descalarization was also shown to occur in the nonlinear scalarization
models (see~Sec.~\ref{sec:beyond_scalarization}) by~\textcite{Doneva:2022byd},
that simulated head-on \bh collisions and worked in the decoupling limit.

A natural question that follows is: what would happen when $\varepsilon = -1$, i.e., the case in which
spin-induced scalarization occurs? This is relevant from an observational point of view, because \bh remnant of binary \bh coalescences have typical dimensionless spins of $j \sim 0.7$. This value comfortably meets the criteria for spin-induced scalarization to happen.
This question was explored in~\textcite{Elley:2022ept}, which found through a suite of numerical relativity simulations that:
(i) spin-induced dynamical descalarization can happen, when the remnant \bh spin is smaller
than that of the initially spinning (and scalarized) binary components and
(ii) spin-induced dynamical scalarization in an initially nonspinning \bh binary system
can result in a scalarized rotating \bh.
In the latter case, for values of the coupling constant near the scalarization threshold for the remnant \bh,
the tachyonic instability can be delayed to some time $\approx 100 M$ after the binary merger, and thus, hide the scalar field
throughout the binary's inspiral.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{"figures/spin-induced-bet-143_Phi_xz.pdf"}
\caption{Snapshots of the scalar field, here labeled $\Phi$, and the $\mathscr{G}$ invariant in the $xz$-plane llustrating the phenomenon of spin-induced dynamical scalarization.
The color map indicates the amplitude of the scalar field.
The isocurvature contours are the same as Fig.~\ref{fig:descalarization_snapshot}.
of the $\mathscr{G}$ invariant correspond to
$1 M^{-4}$ (solid line),
$10^{-1} M^{-4}$ (dashed line),
$10^{-2} M^{-4}$ (dot-dashed line), and
$10^{-3} M^{-4}$ (dotted line).
Black (red) lines correspond to positive (negative) values of $\mathscr{G}$.
We show the inspiral (top left), half an orbit before merger (top right),
formation of the first common apparent horizon (bottom left) and about $200M$
after the merger.
Credit:~\textcite{Elley:2022ept}.
}
\label{fig:spin_induced_scalarization}
\end{figure}

A restriction of the decoupling limit is that it does not capture the
backreaction of the scalar field onto the spacetime metric, thus
not allowing one to derive the modifications to the gravitational-waves in the
binary set-ups mentioned previously.
These difficulties can be partially overcome by employing
a formulation of the theory's field equations in a generalized harmonic
gauge, first developed by~\textcite{Kovacs:2020pns,Kovacs:2020ywu}
[see also~\textcite{Ripley:2022cdh}].
This formulation in principle allows one to evolve the complete scalar field and metric dynamics
in scalar-Gauss-Bonnet gravity.
However, the first binary \bh coalescence simulations using this formalism  in scalar-Gauss-Bonnet
theories that allow for \bh scalarization~\cite{East:2020hgw,East:2021bqk} demonstrated that evolution ceases to be well-posed for larger values of the coupling. How different gauge choices and different models perform in this respect is not well explored yet.

\textcite{Franchini:2022ukz}, motivated by going beyond the small-coupling approximation, explored
\emph{fixing-the-equations} approach to perform numerical simulations,
an approach inspired by works in dissipative relativistic hydrodynamics~\cite{Cayuso:2017iqc}.
The main idea is to modify the evolution equations of the theory, such that
the short wavelength modes (with respect to some chosen scale) in the solution are somewhat controlled,
while longer wavelength modes (associated to strong-coupling and the breakdown of well-posedness) become tractable.
\textcite{Franchini:2022ukz} were able to use the formalism to study scalar field collapse in spherical symmetry in the scalarization model of~\textcite{Silva:2018qhn}.
This allowed them to evolve the system past situations were the original equations fail, due to appearance of regions in spacetime
where the evolution changes character from hyperbolic to elliptic, which prevents further evolution of the system in time~\cite{Ripley:2020vpk}.
Moreover, the ``fixed-equations''  system evolves towards static \bh{s} predicted by the original system of equations, thus indicating
that it may provide a suitable ``completion'' of the original theory.
Whether this approach remains valid in a situation with less symmetry, for example the case of \bh binaries,
remains to be explored.

Progress has also been made in modeling compact binary dynamics in the wider separation and lower velocities, i.e., in the domain of validity of \pn theory.
In particular, \textcite{Shiralilou:2020gah,Shiralilou:2021mfl} extended to a general $f(\varphi)$ coupling the calculations of~\textcite{Yagi:2011xp}, that applied for the shift-symmetric scalar-Gauss-Bonnet theory, i.e., $f(\varphi) \sim \varphi$.
These works calculated the leading-order corrections due to curvature nonlinearities in the \gw and scalar waveforms,
finding that corrections due to the Gauss-Bonnet term appear at 1\pn order in \gw{s}.
They also obtained the \gw polarization and phasing.
In addition,~\textcite{Julie:2019sab} argued that during the adiabatic inspiral of two \bh{s}, the Wald entropy
of each \bh is constant. This allows a precise definition of the ``sensitivities'' of \bh{s} in scalar-Gauss-Bonnet gravity in terms of the variation of the \ADM mass as a function of the ambient scalar field the \bh is embedded in, at fixed Wald entropy.
This parallels the notion that the baryonic mass of \ns{s} is constant during
the inspiral~\cite{Damour:1996ke}. They also derived the two-body Lagrangian at 1\pn order.
\textcite{Julie:2022huo} used this result in conjunction with the numerical calculation of sensitivities of spontaneously scalarized \bh{s} to study binaries. They found that, in principle, \bh{s} can evolve towards a situation in which the inner singularity approaches the event horizon of the \bh before they merge.
This suggest the possibility that, dependent on the model of scalar-Gauss-Bonnet gravity and the binary parameters, the field equations might simply not be able to predict a full binary evolution from inspiral to merger.

\subsection{Black hole scalarization in the presence of matter}
\label{sec:bhs_nvac}

\textcite{Stefanov:2007eq} suggested the earliest model of \bh scalarization which
consisted of the \DEF model coupled to nonlinear electromagnetism~\cite{Stefanov:2007bn,Stefanov:2007qw}.
They introduced a coupling between the scalar field and the Born-Infeld Lagrangian, and found
that \bh solutions with scalar hair branch-off from the general
relativistic sequence of solutions. These \bh{s} violate the no-scalar-hair theorems by~\textcite{Mayo:1996mv,Sen:1998ftn} which are valid for
charged, nonrotating \bh{s} in scalar-tensor theory.
The bifurcation point occurs where the \gr sequence of solutions becomes unstable to scalar perturbations~\cite{Doneva:2010ke}.

A variety of scalarization models have been studied in Einstein-Maxwell theory following~\textcite{Herdeiro:2018wub},
which introduced a coupling between the scalar field and the Lorentz invariant ${\cal F}^2 = F_{\mu\nu}F^{\mu\nu}$ of the form $\exp(- \alpha \varphi^2) {\cal F}^2$,
where $\alpha$ is a dimensionless constant.
In this model, \bh{s} can scalarize, with ${\cal F}^2$ playing the role
of the trace of the energy-momentum tensor in the \ns scalarization in \DEF-like models, or the Gauss-Bonnet invariant in \bh scalarization models.

The simpler nature of the model allowed~\textcite{Herdeiro:2018wub} to perform numerical relativity simulations to
study the scalarization process in the time-domain by expanding upon already available code bases~\cite{Sanchis-Gual:2015lje,Sanchis-Gual:2016tcm}.
This allowed~\textcite{Herdeiro:2018wub} to show that the endpoint of the instability
of the Reissner-Nordstr\"om \bh is indeed a \bh
with scalar hair.
Additional studies of this model were done by~\textcite{Fernandes:2019rez,Astefanesei:2019pfq,Herdeiro:2020htm}.
Within this model,~\textcite{Myung:2018vug} studied the instability of the Reissner-Nordstr\"om solution, while~\textcite{Myung:2018iyq,Myung:2019oua,Blazquez-Salcedo:2020jee} analyzed the stability of the scalarized solutions.
\textcite{Fernandes:2019kmh,Herdeiro:2020iyi} explored other variations of this theme, including
scalar-axion couplings, i.e., $\propto h(\varphi) F^{\mu\nu} \tilde{F}_{\mu\nu}$, where
$\tilde{F}_{\mu\nu} = (1/2) \epsilon\indices{_{\mu\nu}^{\alpha\beta}} F_{\alpha\beta}$
In addition, \textcite{Ikeda:2019okp} considered the coupling between the double-dual of the Riemann tensor and $F_{\mu\nu}$.

\textcite{Cardoso:2013fwa,Cardoso:2013opa} explored a relevant scenario from an astrophysical perspective: matter in the vicinity of a \bh could trigger scalarization and, in turn, endow the \bh with scalar hair. They showed that this ``matter-induced'' scalarization (see Sec.~\ref{sec:matter_induced_other_field}) is, in principle, possible for the \DEF model.
The viablility of this proposal in realistic matter configurations, such as an accretion disk or a dark matter halo, has not been explored yet. Matter-induced scalarization seems unlike to happen in light of the severe constraints on the \DEF model. However, this might not be the case for other scalarization models, and the topic deserves further work.

After it was understood that in scalar-Gauss-Bonnet theories the Schwarzschild solution can scalarize, it was natural to ask if the same can happen to its charged counterpart, i.e., the Reissner-Nordstr{\"o}m solution.
\textcite{Doneva:2018rou} analyzed this and found the existence of
two bifurcation points, one at larger masses where the scalarized solutions
bifurcated from the Reissner-Nordstr{\"o}m one, and one at smaller masses where the scalar
charge of the solution decreases again to zero and the branch merges again with the
\gr one.
Scalarized charged \bh{s} (in the fundamental mode) are also thermodynamically favored over the
Reissner-Nordstr{\"o}m solution.
\textcite{Brihaye:2019kvj} showed that for near-extremal Reissner-Nordstr{\"o}m case, scalarization can happen for either sign of the scalar-to-Gauss-Bonnet couplings constant~\eqref{eq:quadratic_coupling}.
Other works varying the choice of the coupling function $f$ between the scalar field (or the axion field) and
the Maxwell invariant, and studying the scalarization of dyonic \bh{s} include~\textcite{Fernandes:2019rez},~\textcite{ Fernandes:2019kmh},~\textcite{Blazquez-Salcedo:2020nhs} and~\textcite{Blazquez-Salcedo:2020jee}.

\textcite{Herdeiro:2019yjy} put forward a motivation to study these models.
They noted that quantum effects can break the scale invariance and vanishing energy-momentum trace properties of electro-vacuum \gr \bh{s}.
As a consequence, even the simplest nonminimally scalar field coupling $\xi \varphi^2 R$ can result in \bh scalarization when these quantum effects are taken into account.
Within these models, they discussed the scalarization of the Reissner-N\"ordstrom solution
and a noncommutative generalization of the Schwarzschild solution.
They found that resultant scalarized \bh{s} are, in general, entropically favored over the \gr solutions.

\subsection{Variation of the curvature-induced scalarization model}
\label{sec:bh_scalarization_others}

In light of the nontrivial effect of rotation on scalarization, it is quite
natural to consider what happens when Gauss-Bonnet invariant is replaced with the
Pontryagin invariant $\pnt$ as the ``curvature source''.
The Pontryagin density is known to vanish in spherically symmetric spacetimes
(such as that of Schwarzschild \bh) and becomes nonzero for stationary
and axisymmetric spacetimes (such as that of a Kerr \bh).
In this sense, all possible scalarized BHs in this theory are necessarily
``spin induced''; a fact naturally shared with \bh solutions in dynamical
Chern-Simons gravity (see e.g.,~\textcite{Yunes:2009hc,Konno:2009kg},
and~\textcite{Alexander:2009tp} for a review) which features a linear coupling
between the scalar field and the Pontryagin invariant.
We remark that the coupling between scalar field and Pontryagin density results in field equations
with higher order derivatives (contrary to the case of Gauss-Bonnet) and, as such, these models are usually treated as \EFT{s},  see e.g.,~\textcite{Motohashi:2011ds,Delsate:2014hba}.
\textcite{Myung:2019wvb} studied the combined effect of Gauss-Bonnet and Pontryagin densities.
In the test field limit, the scalar field dynamics with an effective mass proportional
to $\phi^2 \, \pnt$ was also studied in a Kerr background~\cite{Gao:2018acg,Doneva:2021dcc}
and in the Schwarzschild-NUT background~\cite{Brihaye:2018bgc}.

Another extension of the original ``Gauss-Bonnet'' spontaneous scalarization
model involves considering $n>1$ scalar fields $\varphi^{a}$, whose interaction
is determined by their ``target-space'' $\gamma_{ab}(\varphi^{c})$, an
$n$-dimensional Riemannian manifold~\cite{Damour:1992we,Horbatsch:2015bua}.
These are similar to the tensor-multi-scalar \DEF models for \ns{s} scalarization presented in Sec.~\ref{sec:ns_extended}.
In these models, the scalar field dynamics is determined by the quantity $\gamma_{ab}(\varphi) g^{\rho\sigma}\partial_{\rho}\varphi^{a}\partial_{\sigma}\varphi^{b}$.
\textcite{Doneva:2020qww} numerically obtained scalarized \bh for the case $n=3$ and maximally symmetric target-space geometries.
An important feature of these solution is that the scalar fields do not decay asymptotically as $1 / r$ (i.e., the \bh{s} do not have a monopole scalar charge). This shows that \bh{s} in these theories will not emit dipole scalar radiation when placed in binaries.

\subsection{Scalarization and cosmology}

As already discussed in Sec.~\ref{sec:mixed_model}, one of the key challenges for scalarization is understanding whether it is compatible with cosmology. Recall that the main premise of scalarization is that there exists some constant value of the scalar field, $\varphi_0$, for which spacetimes of \gr become admissible solutions to the field equations and it is these solutions that describe stationary objects that we expect (from an observational perspective) to carry no scalar charge. For this to be true, cosmic evolution needs to comply and drive $\varphi$ to $\varphi_0$ in the late universe, or else the whole universe will in effect be ``scalarized.'' It has been pointed out early on by~\textcite{Damour:1992kf} and more recently by~\textcite{Anderson:2016aoi, Franchini:2019npi} that reaching this preferred value for the scalar fields in the late universe is not generic for simple models of scalarization and instead requires severe fine-tuning of initial conditions.

This could well be an artefact of not having the complete theory and, ideally, one would expect the need to fine-tuned initial conditions to disappear by addressing further corrections to the model. A first step in this direction has been made by~\textcite{Antoniou:2020nax}, inspired by~\textcite{Damour:1992kf}, by showing that the ``mixed'' model of scalarization discussed in Sec.~\ref{sec:mixed_model} has \gr with a constant scalar as a late-time cosmic attractor for the right sign of the coupling between the scalar and the Ricci scalar. This demonstrated that, in principle, scalarization models can be made compatible with late-time cosmology and provided a recipe for doing so: include in the action terms that will dominate in late cosmology (e.g.,~Ricci coupling) over the terms that control the onset of scalarization (e.g.,~Gauss-Bonnet invariant) and hence impose the desired cosmological behavior. However, the specific model discussed as an example in \textcite{Antoniou:2020nax} is by no means unique, and there are many other attempts to address similar concerns~\cite{Chen:2015zmx,Minamitsuji:2022qku,Erices:2022bws}. Perhaps more importantly, understanding the efficiency of the attractor behavior in reproducing the behavior of \gr (with a cosmological constant) quantitatively, down to evolution of perturbations, structure formations, etc. certainly deserves further investigations and could lead to constraints on realistic models of scalarization.

A second point of friction between scalarization and cosmology relates to the early universe. Broadly speaking, scalarization is controlled by the coupling between the scalar field and curvature invariants, so it is inevitable that these couplings will become increasingly important as one runs cosmic evolution backwards and move to higher and higher curvatures. In particular, they will tend to dominate when the size of the universe will become comparable to the size of the compact objects we want scalarization to affect today (see e.g.~\textcite{Antoniou:2020nax}). It has been pointed out by~\textcite{Anson:2019uto,Anson:2019ebp} that quantum fluctuations could then seed a scalarization-like instability during inflation and that it would be hard to prevent this (linear) instability by adding corrections to the action. It should however be stressed that one does not need the scalar to necessarily remain constant or have its current value through the evolution of the universe. Hence one might not need to prevent such an instability, but instead just quench it nonlinearly, exactly as it happens in scalarization itself. The scalar could then smoothly evolve away from the vacuum that leads to this instability before it becomes relevant, as one moves backwards in cosmic time. The key question here is whether there exists an \EFT applicable to the early universe, which contains the scalarization models as late universe limits, but can also host an inflationary scenario compatible with observations.

As a final remark in relation to cosmology, we should emphasize that we have discussed here only the cosmological implications of known scalarization models and only under the assumption that the scalar field is cosmologically subdominant at late times. There is clearly a plethora of generalized scalar-tensor theories that have been used in the context of inflation or dark energy, which exhibit couplings between the scalar and curvature invariants similar to those employed for scalarization models, but in most cases the nature of the couplings differ quite significantly from those of scalarization models. Reviewing such attempts goes well beyond the scope of this manuscript. Exploring whether a scalar field that exhibits scalarization around compact objects could also play a role on inflation or account for dark energy would certainly be interesting, but putting together such a model (with a single scalar) seems particularly challenging. It is worth mentioning in this context that evading constraints on the coupling between a scalar and the Gauss-Bonnet invariant coming from the speed of \gw{s} does require the scalar to be cosmologically subdominant \cite{Franchini:2019npi}.

\subsection{Dynamical evolution}

The study of the dynamical nonlinear regime of scalarization requires
nonperturbative time-evolution schemes, whose availability varies between
models. Finding formulations of gravity theories that are amenable to numerical
time evolution, i.e.,~numerical relativity, is nontrivial even in
\gr~\cite{Pretorius:2005gq,Campanelli:2005dd,Baker:2005vv}, and the issue can
become even more complicated for alternative theories.

Time-domain, nonlinear evolution of DEF-like models have been done since
the early days by~\textcite{Novak:1997hw,Novak:1998rk,Novak:1999jg}, and binary inspirals have been performed
for nearly a decade~\cite{Barausse:2012da}.
This was possible thanks to the relative simplicity of the field equations that allows one
to establish their well-posedness~\cite{Salgado:2005hx,Salgado:2008xh}.
However, detailed and fully nonlinear numerical results in relation to
\gw{s}, a primary tool to test scalarization, are still few, such as the study by~\textcite{Shibata:2013pra}.

The picture is quite different for \bh scalarization where the coupling to the Gauss-Bonnet term can result in a complicated causal structure of the spacetime in dynamical situations. It was demonstrated both at linear~\cite{Blazquez-Salcedo:2018jnn,Blazquez-Salcedo:2020rhf,Blazquez-Salcedo:2020caw} and nonlinear~\cite{East:2021bqk,East:2020hgw,Ripley:2020vpk} levels that, part, but not all, of the parameter space of spontaneous scalarization where scalarized \bh{s} can be found loses hyperbolicity when the theory is evolved in time in a certain set of gauges, which means that no predictions can be obtained.
For the moment, there is no consensus on whether this loss of hyperbolicity can be cured by a better gauge choice or if it is intrinsic to the evolution equations, which calls for further work on the issue. Furthermore, there are no detailed studies of binary evolution and merger even for the part of the parameter space for which hyperbolic time evolution has been shown to exist.
The dynamics of more general models where various coupling terms are present as in Eq.~\eqref{eq:ActionCaseI} are not available either. The lack of results on all these fronts presents important future research directions.

It can be the case that an \EFT, obtained as a certain limit or truncation of a more fundamental theory, is ill-posed while the latter is or is expected to be well-posed. A typical example is relativistic hydrodynamics once viscocity has been taken into account. Indeed, it has been suggested to introduce techniques used in hydrodynamics to gravity in order to render the problematic time evolution of some theories hyperbolic by suitably modifying the equations, while keeping the end point of evolution the same~\cite{Cayuso:2017iqc}.
This approach has seen increasing use in alternatives to \gr, including the specific case of scalar-Gauss-Bonnet theories under certain symmetries~\cite{Franchini:2022ukz}, and it is another avenue to explore the dynamics of scalarization in cases where current methods are inadequate. The applicability and power of such methods to general cases, e.g., fully nonlinear 3+1 evolution of a highly dynamical system, are yet to be confirmed.

\subsection{Model building in and beyond scalarization}

We have covered several models of spontaneous scalarization and their phenomenology for \ns{s} and \bh{s}. What they all have in common is that scalarization is triggered by a tachyonic instability, whose threshold is controlled by less than a handful of terms in the action \cite{Andreou:2019ikc}. Yet, the properties of the final scalarized configuration depend quite delicately on the nonlinearities and coupling terms that contribute beyond the linear level. The specific choices of these nonlinear interactions have been shown to crucially affect the stability  \cite{Blazquez-Salcedo:2018jnn,Silva:2018qhn,Antoniou:2022agj} and the scalar charge \cite{Silva:2018qhn,Macedo:2019sem,Antoniou:2021zoy,Doneva:2018rou} of scalarized solutions, as well as the cosmological behavior of the corresponding theory \cite{Antoniou:2020nax,Anson:2019uto,Erices:2022bws}. Thus it is clear that such choices affect viability and observability of the models. As a result, although certain models of scalarization, such as the original DEF model, have been studied fairly exhaustively, the exploration of the broader class of theories that exhibit scalarization, and their phenomenology, has only started.

Furthermore, one could consider introducing new classes of scalarization based on instabilities that are not tachyonic in nature. As discussed in Sec.~\ref{sec:other}, so far success in this direction is limited but it is nevertheless an existing prospect. Alternatively, one can abandon the idea of a linear instability quenched by nonlinearities altogether as we discussed in Sec.~\ref{sec:beyond_scalarization}. It was recently demonstrated in \textcite{Doneva:2021tvn} that a theory that features a $\varphi^4$ coupling to the Gauss-Bonnet invariant can still have scalarized \bh{s} below a certain mass threshold, although the Kerr solution is linearly stable with respect to (massless) scalar perturbations~\cite{Blazquez-Salcedo:2022omw}. Such observation remains true in more general theories, such as Gauss-Bonnet gravity with multiple scalar fields \cite{Staykov:2022uwq}. This strongly suggests the existence of a class of theories in which scalarization is a purely nonlinear phenomenon.

Scalarization models that have been considered so far typically assume that the scalar field does not couple to matter directly in a suitable choice of variables, usually called the Jordan frame, which suffices for the theory to satisfy the \WEP. However, as pointed out in \textcite{Coates:2016ktu}, this assumption might not be necessary because scalarization itself forces the scalar field into a trivial configuration away from a specific compact object. Therefore, scalarization models can be indistinguishable from \gr for what regards the \WEP for all current observations. It was further argued in \textcite{Coates:2016ktu} and \textcite{Franchini:2017zzx} that particular couplings to matter that do not disrupt scalarization could introduce a Higgs-like mechanism in gravity: the scalar field changes value only near compact objects, and its coupling to matter changes the properties of the Standard Model in these regions. This is a largely unexplored possibility.

As was discussed in detail in Sec.~\ref{sec:other}, another important open question is whether one can generalize the scalarization mechanism to other fields. Doing so in a context of a theory that is free of pathologies has so far proven to be quite challenging, except perhaps for the case of spinors~\cite{Minamitsuji:2020hpl}. This is in part because controlling the dynamics of the extra degree of freedom is notoriously difficult in general in gravity theories with additional vector or tensor fields~\cite{deRham:2014zqa}, and vectorization or tensorization models are no exception.

We should stress that the work on scalarization, or its extensions, has focused on isolated compact objects. In fact, theoretical considerations for putting together scalarization models have also been heavily influence by studies of isolated \bh{s} and \ns{s}. However, as we discussed already above, scalarization could happen dynamically in a binary \cite{Barausse:2012da,Palenzuela:2013hsa,Shibata:2013pra,Taniguchi:2014fqa,Elley:2022ept,Silva:2020omi,Doneva:2022byd}. Apart from finding ways to model this effect~\cite{Khalil:2022sii} and potentially constraining it with \gw observations (see also below), it would be particularly interesting to explore further what properties of a binary control dynamical scalarization, and how this differs between different scalarization models.

\subsection{Observational prospects}

Perhaps the most attractive feature of scalarization is that the tachyonic instability typically leads to large scalar field amplitudes before it is quenched, which leads to nonperturbative deviations from \gr in high curvature regions.
At the same time, the theory mimics \gr closely for weak gravity, easily satisfying existing tests.
This suggests that new fundamental scalar that exhibit this behavior could be discovered with strong gravity observations.

Current observations put stringent bounds on scalar radiation from pulsars, which have already almost completely ruled out the original (massless) DEF model~\cite{Zhao:2022vig}. Electromagnetic observations of compact objects in binaries or surrounded by some form of matter also have the potential to extend such results to more general models of scalarization.
For example, x-ray data from \ns surfaces have been recently considered as a possible way to detect scalarized objects \cite{Silva:2018yxz,Silva:2019leq}.
The masses, radii and moment of inertia of \ns{s} can be inferred through such observations~\cite{Bogdanov:2019ixe,Bogdanov:2019qjb,Bogdanov:2021yip}. Future observation
could be used to distinguish scalarized \ns{s} from nonscalarized ones, for instance, using
\EOS-independent relations. The astrophysical signatures of most extensions of the \DEF model, especially those that allow \bh scalarization have not been studied in details yet, because
they are recent.

The most exciting current development in strong gravity is the advent of \gw astronomy, which is an arena for exploring spontaneous scalarization. Effects of scalarization on \gw{s} of coalescing \ns{s} and/or \bh{s} have been investigated to some extend, mostly for the \DEF model \cite{Sennett:2016rwa,Khalil:2019wyy,Khalil:2022sii,Niu:2021nic} and scalar-Gauss-Bonnet gravity~\cite{Wong:2022wni}.
The most pressing issue in this front is the development of accurate waveform models, that
ideally would cover the inspiral, merger and ringdown stages of coalescing binaries.
Developing such waveform models in \gr has been (and remains) a formidable task that combines tools from post-Newtonian theory, numerical relativity, black-hole perturbation theory, and gravitational self-force. Progress in all these fronts are necessary in order for one to have accurate  ``scalarized'' waveform models.
Stars that undergo core collapse are another interesting source of \gw{s}. Albeit expected to
be rare, such events may produce large bursts of scalar radiation~\cite{Sperhake:2017itk,Rosca-Mead:2019seq,Kuan:2022oxs} and may contribute to the stochastic \gw background which could, in principle, be probed by pulsar timing arrays.
In addition, some model of scalarization (but not exclusively) predict the existence of extra polarizations in \gw{s}, which provides another way to test them.

Last but not the least, future space-borne \gw detectors such as LISA will bring another channel for testing modified gravity. Especially interesting with respect to scalarization is the case of extreme mass-ratio inspirals when the secondary object is scalarized~\cite{Maselli:2020zgv,Maselli:2021men}.
Data from future \gw detectors, such as LISA, or third generation ground-based detectors, such as the Einstein Telescope and the Cosmic Explorer, combined with more precise electromagnetic observations, such as those from improved x-ray observatories, will be crucial for detecting or ruling out scalar fields that exhibit scalarization.

\subsection{Equilibrium neutron stars in the \DEF model}

In this subsection we will review the equilibrium properties of \ns{s} in the
\DEF theory, both in its original form and in extensions of the theory.
The strongest constraints on the latter come from binary pulsar constraints,
which we will also review in detail.

\subsubsection{The original Damour-Esposito-Far\`ese model and binary pulsar constraints}
\label{sec:DEF_model_equil}

\paragraph{Static neutron stars:}~As we saw in Sec.~\ref{sec:theory_def_model},
the effective mass squared of scalar field perturbations in the \DEF model is given by
Eq.~\eqref{eq:mueffsqr_def}, namely,
\begin{equation*}
    \mu_{\rm eff}^2 = -4\pi G_{*} \beta(\varphi_0) T,
\end{equation*}
where we recall that $\beta(\varphi_0)$ is the derivative of the scalar-matter
coupling $\alpha(\varphi)$ evaluated at a constant background scalar field value
$\varphi_0$.
Thus, the condition for scalarization, $\mu_\eff^2 < 0$, can be satisfied when
both the trace of energy-momentum tensor $T$ and $\beta(\varphi_0)$ have the
same sign.
For realistic \ns{s}, which are modeled as a perfect fluid, with pressure $p$ and energy density $\varepsilon$,
we normally have,
\begin{equation}
    T = 3 p - \varepsilon < 0,
\end{equation}
and hence scalarization can happen when $\beta(\varphi_0) < 0$.
In particular, for the coupling function~\eqref{eq:DEF_alpha} studied by~\textcite{Damour:1993hw},
we have that $\beta(\varphi_0) = \beta_0$ is a constant.
Thus most of the studies focus on the situation where both $T$ and $\beta_0$
are negative, but we do remark that at sufficiently high-densities some \EOS{s}
predict a positive sign of $T$, resulting in scalarization when $\beta_0 > 0$;
see e.g.,~\textcite{Mendes:2014ufa,Podkowka:2018gib}.
For most \ns models, $p$ and $\varepsilon$ are related by a barotropic \EOS,
$\varepsilon = \varepsilon(p)$.

In Fig.~\ref{fig:DEF_NSmodels_Nonrot} we show some properties of nonrotating scalarized \ns{s}
for $\alpha_0 = 0$ and some illustrative values of $\beta_0 < 0$. The \EOS is taken to be that
of \APR~\cite{Akmal:1998cf}.
On the left panel we show the \ADM mass of a sequence of \ns solutions,
parametrized by the energy density at the center of the star, $\varepsilon_{{\rm c}}$.
When $\alpha_0 = 0$, \ns solutions of \gr are also solutions of the \DEF model
and we indicate them by the black solid lines. These solutions are
characterized by having zero scalar field.
We see that when a specific critical central energy-density value is reached
(say, $\varepsilon_{{\rm c}, 1}$), the \gr sequence becomes unstable and a new branch of
stable solutions with a nontrivial scalar field (i.e., scalarized stars) bifurcates from it.
In our example, the value of $\varepsilon_{{\rm c}, 1}$ depends only on $\beta_0$
and on the nuclear matter \EOS.
The scalarized branch merges again with
the \gr branch at a second bifurcation point at a larger energy density
$\varepsilon_{{\rm c}, 2}$. Hence, scalarized \ns{s} exist only in the range
$\varepsilon_{\rm c} \in [\varepsilon_{{\rm c}, 1}, \, \varepsilon_{{\rm c}, 2}]$.
We see that the larger $|\beta_0|$ becomes, the more dramatic are deviations
in the \ns mass relative to \gr and also that the range of $\varepsilon_{{\rm c}}$
in which scalarization can happen increases.
This is also shown in the middle panel, where we plot the
\ADM mass as a function of the radius $R$, and in the right panel,
where we show the scalar charge $D$ as a function of the \ADM mass.
The scalar charge is defined in terms of an expansion at spatial infinity
of the scalar field $\varphi$, namely,
\begin{equation} \label{eq:scalar_charge}
    \varphi = \varphi_0 + D/r + {\mathcal O}\left(r^{-2}\right) ,
\end{equation}
where $\varphi_0$ is the cosmological background value of $\varphi$,
often assumed to be zero for simplicity.
The right panel of Fig.~\ref{fig:DEF_NSmodels_Nonrot} shows how $D$ has
a small magnitude near the bifurcation point where scalarization kicks in, it
grows monotonically with $M$, and then it approaches zero again once
$M$ approaches the mass corresponding to the second bifurcation point at $\varepsilon_{{\rm c}, 2}$.

\begin{figure*}[t]
\includegraphics[width=0.315\textwidth]{figures/M_eps_c_APR.pdf}
\includegraphics[width=0.3\textwidth]{figures/M_R_APR.pdf}
\includegraphics[width=0.3\textwidth]{figures/D_M_APR.pdf}
\caption{
Some properties of \ns{s} in the \DEF model with $\alpha_0 = 0$, $\beta_0 < 0$,
and using the \APR \EOS.
We show the mass $M$ as a function of the central energy density $\varepsilon_{\rm c}$ (left panel),
the mass $M$ as a function of the radius $R$ (middle panel), and
the scalar charge $D$ as a function of $M$ (right panel).
In all panels, the solid black lines correspond to the \gr solutions, which are
also solutions in the \DEF model. We use different color lines to represent the
scalarized branches for $\beta_0 = -5.0,\, -4.8$, and $-4.5$.
\ns{s} past the maximum mass in the $M$-$\varepsilon_{\rm c}$ plane are unstable to radial oscillations.
These sequences of stars are shown by the dotted lines.
}
\label{fig:DEF_NSmodels_Nonrot}
\end{figure*}

Recall that the effective potential $V_\eff$ needs to be sufficiently negative in order to support at least one tachyonic mode [cf.~Eqs.~\eqref{eq:eff_potential} and~\eqref{eq:bound_state_criteria}].
For typical \ns densities, Eq.~\eqref{eq:mueffsqr_def} implies that scalarization exists for $\beta_0 \lesssim -4$. \textcite{Damour:1993hw} made this estimate using a Newtonian approximation and also confirmed
it by integrating the fully relativistic equations of stellar equilibrium.
Subsequent works refined the threshold for scalarization to $\beta_0 \lesssim -4.35$ and also showed
that this bound is  not very sensitive to the \EOS~\cite{Novak:1997hw,Silva:2014fca,Motahar:2017blm}.
We note that similar to
the \gr branch of \ns{s}, in the original \DEF model the scalarized solutions are
stable up to the maximum mass of the corresponding branch\footnote{For an
exception in the case of massive scalar field see Sec.~\ref{sec:ns_massive}.}
and the stable solutions are generally energetically favorable over the \gr
\ns{s}. This will be discussed in detail in Sec.~\ref{sec:NS_dynamics}.

Let us briefly comment on the exact definition of mass in scalar-tensor
theories (shown in Fig.~\ref{fig:DEF_NSmodels_Nonrot} and elsewhere).
In contrast to \gr, the definition of mass in
scalar-tensor theories is subtle, due to the fact that
these theories violate the strong equivalence principle.
This results
in the appearance of different possible masses as a measure of the total energy
of the star~\cite{Lee:1974pt,Scheel:1994yr,Scheel:1994yn,Whinnett:1999ws,Yazadjiev:1999hy}.
These works showed that only the so-called ``tensor'' mass
has natural energy-like properties.
For example, the tensor mass is positive definite, it decreases monotonically by the
emission of \gw{s} and it is well defined even for dynamical spacetimes~\cite{Lee:1974pt,Scheel:1994yr,Scheel:1994yn}.
In addition, only the tensor mass leads to a physically acceptable picture
since it peaks at the same point as the particle number, a property crucial for
the stability of the static stars~\cite{Whinnett:1999ws,Yazadjiev:1999hy}.
Therefore, the tensor mass, that is defined as the \ADM mass in the Einstein
frame, should be taken as the physical mass. As a matter of fact though, for most of the commonly used coupling functions the Jordan frame and the Einstein frame \ns masses are identical.

After the discovery of the phenomenon, scalarization was examined for a larger
set of parameters and in more detail by~\textcite{Damour:1996ke}.
They showed that the presence of some externally-imposed scalar field background
$\varphi_0$, as well as considering $\alpha_0\ne 0$, smooths the transition to
a scalarized state. This is what we described as induced scalarization in Sec.~\ref{sec:theory_dyn_sca}.
(Strictly speaking, we do not have pure scalarization in this case since \gr is
not a solution of the field equations anymore).
Spontaneous scalarization was further addressed in
\textcite{Salgado:1998sg}, which considered the problem in the Jordan
frame and performed an approximate Newtonian analysis of the system.
They showed that scalarization can be also associated with the fact that the
effective gravitational constant in scalar-tensor theories decreases for large
scalar fields.
It was further argued by \textcite{Whinnett:2004rr} that scalarization leads to violation of the weak
energy condition in the inner regions of \ns{s}, which can cause instabilities.
It was later argued by~\textcite{Salgado:2004tg}, that this is not a general feature of scalar-tensor theories
and there are subclasses of the theory where the weak energy condition is
easily satisfied.

Another consequence of scalarization is that for sufficiently large $\beta_0$ the maximum allowed mass for \ns{s} increases compared to \gr\footnote{This happens not only in scalar-tensor theories, but also in other modified gravity theories. Upper bounds on the maximum mass of \ns{s} in \gr can be found under minimal assumptions on the \EOS using the approach of~\textcite{Rhoades:1974fn}. See~\textcite{Hartle:1978201} for an early account of applications of this method to other gravity theories.}, which can have various observational consequences. This problem was studied by~\textcite{Sotani:2017pfj} taking into account the effects of various microphysics parameters. Empirical relations were derived for the maximum mass of scalarized \ns{s} that are parametrized with respect to the nuclear saturation parameters and the maximum sound velocity in the core.

Until now we have discussed works that considered negative $\beta_0$.
However, as we mentioned, when the trace of the energy-momentum tensor
is positive, scalarization can also occur for $\beta_0>0$ as well, and this
scenario introduces qualitative differences relative to our story so
far~\cite{Mendes:2014ufa}.
In particular, for the commonly used coupling function $\alpha(\varphi) = \beta_0 \varphi$ of the\DEF theory, scalarized stars are not stable for very high values of $\beta_0$ ($\beta_0 \gg 1$), which will be further discussed in Sec.~\ref{sec:ns_dynamics}.

The calculation of \ns parameters spanning the \DEF theory parameter space
provides a challenging technical task. In particular, tests of this theory
against binary pulsar observations (described in the next section) require the
knowledge of the scalar charges (and its derivatives with respect to the star's
mass) for a large catalog of \EOS{s}.
This calculation has been performed most extensively in the works by~\textcite{Anderson:2019hio,Guo:2021leu},
which provide the results in tabulated form or through surrogate models.
\textcite{Yagi:2021loe} (see also~\textcite{Horbatsch:2010hj}) computed scalar charges in the \DEF model analytically by a combination of perturbative weak-field expansion and Pad\`e resummation,
and found excellent agreement with numerical calculations.

Spontaneous scalarization for the case of several close-by compact objects was considered in~\textcite{Cardoso:2020cwo}. The analytical analysis clearly shows that even though an isolated body might be below the threshold for scalarization, a collection of such bodies can develop a nonzero scalar field while maintaining an average compactness much below the scalarization limit.

We have discussed spontaneous scalarization of \ns{s} so far, but other compact objects can also scalarize.
The case of \bh{s} will be considered in Sec.~\ref{sec:bh_scalarization}, but as an additional nonvacuum example
let us briefly discuss the case of boson stars (cf.~\textcite{Liebling:2012fv} for a review), which were also shown to scalarize~\cite{Whinnett:1999sc}. In such systems, one has a complex scalar field as a matter source for the boson star and an additional real scalar field responsible for scalarization similarly to the \DEF model.
The dynamics of this process was examined by~\textcite{Alcubierre:2010ea} showing that nonlinear development of the scalar field is also observed in the absence of self-interactions in the complex scalar field.
\textcite{Ruiz:2012jt} studied spontaneous and induced scalarization starting with initial data corresponding to stable boson stars in \gr. They showed that a strong emission of scalar radiation occurs during the scalarization process.

\paragraph{Observational constraints from binary pulsars:}
\label{sec:ns_binary_pulsar}
Up to date the binary pulsars set the best constraints on scalarized \ns{s} in the \DEF models~\cite{Kramer:2021jcw,Zhao:2022vig}. The reason for this lies in the remarkable accuracy of pulsar timing observations. One of the important observables is the shrinking of the orbit of close binary systems due to energy loss caused by \gw emission.
In contrast to \gr, \DEF model has an additional scalar degree of freedom that leads to a new channel of energy loss. Thus, the shrinking of the orbit should happen faster. The flux of the scalar dipole radiation, that gives the dominant contribution, is given by~\cite{Damour:1992we},
\begin{equation}\label{eq:NS_scalar_dipole_radiation}
F^{\rm dipole}_{\rm scalar}= A \, (1+ q_1 q_2)^3 (q_2-q_1)^2 ,
\end{equation}
where $A$ is a function depending on the properties of the binary, such as its total and reduced mass, and its eccentricity. $q_2$ and $q_1$ are the scalar charges of each binary component normalized with respect to its mass, $q=D/M$.

Binary pulsar observations were considered in the context of scalarized \ns{s} for the first time by~\textcite{Damour:1996ke}. They calculated the gravitational form factors (also known as ``sensitivities'') of slowly-rotating \ns{s}, that is the  set of coupling constants appearing within scalar-tensor theory in the description of the relativistic motion and timing of pulsars, which led to constraint on the $\beta_0$ parameter. Only a few such binary systems were known at the time, and $\beta_0$ was constrained to be greater than $-5$ using polytropic \EOS{s}. These results were later refined to include realistic \EOS{s} and a limit of $\beta_0>-4.5$ was derived in \textcite{Damour:1998jk}. In addition, an estimate was made that it would be difficult for LIGO and VIRGO to improve $\beta_0$ bounds for most \EOS possibilities (some exceptions are pointed out by \textcite{Sampson:2014qqa} and \textcite{Shao:2017gwu}). The reason is that even though the merger events observed by such \gw detectors can lead to much stronger scalar dipole radiation, they are inferior in accuracy compared to the radio observations of binary pulsars, leading to weaker overall constraints. The next-generation detectors, though, such as Cosmic Explorer and Einstein Telescope, will be able to improve the bounds on scalar-dipole radiation~\cite{Sampson:2014qqa,Shao:2017gwu}.

An example of how constraints can be imposed on \DEF models is presented in Fig.~\ref{fig:BinaryPulsarObservationsMAMBplane}, which illustrates how
the model fails the double pulsar test for specific values of the theory parameters. In this case, the failure is due to the additional energy loss from scalar GWs predicted by the \DEF model, predominantly the dipolar contribution.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{figures/KramerEtAl-fig-mm0737-T1.pdf}
\caption{The mass-mass diagram for the double pulsar PSR~J0737-3039A/B for the
\DEF model with $\alpha_0 = 5 \times 10^{-4}$, $\beta_0 = -4$, and
assuming the \ns{s} are described by \EOS MPA1.
Different measured post-Keplerian parameters are shown by the different curves,
with the width indicating the measurement uncertainty in each parameter.
For the point $\alpha_0 = 5 \times 10^{-4}$ and $\beta_0 = -4$ of the
parameter space of the \DEF model to be consistent with observations, all
curves would have to intersect in a region of the mass-mass plane.
This is not the case, as can be seen with the $\dot{\omega}$ (periastron
advance) and $\dot{P}_{\rm b}$ (change of orbital period) curves.
For further details and for the definitions of the other post-Keplerian parameters
see~\textcite{Kramer:2021jcw}.~Credit:~\textcite{Kramer:2021jcw}.}
\label{fig:BinaryPulsarObservationsMAMBplane}
\end{figure}

With the advances in observational astronomy, more pulsars in binary systems suitable for constraining the scalar dipole radiation have been discovered, a complete and up to date list can be found in~\textcite{FreireWebSite}. Consequently, observational bounds on the scalar-tensor gravity parameter $\beta_0$ for such systems were widely discussed in the literature~\cite{EspositoFarese:2004cc,Antoniadis:2013pzd,Freire:2012mg,Shibata:2013pra,Wex:2014nva,Shao:2017gwu,Voisin:2020lqi,Kramer:2021jcw,Chiba:2021rqa,Zhao:2022vig}.
The strongest current limit comes from \textcite{Zhao:2022vig} that practically closes the scalarizaton window for the original \DEF model, i.e. the possibility for scalarization is ruled out in this specific case. As we will discuss below, though, there are a number of other well-motivated  models where scalarization is still possible, or can not be constrained at all by binary pulsar observations. These include theories with a massive scalar field, Gauss-Bonnet theories of gravity, tensor-multi-scalar theories, or even the standard \DEF model when rapid rotation of \ns{s}, which enhances the effect of scalarization, is considered. Furthermore, the theoretical and numerical approaches developed for the study of the \DEF model are still applicable to these generalized theories in most situations. Thus, we will spend considerable time on the aspects of the \DEF model despite its original form being essentially ruled out.

The constraints on scalarization with $\beta_0>0$ using pulsar-timing observations have been investigated in \textcite{Mendes:2019zpw}. Due to the fact that the scalar charge is suppressed as $\beta_0$ increases, while the range of masses allowing spontaneous scalarization decreases, it turns out that only weak constraints can be imposed by binary pulsar observations in this part of the parameter space.

\paragraph{Rotating scalarized neutron stars:}
So far we have commented only on static \ns models.
All observed \ns{s} are at least slowly rotating, and
some dynamical processes such as \ns mergers or stellar core-collapse can
produce relatively long-lived rapidly rotating (supramassive) protoneutron
stars.
Hence, the inclusion of rotation in \ns physics is an inseparable part of the
goal to explore their astrophysical implications.

\textcite{Damour:1996ke} were the first to study slowly rotating scalarized \ns{s} to leading order in rotation frequency, ${\mathcal O} (\Omega)$, using the formalism of~\textcite{Hartle:1967he,Hartle:1968si},
which allowed them to calculate the \ns moment of inertia (also see \textcite{Sotani:2012eb}).
It is interesting that in this case there is an exact analytical solution for the \ns exterior~\cite{Damour:1996ke}. Static and slowly-rotating \ns{s} for a wide range of realistic \EOS{s}, including examples with hyperons or quark matter, were considered  in \textcite{Motahar:2017blm}.
The extension to second order in the rotational frequency, ${\mathcal O} (\Omega^2)$, was made in \textcite{Pani:2014jra}, which allowed to calculate rotational corrections to star's radius and mass, and also its quadrupole moment.

\begin{figure}
	\centering
	\includegraphics[width=0.45\textwidth]{figures/M_R_APR_RNS}
\caption{The mass-radius relation for \ns{s} in the \DEF-model for $\beta_0 = 0,\, -4.5$~and~$-4.8$
and employing the \APR \EOS. Nonrotating \ns{s} are shown by the solid curves, whereas stars rotating at the mass-shedding limit with the dashed lines. As in Fig.~\ref{fig:DEF_NSmodels_Nonrot}, stars unstable to radial oscillations
are shown with dotted lines. The radius of rotating stars refers to their equatorial radius.}
	\label{fig:DEF_NS_RNS}
\end{figure}

Rapidly uniformly rotating scalarized stars, without
approximation, were obtained by \textcite{Doneva:2013qva}.
They showed that, for a fixed $\beta_0$, the maximum deviation from \gr that is
achieved at the mass-shedding limit is considerably larger compared to the
static case, and the range of central energy densities where scalarization is
possible is significantly broadened.
This can be seen in Fig.~\ref{fig:DEF_NS_RNS} where sequences of static
scalarized \ns{s} are compared to \ns{s} rotating at the Kepler limit, for the
same values of $\beta_0$.

There are a number of factors leading to differences with the static case. The first, more intuitive, one is that the rotational energy of the star also acts as a source for the scalar field, and thus can change the onset and degree of scalarization. Meanwhile, the rapidly rotating models tend to be less compact, which can reduce the degree of scalarization. The large deviations from \gr compared to the static case, on the other hand, are mainly due to the fact that scalarized stars can sustain much larger angular momentum before reaching the Kepler limit. This is a nonlinear effect that could not be normally caught in the slow rotation approximation.

A natural consequence of the rotation effects above is that the minimum $|\beta_0|$ where scalarization is possible changes compared to the static case. Thus, for the same \EOS II~\cite{1985ApJ...291..308D} used by \textcite{Damour:1996ke}, scalarization happens for $\beta_0<-3.9$ in rotating stars \cite{Doneva:2013qva} compared to $\beta_0<-4.35$ in the static case \cite{Damour:1996ke}.
Therefore, one can conclude that although binary pulsar
observations seem to rule-out \DEF scalarization for static or
slowly-rotating \ns{s}~\cite{Zhao:2022vig}, there is still an observationally viable range of $\beta_0$ where only rapidly rotating \ns{s} can scalarize.
This is potentially relevant for binary mergers and stellar core collapse where
such rapidly rotating \ns{s} can form.

A caveat in the previous argument is that, one does not expect the star to rotate uniformly in these extremely dynamical events.
Such \emph{differential} rotation was first studied in~\textcite{Doneva:2018ouu}, which adopted  a simple rotation law that can still capture some of the main properties of the merger remnants, especially a few tens of milliseconds after the merger of the binary~\cite{Bauswein:2017aur}.
When scalarization is considered, larger values of the maximum mass as well as
of the angular momentum can be achieved for supramassive \ns{s} compared to
\gr.
Moreover, the scalar field causes rapidly rotating models to be
less quasitoroidal than their general-relativistic counterparts.
This can have direct astrophysical implications especially for binary \ns mergers where the maximum possible mass and angular momentum of a \ns can sustain are crucial to determine the merger outcome and the lifetime of the merger remnant (in case a supramassive or hypermassive \ns forms\footnote{Supramassive \ns{s} do not have a stable static limit, but are supported against collapse due to rapid rotation. Hypermassive \ns{s} do not have a stable uniformly rotating limit, but are supported against collapse due to differential rotation \cite{Paschalidis:2016vmz}.}).
See e.g.,~\textcite{Bauswein:2013jpa,Takami:2014zpa,Weih:2017mcw,Bauswein:2017aur}.

\subsubsection{Massive scalar field}\label{sec:NS:MassiveSF}
\label{sec:ns_massive}
\begin{figure}
	\includegraphics[width=\columnwidth]{figures/Doneva2016xmf_Fig1.pdf}
\caption{The mass as a function of the radius for \EOS APR4. The results for different values of the coupling constant $\beta_0$ and mass of the scalar field $\mu$ are plotted. Here the potential is assumed to have the form $V(\varphi) = (\mu^2 / 2) \varphi^2$ with a scalar field defined though the action \eqref{actionsgb}. Credit: \textcite{Kruger:2021yay}.}
	\label{fig:DEF_massive_SF}
\end{figure}

A key property of scalar-tensor theories that was neglected in the studies
discussed above is the possibility to have a nonzero scalar field potential.
The simplest case is to take a potential that leads to a nonzero scalar field
mass $\mu$, but more complicated potentials, such as those with
self-interaction terms, can be considered [see Eq.~\eqref{eq:phi4P}].
Although this seems like a simple extension, it has a dramatic effect on the observational properties of \ns{s}, especially on the \gw emission. The reason lies in the different asymptotic behavior of the scalar field.  In the case of zero potential the scalar field decreases as $D/r$ at infinity according to Eq.~\eqref{eq:scalar_charge}. This leads to a nonzero scalar charge $D$ and thus nonzero scalar dipole radiation. In the presence of nonzero scalar field mass $\mu$, though, the scalar field tends exponentially to zero after some characteristic distance related to its Compton wavelength $\lambda_\varphi = 2\pi/\mu$ as discussed by \textcite{Ramazanoglu:2016kul}. Hence, the scalar field is effectively confined to some characteristic radius and its scalar charge is zero. If the orbital separation between the two objects in a binary pulsar system is much larger that $\lambda_\varphi$, the dynamics will not be directly altered by the scalar field and there is no significant emission of scalar dipole radiation~\cite{Alsing:2011er}.
Since $\lambda_\varphi$ is controlled by $\mu$ in the simplest case of scalar
field potential, one can reconcile the \DEF-model (for arbitrarily small
$\beta_0$) with binary pulsar constraints by giving the scalar field a mass
$\mu \gg 10^{-16}$~eV.

Scalarized \ns{s} in massive scalar-tensor theories were first studied in the static case~\cite{Chen:2015zmx,Popchev2015,Ramazanoglu:2016kul}, and later extended to slow~\cite{Yazadjiev:2016pcb} and rapid rotation~\cite{Doneva:2016xmf}.
The inclusion of a quartic self-interaction term to the potential was considered in \textcite{Staykov:2018hhc}.
These works showed that the mass of the scalar field and the self-interaction have similar effects on the scalar field around \ns{s} and they both suppress scalarization.
The quartic interaction by itself cannot affect the range of central energy
densities where scalarized solutions exist, because it is a nonlinear
contribution to the linearized scalar field equation of motion [recall
Eq.~\eqref{genLinEq}]. This self-interaction was found to lead to smaller
deviations from \gr in comparison to the massless case.
In contrast, the mass term shrinks the domain of existence of scalarized \ns{s}
and for large enough masses no scalarization is possible at all.
This is evident in Fig. \ref{fig:DEF_massive_SF}, where \ns mass is plotted as a function of its radius for different combinations of $\beta_0$ and
$\mu$.
The massive scalar field solutions are confined between the zero scalar field
mass models (the original \DEF models) and the \gr ones, corresponding loosely
speaking to $\mu \rightarrow \infty$.
It is interesting that for $\mu \gtrsim 10^{-16}$~eV, the solutions are almost
indistinguishable from the massless \DEF model. For this reason the latter
represents an upper limit on the possible deviations from \gr in massive
scalar-tensor theories.

The exponential asymptotic behavior of the massive scalar field brings
computational challenges to the construction of scalarized \ns{s}, which lead
to new numerical approaches~\cite{Rosca-Mead:2020bzt}, which also facilitated the construction of \ns{s} for very negative $\beta_0$ and large scalar field masses.
\textcite{Rosca-Mead:2020bzt} showed that for sufficiently negative $\beta_0$
qualitative changes in the strongly scalarized branch of solutions are
possible.
For example, the maximum of the scalar field can be located away from the stellar center and in its most extreme form, these solutions are composed of a highly compact \ns model surrounded by a scalar field shell.
Also,~\textcite{Tuna:2022qqr} showed that some scalarized solutions in this
part of the $(\beta_0, \mu)$-parameter space show indications of metastability:
they are stable to small perturbations, but have lower binding energy
compared to their \gr counterparts.

An extension of these results to other forms of coupling functions and
scalar field potentials is the asymmetron model~\cite{Chen:2015zmx,Morisaki:2017nit}, which realizes proper cosmic evolution and can account for the cold dark matter.
These works studied a massive scalar field together with a conformal factor
${\mathcal A}(\infty)>0$ [see Eq.~\eqref{eq:action_st}] that suppresses
deviations from \gr for large scalar fields.
Highly negative values of $\beta_0$ together with a range of scalar
field masses covering several orders of magnitude were explored, including the
case when the scalar field Compton wavelength is much smaller compared to the
size of the \ns.
This can significantly affect the stellar structure, but the interplay of
$\beta_0$, the scalar mass and the form of ${\mathcal A}(\varphi)$ can keep the
deviations from \gr within observational bounds.

\subsubsection{Incorporating further physics}
New aspects of the original \DEF model were recently studied by the inclusion of different physical details.
For instance,~\textcite{Silva:2014fca} studied the presence of anisotropic
pressure of the nuclear matter both for static and slowly-rotating \ns{s}.
The motivation comes from the fact that some theoretical considerations,
e.g., with magnetic fields or within the Skyrme model (a low-energy \EFT of quantum chromodynamics)~\cite{Nelmes:2012uf}
suggest that at high densities the \ns \EOS might have a significant degree of anisotropy~\cite{Herrera:1997plx}.
In such a case, the effects of scalarization increase (decrease) when the
tangential pressure is bigger (smaller) than the radial pressure. The threshold
value of $\beta_0$ for the development of scalarization, which in the
isotropic case is $\beta_0<-4.35$, can be increased due to the
presence of anisotropy, thus widening the range of parameters where
scalarization is possible.

Another astrophysically interesting extension of scalarization is to include
the magnetic fields. According to observations and modeling, \ns magnetic
field values can span from $10^8-10^{12}$~G for standard ``old''  pulsars,
ranging through $10^{16}$~G at the surface of some magnetars, and reaching
hypothetically as high as $10^{17}-10^{18}$~G in the cores of newly formed
protoneutron stars.
Such strong magnetic fields impact the properties of scalarized \ns{s}, including their magnetic deformability, maximum mass, and range of scalarization as studied in~\textcite{Soldateschi:2020hju}.
They also found a magnetically-induced spontaneous scalarization, whose essence is the following:
strong toroidal magnetic fields can support descalarized configurations, and if
the star's magnetic field decreases during some nonideal magnetohydrodynamical process, the star can
undergo a rapid growth of the scalar field, i.e., it scalarizes.
The magnetic quadrupolar deformations of scalarized \ns{s} and the related
\gw{s} produced by rotating magnetars were studied in~\textcite{Soldateschi:2020zxb}

Another interesting extension to the standard \DEF model is related to
challenging the idea that  fundamental physics remains unchanged in the star's
interior, which is a common assumption when constructing a nuclear matter \EOS.
This was studied in \textcite{Coates:2016ktu} where two models were considered
in which the mass of the photon has a different value in the interior and the
vicinity of a compact star compared to the mass measured by experiments
performed in a weak gravity regime. The first model is based on a Proca-like
mass with an effective mass term dependent on $\varphi$.
The second model can be thought of as a gravitational Higgs mechanism where the
Higgs potential is replaced by the scalar-gravity coupling. In both cases the
scalar field undergoes spontaneous scalarization, thus acquiring nonzero value
if the compactness passes a certain threshold, providing a mass to the photon
by coupling to it in an appropriate manner. Although the focus of
\textcite{Coates:2016ktu} is on the electromagnetic field as a proof of
principle, these results can be extended to other fields of the Standard Model.
The signatures of such a gravitational Higgs mechanism on the
behavior of magnetic field of \ns{s} in Einstein-Maxwell theory was studied by~\textcite{Krall:2020kto}.

\subsection{Dynamics of scalarized neutron stars and binary mergers}

The dynamics of isolated \ns{s} can be studied by solving the full nonlinear field equations of scalar-tensor gravity,
which is often a challenging task. Instead, one usually first approaches the problem by linearizing the field equations around a background solution, and then analyze the resulting linearized dynamics. The study of nonlinear dynamics is then done when necessary and feasible.
We follow this sequence in this subsection. We will then review what happens when \ns{s} is scalar-tensor theories are placed in binary systems.

\subsubsection{Linearized dynamics}\label{sec:NS_dynamics}

The studies of the linearized dynamics concern the stability of scalarized stars and
the analysis of their \qnm spectrum. The latter involves the study of different classes
of \ns oscillations modes, and which are tied to emission of \gw{s}.
See~\textcite{Kokkotas:1999bd} for a review.

\paragraph{Stability:}
\label{sec:NS_Stability}

We have seen that scalarized \ns{s} coexist with their nonscalarized counterparts
as solutions the \DEF model. Which of these branches of solutions is the one realizable in Nature?
One way of answering this question consists in calculating the (fractional) binding
energy $M_0 / M - 1$, where $M_0$ is the star's baryonic mass.
Through this calculation,~\textcite{Damour:1993hw} showed
that scalarized \ns{s} are usually energetically favorable over the \gr ones.
We show this in Fig.~\ref{fig:DEF_NSstability_Nonrot}, where we plot the binding energy
as a function of the baryonic mass. At constant baryonic mass, the scalarized solutions (solid curve) have larger binding energies relative to the \gr solutions (dashed curve), and are then energetically favorable.
In addition, we can see a cusp at the maximum of the mass for both for the scalarized and nonscalarized \ns{s}. This suggests that solutions beyond the maximum of the mass (dotted curve in the inset) are unstable.

\begin{figure}
\includegraphics[width=\columnwidth]{figures/M_M0_APR.pdf}
\caption{The fractional binding energy $M_0 / M - 1$ as a function of the baryonic mass $M_0$ for scalarized \ns{s} in the \DEF model with $\beta_0 = - 5$ and using the \APR \EOS. The inset zooms in the scalarized branch and shows the formation of a cusp at the maximum of the mass.
}
\label{fig:DEF_NSstability_Nonrot}
\end{figure}

The stability analysis above relied on the bulk properties of the star. A complementary approach
considers the linear perturbations of the star.
The first step in this direction was taken by~\textcite{Harada:1997mr}, who studied scalar field perturbation in the background of a \ns with zero (or constant) scalar field the \DEF model.
\textcite{Harada:1997mr} studied the perturbation equations in the frequency domain and showed that the \gr solution becomes unstable after a specific critical central energy density.
This is the point where the scalarized solution branch out from the \gr ones; see Fig.~\ref{fig:DEF_NSmodels_Nonrot}.
\textcite{Harada:1998ge} reached similar conclusions, but worked in the context of catastrophe theory.
The radial stability of scalarized \ns{s} was also studied by~\textcite{Mendes:2018qwo}.
They considered metric and scalar field perturbation, and also both signs of $\beta_0$.
They found that scalarized \ns{s} are stable against linear perturbations and
that instability takes place past the point of maximum mass.

\paragraph{Gravitational waves from perturbed \ns:}
\label{sec:NS_QNMs}
\begin{figure*}[th]
\includegraphics[width=0.9\textwidth]{figures/Mendes_2018qwo_Fig3.pdf}
\caption{Frequency (left panel) and inverse damping time (middle panel) of the three lowest-frequency radial modes of stellar models in scalar-tensor theories, for coupling function ${\mathcal A}(\varphi) = \left[ \cosh \left(\sqrt{3} \beta_0 \varphi \right) \right]^{1/(3\beta_0)}$ denoted by M1 and ${\mathcal A}(\varphi) = e^{\beta_0 \varphi^2/2}$ denoted by M2 with $\beta_0 = -5$, as functions of compactness.
Lines (points) represent values computed with frequency (time) domain techniques. In the figure $F$ stands for the fundamental radial oscillation mode (with no nodes) and $H_1$ for its first overtone (with one node), while with $\phi$ is denoted the fundamental scalar mode (with no nodes). The right panel focuses to the details of the fundamental scalar mode, showing its real and imaginary parts. The employed \EOS is a two-phase polytrope consisting of a stiff core and a soft crust.  Credit:~\textcite{Mendes:2018qwo}.}
\label{fig:DEF_Mendes_RadialOscillation}
\end{figure*}

Linearized perturbations are also helpful to studying the \ns oscillation modes directly related to \gw emission.
First, let us note that \gw{s} in scalar-tensor gravity can carry additional polarizations compared to \gr. In particular, one can have breathing modes in addition to the standard ``plus'' and ``cross'' polarizations of \gr~\cite{Will:2018bme}.
Moreover, radial perturbations in scalar-tensor gravity can excite \gw{s}, contrary
to what happens in \gr. These perturbations source monopole scalar waves, that result in a nonvanishing contribution to the perturbed Jordan-frame Riemann tensor (linearized around a Minkowski background) in the transverse-traceless gauge~\cite{Damour:1992we,Novak:1999jg}. For \DEF-like scalar-tensor theories
this requires $\alpha_0 \neq 0$. This contribution is then linked to existence of a breathing polarization mode of the \gw.
We can then conclude that radially oscillating scalarized \ns{s} in a scalar-tensor theories with $\alpha_0 \neq 0$ will emit \gw{s} with amplitude controlled by $\alpha_0$.

The first study of radial oscillations of \ns{s} in the \DEF model was performed in~\textcite{Sotani:2014tua} in the Cowling approximation~\cite{Cowling:1941MNRAS,McDermott:1983ApJ}.
In this approximation, the spacetime is held fixed while only the fluid and the scalar field are perturbed.
Despite its limitations, the Cowling approximation allows one to study the influence of the scalar field on the fluid's radial oscillation frequencies and the total energy radiated by scalar waves.
The full problem, i.e., when both the metric and the scalar field perturbations are taken into account, was addressed by~\textcite{Mendes:2018qwo}. This work found a new family of modes, which were named scalar modes, that have no counterpart in \gr.
The results show that they have quite distinct frequencies and damping times compared to the fluid radial oscillation modes.
In Fig.~\ref{fig:DEF_Mendes_RadialOscillation} we show the frequencies and damping times of the fundamental fluid radial oscillation mode and its first overtone, as well as the fundamental scalar radial mode, as functions of the stellar compactness. We see in the right panel of the figure that scalarized \ns{s} become unstable at the maximum mass and it is the scalar mode which is responsible for the instability. This contrasts with \gr where it is the fundamental radial fluid mode that becomes unstable at the maximum of the mass.

Not only the radial oscillations of scalarized \ns{s} will differ significantly from the \gr case, but the nonradial modes, that are related to the tensorial gravitational wave emission, can be strongly influenced by the scalar field.
The behavior of nonradial oscillations modes, which are sources of the usual tensor polarizations of \gw{s}, can also be strongly affected by the scalar field.
These perturbations can be classified as ``polar'' or ``axial'', depending on how
they behave under parity transformations~\cite{Regge:1957td,Thorne:1967ApJ}.

In particular, scalar field perturbations are of the polar type, meaning that
they only couple to polar perturbations of the fluid and metric. In \gr, these
perturbations are the most efficient \gw{s} sources.
\textcite{Sotani:2004rq} were the first to study the nonradial oscillations of scalarized \ns{s}. They used the Cowling approximation, in which only polar-parity fluid perturbations are dynamical. This simplifies the problem considerably and \textcite{Sotani:2004rq} calculated the (fundamental) $f$- and the (pressure) $p$-modes frequencies.
In a follow-up work,~\textcite{Sotani:2005qx} went beyond the Cowling approximation and
derived equations, for both axial and polar perturbations, that included metric perturbations. They analyzed only the simpler axial perturbations (discussed below).
A study of the polar perturbations only became possible after new techniques were
developed in \gr by~\textcite{Kruger:2019zuz,Kruger:2020ykw}. These techniques
were then used in scalar-tensor theory by~\textcite{Kruger:2021yay}, who also considered
self-interacting massive scalar fields.\footnote{We refer the reader to~\textcite{Staykov:2015cfa,Blazquez-Salcedo:2020ibb} for the case of massive scalar-tensor theory that is mathematically equivalent to $R^2$-gravity.}
\textcite{Kruger:2021yay} found that the scalar field leaves very ``clear'' imprints on the oscillation frequencies of \ns{s}. The extension of these results for rapidly-rotating \ns{s}, but now back to the Cowling approximation, was made by \textcite{Yazadjiev:2017vpg}. They also studied in details the Chandrasekhar-Friedman-Schutz instability driven by rotation~\cite{Chandrasekhar:1992pr,Friedman:1978hf}.

\begin{figure}
\includegraphics[width=\columnwidth]{figures/fmode_l2_MR3}
\caption{The $f$-mode frequency (for $\ell = 2$) as a function of the average \ns density $\sqrt{M/R^3}$ for several values of $\beta$. The SLy \EOS is used~\cite{Douchin:2001sv}. Credit:~\textcite{Kruger:2021yay}.}
\label{fig:DEF_fmodes_l2}
\end{figure}

An interesting property of \ns{s} oscillations in \gr, is the almost linear
and \EOS-independent relation between the quadrupole ($\ell = 2$) $f$-mode frequency
and the average density of the \ns~\cite{Andersson:1997rn}.
\textcite{Sotani:2004rq} and~\textcite{Kruger:2021yay} showed that this linear scaling
is dramatically broken by scalarization.
We show this in Fig. \ref{fig:DEF_fmodes_l2}, where we also see that the deviations
from \gr increase with the increase of $|\beta_0|$.

Another class of modes that can be loosely speaking attributed to the ``oscillations'' of the spacetime itself are the axial spacetime $w$-modes. These modes are somehow easier to calculate (without approximations) because in this case the perturbations of the fluid and the scalar field are zero.
The axial modes of scalarized \ns{s} were considered for the first time in \textcite{Sotani:2005qx}, who the frequencies and the damping times of different classes of $w$-modes were calculated.
The extension of these results to a variety of realistic \EOS{s} including nuclear, hyperonic and hybrid matter, was done by~\textcite{AltahaMotahar:2018djk}, while the case of massive self-interacting scalar fields was studied by~\textcite{AltahaMotahar:2019ekm}.
The effect of scalarization in stronger on the damping times compared to the frequencies, and, in general, their values are reduced compared to \gr.
In addition, \EOS-independent relations between $w$-mode properties known to exist in \gr can also be obtained in scalar-tensor theory.

A class of \ns modes related to the crustal torsional oscillation was studied in \textcite{Silva:2014ora} for the \DEF model in the Cowling approximation. These oscillations most probably follow the giant flares in soft gamma-ray repeaters~\cite{Israel:2005av,Strohmayer:2005ks,Strohmayer:2006py} and are associated to motions in the \ns crust. \textcite{Silva:2014ora} found that for values of $\beta_0$ consistent with binary pulsar constraints, the effect of scalarization on the torsional oscillation frequencies is smaller than the uncertainties in the microphysics modeling of the crust.

\subsubsection{Nonlinear stability and collapse to a black hole}
\label{sec:ns_dynamics}
While the stability and the oscillations of a scalarized \ns can be studied perturbatively, the formation of the scalar hair starting from a \gr \ns solution is a full nonlinear process.
More precisely, the initial exponential growth of the scalar field, starting from unstable \gr solution, can be modelled by linearized dynamics, but the subsequent saturation of the scalar field to an equilibrium value is a nonlinear phenomenon.

The transition from a nonscalarized to a scalarized \ns was first considered by~\textcite{Novak:1998rk} (also see~\textcite{Degollado:2020lsa}), who studied the nonlinear evolution in spherical symmetry.
The numerical simulations of \textcite{Novak:1998rk} showed that the scalar field for an unstable \gr \ns first grows exponentially, and then saturates, with the system saturating to an equilibrium scalarized endstate\footnote{A realistic astrophysical scenario for that of a low mass \ns that can not scalarize on its own. If this star accretes matter, its mass will gradually increase, eventually crossing the point of instability, and then scalarizing}.
A consequence of the work by~\textcite{Novak:1998rk} was the proof of nonlinear stability of scalarized \ns{s}. That is, the numerical simulations showed that the scalarized \ns{s} have stable evolution if the full nonlinear system of field equations (in spherical symmetry) is considered.
Although not strictly dedicated to the stability analysis, fully nonlinear evolution of scalarized \ns{s} were performed in a series of papers that we discuss later.
This provides strong support for the stability of these objects under the most general circumstances.

\textcite{Novak:1997hw} studied the collapse of a scalarized \ns to a \bh in spherical symmetry.
Since \bh no-hair theorems include \DEF models (see~e.g.,~\textcite{Herdeiro:2015waa} for a review), the resulting \bh will be bald and the scalar field has to be radiated away
during collapse.
This takes places in the form of scalar waves, similarly to what we discussed in the case of radial \ns oscillations in Sec.~\ref{sec:NS_dynamics}.

Our discussion so far covered the case $\beta_0<0$ case.
\textcite{Palenzuela:2015ima} and \textcite{Mendes:2016fby} studied the end state of the tachyonic instability in scalar-tensor theories for representative coupling functions with $\beta_0>0$ and realistic \EOS{s}.
This was done both through an energy balance analysis of the existing equilibrium configurations, and by nonlinear dynamical simulations.
They found (contrary to the $\beta_0<0$ case) that the final state of the instability is highly sensitive to the details of the coupling function, varying from gravitational collapse to spontaneous scalarization.
They also found that in the original \DEF model (remember, $\alpha(\varphi)=\beta_0 \varphi$) that scalarized solutions can become unstable compared to the \gr ones when $\beta_0 \gg 1$.
However, stability can be recovered for all values of $\beta_0$, by considering different coupling functions. This is the case for coupling functions with bounded values, such as $\alpha(\varphi \to \infty) = \alpha_\infty$ for some constant $\alpha_\infty$.
This distinction could give rise to novel astrophysical tests for determining the detailed form of the coupling.

\subsubsection{Stellar core-collapse}
Once we know that equilibrium scalarized \ns solutions exist in the \DEF model and they are stable against linear perturbations, the next step is to study how they are formed. Isolated \ns{s} can form after the collapse of the core of a massive star during a supernova explosion. A \ns on the other hand can collapse to a \bh if a threshold for the mass and the angular momentum is reached.
All these are highly dynamical nonlinear processes that can have strong observational signatures both in the electromagnetic and the \gw signals.

The first study of a degenerate stellar core collapse (more specifically with white dwarf initial data) to a scalarized \ns through a bounce and the formation of a shock was done by \textcite{Novak:1999jg}. The simulations were done in spherical symmetry, which allowed them to calculate only the resulting gravitational monopolar radiation. They found that the emitted breathing modes can be potentially detected by LIGO or VIRGO. More interestingly, the emitted signal will be substantially different than the collapse of a \ns to a \bh, allowing one to distinguish between them.

\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{figures/Sperhake_2017itk_fig1.pdf}
\caption{The scalar wave $r \varphi$ during stellar core-collapse extracted at $5\times 10^4\,{\rm km}$. The legend lists deviations from the fiducial parameters  $\mu=10^{-14}\,{\rm eV},~\alpha_0=10^{-2},~\beta_0=-20,~ \Gamma_1=1.3,~\Gamma_2=2.5,~\Gamma_{\rm th}=1.35$, where $\mu$ is the scalar field mass, and $\Gamma_1$, $\Gamma_2$ and $\Gamma_{\rm th}$ are parameters of the polytropic \EOS and the thermal contribution. See \textcite{Sperhake:2017itk} for more details. Credit:~\textcite{Sperhake:2017itk}.}
\label{fig:NS_core_collapse}
\end{figure}

\citet{Gerosa:2016fri} studied the problem of spherically symmetric core collapse in the \DEF model in further details. Two types of initial data were used: collapse of a stellar iron core and collapse of ``realistic'' \ns progenitors that were obtained from computations of stellar evolution. Depending on the theory parameters three possible outcomes of the core collapse are possible -- collapse to a \gr \ns, collapse to a scalarized \ns{s}, collapse to a short-lived protoneutron star followed by a (nonscalarized) \bh formation. It was in the latter case that the most prominent \gw signal with a clear signature from the presence of nontrivial scalar fields was observed. While the fluid dynamics during the collapse is only weakly affected by the scalar field, the contrary is not true -- the scalar radiation depends strongly on the specifics of the matter collapse as well as on the choice of the coupling parameters $\alpha_0$ and $\beta_0$.

The inclusion of a scalar field mass to the core collapse simulations has  some  interesting consequences~\cite{Sperhake:2017itk}. The waveform $r\varphi$ during such events is shown in Fig. \ref{fig:NS_core_collapse}. The mass itself does not have a large influence on the dynamics, but it allows to reconcile the theory with binary pulsar observations for a much larger range of $\alpha_0$ and $\beta_0$, as discussed in Sec.~\ref{sec:NS:MassiveSF}. Thus, a dramatic increases in the radiated \gw signal observable even with the existing \LVK detectors was predicted. A prominent feature is that we expect to receive an inverse chirp signal from such events that will last for years, with a near monochromatic signature on time scales of $\sim 1$ month. The extension of these results to the case of self-interacting scalar field potential was made by~\textcite{Cheong:2018gzn} and~\textcite{Rosca-Mead:2019seq}. Constraints on the theory based on the scalar field evolution in the Einstein frame were imposed by~\textcite{Geng:2020slq}. The problem was considered in greater detail in~\textcite{Rosca-Mead:2020ehn} where the three possible scenarios of the collapse outcome in scalar-tensor theory with sufficiently negative $\beta_0$ were found. These constitute the formation of a \bh following multiple \ns stages, the multi-stage formation of a strongly scalarized \ns, and the single-stage formation of a strongly scalarized NS.~\textcite{Rosca-Mead:2019seq} found that the resulting \gw signal can reach a signal-to-noise ratio of over $20$ for the existing \gw detectors, which has the potential to put strong constraints on the theory.

\subsubsection{Dynamical scalarization and neutron star mergers}

Another highly dynamical and nonlinear process that has important astrophysical, and especially \gw implications is the binary \ns mergers.
This problem can be more challenging to solve compared to the stellar core-collapse because it is not possible by construction to apply certain approximations, such as spherical symmetry.
That is why the binary merger dynamics in the \DEF model was addressed only years later, first by~\textcite{Barausse:2012da} and~\textcite{Shibata:2013pra}. The overall conclusions of both works are similar and they reside in the fact that even if one or both of the \ns{s} are not scalarized before the merger, they can develop nonzero scalar fields during the inspiral. This was called \emph{dynamical scalarization}.
The main advances in~\textcite{Shibata:2013pra} constitute in using several realistic \EOS with consistently derived bounds on the parameter $\beta_0$ as well as developing an initial data solver for scalarized binary \ns{s}. The overall results demonstrated the significant change of the inspiral \gw signal at the moment of dynamical scalarization and afterwards. The reason is that the inspiral is accelerated due to the scalar dipole radiation and the total number of \gw cycles is significantly decreased with respect to the general relativistic case.

\begin{figure}
	\centering
	\includegraphics[width=0.9\columnwidth]{figures/Shibata_2013pra_fig3a.pdf}
	\includegraphics[width=0.9\columnwidth]{figures/Shibata_2013pra_fig3b.pdf}
\caption{Evolution of the maximum values of the rest-mass density $\rho_{\rm max}$ and scalar field $\varphi_{\rm max}$ during \ns merger for several models with total mass $M = 2.7~{\rm M}_{\odot}$ with the APR4 \EOS (top panel) and the H4 \EOS (bottom panel). The merger sets in at the time as the maximum density steeply increases. We note that for $\beta_0 \ge -4.2$ with APR4 \EOS and for $\beta_0=-4.75$ with H4 \EOS, the scalarization already occurred at $t=0$. Credit:~\textcite{Shibata:2013pra}.}
	\label{fig:NS_merger}
\end{figure}

The actual merger and the post-merger phase were addressed only in~\textcite{Shibata:2013pra} which showed that dynamical scalarization can happen not only after the inspiral, but also during the merger since a very massive compact star is formed in this process. The evolution during \ns merger of the maximum values of the rest-mass density $\rho_{\rm max}$ and scalar field $\varphi_{\rm max}$ for several binary \ns{s} with total mass  $M = 2.7~{\rm M}_{\odot}$ and several values of $\beta_0$ is shown in Fig. \ref{fig:NS_merger}. Two \EOS{s} are considered in the two panels and the cosmological value of the scalar field is taken to be $\varphi_0=10^{-5}$. Some of the binary \ns{s} are not scalarized at the beginning of the evolution (the models with  $\varphi_{\rm max}=0$ at the beginning), but the values of $\beta_0$ are chosen such that all of them develop scalar fields at a certain point of the evolution. The actual merger of the two \ns{s} is marked by the rapid increase of $\varphi_{\rm max}$. Afterwards, either a hyper- or supramassive \ns is formed, or the merger remnant collapses to a bald \bh and the scalar field is radiated away.
In agreement with the studies of equilibrium differentially rotating \ns{s} discussed in Sec.~\ref{sec:DEF_model_equil}, the scalarized merger remnant can sustain a larger mass without collapsing to a \bh~\cite{Doneva:2018ouu}. This is evident in the right panel of Fig.~\ref{fig:NS_merger} where collapse to a \bh is observed in pure \gr, while a supramassive \ns is formed after the merger scalar-tensor gravity for small enough $\beta_0$.

The quasiperiodic oscillations of the merger remnant were examined by \textcite{Shibata:2013pra} showing a clear distinction compared to the \gr case. Such oscillations were studied in a series of papers in \gr (see e.g. \textcite{Bauswein:2011tp,Bauswein:2012ya,Hotokezaka:2013iia,Bauswein:2014qla,Takami:2014zpa,Clark:2014wua,Takami:2014tva,Rezzolla:2016nxn,Maione:2016zqz}) mainly as a tool to determine the nuclear matter \EOS from the post-merger GW signal.
The observed differences with the scalarized case can potentially be used to discriminate between \gr and modified gravity theories. This is specially the case, because scalarized merger remnants are produced only if specific initial conditions, related to the mass of the merging compact objects and the specifics of the \EOS, are met~\cite{Shibata:2013pra}.

\textcite{Taniguchi:2014fqa} took a different approach, where quasi-equilibrium sequences of binary \ns{s} at different time instants were calculated instead of performing time evolution.
This approach is well known in \gr~\cite{Gourgoulhon:2000nn,Taniguchi:2003hx,Taniguchi:2010kj}, and it is assumed that the characteristic time of the system to settle to an equilibrium is much smaller than the inspiral time scale. This is supposed to give a relatively accurate picture of the binary evolution even close to the merger. The results are in agreement with \textcite{Shibata:2013pra}, while the small deviations relative to~\textcite{Barausse:2012da} are probably due to the fact that \gr initial data were used in the latter.
\textcite{Taniguchi:2014fqa} showed that the absolute value of the binary binding energy is smaller than in \gr. In addition, the \gw cycles prior to the merger were significantly reduced compared to pure \gr once scalarization kicks in, and the effect is considerably stronger compared to the one due to tidal interactions.

The \pn approximation has also been a valued tool to model \ns inspiral
in the \DEF model. Initial work has focused in tensor-multi-scalar theories
and simple scalar-tensor theories that do not allow scalarization~\cite{Damour:1992we,Damour:1995kt,Mirshekari:2013vb,Lang:2013fna}.
This approach cannot be immediately applied to the case dynamical scalarization since this is nonperturbative effect absent in the weak field regime.
\textcite{Palenzuela:2013hsa} was the first to address this issue. They used the
equations of motion at 2.5PN order, derived in scalar-tensor gravity of~\textcite{Mirshekari:2013vb}, and supplemented them with a ``Newtonian feedback'' mechanism to account scalarization in the inspiral.
The resulting inspiral evolution was found to be in agreement with the numerical simulations of~\textcite{Barausse:2012da}. The approach of \textcite{Palenzuela:2013hsa}
is computationally inexpensive, and this allowed them the study large portions of the parameter space, including (un)equal-mass and eccentric binaries. This approach can also be used to efficiently generate inspiral \gw templates in the \DEF model.
A further step forwards was taken by~\textcite{Sennett:2016rwa}, that introduced a methodology called post-Dickean expansion. Their main improvements with respect to \textcite{Palenzuela:2013hsa} is a 1PN extension of the feedback mechanism.
The post-Dickean expansion was compared against the quasi-equilibrium calculations of~\textcite{Taniguchi:2014fqa}, and it was shown that it can accurately predict the onset and magnitude of dynamical scalarization.

\textcite{Sennett:2017lcx} took a different approach in the perturbative study of the inspiral, that consisted of an analytic model of dynamical scalarization using an effective action. The motivation was to cure two deficiencies in the previous \pn studies.
In effective action approach, the nonlinear scalarization process is reduced to a pair of cubic equations that have a closed-form solution depending on the binary separation when dynamical scalarization happens and on the magnitude of the developed scalar charge. This simplifies the problem relative to \textcite{Palenzuela:2013hsa} and \textcite{Sennett:2016rwa}.
In addition, the effective action approach allows one to construct a simple two-body Hamiltonian, that can be used to compute the binary's binding energy.~\textcite{Sennett:2017lcx} used this stablished dynamical scalarization as a second-order phase transition.
\textcite{Khalil:2019wyy} extended these results to general theories admitting scalarization, either for \bh{s} or \ns{s}, and are valid for adiabatic (quasistationary) and quasicircular orbits.
\textcite{Khalil:2022sii} took a step further, where these approximations were dropped and the dynamical evolution around the phase transition to the scalarized regime was studied. The results showed that in some cases, assuming a quasistationary evolution might not be accurate enough even for quasicircular binaries.

At last, \textcite{Ponce:2014hha} studied the effect of the scalar field on the electromagnetic radiation emitted in the merger of magnetized \ns{s}. They found that deviations in the emitted electromagnetic flux due to scalarization are not negligible yet challenging to measure.
However, if combined with \gw observations, constraints on scalar-tensor theory can, in principle, be placed, showing the usefulness of multi-messenger astronomy.

\subsection{Astrophysical implications of scalarized \ns{s} in the \DEF model}

In this subsection we will discuss in greater detail the astrophysical implications of scalarized \ns{s}. Many aspects have already been covered in the previous sections, such as the binary pulsar observations, \ns oscillations, stellar core-collapse, and binary mergers, as they naturally appeared in the presentation. Here, we will shed further light on possible astrophysical implications of scalarization, trying to be as complete as possible under two main parts. First, we will discuss the astrophysical implications directly related to electromagnetic observables. Afterwards, we will address the problem of universal relations for scalarized \ns models.

\paragraph{Electromagnetic observations:}
Scalarized \ns{s} in the \DEF model and its extensions were studied in a variety of astrophysical scenarios in attempt to probe the existence of the scalar field.
Since \ns{s} are often times surrounded by accretion disks, it is natural to study the effect of a nontrivial scalar field on disk' properties.
The simplest model is called ``thin disk'', where particles are assumed to move on geodesics around the central compact object.
Any small perturbation acting upon these particles will lead to oscillations around their equilibrium orbit with some characteristic (epicyclic) frequencies. Epicyclic and
orbital frequencies are thought be, in one way or another, related to the interpretation
of different accretion disk properties. An example are the \qpos observed in the spectrum of low-mass x-ray binaries.
For this reason, some works calculated these frequencies in the spacetime of scalarized \ns{s}. \textcite{DeDeo:2004kk} did this for nonrotating \ns{s},~\textcite{Staykov:2019pwj} included slow-rotation and also a mass to the scalar field, while \textcite{Doneva:2014uma} studied rapidly-rotating \ns{s}.
The conclusion of these works is that if one takes into account current observational constraints, only the cases of rapid rotation and massive scalar fields leave room
for significant effect of the scalar field.
Yet, a question that deserves further works is how one can disentangle uncertainties in the \ns \EOS and in the disk model from modifications to \gr.

Another promising observational property of accretion disks around a compact objects
is the shape of the Fe~$K_\alpha$ fluorescent line at 6.4~keV. This line is expected to be observed in greater details with future x-ray mission such as ATHENA.
\textcite{Bucciantini:2020owm} argued that both the intensity of the low-energy tails and the position of the high-energy edge of the line change due \gr modification are potentially observable in the future.

Accretion onto \ns{s} can also trigger a process called ``gravitational phase transition'', named in analogy with matter phase transitions from confined nuclear matter to deconfined quark matter~\cite{Kuan:2022oxs}. This process can happen when the maximum mass of the scalarized \ns{s} is smaller than the maximum mass of the zero scalar field (\gr) solution. The idea is then that if a scalarized \ns close to this maximum mass accretes some matter it may pass beyond the stability point. In \gr this would cause a collapse to a \bh. In scalar-tensor theory, the star will radiate its scalar hair and evolve towards the zero scalar field (i.e.,~\gr) branch. A significant amount of \gw{s} can be produced in this process, potentially detectable with the next generation of \gw detectors.

The \ns surface as an emitter of x-ray radiation is also an important probe in strong-field gravity, because the electromagnetic radiation is emitted from a
region with large spacetime curvature.
In this regard, the simplest observable is the gravitational redshift of surface atomic lines. The redshift carries information about the \ns{s} mass, radius, spin, and in scalar-tensor gravity, about the scalar field. \textcite{DeDeo:2003ju} showed that scalarization has significant effect on the redshift, but only for negative enough values of $\beta_0$. Such values are already ruled out by binary pulsar observations. One could entertain the idea that a small scalar field masses could allow for large deviations in the redshift relative to \gr, while reconciling the theory with with binary-pulsar constraints~\cite{Popchev2015,Ramazanoglu:2016kul,Doneva:2016xmf}. A study of this problem has not been performed yet.

The emission of x-ray pulse profiles emitted by hot spots at \ns surfaces have also received considerable attention recently. Observations of these signals allow for a relatively clean inference of \ns masses and radii~\cite{Watts:2016uzu}. This potential was met with observations~\cite{Bogdanov:2019ixe,Miller:2019cac,Riley:2019yda} by the \NICER~\cite{Arzoumanian2014:SPIE}.
\textcite{Sotani:2017rrt} and~\textcite{Silva:2018yxz} developed pulse-profile models for \ns{s} in massless \DEF model.~\textcite{Xu:2020vbs} considered the case of massive scalar fields and~\textcite{Hu:2021tyw}
studied the positive $\beta_0$ case. \textcite{Silva:2019leq} made the first study of what constraints can be placed on the \DEF model with pulse-profile observations and found they can be competitive with binary pulsar observations. However, this work is rather simple and not at the level of realism found in analysis of real data. This remains an important avenue for future work.

\begin{figure}
\includegraphics[width=0.95\columnwidth]{figures/constant_mphi_likelihoods_review_paper_4.pdf}
\caption{The posterior probability density for  $\beta_0 \in [-100,\ 0]$ for various assumed values of the scalar field mass $\mu$ in the massive version of the \DEF model introduced by \textcite{Ramazanoglu:2016kul} (marginalized over various \EOS). One can obtain the bound $\beta_0 \gtrsim -20$ for $\mu \lesssim 2\times 10^{-11}$~eV (left panel), but there is no effective bound for higher masses, at least in the interval $\beta_0 \in [-100,\ 0]$ (right panel). No bound is possible for $\mu \gtrsim 10^{-10}$~eV using this data, since scalarization does not occur for such high $\mu$ values within the considered $\beta_0$ interval. Credit:~\textcite{Tuna:2022qqr}.
}
\label{fig:mass_radius_constraint}
\end{figure}

At last, \textcite{Tuna:2022qqr} used the \ns mass and radius measurements~\cite{Ozel:2015fia,Bogdanov:2016nle} to constrain the massive extensions of the \DEF model. They obtained a weak lower bound $\beta_0 \gtrsim -20$ for scalar field masses $\mu \lesssim 2 \, \times10^{-11}$~eV; see Fig.~\ref{fig:mass_radius_constraint}. This is significant since no other bound is known for $\mu \gg 10^{-16}$~eV. These results show that large scalar masses enable agreement with observations, even for extremely negative $\beta_0$, and demonstrate the difficulty to constrain scalarization when the scalar field is massive.

\paragraph{Universal relations:}

One of the biggest obstacles in using \ns observations to test modified theories of gravity is the uncertainty in the \ns \EOS, which remains unknown at high-densities~\cite{Lattimer:2015nhk,Baym:2017whm}. In general, modifications to \ns properties predicted by different \EOS{s} are degenerate with changes to the underlying gravity theory used to model these stars. One way to break this degeneracy is to consider relations between different \ns properties that depend weakly on the \EOS~\cite{Yagi:2016bkt,Doneva:2017jop}. In fact, in Sec.~\ref{sec:NS_QNMs}, we had already met one such example relating the \qnm{s}' frequencies and damping times and the mean density of \ns{s}. Here we will focus on other \EOS-independent relations connecting various \emph{equilibrium properties} of scalarized \ns{s}.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{figures/Doneva_2014faa_Fig2.pdf}
\caption{Comparison between ${\bar I}-{\bar Q}$ relations for \gr and scalar-tensor theory, for several values of $\beta_0$, where ${\bar I}=I/M^3$ and ${\bar Q}=QM/J^2$. The presented data are only for APR \EOS in order to have better visibility, but the deviations for other \EOS are very small, typically below roughly $2\%$. Sequences for different values of the normalized rotational parameter $\alpha=M\Omega/2\pi$ are given, where $\Omega$ is the angular velocity of the star. Credit:~\textcite{Doneva:2014faa}.}
\label{fig:NS_STT_ILoveQ}
\end{figure}

A class of universal relations that attracted considerable attention is the I-Love-Q relations, connecting the normalized moment of inertia, tidal Love number and quadrupole moment of \ns{s}~\cite{Yagi:2013bca,Yagi:2013awa}. If any two quantities in the I-Love-Q trio are measured independently, one can constrain deviations from \gr by checking whether these observations agree with the \gr I-Love-Q relation within the observational accuracy. \textcite{Silva:2020acr} applied this strategy to constrain dynamical Chern-Simons gravity~\cite{Jackiw:2003pm,Alexander:2009tp} using observational data from \LVK and \NICER. The application to scalarization models has not been done yet.

The I-Love-Q relations for scalarized \ns{s} were studied in slow~\cite{Pani:2014jra} and rapid~\cite{Doneva:2014faa} rotation, including the case of massive scalar fields~\cite{Hu:2021tyw,Doneva:2016xmf}. In Fig.~\ref{fig:NS_STT_ILoveQ} we show the relation between the normalized moment of inertia $\bar{I}$ and quadrupole moment $\bar{Q}$ for sequences of scalarized \ns{s} with fixed rotational parameter $\alpha=M\Omega/2\pi$ is given in.
We see that for small absolute values of $\beta$, e.g., $\beta_0=-4.5$, the deviations from \gr are practically negligible for slow rotation (e.g.,~$\alpha=0.32$) while they are enhanced by rapid rotation. As we already saw, massive scalar fields allow much smaller values $\beta_0$, while remaining consistent with observations. Hence, these scalarization models allow for larger differences between the scalarized \ns universal relations and the \gr ones~\cite{Doneva:2016xmf}.

The moment of inertia $I$ and the quadrupole moment $Q$ are among the leading order multipole moments in the asymptotic expansion of the metric functions at infinity. There are infinitely many other higher order multipole moments, though, and it is interesting to explore whether similar universal relations hold for them. This was first addressed in \gr~\cite{Pappas:2013naa,Yagi:2014bxa} where such universal relations were derived for the higher multipole moments when proper normalization is applied. \textcite{Pappas:2014gca} and \textcite{Pappas:2018csu} extended these results to scalarized \ns{s}. This required the generalization of the multipole moment formalism of~\textcite{Geroch:1970cd,Hansen:1974zz} to scalar-tensor gravity.
\textcite{Pappas:2018csu} showed that future observations of \qpos of low-mass x-ray binaries can,
in principle, be used can be measure different \ns properties and tell apart different gravity.

Another class of universal relations studied in the context of the \DEF model connects the \ns moment of inertia and compactness. This is motivated by \textcite{Lattimer:2004nj} (see also \textcite{Breu:2016ufb}), that used this relation to argue that pulsar timing
observations could lead to a measurement of the moment of inertia to within $10\%$.
These relations were generalized to scalarized \ns{s} with massless scalar field potential by~\textcite{Motahar:2017blm}, and later to massive self-interacting scalar fields by~\textcite{Popchev:2018fwu}. Similar to the I-Love-Q relations, only the massive scalar field case leads to large deviations from \gr when observational constraints are taken into account.

\textcite{Ofengeim:2020zuc} proposed to connect the physical parameters of static \ns{s}, such as mass, radius, central energy density, pressure, and sound speed, at the maximum-mass point for a given \EOS. This is based on the observations that the nuclear matter \EOS (without phase transitions) can be well parameterized by only two parameters~\cite{Lindblom:2010bb}.  As a result, \textcite{Ofengeim:2020zuc} derived multiple constraints on the nuclear matter \EOS. \textcite{Danchev:2020zwn} generalized these relations to scalarized \ns{s}. They showed how sensitive these relations can to the underlying theory of gravity. Thus, all the \EOS restrictions derived in \gr should be taken with care, since the possibility for \gr modifications are rarely taken into account when interpreting observational data.

\subsection{Extended scalar-tensor theories beyond the Damour-Esposito-Far\`ese model}
\label{sec:ns_extended}

Until now we discussed scalarization in the \DEF model or other \DEF-inspired models. However, s we saw in Sec.~\ref{models}, spontaneous scalarization is also possible for other modified theories of gravity. Here, we will discuss such theories where scalarized \ns solutions were obtained.

\paragraph{Scalar-tensor theories with disformal coupling}

\textcite{Minamitsuji:2016hkk} studied \ns scalarization in scalar-tensor theories with disformal coupling of the form \eqref{eq:disformal}. More specifically, instead of the standard conformal transformation between Einstein frame metric ($g_{\mu\nu}$) and Jordan frame metric (${\tilde g}_{\mu\nu}$)
\begin{equation}
   g_{\mu\nu}={\cal A}^2(\varphi){\tilde g}_{\mu\nu},
\end{equation}
where ${\mathcal A}(\varphi)$ is a function of the scalar field related to the Jordan frame coupling between the scalar field and the Ricci scalar, we have a more general and complicated transformation that also involves scalar field derivatives:
\begin{equation}
    g_{\mu\nu}=A^2(\varphi) \left[{\tilde g}_{\mu\nu} + \Lambda B^2(\varphi)\nabla_\mu \varphi \nabla_\nu \varphi \right],
\end{equation}
where $B(\varphi)$ is a function of the scalar field and $\Lambda$ is a parameter of dimensions of [length]$^{2}$. The motivation behind this modification of the conformal factor comes from \textcite{Bekenstein:1992pj} that aimed to find the most general coupling constructed from the metric $g_{\mu\nu}$ and the scalar field $\varphi$ that respected causality and the \WEP. More recently, \textcite{Andreou:2019ikc} showed that such disformal coupling is actually equivalent to a kinetic coupling with the Ricci scalar and thus is contained in the minimal action \eqref{eq:ActionCaseI} (see Sec. \ref{sec:minimal_action} for a detailed discussion).

\textcite{Minamitsuji:2016hkk} showed that for negative values of the disformal coupling parameter $\Lambda$, scalarization can be suppressed, while for large positive values of $\Lambda$ stellar models cannot be obtained. This study also explored the universal relations between the moment of inertia and the compactness of \ns{s}, and they determined the range of parameters where these relations can deviate significantly from \gr.

\paragraph{Scalar-Gauss-Bonnet theory:}

Another well-studied modified gravity theory that allows for \ns scalarization is scalar-Gauss-Bonnet gravity discussed in Sec. \ref{models}. Interestingly, spacetime curvature itself as well as matter can be the source of the scalar field in this theory. Scalarized \ns solutions in scalar-Gauss-Bonnet gravity were studied for the first time in \textcite{Doneva:2017duq,Silva:2017uqg}. They exhibit quite different features compared to the \DEF model while being qualitatively similar to \ns solutions in other classes of Gauss-Bonnet theories where scalarization is not possible~\cite{Pani:2011xm,Kleihaus:2016dui,Saffer:2019hqn}. For example, for a given central energy density the stellar mass is always smaller compared to \gr. Moreover, normally there is only one bifurcation point at small central energy densities and afterwards the branches of solutions are terminated because of violation of the regularity conditions.

In Fig.~\ref{fig:NS_MR_GB}, we show the branches of scalarized \ns solutions for a coupling function $f(\varphi) = - \lambda^2 (2 \beta)^{-1} [1 - \exp(-\beta \varphi^2/4)]$. Here, $\lambda$ is a parameter with dimensions of [length]$^{-1}$ that controls the coupling strength in the action \eqref{actionsgb}. Note that the sign of $f(\varphi)$ is opposite to the one discussed in Sec. \ref{sec:theory_sgb} and thus can not lead to \bh spontaneous scalarization since in vacuum $\mu_{\eff}^2>0$ as defined in \eqref{eq:mueffsqr_sgb}. However, due to the presence of matter \ns{s} can scalarize. Scalarized \ns solutions were found also for the more ``conventional'' sign of $f(\varphi)$, namely $f(\varphi) = \lambda^2 (2 \beta)^{-1} [1 - \exp(-\beta \varphi^2/4)]$ \cite{Doneva:2017duq}.

\begin{figure}
\includegraphics[width=0.95\columnwidth]{figures/Doneva_2017duq_fig2.pdf}
\caption{Scalarized \ns{s} in scalar-Gauss-Bonnet gravity. We show the mass as a function of the radius for \gr (trivial branch with $\varphi=0$) as well as several scalarized branches with different values of the scalar field coupling parameters $\lambda$ and $\beta$. Credit:~\textcite{Doneva:2017duq}.}
\label{fig:NS_MR_GB}
\end{figure}

Little is currently known about the astrophysical implications of these stars. \textcite{Kuan:2021lol} considered the spherically symmetric core collapse of a non-compact star either to a protoneutron star or to a \bh in Gauss-Bonnet theory. They also proposed a realistic physical mechanism for the formation of isolated scalarized \bh{s} and \ns{s} The complexity of the problem is greatly increased with respect to \gr, though, and there are fundamental difficulties such as the loss of hyperbolicity of the evolution equations for certain ranges of parameters, effectively limiting the maximum possible scalar field to relatively low values (see also \textcite{East:2021bqk,East:2020hgw,Ripley:2020vpk} for the case of \bh{s}). Core-collapse simulations show that the remnant in scalar-Gauss-Bonnet theory can be very rich with (de)scalarization happening at intermediate or final stages of the collapse depending on the properties of the progenitor and the theory parameters. Since breathing modes are absent in this theory, the effect on \gw emission can be estimated only if one drops the assumption of spherical symmetry. This is something that has not been done yet in any modified gravity theory due to problem's complexity.

\paragraph{The Ricci-Gauss-Bonnet model:}

One can further modify the action in scalar-Gauss-Bonnet gravity to include additional terms. This is the case of Ricci-Gauss-Bonnet model discussed in Sec.~\ref{sec:mixed_model}. This model has interesting advantages such as the possibility to reconcile scalarization with cosmology. Here we will focus on the \ns solutions this theory that were considered by \textcite{Ventagli:2021ubn}. This work made thorough exploration of the theory parameter space in order to find the sectors where \ns solutions exist. Since their scalar charge is nonzero, in these sectors one can put severe constraints on the theory based on the binary pulsar observations, similar to the \DEF model (see Sec. \ref{sec:ns_binary_pulsar}).

The free parameters in the theory, as evident from the action \eqref{action:RGmixed}, are $\beta$ and $\alpha$ that control the coupling strength to the Ricci scalar and the Gauss-Bonnet invariant respectively. In Fig. \ref{fig:NS_MixedModel_MPA1} the existence of \ns solutions is shown in a two-dimensional plot spanning the parameters $\alpha$ and $\beta$ in a broad range, for a fixed central energy density and the MPA1 \EOS~\cite{Muther:1987xaa}. The white area corresponds to the region of the parameter space where the \gr solution is stable against scalar perturbations. Throughout the parameter space that is spanned by $\alpha$ and $\beta$, a new unstable mode appears every time one crosses a black line. These unstable modes can be labeled by the number of scalar field nodes (denoted by $n=0,1,2,3$ in the figure).  Any point in the parameter space that lies within some grey region corresponds to a configuration where the \gr neutron star solution is unstable. The red (respectively
blue) area corresponds to the region where scalarized solutions with $n=0$ (respectively $n=1$) nodes exist.

As one can see, scalarized solutions exist in only part of the parameter space. In the $\beta>0$ and small $\alpha$ range (in white in the figure) the \gr neutron starts are stable. This is exactly the region most interesting from the point of view of \bh scalarization and connection to cosmology \textcite{Antoniou:2020nax,Antoniou:2021zoy}, see also Sec.~\ref{sec:bhs_stability_nr}. Thus, including a Ricci coupling to the scalar-Gauss-Bonnet action seems to allow for \bh scalarization while evading the binary pulsar constraints.

Fig. \ref{fig:NS_MixedModel_MPA1} is only for one \EOS and a specific central energy density. Different \EOS{s} and central densities can lead to significant deformations of the existence and stability regions. The general point made in~\textcite{Ventagli:2021ubn}, though, is that a stability (white) region always survives.

\begin{figure}[t]
	\includegraphics[width=\columnwidth]{figures/MPA1n1.pdf}%
	\caption{ Regions of existence of scalarized solutions in the $(\alpha,\beta)$ plane for the MPA1 \EOS with central energy density $\varepsilon_c=6.3\times 10^{17}~\text{kg}/\text{m}^3$. In \gr, a star with this choice of $\varepsilon_c$ and \EOS is light, with $M_\text{GR}=1.12~{\rm M}_{\odot}$. The red (respectively blue) region is the region where scalarized solutions with $0$ (respectively $1$) scalar field node exist. The grey contours obtained by~\textcite{Ventagli:2020rnx} are superimposed, which represent the lines beyond which \gr solutions with the same density are unstable to scalar perturbations with $0, 1, 2,$ etc nodes. We see that the region where there exist scalarized solutions with $n$ nodes is included in the region where the \gr solutions are unstable to scalar perturbations with $n$ nodes, but much smaller. The dashed boundary for the blue region corresponds to a breakdown of the integration inside the star. Credit:~\textcite{Ventagli:2021ubn}.}
	\label{fig:NS_MixedModel_MPA1}
\end{figure}

\paragraph{Tensor-multi-scalar theories:}
\label{sec:TMST}

Another class of gravity theories where \ns scalarization was considered, is the \TMST~\cite{Damour:1992we}, whose basics we discussed in Sec.~\ref{models}. In this theory the gravitational interaction is mediated  by the spacetime metric $g_{\mu\nu}$ and $N$ scalar fields $\varphi^{a}$  which take values in a coordinate patch of an $N$-dimensional Riemannian (target) manifold ${\cal E}_{N}$  with a positive definite metric $\gamma_{ab}(\varphi)$ defined on it \cite{Damour:1992we,Horbatsch:2015bua}.

The main features of \TMST is the inclusion of more than one scalar field and a structure called the target space metric. In that sense, there are two directions to go in order to obtain scalarized \ns{s}. The first one is to consider a mixture of several, more or less equivalent, scalar fields that was the approach in \textcite{Horbatsch:2015bua}. There, the authors examined the case of two scalar fields in the form of a complex scalar with a maximally symmetric target space metric $\gamma_{ab}$.
\textcite{Doneva:2020afj} studied scalarization in tensor-multi-scalar gravity in a more complicated setup, namely when $\gamma_{ab}$ is a three-dimensional maximally symmetric space together with a nontrivial map $\varphi:$~spacetime~$\rightarrow$~target~space.
While the solutions in \textcite{Horbatsch:2015bua} can be viewed as generalization of the \DEF model to multiple scalar fields, the scalarized \ns{s} in \textcite{Doneva:2020afj} have some quite distinct properties and they can be considered more as a limiting case of the topological \ns{s} discovered by \textcite{Doneva:2019ltb}. More specifically, the scalar field is zero in the stellar center, while only its first derivative is zero in the \DEF model. The scalar charge for the scalarized stars in tensor-multi-scalar theories is zero as well. This automatically reconciles this theory with the binary pulsar observations (due to the absence of scalar dipole radiation), while still allowing for large deviations from \gr. Another interesting property is the fact that there is non-uniqueness of the scalarized solutions with respect to the central energy density -- something that was never observed before at least for scalarized \ns{s}. Fig.~\ref{fig:TMST_Scalarization_NS} illustrates this property. Interestingly, \textcite{Kuan:2021yih} showed that all these solutions are stable until the maximum \ns mass regardless of the non-uniqueness.

\begin{figure}
	\includegraphics[width=0.94\columnwidth]{figures/Doneva_2019ltb_Fig4.pdf}
\caption{Tensor-multi-scalar theories: the mass as a function of the central energy density (left panel) and as a function of the stellar radius (right panel) for scalarized \ns{s} in tensor-multi-scalar theories where the target space metric represents a three-dimensional maximally symmetric space and we have a nontrivial map $\varphi:$~spacetime~$\rightarrow$~target~space.  Sequences of scalarized solutions for different combinations of parameters are given. Observe the appearance of non-uniqueness of the solutions with nonzero scalar field for certain combinations of the parameters (e.g., the small central energy density region for the blue dotted line). Credit:~\textcite{Doneva:2020afj}. }
	\label{fig:TMST_Scalarization_NS}
\end{figure}

\paragraph{Other scalarization models:}
Apart from models of scalarized \ns{s} already discussed, other works have attempted to consider scalarization in other more exotic theories of gravity or even different types of scalarization.

For example,~\textcite{Azri:2020agc} considered \ns scalarization in what is called scalar-connection gravity, where gravity is mediated by a scalar field and the connection.
Other types of scalarization in standard scalar-tensor theories, but if we allow for different scalar field potentials. \textcite{Minamitsuji:2022qku} showed this by
considering the potential $V(\varphi)=m_\phi^2 f_B^2 [1 + \cos(\phi/f_B)]$ of a pseudo-Nambu-Goldstone boson. Here $(\mu_\phi$, $f_B)$ are constants with dimensions of length and the scalar field has negative mass, $\mu^2 = - m_\phi^2$. \textcite{Minamitsuji:2022qku} showed that in this case the scalar field sits at its vacuum expectation value far away from the source, while inside a \ns, a symmetry restoration can take place resulting in a new type of scalarization. This theory also has the advantage that it avoids cosmological instabilities present in some other scalarization models.

\subsection{Spontaneous vectorization}
\label{sec:vectorization}
What happens if we replace the scalar field of spontaneous scalarization with a vector? In scalarization, we started with the most general context in Sec.~\ref{sec:minimal_action}, considering all the allowed couplings to the metric. Here we will follow the reverse path, starting with individual models of vector-tensor theories, and considering the more encompassing theory later. This exposition is preferable because the study of specific models, that are straightforward generalizations of known scalarization models, reveal some pathologies, and provides guidance for further model building.

The idea and the first specific model of vectorization was introduced in analogy with a massive version of the \DEF model in Eq.~\eqref{eq:action_st}.
Consider the vector-tensor theory action~\cite{BeltranJimenez:2012kby,Ramazanoglu:2017xbl}
\begin{align}\label{action_vt}
S &= \frac{1}{16\pi G_{*}} \int {\rm d}^{4}x \, \sqrt{-g} \, \left[R - F^{\mu\nu} F_{\mu\nu} - 2 \mu_X^2 x \right] \nonumber \\
&\quad + S_{{\rm m}}\left[\Psi_{{\rm m}};\, {\cal A}^{2}(x) g_{\mu \nu} \right],
\end{align}
where $F_{\mu\nu} = \nabla_\mu X_\nu -\nabla_\nu X_\mu$, $x =g^{\mu\nu}X_\mu X_\nu$. The vector field $X_\mu$ has the canonical kinetic term in this frame, and the matter degrees of freedom couple to the metric $\tilde{g}_{\mu\nu}={\cal A}^{2}(x) g_{\mu \nu}$ which is conformally scaled with respect to $g_{\mu\nu}$. The vector field equation is
\begin{align} \label{eom_vt1}
\nabla_\rho F^{\rho \mu}& = [-8\pi G_{*}{\cal A}^4 \Lambda \tilde{T} + \mu_X^2 ] \, X^\mu,
\end{align}
where $\Lambda \equiv {\rm d} \ln {\cal A} / {\rm d} x$, $T_{\mu \nu}$ is the energy-momentum tensor in the Einstein frame, and $\tilde{T}$ is the trace of the stress-energy tensor in the Jordan frame, i.e., with respect to the metric $\tilde{g}_{\mu\nu}$.

Equation~\eqref{eom_vt1} is that of a massive vector (Proca) field where the expression inside the square brackets acts as an effective mass $\mu_{\rm eff}^2$. In parallel to the \DEF model, when $\cal A$ has an appropriate form, e.g. ${\cal A}= \exp(\beta_0 x/2)$ with sufficiently negative $\beta_0$, $\mu_{\rm eff}^2$ becomes negative. Furthermore, for dense enough objects such as \ns{s}, this occurs for order-of-unity values of $\beta_0$. The expectation is that the vector field grows from arbitrarily small perturbations around $X_\mu = 0$ due to this tachyonic behavior, which can be called spontaneous vectorization, in exact analogy to spontaneous scalarization. However, we will later see that there are many subtle points in this narrative.

Even though we presented spontaneous vectorization in the Einstein frame with a nonminimal coupling to matter, it can be converted to a theory of minimal matter coupling and vector-curvature couplings as in spontaneous scalarization~\cite{Ramazanoglu:2019tyi}. Even more directly, there are theories of spontaneous vectorization that are purely conceived through curvature coupling, the first example being the Hellings-Nordtvedt theory~\cite{Hellings:1973zz} studied by~\textcite{Annulli:2019fzq}
\begin{align}\label{action_vec_hl}
    S &= \frac{1}{16 \pi G} \int {\rm d}^{4}x \, \sqrt{-g} \,
    \left[ R - F^{\mu\nu} F_{\mu\nu} -\Omega x R \right.
    \nonumber \\
    &\quad \left. - \, \eta X^\mu X^\nu R_{\mu\nu} \right] + S_{{\rm m}} [\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{align}
where $\Omega$ and $\eta$ are coupling constants. Scalarization through coupling to the Gauss-Bonnet term also has a vector analog as in the action~\cite{Ramazanoglu:2019gbz}
\begin{align}\label{action_vec_gb}
    S &= \frac{1}{16 \pi G} \int d^4x \sqrt{-g} \left[ R - F^{\mu\nu} F_{\mu\nu} + f(x) \mathscr{G} \right]
    \nonumber \\
    &\quad + S_{{\rm m}} [\Psi_{{\rm m}};\, g_{\mu\nu}],
\end{align}
which has a similar structure to Eq.~\eqref{actionsgb}. Action~\eqref{action_vec_gb} leads to \bh vectorization as well, unlike the other models we have seen so far.

The models of vectorization we have examined are only specific examples, and any theory of scalarization can be generalized to vector fields in principle. In this regard, \textcite{Brihaye:2020oxh} studied the coupling of a vector field to a scalar, and \textcite{Oliveira:2020dru} and \textcite{Brihaye:2021qvc} studied the coupling between two vector fields to obtain vectorized compact objects. \textcite{Kase:2020yhw} considered generalized Proca theories~\cite{Heisenberg:2014rta,Tasinato:2014eka} with various different couplings of a massive vector field, and showed that spontaneous vectorization occurs in these models as well. \textcite{Ramazanoglu:2018tig} investigated models where the effective mass of the vector field is generated by the Higgs-mechanism, preserving the gauge symmetry. On a separate front, conformal scaling of the metric in the matter action~\eqref{action_vt} can be generalized to disformal transformations, and spontaneous vectorization still occurs~\cite{Ramazanoglu:2019jrr,Minamitsuji:2021rtw}. In terms of approximate solutions, \textcite{Garcia-Saenz:2022wsl} recently calculated the Schwarzschild \qnm{s} of nonminimally coupled vector fields. There have also been efforts to study all these phenomena in a unified manner in compact binaries, using more generic \EFT tools~\cite{Khalil:2019wyy}.

To summarize, coupling of the vector field to any non-vanishing term in the action can be considered for spontaneous vectorization~\cite{Ramazanoglu:2019gbz}, and most options have been considered at least in a preliminary sense. The case of all possible couplings, an analog of the minimal action of scalarization in Sec.~\ref{sec:minimal_action}, has also been studied by~\textcite{Garcia-Saenz:2021uyv}. This work also revealed some of the fundamental problems of vectorization, as we will discuss more broadly later.

Vectorizations models strongly resemble scalarization models in terms of their actions, but this can be misleading. One major difference is in the number of additional degrees of freedom. Scalar fields contribute a single new degree of freedom irrespective of whether they are part of spontaneous scalarization or not, provided that the field equations are second order in derivatives. For vectorization, this is not the case. The action~\eqref{eom_vt1} and all the vector-tensor theory actions we have discussed break the gauge freedom found in massless vector fields.\footnote{Even though we have an intrinsic mass $\mu_X$ for the vector field in our discussion, vectorization can occur without this term, as is the case for scalarization.} As a result, the vectors of vectorization models carry three degrees of freedom instead of the two found in electromagnetism. This is not an immediate reason of concern, since minimally coupled  massive vectors, known as Proca fields~\cite{Proca:1936fbw}, also break the gauge symmetry, and they still provide a classical field theory with well-behaved dynamics. However, as we will see, the extra degree of freedom is far more problematic in vectorization.

Although the vectorization models are designed to have a linear tachyonic instability in analogy with scalarization, the presence of the third degree of freedom, whose dynamics is not obvious, casts doubt on the success of this design. Namely, whether vectorization can indeed be described as a tachyonic instability of the vector modes, that is then quenched nonlinearly is not apparent. This concern is amplified by the fact that vectorized compact objects in the models studied so far seem to have some striking differences with respect to their scalarized counterparts.

As an example, consider vectorized \bh{s} in the theory described by action~\eqref{action_vec_gb}. Entropy can be a measure of stability for \bh{s}, where higher entropy solutions are favored in terms of stability. Entropy of spherically symmetric vectorized \bh{s} were numerically computed by \textcite{Barton:2021wfj}, and \gr \bh{s} were shown to be entropically favored over vectorized ones in all cases. The verdict of stability is not definite without time evolution, but this is a strong indication for the instability of vectorized \bh{s}. In contrast, scalarized \bh{s} in scalar-Gauss-Bonnet theories, e.g.,~Eq.~\eqref{actionsgb}, can be entropically favored to the ones in \gr for appropriate coupling functions~\cite{Doneva:2017bvd}, and their stability is also indicated by other methods (see Sec.~\ref{sec:bhs_stability_nr}), which makes the stability of vectorized objects even more suspect.

\begin{figure}
\includegraphics[width=\columnwidth]{figures/MassiveVT_3by3_v3.pdf}
\caption{
The strength of vectorization measured by ${\cal A}(r=0)-1$ as a function the \ns mass $M$ for various values of $\beta_0$ and $\mu_X$ in Eq.~\eqref{eom_vt1}. Rows: upper $\mu_{X}= 1.6 \times 10^{-11}$~eV, middle $\mu_{X}= 8.0 \times 10^{-12}$~eV, lower $\mu_{X}= 4.8 \times 10^{-12}$~eV. Columns: left $\beta_0 =-4$, middle $\beta_0=-5$, right $\beta_0=-6$. Dashed and dot-dashed lines indicate solutions where the energy density $\varepsilon$ or $\Phi=-n^\mu X_\mu$ ($n^\mu$ being the normal vector field to spatial hypersurfaces) do not monotonically decrease with radius, respectively.
Credit: \textcite{Ramazanoglu:2017xbl}.
}
\label{vec_parameter_space}
\end{figure}
Scalarized and vectorized \ns{s} also show major differences. For example, the dependence of the strength of vectorization on the \ns mass has no discernible pattern as can be seen in Fig.~\ref{vec_parameter_space}, whereas scalarization occurs in a finite \ns mass interval, and vanishes at the boundary mass values (see Fig.~\ref{fig:DEF_NSmodels_Nonrot}). It is also known that while scalarized \ns{s} are continuously connected to those of \gr as the theory parameters are smoothly changed, e.g. in the limit $\beta \to 0$ in Fig.~\ref{fig:DEF_NSmodels_Nonrot}, this is not the case in vectorization~\cite{Ramazanoglu:2017xbl,Minamitsuji:2020pak}. The latter is known to be an indication of instability in certain scalarization theories, see~\textcite{Mendes:2016fby}. Thus, qualitative differences between scalarized and vectorized objects are also apparent for \ns{s}, and moreover, point to the instability of the vectorized ones.

These observations suggest that the vectorized compact object solutions obtained so far may not be stable, hence they are not astrophysically relevant, at least for the case of spherically symmetric spacetimes. A close look at the instability mechanism of vectorization is going to reveal the underlying reasons for this in the following section.

\subsection{Vectorization and ghosts}
\label{sec:vectorization_ghosts}
Despite the fact that the vectorization phenomena we have presented are directly inspired by scalarization, we have seen noteworthy differences between the two mechanisms. Even though it was puzzling for a time, this is now strongly suspected to be related to the fact that vectorization models suffer from ghost and/or gradient instabilities in addition to tachyonic ones, as demonstrated by~\textcite{Garcia-Saenz:2021uyv} and~\textcite{Silva:2021jya}.
Here, we will first follow the latter, and later outline the methodology of the former.\footnote{\textcite{Esposito-Farese:2009wbc} is an earlier similar study that specifically investigates cosmological scenarios.}

Let us first review the basic aspects of different types of instabilities~\cite{Demirboga:2021nrc}.
Recalling our discussion of the tachyon in Secs.~\ref{tachyonquenched} and~\ref{sec:tachcurved}, consider the linearized scalar field equation~\eqref{phi4lincurv} in $1+1$ dimensions
\begin{align}
\label{eq:instability_eom}
g^{tt} \partial_t^2 \delta\varphi
+ g^{xx} \partial_x^2 \delta\varphi
=\mu^2 \varphi,
\end{align}
where we assume the metric is diagonal with constant components for simplicity. The common case without instabilities is when $g^{tt}<0$, $g^{xx}>0$, and $\mu^2>0$. What happens when each of these three terms changes its sign, while the other two stay the same?

For a plain wave mode $\delta\varphi(t,x) \propto \exp[{\rm i}\left(\omega t - kx\right)]$, the dispersion relation is given by
\begin{align}
\omega(k) = \sqrt{ (\mu^2+g^{xx} k^2) / (-g^{tt}) },
\end{align}
We saw that for $\mu^2<0$ the mode behaves as a tachyon, where $\omega(k)$ becomes imaginary for sufficiently small $|k|$. This leads to exponential growth as discussed, and the fastest growing mode behaves as $\sim \exp[\sqrt{\mu^2/g^{tt}} \, t]$, i.e., the growth rate has an upper limit.

A \emph{ghost} occurs when $g^{tt}>0$. There is also exponential growth in this case, however the growth rate diverges with increasing wave number as $\sim \exp[\sqrt{g^{xx}/g^{tt}}k t]$, hence has no upper bound.
The case $g^{xx}<0$ also gives the same asymptotic growth rate as the ghost, and is called a \emph{gradient instability}.\footnote{For the simple equation~\eqref{eq:instability_eom}, a simultaneous change of sign of two of the coefficients is equivalent to changing the sign of the third one, i.e. $g^{tt}>0$, $g^{xx}<0$ and $\mu^2>0$ is also a tachyon. Similarly, all three coefficients changing signs gives a stable theory. This is due to the fact that the equation as a whole or the action that leads to it can be multiplied by an overall factor of $-1$ without changing physics. However, the scalar, or any field we will investigate which carries instabilities, is always coupled to the metric and other terms in the action, which means the overall sign is meaningful. Hence, we will use the aforementioned classification, and say that e.g. a theory with $g^{tt}>0$, $g^{xx}<0$ and $\mu^2>0$ carries both ghost and gradient instabilities following the literature~\cite{Demirboga:2021nrc}.} In summary, ghost and gradient instabilities as we presented them here have qualitative differences from tachyons, a main one being the arbitrarily fast growth rates. This means a well-defined time evolution can be problematic in such theories, and it is harder to devise nonlinear quenching mechanisms.

The spontaneous vectorization model of Eq.~\eqref{action_vt} modifies the ``mass term'' in the field equation for the vector field. Hence, the first impression might be that it leads to a tachyonic instability as in scalarization. However, this approach is naive, and a thorough analysis reveals that theories of vectorization necessarily carry ghost instabilities in addition to tachyonic ones. Following \textcite{Silva:2021jya}, one way to understand this is using the Stueckelberg trick~\cite{Ruegg:2003ps} where we introduce a scalar field $\psi$ in the action~\eqref{action_vt} by the substitution
\begin{align}
X_{\alpha} \to X_{\alpha} + \mu_X^{-1} \nabla_{\alpha} \psi.
\label{eq:trick}
\end{align}
This leads to the scalar equation of motion
\begin{align} \label{eq:box_psi_eff_metric}
\overline{\Box}\psi = -\bar{g}^{\mu\nu} \,
[\, \mu_X \overline{\nabla}_\mu X_\nu +\tfrac{1}{2}
( \overline{\nabla}_\mu \log \hat{z} )
( \overline{\nabla}_\nu \psi+\mu_X X_\nu )
\,].
\nonumber \\
\end{align}
where the covariant derivative $\overline{\nabla}$ is compatible with the new metric
\begin{equation}\label{eq:geff}
\bar{g}^{\mu\nu}= \hat{z} \, g^{\mu\nu}, \qquad
\hat{z} \equiv \frac{\mu^2_{\rm eff}}{\mu_X^2} = 1-\frac{8\pi G_{*} {\cal A}^4 \Lambda \tilde{T}}{\mu_X^2}.
\end{equation}
The signature of $\bar{g}_{\alpha\beta}$ changes with the sign of $\hat{z}$, which controls the change of sign of the effective mass of $X_\mu$ in Eq.~\eqref{eom_vt1}. In other words, when we try to instigate a tachyon on $X_\mu$ by having negative $\hat{z}$ in parts of the spacetime, $\psi$ carries ghost and gradient instabilities [cf. Eq.~\eqref{eq:instability_eom}]. This is a degree of freedom whose dynamics is governed by the \emph{effective metric} $\bar{g}_{\mu\nu}$, not the spacetime metric $g_{\mu\nu}$.

It is curious to see that the Stueckelberg trick reveals the ghost, even though it is not apparent in the vector field equation~\eqref{eom_vt1}. However, a more careful analysis can clarify this issue~\cite{Silva:2021jya}. First note that the vector field equation~\eqref{eom_vt1} is not manifestly hyperbolic, i.e. its principal part, the term with the highest derivatives, is not the wave operator, since
\begin{equation}
\nabla_\mu F^{\mu\nu} = \nabla_\mu \nabla^\mu X^\nu -\nabla_\mu \nabla^\nu X^\mu.
\end{equation}
The second term on the right hand side of Eq.~\eqref{eom_vt1} vanishes for the Proca theory, ${\cal A}=1$, making the principal part the wave operator. This is not the case in general due to the so-called (generalized) \emph{Lorenz condition}
\begin{equation} \label{lorenz_vt}
\nabla_\mu \left(\hat{z}  X^\mu \right) = 0,
\end{equation}
which is obtained by acting on both sides of Eq.~\eqref{eom_vt1} by $\nabla_\mu$, and recalling the antisymmetry of $F^{\rho\mu}$. We can use Eq.~\eqref{lorenz_vt} to obtain a manifestly hyperbolic form of the linearized field equations on a fixed background of vanishing vector fields as
\begin{align}
\nabla_\mu \nabla^\mu X_\nu
+ ( \nabla_\mu \ln \hat{z} ) \nabla_\nu X^\mu
= \mathcal{M}_{\nu\mu} X^\mu ,
\label{eq:vec_component}
\end{align}
where we defined the mass-squared tensor
\begin{align}
\mathcal{M}_{\nu\mu} = \hat{z} \mu_X^2 g_{\nu\mu}  + R_{\nu\mu}
- \nabla_\nu \nabla_\mu \ln |\hat{z}| .
\label{eq:mass-tensor}
\end{align}
Here, the metric and $\hat{z}$ have their fixed background \gr values. This is a generalized massive wave equation.

Since $\hat{z}=1$ in the absence of matter, and it becomes negative in a continuous manner inside matter for an instability to exist at all, $\hat{z}$ necessarily vanishes at some points of a \ns spacetime. The $\nabla_\nu \nabla_\mu \ln \hat{z}$ term introduced by the vectorization contains powers of $\hat{z}^{-1}$, which means the effective mass-squared tensor diverges at such points. Alternatively, if we move the $\hat{z}^{-1}$ terms to the other side of the equation, the principle part becomes $\hat{z}\Box X_\nu$, i.e. it changes its sign with $\hat{z}$ the same way as the $\psi$ field in equation~\eqref{eq:box_psi_eff_metric}. Hence, if we properly rewrite the vector field equation into a manifestly hyperbolic form, the ghost is apparent.\footnote{There is no change of sign in the case of intrinsically massless vector fields, $\mu_X=0$, but it can be shown that the coefficients in the field equations still diverge~\cite{Silva:2021jya}.}

Equations~\eqref{eq:vec_component} and~\eqref{eq:mass-tensor} show that calling $\mu^2_\eff = \hat{z} \mu_X^2$ the effective mass squared of the theory is misleading, unlike the case of scalarization. Even though $\mu^2_\eff$ replaces the mass term in the vector field equation of a minimally coupled Proca field, it is not actually the effective mass of physically propagating degrees of freedom. The \emph{physical mass} is determined by Eq.~\eqref{eq:mass-tensor}, which diverges at points where $\hat{z}=0$.

One might at first think that the ghost and gradient instabilities can be regularized by nonlinearities. After all, the tachyon of scalarization is also an instability, but it is eventually quenched by nonlinear effects as explained in Sec.~\ref{sec:scalarization_mechanism}. In such a scenario, one would expect the vectorized objects to be free of instabilities, even though the \gr solutions are unstable. However, this is not the case.~\textcite{Demirboga:2021nrc} showed
that spherically symmetric vectorized \ns{s} in the theory~\eqref{action_vt} are ghost and gradient
unstable at the perturbative level.

Even if stable solutions can be found in other vectorization models, there is a more fundamental problem due to the modification to the principal part of the vector equations in Eq.~\eqref{eq:box_psi_eff_metric} or~\eqref{eq:vec_component}. The coefficient of the wave operator vanishes at certain points in spacetime, or equivalently, the physical effective mass diverges if we move these coefficients to the other side of the equation. In  either case, this leads to divergent time derivatives for arbitrarily small perturbative vector field values, which renders the initial assumptions of a perturbative approach invalid. Overall, it is not clear whether one can even define a well-posed time evolution, let alone one that can somehow suppress the exponential growth of the vector field~\cite{Silva:2021jya}. Further elucidation of these issues requires a more rigorous mathematical investigation, which has not been attempted yet.

We  have so far investigated the ghost and gradient instabilities of the specific vectorization theory in action~\eqref{action_vt}, but similar results are known to exist for all theories of vectorization~\cite{Garcia-Saenz:2021uyv,Silva:2021jya} due to similar mechanisms. To summarize, replacing the scalar of scalarization with a vector field results in theories  where the central issue is not the nature of the field the instability lives on, but the type of instabilities that arise. We cannot construct vectorization models that only exhibit the tachyonic instability commonly associated with the scalarization mechanism, which can be benign and tamed by nonlinear effects. Rather, vectorization also exhibits far more threatening instabilities.

Above, we closely followed \textcite{Silva:2021jya}, but \textcite{Garcia-Saenz:2021uyv} reach similar conclusions via an alternative approach. The authors consider the most general action that contains a vector field and the metric~\cite{Heisenberg:2014rta,Tasinato:2014eka}, and truncate its action at the quadratic order around \gr to obtain the linearized field equations. The main qualitative difference from the scalar minimal action of Sec.~\ref{sec:minimal_action} is that one can demonstrate that if any of the physical degrees of freedom has a tachyonic instability, then there are also degrees of freedom with ghost or gradient instabilities, which are manifested by change of sign of coefficients in the dispersion relationship. The conclusion is the same: vectorization based on tachyonic instabilities necessarily has ghost or gradient instabilities as well.

Recent works showed that the issues with vectorization also happen in simpler theories, an example being minimally coupled self-interacting vector field theories. For instance, for $\mathcal{A}(x)=1$ the action~\eqref{action_vt} becomes the minimally coupled Proca theory, which is well-posed. However, generalizations of the potential beyond a mass term, that is $2 \mu_X^2 x\to 4V(x)$ for some generic function $V$, results in ill-posedness~\cite{Mou:2022hqb,Clough:2022ygm,Coates:2022qia}. The underlying reason is the fact that there is an analog of the generalized Lorenz condition~\eqref{lorenz_vt} in this case as well, which ultimately leads to the same result: there is a degree of freedom in the theory whose dynamics is governed by an effective metric, an analog of $\bar{g}_{\mu\nu}$ in Eq.~\eqref{eq:geff}. This new metric depends on the vector field itself, and can lose its Lorentzian nature depending on how $X_\mu$ evolves in time~\cite{Coates:2022qia}.

\textcite{Esposito-Farese:2009wbc} provided some of the earliest examples of the ill-posedness of self-interacting vector field theories we mentioned above in a cosmological context,  and it was recently demonstrated that initially healthy configurations of self-interacting vectors naturally evolve to a point where hyperbolic evolution becomes impossible~\cite{Mou:2022hqb,Clough:2022ygm,Coates:2022qia}. This means, nonlinear extensions of the Proca theory are more delicate than their scalar counterparts in terms of providing physical theories. It seems that the underlying reason for the ``failure'' of vectorization is a much more general phenomenon in theories with conditions such as Eq.~\eqref{lorenz_vt}, at least for the models conceived thus far. In this sense, exploration of vectorization has been leading to a deeper understanding of vector fields in general, which has implications beyond theories of gravitation.

We close this section with a pertinent comment on ghost instabilities. We encountered them as an unappealing artefact in vectorization theories, but one could consider making them the driver of a spontaneous growth mechanism. This option has indeed been considered prior to their discovery in vectorization models. The simplest example is another analog of the DEF model where the conformal scaling function ${\cal A}$ depends on derivative terms~\cite{Ramazanoglu:2017yun}
\begin{align}\label{eq:action_scalar_ ghost}
S&= \frac{1}{16\pi G_{*}} \int d^4x \sqrt{-g} \left[R -
2g^{\mu\nu}\partial_{\mu}\varphi \partial_{\nu}\varphi -
2\mu^2 \varphi^2\right] \nonumber \\
&\quad + S_{{\rm m}}[\Psi_{{\rm m}}; {\cal A}^{2}(K)g_{\mu\nu}],
\end{align}
with $K=-g^{\mu\nu}\partial_\mu\varphi\partial_\nu\varphi/2$. The scalar equation of motion is
\begin{align}\label{eq:ghost_EOM}
\nabla_\mu \left(\hat{z} \nabla^\mu \varphi\right) = \mu^2 \varphi,
\end{align}
where $\hat{z} \equiv 1+4\pi G_* \tilde{T}{\cal A}^4 \left({\rm d} \ln{\cal A}/{\rm d} K\right)$. Unlike the tachyonic instability where the mass term is modified, the derivative terms are modified in this case.

The nature of Eq.~\eqref{eq:ghost_EOM} can be better understood when we consider the linearized equation for perturbative values of the scalar around a \gr background, and highlight the principal part as
\begin{align}\label{eq:ghost_EOM2}
\hat{z} \, \Box \delta\varphi + \dots = 0.
\end{align}
For an appropriate choice such as ${\cal A} = \exp(\beta_0 K)$ with large enough $\beta_0>0$, the sign of $\hat{z}$ can become negative in the presence of matter. In all spontaneous scalarization theories considered so far, the source of the instability was the change of sign of the effective mass term, which resulted in a tachyon. Here, the overall sign of the wave operator changes, hence we have both ghost and gradient instabilities, recalling Eq.~\eqref{eq:instability_eom} and the subsequent discussion. Spontaneous scalarization that might arise from this new mechanism is called \emph{ghost-based spontaneous scalarization}~\cite{Ramazanoglu:2017yun}. It has similar problems as vectorization, in that it is hard to tame the ghost and other instabilities, hence we do not end up with a theory with sensible dynamics.

\subsection{Spontaneous tensorization}
The basic mechanism of spontaneous scalarization can be generalized beyond vector fields, resulting in \emph{spontaneous tensorization}. However, all known such theories also feature ghosts and suffer from similar problems with vectorization, aside from a single possible example of spinor-tensor theory.

The simplest mathematical generalization of scalarization beyond vectors occurs for $p$-form fields, i.e., totally antisymmetric tensors $X_{\mu_1 \dots \mu_p}$ of rank $(0,p)$. This is due to the fact that a vector field $X_\mu$ is a $1$-form field, and the actions we encountered in spontaneous vectorization can be readily generalized to this class of fields as~\cite{Ramazanoglu:2019jfy}
\begin{align}\label{action_pform}
S &= \frac{1}{16\pi} \int {\rm d}^{4}x \, \sqrt{-g} \, \left[R - F_{\mu_1\dots \mu_{p+1}}F^{\mu_1\dots \mu_{p+1}} \right] \nonumber \\
&\quad + S_{{\rm m}} \left[\Psi_{{\rm m}};\, {\cal A}^{2}(x) g_{\mu \nu} \right],
\end{align}
where $F_{\mu_1\dots \mu_{p+1}} = (p+1) \nabla_{[\mu_1} X_{\mu_2\dots \mu_{p+1}]}$, with square brackets denoting antisymmetrization, and $x=X_{\mu_1\dots \mu_{p}} X^{\mu_1\dots \mu_{p}}$. It is then straightforward to show that the analog of $\mu^2_\eff$ becomes negative inside \ns{s} for an appropriate choice of ${\cal A}(x)$. However, these theories carry ghost instabilities similar to those of vectors when we try to establish spontaneous growth from tachyonic instabilities, since they possess an analog of the Lorenz condition~\eqref{lorenz_vt}~\cite{Silva:2021jya}.

Spontaneous tensorization of a symmetric rank $(0,2)$ tensor field $f_{\mu\nu}$, which is spin-$2$, is a natural avenue to investigate after spin-$0$ (scalars) and spin-$1$ (vectors). Ghosts appear in this case once more since the metric in our gravity theories, $g_{\mu\nu}$, is also a spin-$2$ field. Note that spontaneous tensorization theories include terms where the field that tensorizes, in this case $f_{\mu\nu}$, is coupled to the metric $g_{\mu\nu}$. However, interacting spin-$2$ fields are known to generically lead to ghost degrees of freedom, rendering most such theories unphysical~\cite{Boulware:1972yco,deRham:2014zqa}, and ghost-free theories being discovered more recently in the form of massive gravity and bigravity~\cite{deRham:2010kj, Hassan:2011zd, deRham:2014zqa}. Matter can couple to one of the spin-$2$ fields in the novel ghost-free theories, but coupling to both metrics generically invoke a ghost again~\cite{Yamashita:2014fga}. For example, trying to mimic the \DEF model by having a conformal scaling function that depends on $f_{\mu\nu}$, i.e.,~changing the matter action as $g_{\mu \nu} \to {\cal A}(f_{\mu\nu},g_{\mu\nu})\ g_{\mu \nu}$, leads to ghosts. Overall, there is no known form of bigravity theory that features an analog of spontaneous scalarization.

Spinor fields present the most interesting case of generalizing scalarization, perhaps aside from vectors.\footnote{Although a spinor is not a tensor in the technical sense, we will classify spinorization as an example of tensorization.}
The classical Lagrangian for a Dirac spinor is
\begin{align}\label{eq:spinor_action}
\mathcal{L}_\psi = \tfrac{1}{2} \left[ \bar{\psi} \gamma^\mu (\nabla_\mu \psi)
- \nabla_\mu \hat{\gamma}^5 \gamma^\mu \psi \right] -\mu \bar{\psi} \psi,
\end{align}
where the conventions for gamma matrices and the effect of covariant derivatives on spinors can be found in~\textcite{Ramazanoglu:2018hwk}.
This action provides the usual dispersion relation $\omega^2 = k_ik^i + \mu^2$ for a plane wave of the form $\psi \propto \exp[{\rm i} (\omega t-k_i x^i)]$. Since the mass term appears linearly, not quadratically, changing its sign results in the same dispersion relation. However, a tachyonic spinor action is still possible in the form of
\begin{align}\label{eq:spinor_action_tachyon}
\mathcal{L}_\psi^{5} = \tfrac{1}{2} \left[ \bar{\psi} \hat{\gamma}^5 \gamma^\mu (\nabla_\mu \psi)
- (\nabla_\mu \bar{\psi}) \hat{\gamma}^5 \gamma^\mu \psi \right] -\mu \bar{\psi} \psi.
\end{align}
with the field equation
\begin{align}\label{tachtonic_spin_feqn}
 \gamma^\mu \nabla_\mu \psi -   \hat{\gamma}^5\mu \psi= 0.
\end{align}
Here, $\omega^2 = k_ik^i - \mu^2$, that is, the dispersion relation is tachyonic~\cite{Chodos:1984cy,Jentschura:2011ga}.

\textcite{Ramazanoglu:2018hwk} used Eq.~\eqref{eq:spinor_action_tachyon} to obtain the first \emph{spontaneous spinorization} theory where spinor fields are unstable to growth around \gr backgrounds of \ns{s}. Even though the meaning of a ghost as opposed to a tachyon is a subtle issue for spinors, the equation of motion in this theory has divergent coefficients as in the vector field case~\eqref{eq:vec_component}, hence this form of spinorization suffers from similar problems to vectorization.

\textcite{Minamitsuji:2020hpl} proposed an alternative model of spinorization given by the action
\begin{align}\label{action_spinor2}
S &= \frac{1}{16\pi G_*} \int {\rm d}^{4}x \, \sqrt{-g} \,
\left\{ R +\tfrac{1}{2} \left[\bar{\psi} \hat{\gamma}^5 \gamma^\mu (\nabla_\mu \psi) \right. \right. \nonumber \\
&\quad \left.\left. - \, (\nabla_\mu \bar{\psi}) \hat{\gamma}^5 \gamma^\mu \psi \right]  \right\}
+ S_{\rm m} \left[\Psi_{{\rm m}},\, {\cal A}^{2}(\bar{\psi}\psi) g_{\mu \nu} \right].
\end{align}
This model has the unusual feature that the derivative part of the spinor action is not the canonical one in Eq.~\eqref{eq:spinor_action}, but rather the ``tachyonic'' one in Eq.~\eqref{eq:spinor_action_tachyon}. The resulting tachyonic equation of motion is
\begin{align}\label{dirac_eqn_spinorization2}
 \gamma^\mu \nabla_\mu \psi - \hat{\gamma}^5 ( 4\pi {\cal A}^4 \Lambda \tilde{T} ) \psi= 0,
\end{align}
where $\Lambda = {\rm d} \ln{\cal A}/{\rm d}(\bar{\psi}\psi)$. The interpretation is straightforward, we have a tachyonic field with an effective mass $\mu_{\rm eff}=4\pi {\cal A}^4 \Lambda \tilde{T}$ [cf.~Eq.~\eqref{tachtonic_spin_feqn}]. The most important aspect of this theory is the fact that the effective mass term and all coefficients of the field equation are regular everywhere, hence, the dynamics does not have any signs of ill-posedness that we have found in all generalizations of scalarization so far. In other words, action~\eqref{dirac_eqn_spinorization2} is the only known analog of spontaneous scalarization for non-scalar fields that does not suffer from ghost or gradient instabilities.

Our results on spontaneous tensorization suggest a no-go theorem for generalizing scalarization beyond scalar fields, at least if we want to avoid ghosts. It is curious to find the only exception to this trend in the relatively exotic case of spinorization, which invites studies of the deeper reasons that make spontaneous scalarization hard to replicate for other types of fields.
\section{Introduction}
\label{sec:intro}
\input{chapters/intro}

\section{Theoretical background: mechanism and theories}
\label{sec:theory_bg}
\input{chapters/theory}

\section{Neutron star scalarization}
\label{sec:ns_scalarization}
\input{chapters/NS}

\section{Black hole scalarization}
\label{sec:bh_scalarization}
\input{chapters/BH}

\section{Generalizations of scalarization to other fields and instabilities }
\label{sec:other}
\input{chapters/other}

\section{Open problems and future perspectives}
\label{sec:conclusion}
\input{chapters/conclusions}

\begin{acknowledgments}
D.D acknowledges financial support via an Emmy Noether Research Group
funded by the German Research Foundation (DFG) under grant no. DO 1771/1-1.
S.Y would like to thank the University of T\"ubingen for the financial support
and acknowledges the partial support by the Bulgarian NSF Grant KP-06-H28/7.
T.S. acknowledges partial support from the STFC Consolidated Grants No. ST/T000732/1 and No.~ST/V005596/1.
F.M.R acknowledges support from 2020 Bilim Akademisi GEB\.IP award and T\"UB\.ITAK Project No. 122F097.
H.O.S. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) - project number: 386119226.
The authors would also like to acknowledge networking support by the
COST Actions CA16104 and CA15117.
\end{acknowledgments}

\bibliographystyle{apsrmp4-2}

redpajama_set_name: RedPajamaArXiv