\section{Introduction}

Recent research on the lottery ticket hypothesis \citep[LTH,][]{frankle_2019, frankle_2019b} in deep learning has demonstrated the existence of very sparse neural networks that train to performance levels comparable to those of their dense counterparts.
These results challenge the role of overparameterization in supervised learning and provide a new perspective on the emergence of stable learning dynamics \citep{ frankle_2020b, frankle_2020a}.
Recently these results have been extended to various domains beyond supervised image classification. These include self-supervised learning \citep{chen_2020}, natural language processing \citep{yu_2019, chen_2020b} and semantic segmentation \citep{girish_2020}.
But how does the lottery ticket ticket phenomenon transfer to reinforcement learning agents?
One key challenge may be the inherent non-stationarity of the optimization problem in deep reinforcement learning (DRL): The data-generation process is not static, but depends on the changing state of the neural network.
Furthermore, a weight may serve different roles at different stages of learning (e.g.\ during exploration and exploitation). It is not obvious how a simple weight magnitude-based pruning heuristic acting on a well-performing policy shapes the learning process of the agent.
In this work, we therefore investigate winning tickets in reinforcement learning and their underlying contributing factors. We compare supervised behavioral cloning with DRL, putting a special emphasis on the resulting input representations used for prediction and control.
Thereby, we connect the statistical perspective of sparse structure discovery \citep[e.g.][]{hastie_2019} with the iterative magnitude pruning \citep[IMP,][]{han_2015} procedure in the context of Markov decision processes (MDPs). The contributions of this work are summarized as follows:

\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{figures/figure_1_new.png}%
\caption{Representation compression, disentangling baselines and lottery ticket effects in continuous control tasks.
\textbf{Left.} IMP successively prunes task-irrelevant outer rim pixels in a MNIST digit-classification task. A similar representation compression emerges for an IMP-masked agent solving a visual navigation task in DRL. The channel that encodes the patrolling enemy is pruned up to the point where only potential enemy locations are considered.
\textbf{Middle.} To disentangle the contributions of mask, initialization and layer-wise pruning ratio to the winning lottery ticket (corresponding to condition \textit{mask/weights}), we compare three baselines: After each IMP iteration, we permute either only the remaining initial weights (\textit{mask/permuted}) or also the sparsity mask (\textit{permuted/permuted}). The third baseline is created by randomly sampling a sparse mask and random re-initialization of the weights (\textit{random/re-init}). It no longer preserves the layer-wise pruning ratio of the lottery ticket.
\textbf{Right.} Average normalized performance of policies at different sparsity levels and for different baseline conditions as well as four PyBullet \citep{benelot2018} environments. For both behavioral cloning and PPO agents most of the ticket effect can be attributed to the IMP-derived mask as compared to the weight initialization. Furthermore, agents trained with explicit supervision can be pruned to higher sparsity levels before performance deteriorates.}
\label{fig1:conceptual}
\end{figure}

\begin{enumerate}
    \item We show that winning tickets exist in both high-dimensional visual and control tasks (continuous/discrete). A positive lottery ticket effect is robustly observed for both off-policy DRL algorithms, including Deep-Q-Networks \citep[DQN,][]{mnih_2015} and on-policy policy-gradient methods \citep[PPO,][]{schulman2015trust, schulman_2017}, providing evidence that the lottery ticket effect is a universal phenomenon across optimization formulations in DRL.
    \item By comparing reinforcement learning to supervised behavioral cloning (BC), we show that networks trained with explicit supervision can be pruned to higher sparsity levels before performance starts to degrade, indicating that the RL-specific distributional shift requires a larger amount of parameters in order to be addressed (section \ref{sec3a:distillrl}).

    \item By introducing a set of lottery ticket baselines (section \ref{sec3:baselines}; figure \ref{fig1:conceptual}, middle column), we disentangle the contributions of the mask, weight initialization and layer-wise pruning ratio. We demonstrate that the mask explains most of the ticket effect in behavioral cloning and reinforcement learning for MLP-based agents, whereas the associated weight initialization is less important (section \ref{sec3a:distillation} and \ref{sec3b:deeprl}). For CNN-based agents the weight initialization contributes more.

	\item By visualizing the sparsified weights for each layer, we find that early network layers are pruned more. Entire input dimensions can be rendered invisible to MLP-based agents by the pruning procedure. By this mechanism, IMP compresses the input representation of the MDP (e.g. figure \ref{fig1:conceptual}, left column, bottom row) and reveals a minimal task representation for the underlying control problems (section \ref{sec4:representations}).
	\item The IMP input layer mask not only eliminates obviously redundant dimensions, but also identifies complex relationships between input features and the task of the agent (section \ref{sec4a:repr_vision}), e.g. the proximity of an enemy or the speed of an approaching object. We show that the input masking can be transferred to train dense agents with sparse inputs at lower costs.
	\item We show that the weight initialization scheme is important for discovering minimal representations.
	Depending on the input size of different layers of the network, global magnitude-based pruning can introduce a strong layer-specific pruning bias. We compare initializations and show that a suitable initialization scheme enables the removal of task-irrelevant dimensions (section \ref{sec4b:repr_control}).
\end{enumerate}

\section{Background and Related Work}
\label{sec2:background}

\textbf{Iterative Magnitude Pruning.} We use the iterative pruning procedure outlined in \citet{frankle_2019} to identify winning tickets. We train DRL agents for a previously calibrated number of transitions and track the best performing network checkpoint throughout. Performance is measured by the average return on a set of evaluation episodes. Afterwards, we prune $20\%$ of the weights with smallest magnitude globally (across all layers). The remaining weights are reset to their initial values and we iterate this procedure (train $\to$ prune $\to$ reset).\footnote{In supervised learning, the pruning mask is often constructed based on an early stopping criterion and the final network. We instead track the best performing agent.
Thereby, we reduce noise introduced by unstable learning dynamics and exploit that the agent is trained and evaluated on the same environment.
We found that late rewinding to a later checkpoint \citep{frankle_2019b} is not necessary for obtaining tickets (SI figure \ref{figsupp:regularization}).}
The \textit{lottery ticket effect} refers to the performance gap between the sparse network obtained via IMP and a randomly initialized network with sparsity-matched random pruning mask.

\textbf{Lottery Tickets in Deep Reinforcement Learning.}
\citet{yu_2019} previously demonstrated the existence of tickets in DRL that outperform parameter-matched random initializations.
They obtained tickets for a distributed on-policy actor-critic agent on a subset of environments in the ALE benchmark \citep{bellemare_2013} as well as a set of discrete control tasks.
While they provide empirical evidence for the existence of lottery tickets in DRL, they did not investigate the underlying mechanisms. Here, we aim to unravel these mechanisms. To this end, we focus on a diverse set of environments and provide a detailed comparison between supervised behavioral cloning and on-/off-policy Deep RL with a set of carefully designed ticket baselines.
We analyze the resulting masked representations that the agent learns to act upon and the impact of specific weight initializations on the resulting sparse networks.

\textbf{Lottery Tickets with Non-Stationary Data Distributions.} \citet{desai_2019} investigated whether trained lottery tickets overfit the training data distribution under which they were obtained. Using transfer learning tasks on natural language data, they showed that lottery tickets provide general inductive biases. Similar ticket transfer results were reported by \citet{morcos_2019} and \citet{mehta_2019} in the context of optimizers and vision datasets. Unlike our work, these studies do not investigate within-training covariate shift, but instead focus on transferring ticket initializations after a full IMP run.
\citet{chen_2021}, on the other hand, investigate the ticket phenomenon in the context of lifelong learning and class-incremental image classification. They propose new pruning strategies to overcome the sequential nature of tasks and need for increased model capacity. Compared to the DRL setting, the covariate shift is here determined by the curriculum schedule of tasks and not the exploration behaviour of the network-parameterized agent.

\textbf{Deep Reinforcement Learning Background.} In our off-policy DRL experiments, we train Deep-Q-Networks \citep[DQN,][]{mnih_2015} with double Q-learning loss \citep{van_hasselt_2016} and prioritized experience replay \citep{schaul_2015}.
As a representative on-policy algorithm, we chose Proximal Policy Optimization \citep[PPO,][]{schulman2015trust, schulman_2017}.
PPO is a baseline-corrected policy gradient algorithm which uses a clipping strategy to approximate a computationally expensive trust-region optimization method.
For illustrative purposes, we train DQN agents on a visual navigation task, in which an agent has to collect coins in a grid while avoiding poison and two patrollers that are moving in restricted parts of the grid (figure \ref{fig1:conceptual}, left column, bottom row; SI \ref{SI:environments}).
We scale our results to four PyBullet \citep{benelot2018} continuous control and a subset of ALE benchmark \citep{bellemare_2013} environments. Due to computational considerations we limit each individual IMP iteration for the ATARI environments to 2.5 million frames. All other tasks were trained for a pre-calibrated generous amount of transitions.
We focus on feedforward value estimators and policies (MLP \& CNN) and used default hyperparameters with little tuning (SI \ref{SI:hyperparameters}).

\textbf{Supervised Behavioral Cloning.}
While most supervised learning relies on a stationary data distribution provided by a static dataset, reinforcement learning agents have to acquire their training data in an action-perception loop.
Since the agent's behavioural policy is learned over time, the data distribution used in optimization undergoes covariate shift. To study how the covariate shift influences winning tickets, we mimic the supervised learning case by training agents via supervised policy distillation \citep{rusu_2015, schmitt_2018}. We roll out a pre-trained expert policy and train the student agent by minimizing the KL divergence between the student's and teacher's policies.

\section{Disentangling Ticket Contributions in BC and Deep RL}
\label{sec3:baselines}

\begin{wraptable}{r}{10cm}
\begin{tabular}{ |p{2.8cm}||p{1.5cm}|p{1.5cm}| p{2cm}| }
 \hline \hline
 \multicolumn{4}{|c|}{Examined Sparsity-Generating IMP Variants} \\
 \hline
 & Retain weights & Retain \newline mask &  Retain layer \newline pruning ratio  \\
 \hline
 \textit{mask/weights}  & \cmark    & \cmark & \cmark \\
\textit{mask/permuted}  & \xmark    & \cmark & \cmark\\
\textit{permuted/permuted}    & \xmark & \xmark & \cmark \\
\textit{random/re-init} & \xmark & \xmark & \xmark\\
 \hline \hline
\end{tabular}
\caption{Proposed Baselines for Disentangling Ticket Contributions}
\label{table1:baselines}
\end{wraptable}
There are two contributing factors to the lottery ticket effect: The IMP-identified binary mask and the preserved initialized weights that remain after pruning (\textit{mask/weights}).
We aim to disentangle the contributions by introducing a set of counterfactual baselines, which modify the original IMP procedure (figure \ref{fig1:conceptual}, middle column; table \ref{table1:baselines}).
A first baseline estimates how much of the performance of the ticket can be attributed to the initial weights, by means of a layer-specific permutation of the weights that remain after masking (\textit{mask/permuted}).
A second, weaker baseline estimates the contribution of the mask, by also permuting the layer-specific masks (\textit{permuted/permuted}).
Finally, we consider the standard \textit{random/re-init} baseline, which samples random binary masks -- discarding layer-specific pruning ratios -- and re-initializes all weights at each IMP iteration.
Throughout the next sections we use these baselines to analyze and compare the factors that give rise to the lottery ticket effect in different control settings.

\vspace{-0.1cm}
\subsection{Comparing Winning Tickets in Supervised Behavioral Cloning and Deep RL}
\label{sec3a:distillrl}

Does the covariate shift in DRL affect the existence and nature of lottery tickets?
Weights pruned for their small magnitude at the end of the learning process might be needed at earlier stages, e.g., during exploration. If this were the case, the performance of the ticket in DRL should degrade for lower levels of sparsity than in a corresponding supervised task. To investigate this question, we turn to a behavioral cloning setting, in which a student agent is trained to imitate the stochastic policy of a pre-trained expert and apply magnitude pruning to the student network after each training run iteration. By collecting transitions based on the static behavioral policy of the expert, we avoid the need for exploration and the effect of an otherwise non-stationary data distribution.
For both the discrete action space and all continuous control tasks we find that agents trained with RL and the supervised settings start to degrade in performance at different sparsity levels (figure \ref{fig2:distillation_a}). More specifically, the maximal return of agents that were trained with reinforcement learning starts to drop at significantly earlier stages of the IMP procedure. Hence, the covariate shift inherent to the RL problem increases the minimal required size of a winning ticket in the tasks studied here: More parameters are required for an agent that has to solve the additional exploration problem.

\vspace{-0.1cm}
\subsection{Disentangling Ticket Contributions in Supervised Behavioral Cloning}
\label{sec3a:distillation}
To disentangle the contributions of the initial weights and the weight mask to the ticket effect in supervised behavioral cloning, we next compared the three baselines to the full non-randomized ticket (figure \ref{fig1:conceptual}, top, right and \ref{fig3:distillation_b}).
For most agents, training performance does not degrade substantially when the weights are permuted but the mask is kept intact (no significant gap between red and green curves). But we observe a strong decrease in the ability to prune the policy to high sparsity levels if one additionally permutes the IMP-derived weight mask (gap between green and blue curves).
This observation holds for the cloning of expert policies in both discrete and continuous control tasks. Hence, the ticket effect can be mainly attributed to the mask rather than the weights.
Finally, the traditional \textit{random/re-init} baseline performs worse already for moderate levels of sparsity. The resulting performance gap (blue and yellow curves) indicates the contribution of the layer-wise pruning ratio to the mask effect.
These insights emphasize the importance of strong and nuanced baselines to understand the contributions to the full lottery ticket effect.

\vspace{-0.1cm}
\subsection{Disentangling Ticket Contributions in Deep Reinforcement Learning}
\label{sec3b:deeprl}

The observation that ticket information is mostly carried by the mask rather than the initial weights carries over to the full DRL setups (figure \ref{fig1:conceptual}, bottom, right and figure \ref{fig4:learning_algos}).

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{figures/figure_2_new.png}%
\caption{Comparing lottery tickets in DRL and supervised behavioral cloning. Networks trained with explicit supervision can be pruned to higher sparsity levels before performance starts to degrade. Results are averaged over 15 runs for both the Cart-Pole and Acrobot and 10 runs for PyBullet environments. We plot mean best performance and one standard deviation.}
\label{fig2:distillation_a}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{figures/figure_3_new.png}%
\caption{Disentangling baselines for lottery tickets in supervised behavioral cloning. The gap between the ticket (\textit{mask/weights}) and weight-permuted baseline (\textit{mask/permuted}) is small, indicating a strong contribution of the mask. Results are averaged over 15 runs for both Cart-Pole and Acrobot and 10 runs for PyBullet environments. We plot mean best performance and one standard deviation.}
\label{fig3:distillation_b}
\end{figure}

When preserving the information contained in the mask, agents can be pruned to much higher levels of sparsity as compared to the case in which the mask information is distorted (permuted or resampled).
This effect holds for both PPO- and DQN-based agents and across a wide range of qualitatively different tasks.\footnote{For further experiments investigating the role of different network capacities (SI figures \ref{figsupp:network_capacity}, \ref{figsupp:network_capacity_2}), state encoding (SI figure \ref{figsupp:representations}) and regularization (SI figure \ref{figsupp:regularization}) we refer the interested reader to the SI.}
Interestingly, ticket weights are more important for CNN-based agents (e.g., see ATARI plots). We provide further experimental evidence for this observation on four MinAtar games \citep{young_2019} in the supplementary material (SI figure \ref{figsupp:learning_algos_b}), where we compare MLP- and CNN-based agents.
Unlike the behavioral cloning experiments, for some tasks the \textit{permuted/permuted} baseline performs worse than the \textit{random/reinit} configuration. This observation only occurs for the PyBullet control tasks in which both the separate critic and actor networks were pruned. If instead, one only prunes the actor network, as was done for the PPO agents trained on Cart-Pole, Acrobot and ATARI, the \textit{permuted/permuted} baseline performs stronger (see SI figure \ref{figsupp:no_critic_prune}). This result highlights the importance of differences in IMP-derived pruning ratios between agent modules, i.e. separate actor and critic networks.

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{figures/figure_4_new.png}%
\caption{Tickets in on- and off-policy deep reinforcement learning.  The disentangling baselines for tickets in on-policy (PPO) and off-policy (DQN) DRL for a set of continuous control, a visual navigation task and a subset of ATARI environments reveal the consistent importance of the IMP-extracted mask. The results are averaged over 5 independent runs on the Gridworld and ATARI environments and 10 runs on the continuous control environments. We plot mean best performance and one standard deviation.}
\label{fig4:learning_algos}
\end{figure}

\section{Minimally Task-Sufficient Representations via IMP}
\label{sec4:representations}

The previous results establish the importance of the information contained in the mask for the success of winning tickets. To better understand this phenomenon and the resulting inductive biases, we analyzed the sparsified weight matrices, specifically between the input and the units in layer one. The next section sets out to answer the following questions: \textbf{1)} The derived mask can be thought of as a pair of goggles which guide the processing of state information. What do the masked MLP-based agents observe? \textbf{2)} Can we transfer the input layer mask and re-use it as an inductive bias for otherwise dense agents? \textbf{3)} How much influence do weight initialization and input dimensionality have on the discovered input layer mask?

\subsection{Minimal Task Representations in High-Dimensional Visual Tasks}
\label{sec4a:repr_vision}

In the visual navigation experiments, the pruning primarily affects the input layer (figure \ref{fig6:representation_grid}, top right column).
When visualizing the cumulative absolute weights of the input layer for an IMP-derived DQN agent, we find that IMP deletes entire dimensions from the observation vector (figure \ref{fig6:representation_grid}, left column) by removing all connections between those dimensions and the first hidden layer.
These eliminated input dimensions are not task-relevant. A fully-connected network trained on the remaining dimensions learns the task as quickly as when trained on the complete set of input dimensions (figure \ref{fig6:representation_grid}, bottom right column).
Hence, the IMP-derived mask provides a compressed representation of a high-dimensional observation space, enabling successful training even at high levels of sparsity.

\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{figures/f6_representation_grid.png}%
\caption{IMP eliminates task-irrelevant observation dimensions for a high-dimensional visual navigation task ($o_t \in \mathbb{R}^{6 \times 10 \times 20}$).
\textbf{Left.} Channel-/pixel-wise cumulative weight magnitudes. IMP successively prunes redundant input pixels which are not necessary to solve the navigation task. All of the pruned enemy channel pixels encode locations which the patrolling enemy cannot access.
\textbf{Right, Top.} Pruning affects the layers and object channels differentially. The input layer is pruned most strongly, while the agent channel is pruned the least.
\textbf{Right, Bottom.} The input layer pruning mask (at moderate sparsity levels) can be used as an inductive bias for training a dense network.}
\label{fig6:representation_grid}
\end{figure}

This result extends to masked input layers of MLP-agents trained on the MinAtar environments (figure \ref{fig7:representation_minatar} and SI).
For the Freeway task we find that the IMP mask encodes the notion of velocity of moving cars. The channel encoding the slowest moving car is pruned the most, maintaining only the pixels needed to react in time.
The same principle applies to the enemy bullet channel in the SpaceInvaders task. The agent only needs to know about a potential collision at the next time step in order to avoid being hit. Therefore, IMP prunes all information about the enemy bullet that is more than one step away from the agent.
For all visualized MinAtar tasks, only the actionable row/column of the agent channel is preserved. In Breakout, additionally the bottom row of the ball and trail channel is pruned since the game terminates once the ball or trail reach it. There is no action-relevant information encoded and it can be discarded.
In summary, we find that IMP-derived input layer mask provides a visually interpretable and physically meaningful compression of the observation space.

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{figures/f7_representation_minatar.png}%
\caption{IMP eliminates task-irrelevant observation dimensions for a selection of MinAtar environments \citep{young_2019}. The environment depicting figures were adapted from \citet{young_2019}.
\textbf{Left.} Freeway. IMP provides an inductive bias by differentially pruning object channels with different velocities (e.g.\ car speeds).
\textbf{Middle.} SpaceInvaders. IMP only preserves pixels of enemy objects which encode actionable proximity information (e.g. bullets being close to the agent).
\textbf{Right.} Breakout. IMP discards pixels which only change when it is too late to act (e.g.\ when the ball has left the display and the episode terminates).}
\label{fig7:representation_minatar}
\end{figure}

\subsection{Minimal Representations in Low-Dimensional Control Tasks and the Issue of Initialization Scales}
\label{sec4b:repr_control}

The continuous and discrete control tasks, on the other hand, rely on low-dimensional state representations. At first sight, they do not contain obviously uninformative input dimensions that lend themselves to be pruned.
We would expect the cumulative weight strength of each input dimension to decrease at equal speed over the course of iterative pruning.
Contrary to this initial intuition, IMP still aggressively prunes core dimensions while yielding trainable agents (figure \ref{fig7:representation_pybullet}, \ref{fig8:representation_cartpole}).

For a subset of continuous control tasks, we find that one can prune 20 and up to 50 percent of the input dimensions, while still being able to train the agent to the performance level of a dense counterpart (figure \ref{fig7:representation_pybullet}). The exact amount depends on the considered environment indicating a varying degree of observation over-specification across task formulations. For example, the Walker2D environment can train to full performance with only 11 out of 22 observation dimensions, while the Ant environment requires 22 out of 28 dimensions. We note an initial positive effect of the weight pruning on the performance of all agents indicating the effectiveness of the implied regularization.

\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{figures/figure_7_new.png}%
\caption{IMP eliminates task-irrelevant observation dimensions for a set of continuous control tasks. For all three tasks there exists a moderate sparsity level (approximately 10\% non-sparse weights) at which entire input units (columns of the first linear layer) are pruned while the agents still train to the performance level of their dense counterparts. The extent of this observation strongly depends on the considered environment indicating a varying extent of observation over-specification.}
\label{fig7:representation_pybullet}
\end{figure}

In the Cart-Pole task the cart position and its velocity are pruned (figure \ref{fig8:representation_cartpole}, top). Newton's laws are invariant with regard to the choice of the inertial frame of reference, so the input channels encoding the cart position and velocity do not provide additional information over the pole's angle and velocity.
The Acrobot task requires the agent to swing-up a two-link pendulum across a line as fast as possible. Here, IMP eliminates all sine and cosine transformations of the rotational angles of the links (figure \ref{fig8:representation_cartpole}, middle row, right column). Since the swing-up can be achieved by coordinating and increasing the angular velocities of the two links, only these two dimensions have to be preserved from pruning.
These experiments demonstrate that IMP can yield fundamental insights into what physical information is sufficient to solve a task and leverages them to prune the input representation without impairing performance.

During our analysis we discovered that the core dimension discovery of IMP crucially depends on the initialization of the input layer weights. Popular weight initialization heuristics such as the Kaiming family \citep{he_2015} aim to preserve activation magnitudes across network layers. As a result, layers that receive high-dimensional input are initialized at lower values and are therefore prone to more aggressive pruning by IMP. This bias does not occur for Xavier initializations \citep[figure \ref{fig8:representation_cartpole}, top rows, right column]{glorot_2010}. Furthermore, we experimented with different approaches to combat this initialization bias. We find that downscaling the Kaiming initialization of the input layer by a factor of ten leads to a robust discovery after few IMP iterations (figure \ref{fig8:representation_cartpole}, top rows, middle column). While smaller rescaling factors lead to a slower separation of dimensions, larger factors result in overly aggressive pruning of the input layer. In summary, we have shown that the input layer initialization can shape the resulting lottery ticket mask and provide a first initialization heuristic promoting efficient minimal task representation discovery.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{figures/f8_representation_cartpole.png}%
\caption{IMP eliminates task-irrelevant observation dimensions for low-dimensional control tasks. Suitable initialization heuristics can promote and accelerate the discovery of minimal representations.
\textbf{Top.} IMP identifies that the Cart-Pole task can be solved with only two dimensions: The pole angle and angular velocity.
\textbf{Middle.} IMP identifies that the Acrobot task can be solved with only two dimensions: The angular velocities of the two pendulum links.
\textbf{Bottom.} Agents can successfully be trained on the subset of IMP-derived dimensions alone. We plot mean performance and one standard deviation across 10 independent runs.}
\label{fig8:representation_cartpole}
\end{figure}

\section{Conclusion}
\label{section5:conclusion}

\textbf{Summary.} This work has investigated the mechanisms underlying the lottery ticket effect in behavioral cloning and DRL.
We provide evidence that agents which are trained on a stationary data distribution and do not face the exploration problem, can be pruned to higher levels of sparsity while successfully performing their task.
For both setups we find that the effect can be attributed to the identified mask rather than to the corresponding weight initialization.
We have shown that the mask compresses the input representation by removing task-irrelevant information in the form of visual object pixels or entire control observation dimensions.
Finally, we found that the pruning process is sensitive to the layer-specific scale of the weight initialization. The interpretability of the resulting minimal representations can be enhanced by heuristically re-scaling the weights in different layers.

\textbf{Limitations.} While our observations are qualitatively robust, they quantitatively vary across learning paradigms, algorithms, architectures and different tasks. Compared to our results, weight values have been claimed to play a more central role in the context of lottery tickets in supervised image classification \citep{frankle_2020b}. Future work will need to investigate under which conditions the initialization of the weight values significantly contributes to the ticket effect in RL.
Furthermore, whether the interpretation of the mask effect as a representation regularizer also holds for intermediate and higher layer masks is hard to assess and remains an open question. Finally, this study is empirical in its nature and will require further theoretical guidance and foundations.

\textbf{Future Work.} In future work we want to investigate layer-wise effects, e.g., by studying layer-specific pruning ratios, normalization schemes and fine-grained baseline analyses, in the hope of identifying sparser, more efficient and interpretable winning tickets.
Furthermore, our proposed baseline analysis can be applied more generally to all studies considering the lottery ticket procedure.
Finally, we believe that there are many opportunities to modify the original IMP procedure to simplify the discovery of minimal representations. For example, pruning entire units of intermediate layers may improve the interpretability of hidden representations.

\section*{Acknowledgements}

We thank Jonathan Frankle for initial discussions and feedback on the first manuscript draft. Furthermore, we thank Joram Keijser for reviewing the manuscript. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC 2002/1 “Science of Intelligence” - project number 390523135.

\newpage


\section{Introduction}

A \emph{$k$-coloring} of the vertices of a graph $G=(V,E)$ is a map $\phi:V \rightarrow\{1,2,\dots,k\}$. A $k$-coloring $\phi$ is a \emph{proper coloring}, if and only if, for all edge $xy\in E,\phi(x)\neq\phi(y)$. In other words, no two adjacent vertices share the same color. The \emph{chromatic number} of $G$, denoted by $\chi(G)$, is the smallest integer $k$ such that $G$ has a proper $k$-coloring.  A generalization of $k$-coloring is $k$-list-coloring.
A graph $G$ is {\em $L$-list colorable} if for a
given list assignment $L=\{L(v): v\in V(G)\}$ there is a proper
coloring $\phi$ of $G$ such that for all $v \in V(G), \phi(v)\in
L(v)$. If $G$ is $L$-list colorable for every list assignment $L$ with $|L(v)|\ge k$ for all $v\in V(G)$, then $G$ is said to be {\em $k$-choosable} or \emph{$k$-list-colorable}. The \emph{list chromatic number} of a graph $G$ is the smallest integer $k$ such that $G$ is $k$-choosable.

In 1969, Kramer and Kramer introduced the notion of 2-distance coloring \cite{kramer2,kramer1}. This notion generalizes the ``proper'' constraint (that does not allow two adjacent vertices to have the same color) in the following way: a \emph{$2$-distance $k$-coloring} is such that no pair of vertices at distance at most 2 have the same color (similarly to proper $k$-list-coloring, one can also define \emph{$2$-distance $k$-list-coloring}). The \emph{$2$-distance chromatic number} of $G$, denoted by $\chi^2(G)$, is the smallest integer $k$ so that $G$ has a 2-distance $k$-coloring. We denote $\chi^2_l(G)$ the $2$-distance list chromatic number of $G$.

For all $v\in V$, we denote $d_G(v)$ the degree of $v$ in $G$ and by $\Delta(G) = \max_{v\in V}d_G(v)$ the maximum degree of a graph $G$. For brevity, when it is clear from the context, we will use $\Delta$ (resp. $d(v)$) instead of $\Delta(G)$ (resp. $d_G(v)$).
One can observe that, for any graph $G$, $\Delta+1\leq\chi^2(G)\leq \Delta^2+1$. The lower bound is trivial since, in a 2-distance coloring, every neighbor of a vertex $v$ with degree $\Delta$, and $v$ itself must have a different color. As for the upper bound, a greedy algorithm shows that $\chi^2(G)\leq \Delta^2+1$. Moreover, this bound is tight for some graphs (see \cite{lmpv19} for some examples).

By nature, $2$-distance colorings and the $2$-distance chromatic number of a graph depend a lot on the number of vertices in the neighborhood of every vertex. More precisely, the ``sparser'' a graph is, the lower its $2$-distance chromatic number will be. One way to quantify the sparsity of a graph is through its maximum average degree. The \emph{average degree} $\ad$ of a graph $G=(V,E)$ is defined by $\ad(G)=\frac{2|E|}{|V|}$. The \emph{maximum average degree} $\mad(G)$ is the maximum, over all subgraphs $H$ of $G$, of $\ad(H)$. Another way to measure the sparsity is through the girth, i.e. the length of a shortest cycle. We denote $g(G)$ the girth of $G$. Intuitively, the higher the girth of a graph is, the sparser it gets. These two measures can actually be linked directly in the case of planar graphs.

\begin{proposition}[Folklore]\label{maximum average degree and girth proposition}
For every planar graph $G$, $(\mad(G)-2)(g(G)-2)<4$.
\end{proposition}

A graph is \emph{planar} if one can draw its vertices with points on the plane, and edges with curves intersecting only at its endpoints. When $G$ is a planar graph, Wegner conjectured in 1977 that  $\chi^2(G)$ becomes linear in $\Delta(G)$:

\begin{conjecture}[Wegner \cite{wegner}]
\label{conj:Wegner}
Let $G$ be a planar graph with maximum degree $\Delta$. Then,
$$
\chi^2(G) \leq \left\{
    \begin{array}{ll}
        7, & \mbox{if } \Delta\leq 3, \\
        \Delta + 5, & \mbox{if } 4\leq \Delta\leq 7,\\
        \left\lfloor\frac{3\Delta}{2}\right\rfloor + 1, & \mbox{if } \Delta\geq 8.
    \end{array}
\right.
$$
\end{conjecture}

The conjecture was proven for some cases (\cite{wegner}, \cite{tho18}, \cite{har16}, \cite{havet}) and some subfamilies of planar graphs \cite{lwz03}.

Wegner's conjecture motivated extensive researches on $2$-distance chromatic number of sparse graphs, either of planar graphs with high girth or of graphs with upper bounded maximum average degree which are directly linked due to \Cref{maximum average degree and girth proposition}. For a survey of the work done on $2$-distance coloring of planar graphs with high girth, see \cite{lmpv19}.

In this article, we prove the following theorems:

\begin{theorem} \label{main theorem1}
If $G$ is a graph with $\mad(G)<\frac83$ and $\Delta(G)\geq 4$, then $\chi^2_l(G)\leq \Delta(G)+3$.
\end{theorem}

\begin{theorem} \label{main theorem2}
If $G$ is a graph with $\mad(G)<\frac{14}5$ and $\Delta(G)\geq 6$, then $\chi^2_l(G)\leq \Delta(G)+3$.
\end{theorem}

\Cref{main theorem1} improves upon a previous result for graphs with $\Delta(G)=5$ \cite{bu15}.

Due to \Cref{maximum average degree and girth proposition}, we get the following corollaries for planar graphs:

\begin{corollary}
If $G$ is a graph with $g(G)\geq 8$ and $\Delta(G)\geq 4$, then $\chi^2_l(G)\leq \Delta(G)+3$.
\end{corollary}

\begin{corollary}
If $G$ is a graph with $g(G)\geq 7$ and $\Delta(G)\geq 6$, then $\chi^2_l(G)\leq \Delta(G)+3$.
\end{corollary}

\paragraph{Notations and drawing conventions.} For $v\in V(G)$, the \emph{2-distance neighborhood} of $v$, denoted $N^*_G(v)$, is the set of 2-distance neighbors of $v$, which are vertices at distance at most two from $v$, not including $v$. We also denote $d^*_G(v)=|N^*_G(v)|$. We will drop the subscript and the argument when it is clear from the context. Also for conciseness, from now on, when we say ``to color'' a vertex, it means to color such vertex differently from all of its colored neighbors at distance at most two. Similarly, any considered coloring will be a 2-distance coloring. We will also say that a vertex $u$ ``sees'' another vertex $v$ if $u$ and $v$ are at distance at most 2 from each other.

Some more notations:
\begin{itemize}
\item A \emph{$d$-vertex} ($d^+$-vertex, $d^-$-vertex) is a vertex of degree $d$ (at least $d$, at most $d$).
\item A \emph{$k$-path} ($k^+$-path, $k^-$-path) is a path of length $k+1$ (at least $k+1$, at most $k+1$) where the $k$ internal vertices are 2-vertices.
\item A \emph{$(k_1,k_2,\dots,k_d)$-vertex} is a $d$-vertex incident to $d$ different paths, where the $i^{\rm th}$ path is a $k_i$-path for all $1\leq i\leq d$.
\end{itemize}

In both proofs, we will consider $G_1$ (resp. $G_2$) a counter-example to \Cref{main theorem1} (resp. \Cref{main theorem2}) with the fewest number of vertices. The purpose of the proofs is to prove that $G_1$ (resp. $G_2$) cannot exist. We will always start by studying their structural properties, then apply a discharging procedure.

\section{Proof of \Cref{main theorem1}}
\label{sec2}

Since $\mad(G_1)<\frac83$, we have
\begin{equation}\label{equation1}
\sum_{u\in V(G)}(3d(u)-8) < 0
\end{equation}

We assign to each vertex $u$ the charge $\mu(u)=3d(u)-8$. To prove the non-existence of $G_1$, we will redistribute the charges preserving their sum and obtaining a non-negative total charge, which will contradict \Cref{equation1}. Let us study the structural properties of $G_1$ first.

\subsection{Structural properties of $G_1$\label{tutu1}}

\begin{lemma}\label{connected1}
Graph $G_1$ is connected.
\end{lemma}
Otherwise a component of $G_1$ would be a smaller counterexample.

\begin{lemma}\label{minimumDegree1}
The minimum degree of $G_1$ is at least 2.
\end{lemma}

\begin{proof}
By \Cref{connected1}, the minimum degree is at least 1 or $G_1$ would be a single isolated vertex which is $(\Delta(G_1)+3)$-colorable. If $G_1$ contains a degree 1 vertex $v$, then we can simply remove such vertex and 2-distance color the resulting graph, which is possible by minimality of $G_1$. Then, we add $v$ back and color $v$ (at most $\Delta(G_1)$ constraints and $\Delta(G_1)+3$ colors).
\end{proof}

\begin{lemma}\label{2-path1}
Graph $G_1$ has no $2^+$-paths.
\end{lemma}

\begin{proof}
Suppose that $G_1$ contains a $2^+$-path $uvwx$. We color $G_1-\{v,w\}$ by minimality of $G_1$. Observe that $v$ and $w$ each sees at most $\Delta+1$ colors so they have at least two available colors left each. Thus, we can easily extend the coloring to $v$ and $w$.
\end{proof}

\begin{lemma}\label{(111)1}
Graph $G_1$ has no $(1,1,1)$-vertex.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,1,1)$-vertex $u$ with three $2$-neighbors $u_1$, $u_2$, and $u_3$. We color $G_1-\{u,u_1,u_2,u_3\}$ by minimality of $G_1$, then we extend this coloring to the remaining vertices by coloring $u_1$, $u_2$, $u_3$, and $u$ in this order. Observe that this possible since we have $\Delta+3\geq 7$ colors as $\Delta\geq 4$.
\end{proof}

\begin{lemma}\label{110-001}
Graph $G_1$ has no $3$-vertex with a $2$-neighbor and a $(1,1,0)$-neighbor.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a 3-vertex $u$ with a $2$-neighbor $v$ and a $(1,1,0)$-neighbor $w$. Let the $2$-neighbors of $w$ be $w_1$ and $w_2$.

First, observe that if two adjacent $3$-vertices share a common $2$-neighbor, for example, if $u$ is also adjacent to $w_1$, then we color $G_1-\{u,w,w_1\}$ by minimality of $G_1$ and finish by coloring $u$, $w$, and $w_1$ in this order. This is possible since we have $\Delta+3$ colors and $\Delta\geq 4$. Hence, all named vertices are distinct.

Now, we color $G_1-\{u,v,w,w_1,w_2\}$ by minimality. Let $L(x)$ be the list of available colors left for a vertex $x\in\{u,v,w,w_1,w_2\}$. Since we have $\Delta+3$ colors and $\Delta\geq 4$, $|L(v)|\geq 2$, $|L(u)|\geq 2$, $|L(w)|\geq 4$, $|L(w_1)|\geq 3$, and $|L(w_2)|\geq 3$. We remove the extra colors so that $|L(x)|$ reaches the lower bound for each $x\in\{u,v,w,w_1,w_2\}$. Consider the two following cases.
\begin{itemize}
\item If $L(u)\neq L(v)$, then we color $u$ with $c\in L(u)\setminus L(v)$. We finish by coloring $w_1$, $w_2$, $w$, and $v$ in this order.
\item If $L(u)=L(v)$, we color $w_1$ with $c\in L(w_1)\setminus L(u)$ (which is possible since $|L(w_1)|=3$ and $|L(u)|=2$). Then, we color $w$ with $d\in L(w)\setminus (L(u)\cup\{c\})$ (which is possible as $|L(w)|=4$). Finally, we finish by coloring $w_2$, $u$, and $v$ in this order.
\end{itemize}
We thus obtain a valid coloring of $G_1$, which is a contradiction.
\end{proof}

\begin{lemma}\label{110-00(-3)0-011}
Graph $G_1$ has no $3$-vertex with two $(1,1,0)$-neighbors and another $3$-neighbor.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $3$-vertex $u$ with two $(1,1,0)$-neighbors $v$ and $w$ and another $3$-neighbor $t$. Let $v_1$ and $v_2$ (resp. $w_1$ and $w_2$) be $v$'s (resp. $w$'s) $2$-neighbors.

If $v$ and $w$ share a common $2$-neighbor, say $v_1=w_1$, then we color $G_1-\{u,v,w,v_1,v_2,w_2\}$ by minimality of $G_1$ and finish by coloring $u$, $v$, $w$, $v_2$, $w_2$, and $v_1$ in this order. This is possible since we have $\Delta+3$ colors and $\Delta\geq 4$. Note that this coloring also work when $v_2=w_2$. Hence, all named vertices are distinct.

Now, we color $G_1-\{u,v,w,w_1,w_2\}$ by minimality. Let $L(x)$ be the list of available colors left for a vertex $x\in\{u,v,w,w_1,w_2\}$. Since we have $\Delta+3$ colors and $\Delta\geq 4$, $|L(u)|\geq 2$ (as $d(t)=3$), $|L(v)|\geq 2$, $|L(w)|\geq 4$, $|L(w_1)|\geq 3$, and $|L(w_2)|\geq 3$. Note that we obtain the same lower bounds on the lists of colors as in \Cref{110-001}. Thus, the exact same proof holds and we have a valid coloring of $G_1$, which is a contradiction.
\end{proof}

\subsection{Discharging rules \label{tonton1}}

In this section, we will define a discharging procedure that contradicts the structural properties of $G_1$ (see Lemmas \ref{connected1} to \ref{110-00(-3)0-011}) showing that $G_1$ does not exist. We assign to each vertex $u$ the charge $\mu(u)=3d(u)-8$. By \Cref{equation1}, the total sum of the charges is negative. We then apply the following discharging rules.

\medskip

\begin{itemize}
\item[\ru0] Every $3^+$-vertex gives 1 to each of its $2$-neighbors.
\item[\ru1] Every $4^+$-vertex gives 1 to each of its $3$-neighbors.
\item[\ru2] Every $(0,0,0)$-vertex gives 1 to each of its $(1,1,0)$-neighbors.
\end{itemize}

\subsection{Verifying that charges on each vertex are non-negative} \label{verification}

Let $\mu^*$ be the assigned charges after the discharging procedure. In what follows, we will prove that: $$\forall u \in V(G_1), \mu^*(u)\ge 0.$$

Let $u\in V(G_1)$.

\textbf{Case 1:} If $d(u) = 2$, then $u$ receives charge 1 from each endvertex of the 1-path it lies on by \ru0 (as there are no $2^+$-path by \Cref{2-path1}); Thus we get $\mu^*(u) = \mu(u) + 2\cdot 1 = 3\cdot 2 - 8 + 2 = 0$.

\textbf{Case 2:} If $d(u)=3$, then $\mu(u)=3\cdot 3 - 8 = 1$. Since there are no $2^+$-paths due to \Cref{2-path1} and no $(1,1,1)$-vertices due to \Cref{(111)1}, we have the following cases.
\begin{itemize}
\item If $u$ is a $(1,1,0)$-vertex, then $u$ gives 1 to each of its two $2$-neighbors by \ru0. At the same time, $u$ also receives 1 from its $3^+$-neighbor $v$ by \ru1 or \ru2 as $v$ is either a $4^+$-vertex or a $(0,0,0)$-vertex by \Cref{110-001}. To sum up,
$$ \mu^*(u)\geq 1-2\cdot 1 + 1 = 0.$$
\item If $u$ is a $(1,0,0)$-vertex, then $u$ only gives 1 to its $2$-neighbor by \ru0. Hence,
$$ \mu^*(u)\geq 1-1 = 0.$$
\item If $u$ is a $(0,0,0)$-vertex, then $u$ only gives charge to $(1,1,0)$-vertices by \ru2. Let $t$, $v$, and $w$ be $u$'s $3^+$-neighbors.

If $u$ is adjacent to a $4^+$-neighbor then it receives 1 by \ru1 and at worst, it gives 1 to each of the two other neighbors by \ru2. As a result,
$$ \mu^*(u)\geq 1+1-2\cdot 1 = 0. $$

If $u$ is adjacent to three $3$-vertices, then at most one of them can be a $(1,1,0)$-vertex due to \Cref{110-00(-3)0-011}. So, $u$ only gives at most 1 to a $(1,1,0)$-neighbor by \ru2. Consequently,
$$ \mu^*(u)\geq 1-1 = 0. $$
\end{itemize}

\textbf{Case 3:} If $4\leq d(u)\leq \Delta$, then, at worst, $u$ gives 1 to each of its neighbors by \ru0 and \ru1. As a result,
$$ \mu^*(u)\geq 3d(u)-8-d(u) \geq 2\cdot 4 - 8 = 0.$$

To conclude, we started with a charge assignment with a negative total sum, but after the discharging procedure, which preserved this sum, we end up with a non-negative one, which is a contradiction. In other words, there exists no counter-example $G_1$ to \Cref{main theorem1}.

\section{Proof of \Cref{main theorem2}}
\label{sec3}

Since $\mad(G_2)<\frac{14}5$, we have
\begin{equation}\label{equation2}
\sum_{u\in V(G)}(5d(u)-14) < 0
\end{equation}

We assign to each vertex $u$ the charge $\mu(u)=5d(u)-14$. To prove the non-existence of $G_2$, we will redistribute the charges preserving their sum and obtaining a non-negative total charge, which will contradict \Cref{equation2}. Let us study the structural properties of $G_2$ first.

\subsection{Structural properties of $G_2$\label{tutu2}}

Observe that in Lemmas \ref{connected1} to \ref{(111)1}, we only use the fact that $G_1$ is a minimal counter-example to \Cref{main theorem1} and that we have $\Delta(G_1)+3$ colors, which is at least 7. Consequently, since $G_2$ is a minimal counter-example to \Cref{main theorem2} and we have $\Delta(G_2)+3\geq 9$ colors as $\Delta(G_2)\geq 6$, the same arguments allow us to obtain Lemmas \ref{connected2} to \ref{(111)2}:

\begin{lemma}\label{connected2}
Graph $G_2$ is connected.
\end{lemma}

\begin{lemma}\label{minimumDegree2}
The minimum degree of $G_2$ is at least 2.
\end{lemma}

\begin{lemma}\label{2-path2}
Graph $G_2$ has no $2^+$-paths.
\end{lemma}

\begin{lemma}\label{(111)2}
Graph $G_2$ has no $(1,1,1)$-vertex.
\end{lemma}

\begin{lemma}\label{011-5}
A $(1,1,0)$-vertex can only share a common $2$-neighbor with a $\Delta$-vertex.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,1,0)$-vertex $u$, with two $2$-neighbors $v$ and $w$, and let $x$ be the other endvertex of the $1$-path $uwx$ with $d(x)\leq \Delta-1$. We color $G_2-\{u,v,w\}$ by minimality of $G_2$. Then, it suffices to finish coloring $u$, $v$, and $w$ in this order as we have $\Delta+3$ colors.
\end{proof}

\begin{lemma}\label{110-4}
A $(1,1,0)$-vertex has a $(\Delta-1)^+$-neighbor.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,1,0)$-vertex $u$ with two $2$-neighbors $v$ and $w$ and a $(\Delta-2)^-$-neighbor $t$. We color $G_2-\{u,v,w\}$ by minimality of $G_2$. Then, it suffices to finish coloring $v$, $w$, and $u$ in this order as we have $\Delta+3$ colors.
\end{proof}

\begin{lemma}\label{3-010-3}
A $(1,0,0)$-vertex cannot have two $3$-neighbors.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,0,0)$-vertex $u$ with a $2$-neighbor $v$ and two $3$-neighbors. We color $G_2-\{v\}$ by minimality of $G_2$. We uncolor $u$ then it suffices to finish coloring $v$ and $u$ in this order as we have $\Delta+3$ colors and $\Delta\geq 6$.
\end{proof}

\begin{lemma}\label{4-01(-5)0-3}
A $(1,0,0)$-vertex with a $3$-neighbor and a $4$-neighbor can only share a common $2$-neighbor with a $\Delta$-vertex.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,0,0)$-vertex $u$, with a $2$-neighbor $v$, a $3$-neighbor, and a $4$-neighbor, and let $x$ be the other endvertex of the $1$-path $uvx$ with $d(x)\leq \Delta-1$. We color $G_2-\{v\}$ by minimality of $G_2$. We uncolor $u$ then it suffices to finish coloring $u$ and $v$ in this order as we have $\Delta+3$ colors and $\Delta\geq 6$.
\end{proof}

\begin{lemma}\label{1111-5}
A $(1,1,1,1)$-vertex can only share a common $2$-neighbor with a $\Delta$-vertex.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,1,1,1)$-vertex $u$, with four $2$-neighbors $v_1$, $v_2$, $v_3$, and $v_4$, and let $x$ be the other endvertex of the $1$-path $uv_4x$ with $d(x)\leq \Delta-1$. We color $G_2-\{u,v_1,v_2,v_3,v_4\}$ by minimality of $G_2$. Then, it suffices to finish coloring $v_1$, $v_2$, $v_3$, $v_4$, and $u$ in this order as we have $\Delta+3$ colors and $\Delta\geq 6$.
\end{proof}

\begin{lemma}\label{4-1110-3}
A $(1,1,1,0)$-vertex with a $3$-neighbor can only share a common $2$-neighbor with a $(\Delta-1)^+$-vertex.
\end{lemma}

\begin{proof}
Suppose by contradiction that there exists a $(1,1,1,0)$-vertex $u$, with three $2$-neighbors $v_1$, $v_2$, and $v_3$ and a $3$-neighbor, and let $x$ be the other endvertex of the $1$-path $uv_3x$ with $d(x)\leq \Delta-2$. We color $G_2-\{u,v_1,v_2,v_3\}$ by minimality of $G_2$. Then, it suffices to finish coloring $v_1$, $v_2$, $u$, and $v_3$ in this order as we have $\Delta+3$ colors and $\Delta\geq 6$.
\end{proof}

\subsection{Discharging rules \label{tonton1}}

In this section, we will define a discharging procedure that contradicts the structural properties of $G_2$ (see Lemmas \ref{connected2} to \ref{4-1110-3}) showing that $G_2$ does not exist. We assign to each vertex $u$ the charge $\mu(u)=5d(u)-14$. By \Cref{equation2}, the total sum of the charges is negative. We then apply the following discharging rules:

\medskip

\begin{itemize}
\item[\ru0] Every $3^+$-vertex gives 2 to each of its $2$-neighbors.
\item[\ru1]
\begin{itemize}
\item[(i)] Every $4$-vertex gives $\frac12$ to each of its $3$-neighbors.
\item[(ii)] Every $5^+$-vertex gives 2 to each of its $3$-neighbors.
\end{itemize}
\item[\ru2] Let $uvw$ be a $1$-path:
\begin{itemize}
\item[(i)] If $d(u)=5$ and $d(w)\leq 4$, then $u$ gives $\frac15$ to $w$.
\item[(ii)] If $d(u)\geq 6$ and $d(w)\leq 4$, then $u$ gives $\frac23$ to $w$.
\end{itemize}
\end{itemize}

\subsection{Verifying that charges on each vertex are non-negative} \label{verification}

Let $\mu^*$ be the assigned charges after the discharging procedure. In what follows, we will prove that: $$\forall u \in V(G_2), \mu^*(u)\ge 0.$$

Let $u\in V(G_2)$.

\textbf{Case 1:} If $d(u) = 2$, then $u$ receives charge 2 from each endvertex of the path it lies on by \ru0 (as there are no $2^+$-path by \Cref{2-path2}); Thus we get $\mu^*(u) = \mu(u) + 2\cdot 2 = 5\cdot 2 - 14 + 2\cdot 2 = 0$.

\textbf{Case 2:} If $d(u)=3$, then recall $\mu(u) = 5\cdot 3 - 14 = 1$. Since there are no $2^+$-paths due to \Cref{2-path2} and no $(1,1,1)$-vertices due to \Cref{(111)2}, we have the following cases.
\begin{itemize}
\item If $u$ is a $(1,1,0)$-vertex, then $u$ gives 2 to each of its two $2$-neighbors. At the same time, the other endvertices of $u$'s incident 1-paths must be $\Delta$-vertices due to \Cref{011-5}. So, $u$ receives $\frac23$ from each of these $\Delta$-endvertices by \ru2(ii) as $\Delta\geq 6$. Moreover, $u$'s $3^+$-neighbor must be a $(\Delta-1)^+$-vertex due to \Cref{110-4}. As a result, $u$ also receives $2$ from its $(\Delta-1)^+$-neighbor as $\Delta\geq 6$. To sum up,
$$ \mu^*(u)\geq 1 - 2\cdot 2 + 2\cdot \frac23 + 2 = \frac13. $$

\item If $u$ is a $(1,0,0)$-vertex, then $u$ only gives 2 to its $2$-neighbor by \ru0. We distinguish the two following cases.
\begin{itemize}
\item If $u$ has a $5^+$-neighbor, then $u$ receives 2 by \ru1(ii). Thus,
$$ \mu^*(u)\geq 1 - 2 + 2 = 1.$$
\item If $u$ only has $4^-$-neighbors, then $u$ must have at least one $4$-neighbor due to \Cref{3-010-3}.

Now, if $u$ has two $4$-neighbors, then it receives $\frac12$ twice by \ru1(i). Hence,
$$ \mu^*(u)\geq 1- 2 + 2\cdot\frac12 = 0.$$

If $u$ has exactly one $4$-neighbor and the other one is a $3$-neighbor, then the other endvertex of the 1-path incident to $u$ must be a $\Delta$-vertex due to \Cref{4-01(-5)0-3}. So, $u$ receives $\frac12$ from its $4$-neighbor by \ru1(i) and $\frac23$ from the $\Delta$-endvertex by \ru2(ii) as $\Delta\geq 6$. To sum up,
$$ \mu^*(u)\geq 1 - 2 + \frac12 + \frac23 = \frac16. $$
\end{itemize}

\item If $u$ is a $(0,0,0)$-vertex, then $u$ does not gives any charge away. Thus, $$ \mu^*(u)=\mu(u)=1.$$
\end{itemize}

\textbf{Case 3:} If $d(u)=4$, then recall $\mu(u)=5\cdot 4 - 14 = 6$ and observe that $u$ only gives charge 2 or $\frac12$ away respectively by \ru0 or \ru1(i). We have the following cases.
\begin{itemize}
\item If $u$ is a $(1,1,1,1)$-vertex, then $u$ gives $2$ to each of its four $2$-neighbors by \ru0. At the same time, the other endvertices of the 1-paths incident to $u$ are all $\Delta$-vertices due to \Cref{1111-5}. As a result, $u$ also receives $\frac23$ from each of the four $\Delta$-endvertices by \ru2(ii). To sum up,
$$ \mu^*(u)\geq 6 - 4\cdot 2 + 4\cdot \frac23 = \frac23.$$

\item If $u$ is a $(1,1,1,0)$-vertex, then $u$ gives $2$ to each of its three $2$-neighbors by \ru0. Let $v$ be the $3^+$-neighbor.

If $v$ is a $4^+$-vertex, then $u$ does not give anything to $v$. Thus,
$$ \mu^*(u)\geq 6-3\cdot 2 = 0.$$

If $v$ is a $3$-vertex, then $u$ gives $\frac12$ to $v$ by \ru1(i). Due to \Cref{4-1110-3}, the other endvertices of the 1-paths incident to $v$ must be $(\Delta-1)^+$-vertices. As a result, $v$ receives at least $\frac15$ from each of the three $(\Delta-1)^+$-endvertices as $\Delta\geq 6$ by \ru2. To sum up,
$$ \mu^*(u)\geq 6 - 3\cdot 2 - \frac12 + 3\cdot\frac15 = \frac1{10}.$$

\item If $u$ is a $(1^-,1^-,0,0)$-vertex, then at worst $u$ gives 2 twice by \ru0 and $\frac12$ twice by \ru1(i). Thus,
$$ \mu^*(u)\geq 6 - 2\cdot 2 - 2\cdot\frac12 = 1.$$
\end{itemize}

\textbf{Case 4:} If $d(u)=5$, then, at worst, $u$ gives $2+\frac15$ away along each incident edge by \ru0 (or \ru1(ii)) and \ru2(i). As a result,
$$ \mu^*(u)\geq 5d(u)-14- \left(2+\frac15\right)d(u) = \left(2+\frac45\right)\cdot 5 - 14 = 0.$$

\textbf{Case 5:} If $6\leq d(u)\leq \Delta$, then, at worst, $u$ gives $2+\frac23$ away along each incident edge by \ru0 (or \ru1(ii)) and \ru2(ii). As a result,
$$ \mu^*(u)\geq 5d(u)-14-\left(2+\frac23\right)d(u) \geq \left(2+\frac13\right)\cdot 6 - 14 = 0.$$

To conclude, we started with a charge assignment with a negative total sum, but after the discharging procedure, which preserved this sum, we end up with a non-negative one, which is a contradiction. In other words, there exists no counter-example $G_2$ to \Cref{main theorem2}.

\bibliographystyle{plain}


\section{Introduction}

Neutron stars (NSs) are ideal compact astrophysical objects to probe fundamental physics. Due to their high central density that exceeds the saturation density of nuclear matter by several-fold, NSs  can efficiently test nuclear physics. For example, recent observations of $2M_\odot$ pulsars~\cite{Cromartie:2019kug}, x-ray emissions from hot spots on a rotating NS surface by NICER~\cite{Riley:2019yda,Miller:2019cac}, and gravitational waves from colliding NSs by LIGO and Virgo~\cite{TheLIGOScientific:2017qsa} have constrained properties of nuclear/quark matter and certain nuclear parameters (see e.g.~\cite{Abbott:2018exr,Annala:2017llu,Raithel:2018ncd,Lim:2018bkq,Bauswein:2017vtn,De:2018uhw,Most:2018hfd,Annala:2019puf,Malik2018,Zack:nuclearConstraints,Carson:2019xxz,Raithel:2019ejc,Chatziioannou:2020pqz,Raaijmakers:2019dks,Zimmerman:2020eho,Jiang:2019rcw,Dietrich:2020efo}). Due to their large compactnesses, NSs are also perfect sources to probe strong-field gravity. Indeed, observations of binary pulsars~\cite{Freire:2012mg,Yagi:2013qpa,Yagi:2013ava,Berti:2015itd,Shao:2017gwu,Archibald:2018oxs,Seymour:2019tir,Anderson:2019eay,Gupta:2021vdj} and the binary neutron star merger GW170817~\cite{Abbott:2018lct,Baker:2017hug,Sakstein:2017xjx,Ezquiaga:2017ekz,Zhao:2019suc,Silva:2020acr} have constrained various modifications to GR.

Some of the most well-studied modified theories of gravity are scalar-tensor theories in which one introduces scalar fields (either minimally or non-minimally coupled to the metric) to the action. A simple scalar-tensor theory proposed by  Damour and Esposito-Far\`ese has two theoretical parameters $(\alpha_0,\beta_0)$ and allows NSs to scalarize \emph{spontaneously}~\cite{Damour:1992we,Damour:1993hw,Damour:1996ke}. Namely, when $\beta_0$ is sufficiently negative, NSs can have scalar charges of order unity even when $\alpha_0 \ll 1$ (or equivalently, when the value of the scalar field at infinity is small). Such a phenomenon can be understood as NSs undergoing a tachyonic instability, where the effective mass of the scalar field becomes imaginary~\cite{Ramazanoglu:2017xbl}. A catalog of NS scalar charges is provided in~\cite{Anderson:2019hio} while a surrogate model has recently been constructed in~\cite{Zhao:2019suc} based on \emph{numerical} calculations. The parameters $(\alpha_0,\beta_0)$ have been constrained with solar system experiments, various binary pulsars~\cite{Freire:2012mg,Shao:2017gwu,Archibald:2018oxs,Anderson:2019eay} and the binary NS merger event GW170817~\cite{Zhao:2019suc} to be $|\alpha_0| \lesssim 3 \times 10^{-4}$ and $\beta_0 \gtrsim -4.4$. Future forecasts on probing this scalar-tensor theory with black-hole/NS binaries have been made with pulsar~\cite{Berti:2015itd} and gravitational-wave~\cite{Shao:2017gwu,Carson:2019fxr,Zhao:2019suc} observations.

The goal of this paper is to provide accurate expressions for scalar charges of NSs in the scalar-tensor theory in~\cite{Damour:1993hw} by solving the field equations \emph{analytically}. Such expressions are complementary to e.g. the surrogate model~\cite{Zhao:2019suc} mentioned earlier. To achieve this, instead of using realistic equations of state (EoSs) given in tables, we adopt the Tolman VII model~\cite{tolman} that approximates the energy density profile inside a NS to a quadratic function in a radial coordinate. Such a simple profile allows one to solve the Tolman-Oppenheimer-Volkoff (TOV) equations analytically in GR. Unfortunately, due to the complication of the modified TOV equations in scalar-tensor theories including the coupling with the scalar field, it is challenging to solve them analytically. To overcome this, we work in a weak-field approximation and solve the field equations order by order in the expansion\footnote{Similar calculations were carried out in quadratic gravity~\cite{Yagi:2015oca} and Einstein-\AE ther theory~\cite{Gupta:2021vdj} to respectively find  analytic scalar charges and sensitivities, which were then used in tandem with binary pulsar observations to place bounds on these theories.}. To take into account higher order contributions in the expansion, we apply Pad\'e resummation. We consider both spontaneous scalarization mentioned earlier and ``perturbative'' scalarization where the scalar charge is proportional to $\alpha_0$.

Figure~\ref{fig:charge-C-Tol} summarizes our findings, which compares the analytic scalar charges with the Tolman VII model in terms of the stellar compactness (the ratio between the mass and radius) against those with realistic EoSs (WFF1~\cite{Wiringa:1988tp}, SLy~\cite{SLy}, MPA1~\cite{1987PhLB..199..469M}, MS1~\cite{Mueller:1996pm}) computed numerically. Notice that the former can accurately model the latter, especially for the SLy EoS, for both the spontaneous and perturbative scalarization cases. We also find that when we plot the scalar charges as a function of the stellar binding energy (the difference between the gravitational and baryonic mass), the relation between these quantities becomes quasi-universal and is insensitive to EoSs with an error of $\sim 1\%$ (see~\cite{Pappas:2018csu} for other universal relations involving scalar charges in scalar-tensor theories). We estimate analytically the amount of the quasi-universality by comparing analytic expressions for scalar charges with the Tolman VII model and with constant density stars and find that it is quasi-universal with an error of 0.3\%. The technique developed here should easily be applicable to other theories beyond GR to compute the stellar charges (or sensitivities).

\begin{figure}[h]
\includegraphics[width=8.5cm]{charge-C-Tol_v2.pdf}
\caption{\label{fig:charge-C-Tol}
Comparison of analytic scalar charges with the Tolman VII model against numerical scalar charges with realistic EoSs as a function of the compactness. For the former, we show the results for both perturbative (solid) and spontaneous (dashed) scalarization. We present the results with four different choices of $\beta_0$. $\alpha_0$ is fixed to $\alpha_0=10^{-5}$ except for spontaneous scalarization for the Tolman case which assumes $\alpha_0=0$. Observe that the analytic Tolman results accurately model numerical results, especially with the SLy EoS.
}
\end{figure}

The rest of the paper is organized as follows. In Sec.~\ref{sec:ST}, we review the scalar-tensor theory proposed in~\cite{Damour:1993hw}, explain the modified TOV equations, present the relation between mass and radius, and review spontaneous scalarization following~\cite{Horbatsch:2010hj}. In Sec.~\ref{sec:analytic}, we present the formalism for computing scalar charges analytically in the Tolman VII model by combining weak-field expansions and Pad\'e resummation. We compare such analytic results against numerical ones in Sec.~\ref{sec:results} and also present the quasi-universal relation between the scalar charge and binding energy. We conclude in Sec.~\ref{sec:conclusion} and give several different directions for future work. In Appendix~\ref{sec:CD}, we repeat the analytic calculation for constant density stars and derive scalar charges similar to the Tolman VII case.
We use the geometric units of $c=G=1$. The main expressions for scalar charges are summarized in a Supplemental Mathematica notebook~\cite{mathematica}.

\section{Scalar-tensor Theories and Neutron Stars}
\label{sec:ST}

We begin by reviewing scalar-tensor theories and spontaneous scalarization of NSs.

\subsection{Theory}

Presented here are the action and field equations for scalar-tensor theories in the Einstein frame. The former is given by~\cite{Damour:1992we,Damour:1993hw,Damour:1996ke}
\begin{equation}
S=\int d^{4} x \frac{\sqrt{-g}}{16\pi}\left(R-2 g^{\mu \nu} \partial_{\mu} \varphi \partial_{\nu} \varphi\right)+S_{\mathrm{mat}}\left[\psi, A^{2}(\varphi) g_{\mu \nu}\right]\,,
\end{equation}
where $R$ is the Ricci scalar for the metric $g_{\mu\nu}$ in this frame, $g$ is its determinant, $\varphi$ is the scalar field and $\psi$ is the matter field. $A$ is the conformal factor that relates $g_{\mu\nu}$ and the metric $\tilde g_{\mu\nu}$ in the (physical) Jordan frame as $\tilde g_{\mu\nu} = A(\varphi) g_{\mu\nu}$. Notice that we have set $G=1$ for the bare gravitational constant $G$.
Varying the above action with respect to $g_{\mu\nu}$ and $\varphi$, the field equations are given by~\cite{Damour:1993hw,Damour:1996ke}
\begin{eqnarray}
\label{eq:Ein_eq}
R_{\mu \nu}&=&2 \partial_{\mu} \varphi \partial_{\nu} \varphi+8\pi\left(T_{\mu \nu}-\frac{1}{2} g_{\mu \nu} T\right)\,, \\
\label{eq:scalar_eq}
\square \varphi&=&-4\pi \alpha(\varphi) T\,,
\end{eqnarray}
where $T_{\mu\nu}$ is the matter stress-energy tensor and
\begin{equation}
\alpha (\varphi) \equiv \frac{\partial \ln A(\varphi)}{\partial\varphi}\,.
\end{equation}

In this paper, we consider an example scalar-tensor theory first considered by Damour and Esposito-Far\`ese~\cite{Damour:1992we,Damour:1993hw,Damour:1996ke} with
\begin{equation}
A(\varphi) = \exp\left( \frac{\beta_0}{2} \varphi^2 \right)\,.
\end{equation}
Another parameter of the theory is  $\alpha_0 \equiv \alpha(\varphi_0) = \beta_0 \varphi_0$ with $\varphi_0$ representing the scalar field at infinity.

\subsection{Neutron Stars}

To construct a NS solution, we begin with the metric ansatz given by
\begin{equation}
ds^2 = - e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2 (d\theta^2 + \sin^2 \theta d\phi^2)\,.
\end{equation}
For the matter sector, we consider a perfect fluid whose stress-energy tensor in the Einstein frame is given by~\cite{Pani:2014jra}
\begin{equation}
T_{\mu\nu} = A^4(\varphi) \left[(\widetilde \rho + \widetilde P)u_\mu u_\nu + g_{\mu\nu}\widetilde P \right]\,,
\end{equation}
where $u^\mu$ is the 4-velocity of the fluid. Notice that $\widetilde \rho$ and $\widetilde P$ are the energy density and pressure in the Jordan frame that directly enters in the EoS $\widetilde P (\widetilde \rho)$. Plugging these into the field equations in Eqs.~\eqref{eq:Ein_eq} and~\eqref{eq:scalar_eq}, one finds~\cite{Damour:1993hw,Damour:1996ke}\footnote{There is also an equation for $\nu$ that we do not present here since it is unnecessary for deriving scalar charges.}
\begin{eqnarray}
\label{eq:dmdr}
m^{\prime}&=&4 \pi r^{2} A^{4}(\varphi) \widetilde{\rho}+\frac{1}{2} r(r-2 m) \varphi'^{2}\,, \\
\label{eq:dPdr}
\widetilde{P}^{\prime}&=&-(\widetilde{\rho}+\widetilde{P})\left[\frac{m+4 \pi  A^4 \widetilde{P} r^3}{r(r-2 m)}+\frac{1}{2} r \varphi'^{2} +\alpha(\varphi) \varphi'\right]\,, \nonumber \\
\\
\label{eq:dphidrr}
\varphi'^{\prime}&=&4 \pi  \frac{r }{r-2 m}A^{4}(\varphi) \left[\alpha(\varphi)(\widetilde{\rho}-3 \widetilde{P})+r \varphi'(\widetilde{\rho}-\widetilde{P}) \right] \nonumber \\
&&-\frac{2(r-m)}{r(r-2 m)} \varphi'\,,
\end{eqnarray}
where a prime represents an $r$ derivative and
\begin{equation}
e^{-\lambda(r)} = 1 - \frac{2m(r)}{r}\,.
\end{equation}

The asymptotic behavior of the scalar field at infinity is given by
\begin{equation}
\label{eq:varphi_asympt}
\varphi = \varphi_0 - \alpha_A\frac{M_A}{r} + \mathcal{O}\left(  \frac{M_A^2}{r^2}\right)\,,
\end{equation}
where
\begin{equation}
\alpha_A = \frac{\partial \ln M_A}{\partial \varphi_0}
\end{equation}
is the scalar charge while $M_A$ is the total gravitational mass that can be read off from the asymptotic behavior of $m(r)$ at infinity as
\begin{equation}
\label{eq:m_inf}
m(r) = M_A + \mathcal{O}\left(  \frac{M_A}{r}\right)\,.
\end{equation}

We can construct NS solutions as follows. First, we impose the initial conditions of
\begin{equation}
m(0) = 0\,, \quad \widetilde P(0) = \widetilde P_c\,, \quad \varphi (0) = \varphi_c\,, \quad \varphi'(0) = 0\,,
\end{equation}
where $\widetilde P_c$ is the central pressure while $\varphi_c$ is the central value for the scalar field. Under these initial conditions together with a choice of $(\alpha_0,\beta_0)$ and an EoS, we solve Eqs.~\eqref{eq:dmdr}--\eqref{eq:dphidrr} in the interior region numerically and the stellar radius $R$ is determined by $\widetilde P(R) = 0$. We then solve the equations in the exterior region with $\widetilde P = \widetilde \rho = 0$. Finally, we read off $M_A$, $\varphi_0$ and $\alpha_A$ by comparing the asymptotic behavior of the numerically-computed $\varphi$ and $m$ with Eqs.~\eqref{eq:varphi_asympt} and~\eqref{eq:m_inf}. Figure~\ref{fig:MR} presents the mass-radius relation of NSs in GR and the scalar-tensor theory for four representative EoSs (WFF1~\cite{Wiringa:1988tp}, SLy~\cite{SLy}, MPA1~\cite{1987PhLB..199..469M}, MS1~\cite{Mueller:1996pm}) with different stiffness\footnote{These numerical results are computed with a Mathematica notebook developed in~\cite{Pani:2014jra}.}. ``Humps'' in the relation for the scalar-tensor theory correspond to NSs with spontaneous scalarization.

\begin{figure}[h]
\includegraphics[width=8.5cm]{M-R.pdf}
\caption{\label{fig:MR} Mass-radius relation for NSs with four representative EoSs in GR ($\alpha_0=\beta_0=0$) and the scalar-tensor theory ($\alpha_0=10^{-5}$ and $\beta_0 = -4.4$, $-4.5$). The latter deviates from the former due to spontaneous scalarization. Two cyan dots correspond to the NSs of $1.8M_\odot$ with SLy EoS whose energy density profiles are shown in Fig.~\ref{fig:density_profile}.
}
\end{figure}

\subsection{Spontaneous Scalarization}
\label{sec:spnt}

An analytic attempt of computing scalar charges with spontaneous scalarization was taken in~\cite{Horbatsch:2010hj}. First, when $\alpha$ at the center of a star, $\alpha_c$, is small, one can expand both $\alpha_0$ and $\alpha_A$ in terms of $\alpha_c$ as\footnote{The scalar field (and its derivative) enter in even powers in Eqs.~\eqref{eq:dmdr} and~\eqref{eq:dPdr}, and in odd powers in Eq.~\eqref{eq:dphidrr}. This means that $\alpha_0$ enters in odd powers in $\alpha_c$ and $\alpha_A$. Inverting the former order by order (and substituting it to the latter), one finds that $\alpha_0$ and $\alpha_A$ enter in odd powers in $\alpha_c$.}
\begin{eqnarray}
\label{eq:alpha0_in_alphac}
\alpha_0 &=& d_1 \alpha_c + d_2 \alpha_c^3 + \mathcal{O}(\alpha_c^5)\,, \\
\label{eq:alphaA_in_alphac}
\alpha_A &=& e_1 \alpha_c + e_2 \alpha_c^3 + \mathcal{O}(\alpha_c^5)\,.
\end{eqnarray}
Next, one can solve Eq.~\eqref{eq:alpha0_in_alphac} for $\alpha_c$ as
\begin{equation}
\label{eq:alpha_c}
\alpha_c = \omega C_+ + \bar \omega C_-\,,
\end{equation}
with
\begin{equation}
\omega = 1, \quad -e^{\pm i \pi/3}\,,
\end{equation}
and
\begin{equation}
C_\pm = \left( \frac{\alpha_0}{2d_2} \pm \sqrt D \right)^{1/3}\,, \quad D = \left(\frac{\alpha_0}{2 d_{2}}\right)^{2}+\left(\frac{d_{1}}{3 d_{2}}\right)^{3}\,.
\end{equation}
From Eqs.~\eqref{eq:alphaA_in_alphac} and~\eqref{eq:alpha_c}, the scalar charge is given by
\begin{equation}
\alpha_A =  (\omega C_+ + \bar \omega C_-) e_1 +( \omega C_+ + \bar \omega C_-)^3 e_2  + \mathcal{O}(\alpha_c^5)\,.
\end{equation}
For example, when $\alpha_0 = 0$\footnote{Given current stringent bounds on $\alpha_0$, $\alpha_0 \approx 0$ is a valid approximation when discussing spontaneous scalarization.}, $\alpha_c$ can take non-vanishing values as $\alpha_c = \pm \sqrt{-d_1/d_2}$ and
\begin{equation}
\label{eq:charge_spnt}
\alpha_{A}=\pm\left[\left(-\frac{d_{1}}{d_{2}}\right)^{1 / 2} e_{1}+\left(-\frac{d_{1}}{d_{2}}\right)^{3 / 2} e_{2}\right]\,.
\end{equation}
Thus, although $\alpha_0 = 0$, the scalar charge becomes non-vanishing when $d_1/d_2<0$.

For constant density stars, $d_1$ and $e_1$ are given in a closed analytic form as~\cite{Horbatsch:2010hj}
\begin{eqnarray}
\label{eq:d1_CD}
d_{1} &=&\operatorname{HeunG}\left( \widetilde{a}, \widetilde{q} ; \widetilde{\alpha}_-, \widetilde{\alpha}_+, \frac{3}{2}, \frac{3}{2} ; Z \right) \nonumber \\
&&- \left(\widetilde Z \log \widetilde Z\right)  \, \operatorname{HeunG}'\left(\widetilde{a}, \widetilde{q} ; \widetilde{\alpha}_-, \widetilde{\alpha}_+, \frac{3}{2}, \frac{3}{2} ; Z\right)\,, \nonumber \\
\\
\label{eq:e1_CD}
e_{1} &=&\beta_0 \widetilde Z  \operatorname{HeunG}' \left(\widetilde{a}, \widetilde{q} ; \widetilde{\alpha}_-, \widetilde{\alpha}_+, \frac{3}{2}, \frac{3}{2} ; Z \right)\,,
\end{eqnarray}
with
\allowdisplaybreaks
\begin{eqnarray}
\widetilde{a}&=&-\frac{1}{1+3 \widetilde{p}_{c}}\,, \\
\widetilde{q}&=&\frac{3 \beta_0}{2}\left(\frac{3 \widetilde{p}_{c}-1}{3 \widetilde{p}_{c}+1}\right)\,, \\
\widetilde{\alpha}_\pm &= &\frac{3}{2}\left(1\pm \sqrt{1-\frac{8 \beta_0}{3}} \right)\,, \\
Z &=& \frac{\widetilde p_c}{1+3\widetilde p_c}\,, \\
\widetilde Z &= & \frac{1+\widetilde p_{c}}{1+3 \widetilde p_{c}}\,,
\end{eqnarray}
and $\widetilde p_c = \widetilde P_c/\widetilde \rho_c$ with $\widetilde P_c$ and $\widetilde \rho_c$ being the central pressure and energy density.  HeunG is the general Heun function (a generalization of the hypergeometric function) while its prime refers to a derivative with respect to $Z$.

Figure~\ref{fig:d1} shows $d_1$ for constant density stars as a function of the stellar compactness $C_A$ defined by
\begin{equation}
C_A = \frac{M_A}{R}\,.
\end{equation}
One can check numerically that spontaneous scalarization occurs when $d_1 < 0$~\cite{Horbatsch:2010hj}. This means that $d_2 > 0$ in the relevant parameter region, which we checked with the analytic expression for $d_2$ in Appendix~\ref{sec:CD}.
Notice that
 $\beta_0 = -4.33$~\footnote{See e.g.~\cite{Harada:1997mr} for a related, analytic estimate of the critical $\beta_0$ with constant density stars.
The critical value may vary slightly for neutron stars with other equations of state.} is the critical value below which the spontaneous scalarization is realized for a certain range of $C_A$.
In~\cite{Horbatsch:2010hj}, closed analytic expressions are not found for $d_2$ and $e_2$.

\begin{figure}[h]
\includegraphics[width=8.5cm]{d1_CD.pdf}
\caption{\label{fig:d1} The linear coefficient $d_1$ in Eq.~\eqref{eq:alpha0_in_alphac} as a function of the stellar compactness $C_A$ for constant density stars with $\alpha_0=0$ and various choices of $\beta_0$. Spontaneous scalarization is realized when $d_1<0$.
}
\end{figure}

\section{Analytic Scalar Charges}
\label{sec:analytic}

We now explain how we can use the weak-field expansion and Pad\'e resummation to analytically calculate NS scalar charges with the Tolman VII model.

\subsection{Formalism}

To construct approximate, analytic stellar solutions in the scalar-tensor theory, we work under a weak-field expansion. Namely, we decompose each unknown function $f(r)$ as
\begin{equation}
\label{eq:wf}
f(r) = \sum_{k=0} f_k(r) \epsilon^k\,,
\end{equation}
where $\epsilon$ is a book-keeping parameter that counts the order of the GR compactness $C_0 = M_0/R$ where $M_0$ is the GR mass. We have
\begin{equation}
m_0(r) = \widetilde \rho_0(r) = \widetilde{p}_0(r) =   \widetilde{p}_1(r) = 0\,,
\end{equation}
and thus to leading order, $m(r) = \mathcal{O}(\epsilon)$, $\widetilde \rho(r) = \mathcal{O}(\epsilon)$ and $\widetilde p(r) = \mathcal{O}(\epsilon^2)$. Notice that $\widetilde \rho \gg \widetilde p$ in the weak-field limit. For the scalar field, we find it convenient to introduce
\begin{equation}
\bar \varphi (r) \equiv \frac{\varphi (r)}{\alpha_0}\,,
\end{equation}
so that $\bar \varphi (r) = \mathcal{O}(\alpha_0^0)$ when $\alpha_0 \ll 1$. We decompose $\bar \varphi(r)$ as in Eq.~\eqref{eq:wf} and we have
\begin{equation}
\bar \varphi_0(r) = \frac{1}{\beta_0}\,.
\end{equation}

To see the effect of spontaneous scalarization mentioned in Sec.~\ref{sec:spnt}, we need to derive the scalar charge $\alpha_A$ valid up to $\mathcal{O}(\alpha_0^3)$ (which, in turn, means that it is valid to $\mathcal{O}(\alpha_c^3)$ as in Eq.~\eqref{eq:alpha0_in_alphac}). Namely, we need to find solutions up to to $\mathcal{O}(\alpha_0^2)$ higher than the leading. For this, we seek for a solution $\bar \varphi$ to $\mathcal{O}(\alpha_0^2)$. For $m(r)$, since $m_k = \mathcal{O}(\alpha_0^2)$ when $k \geq 2$ (namely, $m(r) = m_1(r)$ in GR), we need to find a solution for $m_k(r)$ valid to $\mathcal{O}(\alpha_0^4)$ such that $\alpha_A$ is valid up to $\mathcal{O}(\alpha_0^3)$. Since $\widetilde p$ does not enter in the $m'(r)$ equation in Eq.~\eqref{eq:dmdr}, we only need a solution for $\widetilde p$ to  $\mathcal{O}(\alpha_0^2)$. Next, we substitute the weak-field expansion of each unknown function in Eq.~\eqref{eq:wf} to the differential equations in Eqs.~\eqref{eq:dmdr}--\eqref{eq:dphidrr} and solve order by order in $\epsilon$. Below, we present the equations and solutions mainly to leading order in $\epsilon$ for an illustration purpose. In the actual calculation, we derived higher order contributions and applied a Pad\'e resummation that we explain later.

\subsection{Interior Solutions }

\begin{figure}[h]
\includegraphics[width=8.5cm]{density_profile.pdf}
\caption{\label{fig:density_profile}
Energy density profiles (normalized by the central value) for NSs with $1.8M_\odot$ in GR and in the scalar-tensor theory with $\alpha_0=10^{-5}$ and $\beta_0 = -4.5$ (corresponding to cyan dots in Fig.~\ref{fig:MR}). We use the SLy EoS. For comparison, we also show the profile for the Tolman VII model given in Eq.~\eqref{eq:Tolman_density}. Observe that such a model can accurately describe the profile for realistic NSs.
}
\end{figure}

Since realistic EoSs are given in tables, one can only construct stellar solutions numerically. To overcome this, we use the Tolman VII model~\cite{tolman,Jiang:2019vmf} that is known to accurately model realistic NSs in GR.

The density profile is given by a simple quadratic form as
\begin{equation}
\label{eq:Tolman_density}
\widetilde \rho_1 =  \frac{15 M_0}{8 \pi
   R^3} \left(1-\frac{r^2}{R^2}\right)\,, \quad \widetilde \rho_i (r) = 0 \quad (i\geq 2)\,.
\end{equation}
Thus, these stars are parameterized by $M_0$ and $R$, where the former is the stellar mass in GR. To check the validity of the above profile, we present in Fig.~\ref{fig:density_profile} the normalized energy density profile for NSs with $1.8M_\odot$ and the SLy EoS found numerically for both GR and scalar-tensor theory. For the latter, the star is spontaneously scalarized (see Fig.~\ref{fig:MR}). Notice that the energy density profiles in both theories are almost identical and can be approximated by Eq.~\eqref{eq:Tolman_density} shown by the black solid curve. This justifies the use of the Tolman VII model even for NSs in scalar-tensor theory.

We next derive the field equations and solutions in the interior region. The leading differential equations are given by
\begin{eqnarray}
m_1'(r) &=&\frac{15}{2} \left( 1+ 4 \frac{\alpha_0^2}{\beta_0} +8 \frac{\alpha_0^4}{\beta_0^2}  \right)\frac{M_0}{R^5} r^2 (R^2-r^2)+ \mathcal{O}(\alpha_0^6)\,, \nonumber \\
\\
\bar \varphi_1''(r) &=&  -\frac{2 \bar \varphi_1'}{r} +\frac{15}{2} \left( 1+ 4 \frac{\alpha_0^2}{\beta_0} \right)\frac{M_0}{
   R^5}(R^2-r^2)+ \mathcal{O}(\alpha_0^4)\,, \nonumber \\
   \\
\widetilde P_2'(r) &=& -\frac{15 M_0}{8 \pi   R^5}(R^2-r^2)\left( \frac{m_1}{r^2} + \alpha_0^2 \bar\varphi_1'\right)\,.
   \end{eqnarray}
Imposing the boundary conditions
\begin{equation}
\label{eq:BC1}
m_k(0)=0\,, \quad \bar \varphi_k(0)=\bar \varphi_{1c}\,, \quad \widetilde P_k(R)=0\,,
\end{equation}
we find
\begin{eqnarray}
m_1^\mathrm{(int)}(r) &=& \frac{5}{2} \left( 1+ 4 \frac{\alpha_0^2}{\beta_0} +8 \frac{\alpha_0^4}{\beta_0^2}  \right)\frac{M_0}{R^3} r^3 \left( 1  - \frac{3}{5} \frac{r^2}{R^2} \right) \nonumber \\
   &&+ \mathcal{O}(\alpha_0^6)\,, \\
\bar \varphi_1^\mathrm{(int)}(r) &=& \bar \varphi_{1c}+\frac{5}{4}  \left( 1+ 4 \frac{\alpha_0^2}{\beta_0} \right)\frac{M_0 }{
    R^3} r^2  \left( 1  - \frac{3}{10} \frac{r^2}{R^2} \right)\nonumber \\
   & &+ \mathcal{O}(\alpha_0^4)\,, \\
   \widetilde P_2^\mathrm{(int)}(r) &=& \frac{15}{16 \pi  } \left( 1+  \frac{4+\beta_0 }{\beta_0} \alpha_0^2\right)\frac{M_0^2}{R^4} \nonumber \\
   &&\times\left(1-\frac{r^2}{R^2}\right)^2 \left(1-\frac{1}{2}\frac{r^2}{R^2}\right)+ \mathcal{O}(\alpha_0^4)\,.
\end{eqnarray}

\subsection{Exterior Solutions and Perturbative Scalar Charges}

Next, we study the exterior region. By setting $\widetilde P = \widetilde \rho = 0$ in Eqs.~\eqref{eq:dmdr}--\eqref{eq:dphidrr}, we find
\begin{equation}
m_1'(r) = 0\,, \quad \bar \varphi_1''(r)=  -\frac{2 \bar \varphi_1'}{r}\,.
\end{equation}
Imposing regularity at infinity, we can solve the above equations to find the exterior solution as\footnote{An analytic solution for a NS exterior in scalar-tensor theories without the weak-field expansion has been found in a different coordinate system~\cite{Damour:1992we,Damour:1996ke}. We found that working in the same coordinate system as the interior case is easier due to some nonlinearity in the coordinate transformation.}
\begin{equation}
\label{eq:m1_phi1_ext}
m_1^\mathrm{(ext)}(r) = M_1\,, \quad \bar \varphi_1^\mathrm{(ext)}(r) = \frac{\bar \omega_1}{r}\,,
\end{equation}
where $M_1$ and $\bar \omega_1$ are integration constants.
We can determine the integration constants both in the interior and exterior solutions by imposing the boundary conditions at the surface:
\begin{eqnarray}
m_k^\mathrm{(int)}(R) &=& m_k^\mathrm{(ext)}(R)\,, \\
\bar \varphi_k^\mathrm{(int)}(R) &=& \bar \varphi_k^\mathrm{(ext)}(R)\,, \\
\bar \varphi_k'^\mathrm{(int)}(R) &=& \bar \varphi_k'^\mathrm{(ext)}(R)\,.
\end{eqnarray}
Similar to Eq.~\eqref{eq:m1_phi1_ext}, we can introduce integration constants $M_k$ and $\bar \omega_k$ at order $\epsilon^k$ as
\begin{eqnarray}
m_k^\mathrm{(ext)}(r) &=& M_k\left[1 + \mathcal{O}\left( \frac{M_0}{r} \right)\right]\,, \\
\bar \varphi_k^\mathrm{(ext)}(r) &=& \frac{\bar \omega_k}{r}+ \mathcal{O}\left( \frac{M_0^2}{r^2} \right)\,.
\end{eqnarray}

Once the integration constants are determined, we can compute the perturbative scalar charge from Eq.~\eqref{eq:varphi_asympt} as
\begin{equation}
\label{eq:charge}
\alpha_A
=  - \alpha_0 \frac{\sum_{k=1}^N \bar \omega_k \epsilon^k}{
\sum_{k=1}^N M_k\epsilon^k}\,.
\end{equation}
In this paper, we computed up to $N=10$.
Equation~\eqref{eq:charge} allows us to express $\alpha_A$ in a series of the GR compactness $C_0$. It would be more useful to express the scalar charges in terms of the physical compactness in the scalar-tensor theory defined by
\begin{equation}
C_A = \frac{M_A}{R}\,, \quad M_A = \sum_{k=1} M_k\epsilon^k\,.
\end{equation}
We can solve order by order to find $C_0$ in terms of $C_A$ and plug this into Eq.~\eqref{eq:charge}. To have the results consistent up to the order analyzed, we expand Eq.~\eqref{eq:charge} about $\alpha_0=0$ and $C_A=0$ and keep only to $\mathcal{O}(\alpha_0^3,C_A^{N})$. To have the series converge, we then construct a Pad\'e approximant of order $N/2$ (when $N$ is an even number) in $C_A$.

At leading order, the integration constants are derived as
\begin{eqnarray}
\bar \varphi_{1c} &=&- \frac{15}{8}  \left( 1+ 4 \frac{\alpha_0^2}{\beta_0} \right) \frac{M_0}{R}+ \mathcal{O}(\alpha_0^4)\,. \\
\label{eq:M1_Tol}
M_1 &=& \left( 1+ 4 \frac{\alpha_0^2}{\beta_0} +8 \frac{\alpha_0^4}{\beta_0^2}  \right) M_0+ \mathcal{O}(\alpha_0^6)\,, \\
\label{eq:omega1_Tol}
\bar \omega_1 &=& -\left( 1+ 4 \frac{\alpha_0^2}{\beta_0} \right) M_0+ \mathcal{O}(\alpha_0^4)\,.
\end{eqnarray}
\newpage
Using Eq.~\eqref{eq:charge}, the scalar charge to this order is given by
\begin{eqnarray}
\label{eq:charge_leading}
\alpha_A &=& \alpha_0 \frac{\left( 1+ 4 \frac{\alpha_0^2}{\beta_0} \right)}{\left( 1+ 4 \frac{\alpha_0^2}{\beta_0} +8 \frac{\alpha_0^4}{\beta_0^2}  \right)} \nonumber \\
&=& \alpha_0 + \mathcal{O}(\alpha_0^5, C_A)\,.
\end{eqnarray}
We have derived $\alpha_A$ valid to $\mathcal{O}(C_A^{10})$. The first few terms are given by
\begin{eqnarray}
\alpha_A &= &\alpha_0-\frac{10}{7}
   \alpha_0 \left(\beta_0+1\right)  C_A+ \frac{5\alpha_0}{3003} \nonumber \\
  && \times \left(1253 \beta_0^2+1514
   \beta_0-1126\right)  C_A^2 + \mathcal{O}(\alpha_0^3, C_A^3)\,. \nonumber \\
\end{eqnarray}
We then construct a Pad\'e approximant to 5th order in $C_A$, which we provide in a Supplemental Mathematica notebook~\cite{mathematica}.

\subsection{Spontaneous Scalarization}

So far, we have constructed $\alpha_A$ in a series of $\alpha_0$ and $C_A$ in a form
\begin{equation}
\label{eq:alphaA_in_alpha0}
\alpha_A = \bar e_1(C_A) \alpha_0 + \bar e_3(C_A) \alpha_0^3 + \mathcal{O}(\alpha_0^5)\,,
\end{equation}
where $\bar e_1$ and $\bar e_3$ are functions of $C_A$ valid to $\mathcal{O}(C_A^N)$. We also have
\begin{eqnarray}
\label{eq:alphac_in_alpha0}
\alpha_c &=& \alpha_0 \beta_0 \sum_{k=1}^N \bar \varphi_{kc} \epsilon^k \nonumber \\
&=& \bar d_1(C_A) \alpha_0 + \bar d_3(C_A) \alpha_0^3 + \mathcal{O}(\alpha_0^5)\,.
\end{eqnarray}
We can invert Eq.~\eqref{eq:alphac_in_alpha0} order by order in $\alpha_c$ to find $\alpha_0(\alpha_c)$ in a form given by Eq.~\eqref{eq:alpha0_in_alphac}. Then, we substitute this to Eq.~\eqref{eq:alphaA_in_alpha0} and expand about $\alpha_c=0$ to find $\alpha_A(\alpha_c)$ in a form in Eq.~\eqref{eq:alphaA_in_alphac}. Following Sec.~\ref{sec:spnt}, we can then find $\alpha_A$ for NSs under spontaneous scalarization.

In fact, it turns out that the contribution for the second term in Eq.~\eqref{eq:charge_spnt} is negligible and we focus on deriving $d_1$, $d_2$, and $e_1$.
We first find these coefficients in terms of a series expansion of $C_A$ valid to $\mathcal{O}(C_A^{10})$. The first few terms of these functions are given by
\newline
\begin{widetext}
\begin{eqnarray}
d_1 &=& 1+\frac{15 \beta_0 }{8}C_A +\frac{5\beta_0}{128}  (17
   \beta_0-14)C_A^2+ \frac{\beta_0}{7168} \left(695 \beta_0^2-2294
   \beta_0-27720\right)
 C_A^3+ \mathcal{O}(C_A^4)\,,\\
d_2 & =&-\frac{1955
   \beta_0}{448} C_A^2+\frac{(1630328-3417469 \beta_0)
   \beta_0}{512512} C_A^3
+ \mathcal{O}(C_A^4)\,, \\
   e_1 &=& 1+\frac{5}{56} (5
   \beta_0-16) C_A
   +\frac{5}{384384}
   \left(5515 \beta_0^2-54170 \beta_0-144128\right) C_A^2 \nonumber \\
&&   + \frac{1}{156828672}\left(941985 \beta_0^3-19472682
   \beta_0^2-258360664 \beta_0-441114624\right) C_A^3+ \mathcal{O}(C_A^4)\,.
\end{eqnarray}
\end{widetext}
Next, we construct the Pad\'e approximants on these. We first tried 5th order Pad\'e approximants but found that there are some unphysical divergences (Fig.~\ref{fig:charge-C-Tol_conv}). Instead, we use 4th order Pad\'e approximants. We substituted these into the following expression to find scalar charges for NSs under spontaneous scalarization:
\begin{equation}
\label{eq:charge_spnt_simple}
\alpha_{A}=\pm\left(-\frac{d_{1}}{d_{2}}\right)^{1 / 2} e_{1}\,.
\end{equation}
The final expression is available in the Supplemental Mathematica notebook~\cite{mathematica}.

\section{Comparison with Numerical Results}
\label{sec:results}

Let us now compare the analytic scalar charge expression for the Tolman VII model with scalar charges for realistic NSs found numerically. Scalar charges can be both positive and negative but we focus on the former.

\subsection{Scalar Charges vs Compactness}

\begin{figure}[h]
\includegraphics[width=8.5cm]{charge-C-Tol_v2_convergence_with_spnt_Pade5.pdf}
\caption{\label{fig:charge-C-Tol_conv} Convergence of scalar charge calculations with the Tolman model for $(\alpha_0,\beta_0)=(10^{-5},-4.2)$  (top) and $(\alpha_0,\beta_0)=(0,-4.5)$ (bottom) under various approximations. We compare the Taylor series results and various Pad\'e resummation results against numerical results for the SLy EoS. Observe that Pad\'e resummation is crucial to accurately model the numerical result. The expressions for the black solid curves are the ones that we will use in the remaining part of this paper.}

\end{figure}

First, let us check the convergence of the analytic result for scalar charges in terms of the order of the Pad\'e resummation. Figure~\ref{fig:charge-C-Tol_conv} compares the analytic scalar charges at various orders of the Pad\'e approximants. We show the results for the perturbative scalarization case $(\alpha_0,\beta_0)=(10^{-5},-4.2)$ and the spontaneous scalarization case $(\alpha_0,\beta_0)=(0,-4.5)$. We also present Taylor series (without the Pad\'e resummation) up to $\mathcal{O}(C_A^{10})$ for perturbative scalarization and to $\mathcal{O}(C_A^{8})$ for spontaneous scalarization, as well as numerical results computed with the SLy EoS. Observe that for the perturbative scalarization case (top panel), the Pad\'e resummation at 3rd order is already a good approximation and is converging fast to the numerical result\footnote{The 3rd and 5th order Pad\'e results for the \emph{perturbative} scalar charge expressions have a noticeable difference for lower $\beta_0$ where the expressions diverge due to the onset of spontaneous scalarization (not shown in Fig.~\ref{fig:charge-C-Tol_conv}) and thus we use the 5th order result in the remainder of this section.}. On the other hand, the Taylor series result becomes less accurate when the compactness is relatively high. For the spontaneous scalarization case (bottom panel), there is little difference between the 4th and 5th order Pad\'e results, though the latter has some unphysical divergence. We therefore use the former to avoid this divergence. Additionally, the former has the advantages of a simpler functional expression and better agreement with the numerical SLy data.
All in all, these findings suggest that it is crucial to perform the resummation to find an accurate modeling.

Having understood the convergence, let us now carry out the comparison in more detail for various values of $\beta_0$ and EoSs. Figure~\ref{fig:charge-C-Tol} shows such a comparison. Observe that the analytic scalar charges for the Tolman VII case beautifully captures the numerical results for realistic NSs, especially for the SLy EoS. Notice that the agreement between the analytic and numerical results is good even for stars with spontaneous scalarization. This shows that the analytic result serves as an accurate, ready-to-use expression for the NS scalar charge in this scalar-tensor theory.

\subsection{Scalar Charges vs Binding Energy}

We next look at the relation between the scalar charge and the stellar binding energy, where the latter is defined by
\begin{eqnarray}
\Omega_A &=&-\frac{1}{2} \int d^{3} x \widetilde \rho(r) \int d^{3} x^{\prime} \frac{\widetilde \rho\left(r^{\prime}\right)}{\left|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right|}\,, \\
&=& - 16 \pi^2\int_0^R dr  r^2 \widetilde \rho(r) \left(\int_r^R dr' r' \widetilde \rho(r')\right)\,,
\end{eqnarray}
where the second equality is valid only for spherically symmetric systems. Physically, this quantity measures the difference between the gravitational and baryonic mass of a star and can be used to e.g. probe certain formation scenario of NSs~\cite{Newton:2016weo}. For the Tolman VII model, the relation between the binding energy and compactness is given by
\begin{equation}
\label{eq:bind_comp_Tol}
\frac{\Omega_A}{M_A} = -\frac{5}{7}C_A + \mathcal{O}(\alpha_0^2)\,.
\end{equation}
Such a GR relation is sufficient to find the scalar charge expression in terms of the binding energy for perturbative scalar charges since the latter is already proportional to $\alpha_0$.
We can invert this relation and substitute it to the $\alpha_A(C_A)$ expressions.

\begin{figure}[h]
\includegraphics[width=8.5cm]{charge-Bind-Tol_v2.pdf}
\caption{\label{fig:charge-Bind-Tol}
Similar to Fig.~\ref{fig:charge-C-Tol} but as a function of the binding energy. Notice that the results are quasi-universal and insensitive to the EoSs for fixed $\beta_0$.
}
\end{figure}

Figure~\ref{fig:charge-Bind-Tol} presents the scalar charge as a function of the binding energy (normalized by the stellar mass) for the same choices of $\alpha_0$ and $\beta_0$ and EoSs as in Fig.~\ref{fig:charge-C-Tol}. Observe  that numerical results are quasi-universal among different EoSs for a fixed $\beta_0$. Observe also that the analytic result accurately describes the numerical ones for perturbative scalarization. By taking the fractional difference between the numerical and analytic results for $\beta_0 = -4.2$ or $\beta_0 = -4.3$, we found that the quasi-universality holds to $\sim 1\%$ for $\Omega_A/M_A > -0.2$.

Analytic expressions help us to study the quasi-universal relation in more detail. When we expand the perturbative scalar charge expression for the Tolman VII model for small binding energy, we find
\begin{eqnarray}
\label{eq:alpha_bind_Tol}
\alpha_A^{(\mathrm{Tol})} &= & \alpha_0 +2 \alpha_0
   \left(\beta_0+1\right) \frac{\Omega_A}{M_A} \nonumber \\
&&   + \frac{7\alpha_0}{2145
  } \left(1253 \beta_0^2+1514
   \beta_0-1126\right)\frac{\Omega_A^2}{M_A^2} \nonumber \\
   &&+\mathcal{O}\left( \alpha_0^3, \frac{\Omega_A^3}{M_A^3} \right)\,.
\end{eqnarray}
In Appendix~\ref{sec:CD}, we derive analytic scalar charges for constant density stars. We show that although $\alpha_A(C_A)$ is quite different from those for Tolman VII and realistic EoS cases, $\alpha_A(\Omega_A/M_A)$ is similar to the latter two cases, supporting the quasi-universal relation. For the perturbative scalar charge, the Taylor-series expansion in small binding energy is given by
\begin{eqnarray}
\alpha_A^{(\mathrm{CD})} &= & \alpha_0 +2 \alpha_0
   \left(\beta_0+1\right) \frac{\Omega_A}{M_A} \nonumber \\
&&   + \frac{5\alpha_0}{21
   }  \left(17
   \beta_0^2+20 \beta_0-16\right)\frac{\Omega_A^2}{M_A^2} \nonumber \\
   && +\mathcal{O}\left( \alpha_0^3, \frac{\Omega_A^3}{M_A^3} \right)\,.
\end{eqnarray}
Comparing this with Eq.~\eqref{eq:alpha_bind_Tol}, we see that the two expressions are identical up to $\mathcal{O}(\Omega_A/M_A)$ and the difference only appears at $\mathcal{O}(\Omega_A^2/M_A^2)$. When taking the fractional difference between $\alpha_A^{(\mathrm{Tol})}$ and $\alpha_A^{(\mathrm{CD})}$, we find
\begin{eqnarray}
\label{eq:alphaA_frac_diff}
\frac{\alpha_A^{(\mathrm{Tol})} - \alpha_A^{(\mathrm{CD})}}{\alpha_0} &=& \frac{2 \left(311 \beta_0^2+1343
   \beta_0+1013\right)}{15015}\frac{\Omega_A^2}{M_A^2} \nonumber \\
  & & +\mathcal{O}\left( \alpha_0^3, \frac{\Omega_A^3}{M_A^3} \right)\,.
\end{eqnarray}
Figure~\ref{fig:charge-Bind-frac-diff} presents this fractional difference for $\beta_0 = -4.2$ and $-4.3$. Observe that the Tolman and constant density cases agree within an error of $\mathcal{O}(10^{-3})$. This provides strong analytic support for  the quasi-universality of the relation between $\alpha_A$ and $\Omega_A/M_A$.

\begin{figure}[h]
\includegraphics[width=8.5cm]{charge-Bind-frac-diff.pdf}
\caption{\label{fig:charge-Bind-frac-diff} Fractional difference in the scalar charge between Tolman VII and constant density cases in terms of the binding energy (Eq.~\eqref{eq:alphaA_frac_diff}) for two choices of $\beta_0$. Notice that the two agree within an error of $\mathcal{O}(10^{-3})$, indicating the quasi-universality of the relation between the scalar charge and binding energy.
}
\end{figure}

For spontaneous scalarization, although the analytic expressions are valid as an order-of-magnitude estimate, the agreement with numerical results are not as good as the perturbative scalarization case. This is partially because we have used the relation between the binding energy and compactness in GR in Eq.~\eqref{eq:bind_comp_Tol}. Though, due to multi-valued scalar charges for some fixed binding energy with small $\beta_0$ (e.g. $\beta_0 = -4.5$), it would be difficult to find a closed, analytic expression for $\alpha_A(\Omega_A/M_A)$.

\section{Conclusion}
\label{sec:conclusion}

We derived analytically scalar charges for NSs in a scalar-tensor theory proposed in~\cite{Damour:1993hw}. This was achieved by considering the Tolman VII energy density profile. We first worked within the weak-field approximation and then resummed the series through Pad\'e approximants, deriving scalar charge expressions for both perturbative and spontaneous scalarization. We found that the analytic scalar charges are in excellent agreement with those computed numerically from realistic EoSs, especially for the SLy EoS. We also found a quasi-universal relation between the scalar charge and binding energy. The analytic expressions derived here allow one to mathematically support the quasi-universality by comparing the Tolman VII result with the constant density one. The analytic result provides an accurate, ready-to-use, and physically-motivated expression for scalar charges.

A similar quasi-universality between the stellar sensitivities and the binding energy were found recently in Einstein-\AE ther theory~\cite{Gupta:2021vdj}. The relation in the weak-field limit was first derived in~\cite{Foster:2007gr}. This was based on the result for weakly-gravitating stars in~\cite{Foster:2006az}, where the dipole moment for the vector field (which depends on the sensitivities) was derived within the parameterized post-Newtonian framework and all the internal structure dependence was found to be encoded in the binding energy. It would be interesting to study if a similar analysis can be carried out in scalar-tensor theories to explain the quasi-universal relation found here.

Various avenues exist for other possible future work. For example, one obvious extension is to apply the analysis presented here to other scalar-tensor theories with spontaneous scalarization, such as the ones with the mass potential~\cite{Chen:2015zmx,Ramazanoglu:2016kul,Yazadjiev:2016pcb,Doneva:2016xmf} or quartic interaction~\cite{Staykov:2018hhc}, scalar-Gauss-Bonnet gravity~\cite{Doneva:2017bvd,Silva:2017uqg,Doneva:2017duq,Silva:2018qhn,Minamitsuji:2018xde,Macedo:2019sem,Dima:2020yac,Doneva:2020nbb,Herdeiro:2020wei,Silva:2020omi}, or Horndeski theories~\cite{Andreou:2019ikc}. One may also apply the calculation to spontaneous vectorization~\cite{Ramazanoglu:2017xbl,Annulli:2019fzq,Kase:2020yhw}, tensorization~\cite{Ramazanoglu:2019gbz}, or spinorization~\cite{Ramazanoglu:2018hwk}. Another possibility might be to consider an analytic representation of the dynamical/induced scalarization~\cite{Barausse:2012da,Palenzuela:2013hsa,Shibata:2013pra,Sampson:2014qqa,Sennett:2016rwa,Khalil:2019wyy} in compact binary mergers. Finally, the parameter region of the coupling parameter $\beta_0$ that gives rise to spontaneous scalarization in the scalar-tensor theories studied here has been shown to be inconsistent with solar system experiments if one includes cosmological evolution of the scalar field~\cite{Damour:1992kf,Damour:1993id,Sampson:2014qqa,Anderson:2016aoi,Anderson:2017phb}. It would be interesting to consider analytic scalar charges in scalar-tensor theories that has a consistent cosmological evolution of the scalar field~\cite{Mendes:2014ufa,Mendes:2016fby,Antoniou:2020nax}.

\acknowledgments
K.Y. acknowledges support from NSF Grant PHY-1806776, NASA Grant 80NSSC20K0523, a Sloan Foundation Research Fellowship and the Owens Family Foundation.
K.Y. would like to also acknowledge support by the COST Action GWverse CA16104 and JSPS KAKENHI Grants No. JP17H06358.