# A universal model for the evolution of tidally stripped systems

Nicole E. Drakos1, James E. Taylor2,3, and Andrew J. Benson4
1Department of Astronomy and Astrophysics, University of California, Santa
Cruz, 1156 High Street, Santa Cruz, CA 95064, USA
2Department of Physics and Astronomy, University of Waterloo, 200 University
Avenue West, Waterloo, ON N2L 3G1, Canada
3Waterloo Centre for Astrophysics, University of Waterloo, 200 University
Avenue West, Waterloo, ON N2L 3G1, Canada
4Carnegie Observatories, 813 Santa Barbara Street, Pasadena, CA 91101, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Accurate models of the structural evolution of dark matter subhaloes, as they
orbit within larger systems, are fundamental to understanding the detailed
distribution of dark matter at the present day. Numerical simulations of
subhalo evolution support the idea that the mass loss associated with tidal
stripping is most naturally understood in energy space, with the particles
that are the least bound being removed first. Starting from this premise, we
recently proposed a zero-parameter “energy-truncation model" for subhalo
evolution. We tested this model with simulations of tidal stripping of
satellites with initial NFW profiles, and showed that the energy-truncation
model accurately predicts both the mass loss and density profiles. In this
work, we apply the model to a variety of Hernquist, Einasto and King profiles.
We show that it matches the simulation results quite closely in all cases,
indicating that it may serve as a universal model to describe tidally stripped
collisionless systems. A key prediction of the energy-truncation model is that
the central density of dark matter subhaloes is conserved as they lose mass;
this has important implications for dark matter annihilation calculations, and
for other observational tests of dark matter.

###### keywords:

dark matter – galaxies: haloes – methods: numerical

††pubyear: 2021††pagerange: A universal model for the evolution of tidally
stripped systems–B

##  1 Introduction

In the $\Lambda$-Cold Dark Matter ($\Lambda$CDM) framework, positive density
fluctuations present in the matter distribution at early times gradually break
away from the cosmological expansion and collapse into roughly spherical dark
matter haloes. These haloes merge hierarchically on progressively larger
scales as the Universe expands, giving rise to the galaxy, group and cluster-
scale haloes known today. Halo mergers are relatively ineffective at mixing
infalling material, and the central parts of smaller haloes merging into a
larger system can survive for many orbits as self-bound subcomponents, or
‘subhaloes’.

It is critical to understand precisely how subhaloes evolve after a merger,
since a number of important tests of dark matter, including dark matter
annihilation signals (e.g. Stref et al., 2019; Delos, 2019), gravitational
lensing (e.g. Limousin et al., 2005; Baltz et al., 2009; Sereno et al., 2016),
and the imprint of substructure on stellar streams (e.g. Carlberg, 2020;
Bonaca et al., 2020) depend sensitively on the properties of the subhalo
population. The exact internal structure and abundance of small subhaloes is
highly uncertain, however; in particular, there are concerns that substructure
seen in cosmological simulations may still be significantly affected by
artificial (numerical) disruption (e.g. van den Bosch et al., 2018; Errani &
Navarro, 2021). These uncertainties, together with baryonic effects,
contribute to the small-scale structure problems often cited as one of the
main challenges to the $\Lambda$CDM model (see, e.g. Bullock & Boylan-Kolchin,
2017, for a review of these problems).

Given the complexity of cosmological structure formation and the relatively
poor resolution of large-volume simulations, subhalo evolution has often been
studied using simplified simulations of a single (sub)halo evolving in a fixed
background potential (e.g. Hayashi et al., 2003; Kazantzidis et al., 2004;
Boylan-Kolchin & Ma, 2007; Kampakoglou & Benson, 2007; Peñarrubia et al.,
2008a, b, 2009; Choi et al., 2009; Peñarrubia et al., 2010; Drakos et al.,
2017; Ogiya et al., 2019; Delos, 2019; Errani & Navarro, 2021). These studies
have produced a number of models for tidal stripping and structural evolution
that describe how mass is lost as a function of radius. These models generally
remain empirical, however, i.e. they contain free parameters or functional
forms that are adjusted to match specific simulation results. Thus, while
accurate for the specific cases considered, these models are restricted to the
particular systems simulated (which are typically isotropic systems with NFW
profiles; e.g. Hayashi et al., 2003; Green & van den Bosch, 2019).

As first pointed out by Choi et al. (2009), tidal stripping of subhaloes is
well approximated as a monotonic, outside-in process in energy space.
Extending this idea, in Drakos et al. (2017)—Paper I hereafter—we showed that
truncation of the particle distribution function (DF) in a dark matter halo,
via a ‘lowering’ operation, produces a set of density profiles very similar to
those measured for tidally stripped systems in numerical simulations. This
approach has several advantages over models based on simulation data alone, as
it predicts the evolution of the central regions of haloes below the
resolution limit of numerical simulations, and it is easily generalized to
other profiles, or to anisotropic systems. It also has a clear physical
interpretation, namely that as mass is lost from a subhalo, the energies
change more-or-less uniformly within the system, pushing particles close to
the threshold for unbinding over this limit, and thus removing them from the
system. In subsequent work (Drakos et al., 2020, Paper II hereafter), we
developed a full model for tidal stripping that combines this lowering
approach with a simple estimate of the mass loss rate. The model has no free
parameters, and agrees well with simulations of mass loss from systems with
NFW profiles. Since then, a number of other authors have successfully studied
the evolution of tidally stripped systems using an energy-based approach (e.g
Stücker et al., 2021; Amorisco, 2021; Errani et al., 2021).

The aim of the current study is to test whether the energy-truncation approach
and the mass-loss model presented in Papers I and II are valid for a broader
range of tidally stripped collisionless systems. The structure of this paper
is as follows: first, in Section 2 we review the energy-based truncation
method introduced in Paper I and the mass-loss model from Paper II. In
Sections 3 and Section 4 we summarize the collisionless systems we consider in
this paper, and the idealized $N$-body simulations we use to test our models,
respectively. In Section 5 we calculate from energy truncation the expected
evolution of the density profiles of these systems, and compare these
predictions to the simulations in Section 6. We show our model predictions for
the central density of subhaloes in Section 7. Finally, we explore the
physical interpretation of this model in Section 8, and discuss our results in
Section 9. In a companion paper, we explore the implications of this work for
dark matter annihilation and galaxy-galaxy lensing signal predictions (Drakos
et al., 2022, in prep).

##  2 Review of energy-truncation model

In this section we briefly review the energy-based description of tidal
truncation introduced in Paper I, and the mass-loss model developed in Paper
II.

###  2.1 Review of distribution functions of isolated spherical systems

Systems of particles can generically be described by a distribution function
(DF) $f(\mathbf{r},t)={\rm d}m/{\rm d}r^{3}{\rm d}v^{3}$ which specifies the
mass per unit volume ${\rm d}r^{3}{\rm d}v^{3}$ at a given location in phase
space $(\mathbf{r},t)$. For isolated, spherically symmetric and isotropic
systems, the DF can be written as a function of a single variable
$f(\mathbf{r},t)=f(r,v)=f(\mathcal{E})$, where $\mathcal{E}=\Psi(r)-v^{2}/2$
is the (conserved) relative energy and $\Psi(r)$ is the relative potential,
defined as $\Psi(r)=-\Phi(r)+\Phi_{0}$. Here $\Phi(r)$ is the usual
gravitational potential, while $\Phi_{0}$ is a reference potential, usually
taken to be the value of $\Phi$ at the outer boundary of the system. Given the
sign convention in the definition of the relative energy, it is positive for
all bound particles and represents the binding energy needed to remove the
particle from the self-bound system. With these definitions,
$f(\mathcal{E})>0$ when $\mathcal{E}>0$, and $f(\mathcal{E})=0$ otherwise.

In terms of these quantities, we can calculate the density profile $\rho(r)$
corresponding to a given distribution function

$\displaystyle\rho(r)$
$\displaystyle=4\pi\int_{0}^{\Psi(r)}f(\mathcal{E})\sqrt{2(\Psi(r)-\mathcal{E})}{\rm
d}\mathcal{E}$ (1)

or we can invert this relationship to determine the isotropic distribution
function corresponding to a given density profile

$\displaystyle f(\mathcal{E})$
$\displaystyle=\dfrac{1}{\sqrt{8}\pi^{2}}\left[\int_{0}^{\mathcal{E}}\dfrac{1}{\sqrt{\mathcal{E}-\Psi}}\dfrac{{\rm
d}^{2}\rho}{{\rm d}\Psi^{2}}{\rm
d}\Psi+\dfrac{1}{\sqrt{\mathcal{E}}}\left(\dfrac{{\rm d}\rho}{{\rm
d}\Psi}\right)_{\Psi=0}\right]\,\,\,,$ (2)

where $\Phi(r)$ (and thus $\Psi=\Psi(r)$) is calculated from $\rho(r)$, using
Poisson’s equation $\nabla^{2}\Phi(r)=4\pi G\rho(r)$ (see Binney & Tremaine,
1987, for the derivation of these results).

###  2.2 Lowering the DF to represent tidal truncation

Subhaloes are clearly not isolated throughout their evolution, and the
relative energy defined above will vary with time for each individual
particle, due to changes in the self-bound potential, and the effects of tidal
heating. The approach of Paper I was to assume that over the course of a full
orbit, the relative energies of all particles change by a constant ‘tidal
energy’ $\mathcal{E}_{T}$. Considering the system near apocentre, where
heating is minimal, it can then be treated as approximately isolated, with
$f\sim f(\mathcal{E})$, where the potential used to calculate $\mathcal{E}$ is
that of the particles remaining bound at that time, and all the energies have
been shifted by $\mathcal{E}_{T}$ with respect to the previous apocentre. By
that point, tidal stripping will have removed from the system all the
particles in a range of $\mathcal{E}$ between zero and $\mathcal{E}_{T}$.

To estimate the resulting adjustment of the distribution function, we ‘lower’
it by a tidal energy $\mathcal{E}_{T}$, and use this modified form to recover
the new, tidally stripped density profile. This method is similar to how the
well-known King model for truncated stellar systems was derived by lowering
the DF of an isothermal sphere (King, 1966). It was originally proposed in
Widrow & Dubinski (2005) as a way to truncate NFW profiles for use as initial
conditions (ICs) in isolated simulations.

In terms of the original distribution function, the lowered version is defined
as

$f_{\mathcal{E}_{T}}(\mathcal{E})=\begin{cases}f_{0}(\mathcal{E}+\mathcal{E}_{\rm
T})-f_{0}(\mathcal{E}_{T})&\mathcal{E}\geq 0\\\ 0&\mathcal{E}\leq
0\,\,\,,\end{cases}$ (3)

where $f_{0}(\mathcal{E})$ is the DF of the original system, and
$\mathcal{E}_{T}$ is the truncation or tidal energy.

Given the lowered distribution function, we could in principle calculate the
density profile from Equation (1), but we would need the relative potential
$\Psi(r)$, which itself depends on $\rho(r)$ through Poisson’s equation.
Substituting the spherical form of Poisson’s equation on the left-hand side of
Equation (1), we can eliminate $\rho$ and write:

$\displaystyle\dfrac{{\rm d}^{2}{\Psi(r)}}{{\rm d}r^{2}}+$
$\displaystyle\dfrac{2}{r}\dfrac{{\rm d}\Psi(r)}{{\rm d}r}$ (4)
$\displaystyle=-16\pi^{2}G\int_{0}^{\Psi(r)}f_{\mathcal{E}_{T}}(\mathcal{E})\sqrt{2(\Psi(r)-\mathcal{E})}{\rm
d}\mathcal{E}\,.$

This differential equation, together with the boundary conditions
$\Psi(0)=\Psi_{0}(0)-\mathcal{E}_{T},\ {\rm d}\Psi(0)/{\rm d}r=0$ (where
$\Psi_{0}(r)$ is the relative potential of the original, un-truncated system)
is easily solved by conventional techniques (see Paper I). Once we have solved
for $\Psi(r)$, the truncation radius $r_{t}$ is given by the condition
$\Psi(r_{t})=0$, and the truncated density profile $\rho(r)$ can be calculated
from Poisson’s equation.

###  2.3 Predicting the mass-loss rate

As explained in Paper II, our proposed mass-loss model is based on the Jacobi
model for tidal truncation on a circular orbit (Binney & Tremaine, 1987;
Taylor & Babul, 2001). A tidal or limiting radius $r_{\lim}$ is defined such
that:

$\displaystyle\bar{\rho}_{\rm sat}(r_{\lim})$ $\displaystyle=\eta_{\rm
eff}\bar{\rho}_{\rm H}(R_{p})$ (5)

where $\bar{\rho}_{\rm sat}(r_{\rm lim})$ is the mean density of the satellite
interior to the tidal radius and $\bar{\rho}_{\rm H}(R_{p})$ is the mean
density of the host halo interior to the pericentre of the satellite’s orbit.
Throughout this paper, we use $R$ when referring to the radius with respect to
the centre of the host halo, and $r$ when referring to the radius within the
satellite.

The tidal radius, $r_{\rm lim}$, and resulting bound mass of the system,
$M_{\rm bnd}$, are then given by the following system of equations:

$\displaystyle M_{\rm bnd}$
$\displaystyle=\dfrac{4}{3}r_{\lim}^{3}\bar{\rho}_{\rm sat}(r_{\lim})$ (6)
$\displaystyle M_{\rm bnd}$ $\displaystyle=\eta_{\rm
eff}M_{H}(<R_{p})\dfrac{r_{\rm lim}^{3}}{R_{p}^{3}}$

We emphasize that the tidal radius $r_{\rm lim}$ is not the same as the
truncation radius $r_{t}$ defined in the previous section.

The constant, $\eta_{\rm eff}$ is defined as the orbital average of the
instantaneous, spherical $\eta$ value defined in King (1962):

$\eta_{\rm eff}=\dfrac{1}{t_{\rm orb}}\int_{0}^{t_{\rm
orb}}\left(\dfrac{\omega^{2}}{\omega_{c}^{2}}-\dfrac{1}{\omega_{c}^{2}}\dfrac{{\rm
d}^{2}\Phi_{H}}{{\rm d}R^{2}}\right){\rm d}t\,\,\,,$ (7)

where $\omega=|\mathbf{V}\times\mathbf{R}|/R^{2}$ is the instantaneous angular
velocity of the satellite, $\omega_{c}=GM_{H}(<R)/R^{3}$ is the angular
velocity of a circular orbit, and $\Phi_{H}(R)$ is the potential of the host.
For a spherically symmetric system, ${\rm d}\Phi_{H}/{\rm
d}R=GM_{H}(<R)/R^{2}$ As demonstrated in Paper II, with this definition of
$\eta$, we obtain a good match to the mass-loss rates measured in our
simulations without needing to add or adjust any free parameters. Recently,
Stücker et al. (2021) proposed that a more natural way to truncate tidally
stripped profiles in energy space is using the “boosted potential". We
consider $\eta$ determined from the boosted potential in Appendix B, and
compare it to the $\eta$ values found from the energy-truncation model.

Given a tidal radius defined as in Equation (5), we assume the mass outside
this radius is lost over the course of an orbit. We adjust the tidal energy
$\mathcal{E}_{T}$ in Equation (3) until it produces a satellite with this
lowered mass. This gives a full model for the stripped system, allowing us to
specify the total mass, distribution function and density profile by the end
of the orbit. We apply this stripping model at successive apocentres, since
this is the point at which the profile is expected to reach equilibrium. As
the system continues to orbit, each successive pericentric passage will
decrease $r_{\rm lim}$ and increase $\mathcal{E}_{T}$.111We take advantage of
the fact that truncating the profile first at energy $\mathcal{E}_{1}$, then
truncating the new truncated profile a second time at $\mathcal{E}_{2}$ is
mathematically equivalent to truncating the original profile at
$\mathcal{E}_{1}+\mathcal{E}_{2}$. Therefore, though $r_{\rm lim}$ is
calculated using the most recent, stripped density profile, the new profile is
always calculated by lowering the original, untrucated DF.

###  2.4 Algorithmic description of method

In practice, given an orbit and a set of initial conditions, the energy-
truncation model can be performed using the following steps:

  1. 1.

Compute $\eta_{\rm eff}$ by integrating over the orbit using Equation (7).

  2. 2.

Solve for the bound mass $M_{\rm bnd}$ and tidal radius, $r_{\rm lim}$, using
Equation (6). Note that $r_{\rm lim}$ is not needed in the following steps.

  3. 3.

Compute the tidal energy that corresponds to $M_{\rm bnd}$. There is a
monotonic relationship between $M_{\rm bnd}$, $r_{t}$ and $\mathcal{E}_{T}$,
and any one of these three variables can be used to uniquely describe the
truncated profile. We create a lookup table to map between $M_{\rm bnd}$ to
$\mathcal{E}_{T}$222For convenience, Paper I provides an empirical
relationship to map between these variables..

  4. 4.

Using $\mathcal{E}_{T}$, compute the lowered DF in Equation (3).

  5. 5.

Calculate $\Psi(r)$ by solving the ODE in Equation (4).

  6. 6.

Compute the corresponding truncated profile $\rho_{\rm sat}$ using Equation
(1).

  7. 7.

Go to ii, and repeat for as many orbits as desired.

##  3 Collisionless satellite models

The main goal of this paper is to examine how well sharp truncation of the
distribution function at a tidal energy, and the mass-loss model reviewed in
the previous section, describe the evolution of a range of different
collisionless systems333A collisionless system is one in which interactions
between individual particles are negligible, and thus the gravitational force
acting on each particle can be treated as a smooth density field rather than a
collection of individual particles., including common models for CDM haloes
such as the NFW, Hernquist, Einasto models, but also models for stellar
systems such as the King model. In this section we review the properties of
each of these model systems.

###  3.1 NFW

In Papers I and II, we restricted our attention to the NFW profile,

$\rho(r)=\dfrac{\rho_{0}r_{\rm s}^{3}}{r(r+r_{\rm s})^{2}}\,\,\,,$ (8)

where $\rho_{0}$ is a characteristic density and $r_{\rm s}$ is the scale
radius, describing the point where the logarithmic slope is ${\rm
d}\log\rho/{\rm d}\log r=-2$.

For convenience, we provide an analytic description of an energy-truncated NFW
profile, using the functional form

$\rho_{\rm NFWT}=\dfrac{e^{-x/y}}{\left[1+(x/y)\right]^{a}}\rho_{\rm NFW}\,,$
(9)

where $x=r/r_{e}$ is the radius normalized by an effective radius, $r_{e}$,
and $y=M_{\rm bnd}/M_{\rm NFW}(<r_{s})$ is the bound mass of the truncated
profile, normalized by the mass of the untruncated profile within the NFW
scale radius, $r_{s}$. This functional form was chosen as we found it
empirically yielded a good fit, and conserves the central density as
$r\rightarrow 0$.

The parameters $r_{e}$ and $a$ can be expressed as follows:

$\displaystyle\log_{10}r_{e}$
$\displaystyle=0.0811y^{3}+0.358y^{2}+0.0781y-0.201$ (10) $\displaystyle a$
$\displaystyle=0.179y^{3}+0.379y^{2}-0.524y-0.952\,\,\,.$

As shown in Fig. 1, this analytic fit generally agrees with the energy-
truncated model to within 5 per cent, except at large radii ($r\gtrsim
0.5r_{t}$), where the density is very low. Our model predicts a conserved
central density, which is captured in Equation (9); this differs from
empirical models calibrated to simulations, which typically include an
explicit drop in the central density (e.g. Hayashi et al., 2003; Green & van
den Bosch, 2019). This point will be discussed further in Section 7.

Figure 1: Top: Energy-truncated NFW profiles, $\rho$ (solid lines) and the
analytic fit, $\rho_{\rm fit}$ (dotted lines). Bottom: Residuals in the fits,
measured as ($\rho_{\rm fit}/\rho-1)$. The 1 per cent errors are shown with
dashed black lines. The analytic fit agrees well with the energy-truncated
model inside $\sim 0.5r_{t}$.

###  3.2 Hernquist

The Hernquist profile (Hernquist, 1990) was originally used to describe
spherical galaxies, but it is also a reasonable approximation for cosmological
dark matter halo profiles. The density profile is given by:

$\rho(r)=\dfrac{M_{\rm tot}}{2\pi}\dfrac{a}{r(r+a)^{3}}\,\,\,,$ (11)

where $M_{\rm tot}$ is the total mass and $a$ is a characteristic radius, that
encloses a mass of $M_{\rm tot}/4$.

The advantage to this model is that it has simple analytic expressions for
many of its properties, including its DF:

$\displaystyle f(\mathcal{E})$
$\displaystyle=\dfrac{M}{8\sqrt{2}\pi^{3}a^{3}v_{g}^{3}}\dfrac{1}{(1-q^{2})^{5/2}}$
(12a) $\displaystyle\hskip
14.22636pt\times\left(3\sin^{-1}q+q(1-q^{2})^{1/2}(1-2q^{2})(8q^{4}-8q^{2}-3)\right)$
$\displaystyle q$ $\displaystyle=\sqrt{\dfrac{a\mathcal{E}}{GM}}$ (12b)
$\displaystyle v_{g}$ $\displaystyle=\left(\dfrac{GM}{a}\right)^{1/2}\,\,\,.$
(12c)

###  3.3 Einasto

The Einasto profile was first used to describe star counts in the Milky Way
(Einasto, 1965). However, this profile is often a better description of
cosmological dark matter halo profiles than the well-known NFW profile (e.g.
Navarro et al., 2004; Gao et al., 2008; Klypin et al., 2016). The Einasto
profile has the following form:

$\displaystyle\rho(r)$
$\displaystyle=\rho_{-2}\exp\left(-\dfrac{2}{\alpha}\left[\left(\dfrac{r}{r_{-2}}\right)^{\alpha}-1\right]\right)\,\,\,,$
(13)

where $\alpha$ is the Einasto shape parameter and $r_{-2}$ is the radius where
the logarithmic slope is $-2$. Compared to the NFW profile, the Einasto
profile has an extra parameter, $\alpha$, that controls the inner slope of the
density profile, and may reflect the mass accretion history of the halo (e.g.
Klypin et al., 2016).

###  3.4 King

The King model resembles an isothermal sphere at small radii, but has a finite
mass within a well-defined tidal radius. This model is typically used to
describe truncated stellar systems such as globular clusters or elliptical
galaxies (King, 1966). The King model is derived by lowering the DF of an
isothermal sphere:

$f(\mathcal{E})=\rho_{1}(2\pi\sigma^{2})^{-3/2}(e^{\mathcal{E}/\sigma^{2}}-1)\,\,\,,$
(14)

where $\sigma$ is the velocity dispersion, $\rho_{1}$ is a characteristic
density.

The density profile of the King model can then be calculated from the DF using
Equation (1), which gives:

$\rho(\Psi)=\rho_{1}\left[e^{\Psi/\sigma^{2}}\mbox{erf}\left(\dfrac{\sqrt{\Psi}}{\sigma}\right)-\sqrt{\dfrac{4\Psi}{\pi\sigma^{2}}}\left(1+\dfrac{2\Psi}{3\sigma^{2}}\right)\right]\,\,\,.$
(15)

To relate the relative potential energy, $\Psi(r)$, to the density of the
profile, $\Psi(r)$ can be solved numerically using Poisson’s equation
(Equation (4)). There are many possible parameterizations of the King model,
but it can be uniquely defined by the total mass, tidal radius, $r_{\rm t}$
and a dimensionless central potential $P_{0}=\Psi(0)/\sigma^{2}$. Alternately,
King models can also be characterized by a total mass, tidal radius and a
concentration parameter, $c_{K}$, which depends on the ‘King radius’ $r_{0}$.
The latter quantities are defined as:

$\displaystyle c_{K}$ $\displaystyle=\log_{10}\left(\dfrac{r_{\rm
t}}{r_{0}}\right)$ (16) $\displaystyle r_{0}$
$\displaystyle=\sqrt{\dfrac{9\sigma^{2}}{4\pi G\rho_{0}}}\,\,\,.$

where $\rho_{0}$ is the central density of the halo.

##  4 Simulations

The simulations were performed using a version of the $N$-body code gadget-2
(Springel, 2005), modified to contain a fixed background potential
corresponding to the host halo.

###  4.1 Initial halo models

We considered four different satellite models for our initial conditions
(ICs), as described in the previous section. The ICs were created using our
public code Icicle (Drakos et al., 2017). For all models, we define the mass
unit to be the initial mass of the satellite, $M_{\rm unit}=M_{\rm sat}$. We
define the length unit to be $r_{\rm unit}=a$ for the Hernquist model, $r_{\rm
unit}=r_{-2}$ for the Einasto models, and $r_{\rm unit}=0.1\,r_{\rm t}$ for
the King models. Then the density, time and energy units are $\rho_{\rm
unit}=M_{\rm unit}r_{\rm unit}^{-3}$, $t_{\rm unit}=\sqrt{r_{\rm
unit}^{3}/GM_{\rm unit}}$ and $E_{\rm unit}=GM_{\rm unit}/r_{\rm unit}$,
respectively. We calculated the softening length for each profile as
$\epsilon=0.5r_{h}N^{-1/3}$,where $r_{h}$ is the radius enclosing half of the
total mass, as in van Kampen (2000).

For the NFW profile, we used the same set of simulations described in Paper I.
Since NFW profiles are infinitely extended with divergent mass as $r$ goes to
infinity, these ICs were truncated as described in Paper I. The resulting
models resemble a truncated NFW profile with an energy truncation of
$\mathcal{E}_{\rm T}\approx 0.27$. Though the simulation results can be scaled
by the length and mass units, for the Einasto models, there is one free
parameter, $\alpha$, to fix; we chose values of $\alpha=0.15$ and $0.3$, as
these are representative of the range found in simulations (Gao et al., 2008).
The King models also have an additional free parameter (which could be
specified as $r_{0}$, $P_{0}$ or $c$); we used $P_{0}=3$ and $P_{0}=12$, which
correspond to King concentrations of $c_{K}=0.7$ and $c_{K}=2.7$. These
concentrations are typical of globular clusters and elliptical galaxies,
respectively (Binney & Tremaine, 1987). A summary of the IC properties are
given in Table 1. To check the stability of the ICs, they were evolved for
$t=1000\,t_{\rm unit}$ in isolation using the N-body code gadget-2 (Springel,
2005), as shown in Fig. 2. All four profiles are extremely stable outside
$r_{\rm unit}=0.1$ at $t=1000\,t_{\rm unit}$.

Table 1: Summary of the parameters used for the ICs; all satellites have an initial mass of $M_{\rm sat}=1$. The columns list (1) the name of the ICs, (2) the number of particles, (3) the radial unit, (4) specified profile parameters, (5) derived profile parameters, and (6) Gadget-2 softening length, roughly equivalent to the Plummer softening length.  Profile Name | $N$ | $r_{\rm unit}$ | IC Parameters | Derived Parameters | $\epsilon/r_{\rm unit}$
---|---|---|---|---|---
NFWT | $\approx 1.3\times 10^{6}$ | $r_{s}$ | $r_{\rm cut}=10\,r_{\rm unit}$ | $\rho_{0}=0.08\,\rho_{\rm unit}$ | $0.01$
Hern | $10^{6}$ | $a$ | —– | —– | $0.01$
EinHigh | $10^{6}$ | $r_{-2}$ | $\alpha_{E}=0.3$ | $\rho_{-2}=0.01\,\rho_{\rm unit}$ | $0.02$
EinLow | $10^{6}$ | $r_{-2}$ | $\alpha_{E}=0.15$ | $\rho_{-2}=0.005\,\rho_{\rm unit}$ | $0.07$
KingHigh | $10^{6}$ | $r_{\rm t}/10$ | $P_{0}=12$ |  $\rho_{1}=0.003\,\rho_{\rm unit}$, $r_{0}=0.02\,r_{\rm unit}$, $c_{K}=0.7$ | $0.01$
KingLow | $10^{6}$ | $r_{\rm t}/10$ | $P_{0}=3$ |  $\rho_{1}=0.001\,\rho_{\rm unit}$, $r_{0}=2.13\,r_{\rm unit}$, $c_{K}=2.7$ | $0.01$
Figure 2: Stability of ICs. Solid lines are the profile at $t=0\,t_{\rm
unit}$; points are the profile at $t=1000\,t_{\rm unit}$.

###  4.2 Satellite orbits

Each satellite model was evolved on two different orbits, corresponding to the
‘Fast’ and ‘Slow’ Simulations in Paper II (Simulations 3 and 4 from Paper I),
which are representative of orbits that lose mass quickly and slowly,
respectively. In Papers I and II we found that the sharp energy truncation
model is a better descriptor of orbits which are losing mass slowly. The
infinitely extended host halo was assumed to have an NFW profile for all
simulations, with a mass $M_{\rm host}=300M_{\rm unit}$ and a scale radius of
$6.69\,r_{\rm unit}$. The orbits in the Slow and Fast Simulations have an
apocentre of $r_{a}=300\,r_{\rm unit}$, and pericentres of $R_{p}=50\,r_{\rm
unit}$ and $10\,r_{\rm unit}$, respectively. The simulations are summarized in
Table 2.

Table 2:  Summary of orbital parameters for the Fast and Slow Simulations.
Columns give (1) the simulation name (2) the apocentric distance (3) the
pericentric distance (4) the tangential velocity at apocentre (5) the
tangential velocity at pericentre (6) the (radial) orbital period (7) mass-
loss factor, $\eta_{\rm eff}$

.  Simulation Name $R_{a}/r_{\rm unit}$ $R_{p}/r_{\rm unit}$ $V_{a}/V_{\rm
unit}$ $V_{p}/v_{\rm unit}$ $t_{\rm orb}/t_{\rm unit}$ $\eta_{\rm eff}$ Slow
Simulation 100 50 1.42 2.84 185.4 2.34 Fast Simulation 100 10 0.51 5.10 129.7
1.72

Subhalo centres and bound particles were identified similarly to the method
described in Papers I and II. Considering only particles that were bound in
the previous snapshot, the centre of the satellite was defined as the densest
point in position and velocity space. The densest point was calculated using
the centre of mass (COM) in progressively smaller spheres, as originally
described in Tormen et al. (1997). Initially, the sphere was centred on the
COM of all particles, and the radius was defined as the distance to the
furthest particles. The radius was decreased by 90 per cent and recentered on
the COM, until there were fewer than 100 particles in the sphere. The velocity
frame was calculated in the same way, by using progressively smaller spheres
to find the densest point in velocity space.

Once the frame of the satellite was defined, we found the self-bound
particles. First, the energy of each particle was calculated in this frame,
assuming a spherical potential (in Paper II we showed the results are
insensitive to this assumption):

$P_{i}\approx-
Gm\left(\dfrac{N(<r_{i})}{r_{i}}+\sum_{j,r_{j}>r_{i}}^{N}\dfrac{1}{r_{j}}\right)\,\,\,,$
(17)

where $r_{i}$ is the distance of particle $i$ from the center of the system.
Particles were iteratively removed, and energies recalculated, until
convergence. The bound satellite mass is the mass of the remaining bound
particles.

##  5 Predicted profile evolution

In the energy truncation approach described in Section 2.1, the evolution of
the stripped system depends on a single parameter, which can be expressed as
the tidal energy $\mathcal{E}_{T}$. First, we explore how the NFW, Hernquist,
Einasto and King models evolve when truncated as function of
$\mathcal{E}_{T}$. We note that applying energy truncation to a King model
simply results in another King model, with parameters modified as follows:

$\displaystyle P_{0,T}$ $\displaystyle=P_{0}-E_{T}$ (18)
$\displaystyle\rho_{1,T}$ $\displaystyle=\rho_{1}\exp(E_{T})$ $\displaystyle
r_{0,T}$
$\displaystyle=r_{0}\sqrt{\dfrac{\rho(P_{0})}{\rho_{T}(P_{0,T})}}\,\,\,.$

For the other profiles, the truncated version no longer has a simple analytic
description, but can be calculated from the lowered distribution function. For
convenience, we provided an analytic fit to our lowered NFW profile in Section
3.1.

Though the model presented in Section 2.1 most naturally expresses the lowered
density profiles in terms of a truncation energy, $\mathcal{E}_{T}$, the model
can alternatively be parameterized by either the truncation radius, $r_{t}$,
or the total bound mass of the satellite, $M_{\rm bnd}$. In Fig. 3 we show the
relationship between these three parameters. Generally, there is a monotonic
relationship between the bound mass and the tidal energy or the tidal radius.
An exception is the tidal radius for the KingLow models; interestingly in this
case $r_{\rm t}$ _increases_ for large values of
$\mathcal{E}_{\rm{}_{T}}/\Psi(0)$. We will show in Section 5 that this
prediction is also consistent with the simulation results.

Figure 3: Relationship between the tidal energy $\mathcal{E}_{T}$, the bound
mass of the satellite (top), and the tidal radius (bottom) for different
initial satellite profiles. The relative energy, $\mathcal{E}$, has been
normalized by the relative potential of the untruncated profile at $r=0$,
$\Psi(0)$. Note that since the NFW profile has a divergent mass profile, we
have chosen $M_{\rm unit}$ to correspond to an NFW profile truncated at
$\mathcal{E}_{T}=0.27$; therefore, $M_{\rm bnd}/M_{\rm unit}$ can be greater
than 1 (dotted line). In general, the bound mass and truncation radius
decrease with tidal energy, with the exception of the KingLow profile. The
KingLow profile is the only cored profile, and its truncation radius increases
as it is tidally stripped.

In Fig. 4 we show the truncated density profiles and various related
quantities, coloured by the tidal energy, $\mathcal{E}_{\rm{}_{T}}$ normalized
by the relative potential of the untruncated profile at $r=0$, $\Psi(0)$. For
most of the halo models, as the satellites are stripped they decrease in
density at all radii, though most of this decrease is at large radii. An
exception is the KingLow models, where $r_{t}$ increases, with a corresponding
increase in density at larger radii. We note that when the density profile is
cored (as in the centre of the King models), the central density decreases
significantly as mass is lost, while for the cuspier profiles the central
density is conserved. This is consistent with results from Peñarrubia et al.
(2010).

We also find that the scale radius of the NFW, Hernquist and Einasto profiles
(the point where the logarithmic slope of the density profile is $-2$,
estimated from the location of the maximum of $\rho r^{2}$) decreases as the
halo becomes tidally stripped. In principle, this could lead to an increase in
the concentration parameter, though the virial radius is also decreasing; we
examine halo concentrations further in a companion paper (Drakos et al., 2022,
in prep). Finally, the location and amplitude of the peak of the circular
velocity curve both decrease with increasing $\mathcal{E}_{T}/\Psi(0)$ for all
profiles except KingLow.

Figure 4: Predicted evolution for different initial density profiles
(columns), assuming sharp truncation of the distribution function in energy.
The EinLow and EinHigh profiles have $\alpha$ values of $\alpha=0.15$ and
$0.3$, respectively, while KingHigh and KingLow profiles have central
potentials of $P_{0}=12$ and $3$, respectively. The curves are coloured by the
value of the tidal energy, $\mathcal{E}_{T}$, normalized by the relative
potential of the untruncated profile at $r=0$, $\Psi(0)$. The density profile
goes to zero at the truncation radius, $r_{t}$. The relation between the tidal
energy, truncation radius and bound mass is given in Fig. 3.

##  6 Model comparison to simulations

We have already demonstrated in Papers I and II that the energy truncation
approach provides an excellent description of the evolution of systems with
NFW profiles. The purpose of this section is to test whether this approach,
and the mass loss model developed in Paper II, works as well for the other
satellite models described in Section 3. We start by comparing the predicted
and observed mass-loss rates, before returning to the profile evolution.

###  6.1 Mass-loss rates

In Fig. 5 we compare the bound mass of the simulated satellites versus time,
calculated as described in Section 4.2, to the predictions of the mass loss
model summarized in Section 2.3. We assume the mass loss predicted by the
model on each orbit is removed by the time the system reaches apocentre. With
this convention for the Slow Simulation (left column), we get excellent
agreement (the bound mass is predicted to 5 per cent or better accuracy over
the first 5 orbits) for the NFW, Hern and EinLow profiles and good agreement
(10 per cent or better accuracy over the first 5 orbits) for the EinHigh and
KingHigh profiles. We get less accurate predictions ($\leq$20 per cent
accuracy over the first 5 orbits) for the KingLow models. In all cases, mass
loss is slightly faster in the simulations than predicted by the model.

For the Fast Simulations (right column), the observed mass-loss rates are
slightly slower than predicted, particularly on the first few orbits, where
the predictions for most models overestimate mass loss by $15$–$20$ per cent.
After a few orbits, however, the simulations have caught up with the model
predictions, and we once again get good agreement ($\leq$10 per cent accuracy)
for all models. Additionally, the model correctly predicts the complete
disruption of the KingLow profile by the second pericentric passage.

Figure 5: Bound satellite mass as a function of time in the Slow and Fast
Simulations, calculated as described in Section 2 (points), compared to
predictions of the mass loss from our model (grey lines). The black crosses
correspond to the model prediction at apocenter. Each row corresponds to a
different density profile, as indicated. In general, the energy-truncation
model agrees with the simulation results to within 20 per cent or better.

While our model does a reasonable job of predicting the evolution of the mass
distribution of the remnants, there are some discrepancies between the total
mass-loss rates predicted by the model and the rates measured in the
simulations. Defining which particles are bound to the simulated subhalo is
not a straightforward procedure, however, and often includes particles that
are temporarily bound but in the process of leaving the system (Peñarrubia et
al., 2009). We show evidence below that marginally bound mass in the outer
parts of the satellites may indeed affect the predicted mass-loss rates (see
e.g. Fig. 6).

The mass-loss picture is further complicated by the fact that it is often
difficult to determine whether mass loss in simulations is due to physical
processes or numerical artefacts (e.g. van den Bosch et al., 2018; van den
Bosch & Ogiya, 2018). Additionally, our mass loss model does not include
dynamical friction, which will cause the pericentric radius to decrease in
time. While dynamical friction can largely be neglected, as the host halo is
modelled as a smooth potential, there is some self dynamical friction, from
the remnant orbiting through tidal debris (e.g. Miller et al., 2020).

###  6.2 Profile evolution

In addition to predicting mass loss, our model also predicts the tidal energy
on each orbit, and thus the full structural evolution of the satellite. In
Fig. 6, we show how the density, mass and circular velocity profiles of the
simulations compare to the model predictions, for the Slow (top panels) and
the Fast (bottom panels) simulations. For the Slow Simulations, the agreement
is excellent, with the only significant discrepancy being at large radii for
the EinLow profile.

In the Fast Simulations, the agreement is also excellent at small radii, with
the exception of a slight disagreement in the circular velocity curves of the
KingHigh profiles at small radii; we suspect this is due to numerical
relaxation in the simulation. There are more significant discrepancies at
large radii for most of the profiles. The density profiles measured in the
simulations drop abruptly close to the predicted tidal radius, but then they
extend well beyond this radius with a shallower slope. As discussed in Paper
II and in Peñarrubia et al. (2009), this part of the profile includes
transitional material that is still bound but moving outwards, and will be
lost mostly on the next orbit. Refining the procedure we use to identify bound
particles in the simulation might resolve these discrepancies, and improve the
agreement between predicted and observed mass-loss rates in the case of orbits
with rapid mass loss. It may also be the case the energy-truncation fails to
accurately predict the evolution of the outskirts of the profile, due to the
equilibrium assumptions inherent in the model. Understanding the detailed
mass-loss and time scales is an interesting question, and one we leave to
future work.

Overall, our model appears to be universally valid for a wide range of density
profiles. The predictions for the mass loss rate are generally very accurate,
and especially so for the models that might describe dark matter structure
(NFW, Hernquist and Einasto), predicting the reminding bound mass to within
5-10 per cent or better in most cases over the first 5 orbits. The predictions
for the structure of the bound remnants are accurate except in the outer
regions close to the tidal radius; the latter discrepancy might be resolved by
adding a timescale for mass loss to our model, or refining our definition of
bound mass in the simulations. The success of the energy-truncation model is
particularly impressive given that it has no free parameters to adjust.

(a)  (b)

Figure 6:  Comparison of simulation results (points) and the energy-truncation
model (lines) for the Slow (top) and Fast (bottom) Simulations. Different
curves correspond to successive apocentric passages. The energy-truncation
model agrees well with simulations, except close to the truncation radius
where some additional mass appears to be temporarily bound in the simulations.
Note that for the cored KingLow profile (last column, bottom panel), our model
correctly predicts the observed drop in the central density. Residuals for the
density and mass profiles are shown in Appendix A.

##  7 Central density predictions

The main difference between the energy-truncation model and previous empirical
models of mass loss in the literature is that it predicts a conserved central
density in the tidally stripped systems, even for small bound mass fractions.
This is illustrated in Fig. 7, where we compare our NFW model predictions to
empirical models of stripped NFW systems from Hayashi et al. (2003),
Peñarrubia et al. (2010) and Green & van den Bosch (2019).444Note that the
bound mass fraction is defined in these papers as $M/M_{\rm NFW}(<r_{\rm
cut})$. In the following discussion, we convert all bound mass fractions to
$M/M_{\rm unit}$. As shown in Paper II, the recent model from Green & van den
Bosch (2019) is likely the most accurate description of tidally stripped NFW
haloes from single-subhalo simulations, while the model from Hayashi et al.
(2003) underestimates the central density.

Figure 7: Our model prediction (DTB20) for the density profile of a tidally
stripped NFW system, compared to previous empirical models of the stripped
profile, including Hayashi et al. (2003) (H03), Peñarrubia et al. (2010) (P10)
and Green & van den Bosch (2019) (G19). Each column is a different bound mass
fraction. The bottom panels show the residuals, calculated as $\rho_{\rm
model}/\rho_{x}-1$, where $\rho_{\rm model}$ is our model from Paper I, and
$\rho_{x}$ are the empirical models. Our model predicts larger central
densities and slightly sharper truncation compared to the others, particularly
at low bound mass fractions.

While the empirical models capture by design the results of the underlying
idealized simulations, it is unclear whether these simulations accurately
describe realistic cosmological situations. Empirical fits to simulations will
also fit numerical artefacts present in the results. For instance, Hayashi et
al. (2003) used an approximation in their initial conditions which caused a
drop in the central density in their simulated haloes that was numerical in
nature; this in turn is reflected in their empirical model for the stripped
profile (Kazantzidis et al., 2004). While the more recent paper by Green & van
den Bosch (2019) uses simulations carefully calibrated to understand the
effects of numerical relaxation on mass loss (Ogiya et al., 2019), it is still
unclear how the very central regions of the haloes are affected. Indeed,
separate work from Errani & Peñarrubia (2019) supports the idea that
centrally-divergent cusps are never disrupted.

Following Binney & Tremaine (1987), given the number of particles within
radius $r$, $N(<r)$, each of mass $m$, we can calculate a relaxation time
scale,

$t_{\rm rel}(r)\approx 0.1\frac{\sqrt{N(<r)}}{\ln
N(<r)}\sqrt{\frac{r^{3}}{Gm}}\,\,\,,$ (19)

which is the time it takes for a typical particle’s velocity to change by an
order of itself, and an evaporation time scale, $t_{\rm evap}(r)\approx
136\,t_{\rm rel}(r)$, which is the time for a typical particle to reach escape
speed. Since both times increase with radius, given a time, $t$, we can then
calculate relaxation and evaporation radii such that $t=t_{\rm rel}(r_{\rm
rel})$ and $t=t_{\rm evap}(r_{\rm evap})$, respectively. This will give us an
indication of the radii inside which resolution effects might dominate the
evolution. Note that in Fig. 7, we cannot determine at which radius numerical
relaxation and evaporation will dominate, since there is no temporal
information; the time for a subhalo to reach the indicated mass fraction would
depend on the specific orbit and background potential considered.

The slow mass loss case is the one in which the assumptions of the energy-
truncation model appear most valid (Paper II), such that we might expect the
central density to be conserved. The relative relaxation rate increases as the
number of particles in a system decreases, however, and we are unable to
evolve subhaloes to very low masses for the Slow Simulation and still have
enough particles to preserve the subhalo against numerical effects. Thus, we
cannot reliably measure the central density of the halo at late times in the
slow mass-loss simulations.

We can, however, use our analytic mass-loss model, as outlined in Section 2,
to predict the density profile after 10, 50, 100, 250 and 500 orbits. Fig. 8
shows this prediction for the orbit of the Slow Simulation. Additionally, we
show the profile predicted from Green & van den Bosch (2019) (using the same
enclosed mass), and the relaxation and evaporation radii (vertical lines). The
simulations used to calibrate the model in Green & van den Bosch (2019) were
run at a similar resolution to ours ($\sim 10^{6}$ particles), and therefore
should have similar resolution effects. With the exception of the density at
large radii (which depends on the details of how bound mass is calculated),
the results of Green & van den Bosch (2019) are nearly indistinguishable from
our model for $r>r_{\rm evap}$. They differ significantly from the model
predictions at smaller radii, but this is precisely where we expect relaxation
to affect the results.

Figure 8:  The evolution of the density profile, as in Fig. 7, but comparing
our model with the model of Green & van den Bosch (2019), after the number of
orbits labelled. The dotted grey lines and dashed black lines show the
relaxation radius and evaporation radius, respectively. The difference in
central density between our model and the model of Green & van den Bosch
(2019) generally lies at radii unresolved in the simulations.

In summary, while we cannot resolve the question of central density
conservation definitively with our own simulations, we have shown the previous
predictions in the innermost parts of the profile need to be treated with
caution, and that the question is not yet resolved by previous work in the
literature. In a companion paper (Drakos et al., 2022, in prep), we explore
the sensitivity of dark matter annihilation rates and galaxy lensing signals
to the assumed behaviour of the innermost part of the density profile.

##  8 Physical interpretation

In the previous sections, we demonstrated the overall accuracy of the energy-
truncation model for a broad range of profiles. In this section, will consider
the physical justification for this model in more detail. As demonstrated in
Choi et al. (2009) and Paper II, during tidal stripping, mass loss appears to
be ordered primarily by the original relative energy of the particles. The
energy-truncation model captures this trend by shifting the relative energy of
all particles by a constant amount over each orbit, and unbinding those under
some threshold binding energy. However it does not explicitly state why the
ordering in relative energy should be conserved, or equivalently, why the
change in relative energy is approximately constant.

In general, over the course of an orbit, a particle can become unbound from a
subhalo because (1) the velocity of the particle changes, e.g. due to a kick
from a rapid perturbation or, (2) the potential of the system changes, due to
a change in the external field, or due to mass loss from the subhalo itself.
When estimating the relative contributions of these effects, we have to be
clear which frame is used to calculate kinetic energies, and which mass is
included in the calculation the potential. In this section, we will always
define particle velocities $v$ and radii $r$ in the frame of the satellite,
and will we calculate the satellite potential from the self-bound mass defined
previously.

Given these conventions, we can write the two contributions to the change in
energy as

$\Delta\mathcal{E}=\Delta\mathcal{E}_{K}+\Delta\mathcal{E}_{P}\,\,\,,$ (20)

where

$\Delta\mathcal{E}_{K}=-\left[\dfrac{1}{2}(\mathbf{v_{0}}+\Delta\mathbf{v})^{2}-\dfrac{1}{2}\mathbf{v}^{2}\right]=-\mathbf{v}\cdot\Delta\mathbf{v}-\dfrac{1}{2}(\Delta\mathbf{v})^{2}\,,$
(21)

assuming the velocity of the particle has changed from $\mathbf{v}$ to
$\mathbf{v}+\Delta\mathbf{v}$.

If we make the impulse approximation — assuming the satellite undergoes a
perturbation or ‘shock’ much shorter than most other timescales in the
problem, so particle positions are approximately constant throughout the event
(e.g. Taylor & Babul, 2001) — we can get an estimate of the energy change from
the shock by treating it as a point-mass perturber $M_{p}$, passing by the
satellite on a linear trajectory with a velocity $v_{p}$ and an impact
parameter $b$. Tidal heating from the shock should accelerate particles,
imparting a change in velocity

$\Delta\mathbf{v}=\dfrac{2GM_{p}}{v_{p}b^{2}}[-x,y,0]\,\,\,,$ (22)

where $(x,y)$ is the position of the particle in the plane of the encounter
(see, e.g. Mo et al., 2010). Since $\Delta\mathbf{v}\sim\mathbf{r}$,
$\mathbf{v}^{2}\sim r^{2}$, and the impulse approximation predicts that
$|\Delta\mathcal{E}_{K}|\propto r^{2}$ at large radii, and that it may also
have some additional energy or velocity dependence from the first term on the
right-hand-side of Equation (21). On the other hand, at small radii, adiabatic
shielding prevents particles from experiencing much net acceleration, since
their internal orbital timescales become shorter than the timescale of the
changing tidal field close to the centre of the satellite (e.g. Gnedin &
Ostriker, 1999). Thus the change in $\mathcal{E}_{K}$ should be approximately
zero in this limit. After this instantaneous change in kinetic energy, the
system re-virializes. The negative heat capacity of the system suggests that
there will be a resulting decrease in kinetic energy (an increase in
$\Delta\mathcal{E}_{K}$).

Fig. 9 shows the change in the specific energy of the bound particles over 5
orbits. The top panels show the change in relative energy as a function of
initial radius (left-hand) and initial energy (right-hand panels), while the
bottom panels show how the mass is distributed in the system, as a function of
the same variables. While there some variation in $|\Delta\mathcal{E}|$ across
the satellite (for particles that remain bound to the system,
$|\Delta\mathcal{E}|$ decreases with radius and also with binding energy),
overall it is fairly constant, varying by $\sim 10$ per cent over the range of
radius or energy that contains most of the mass555Note that in Paper II, we
showed that the scatter in the initial relative energy of particles stripped
on a given orbit was also $\sim 10$ per cent – see Fig. 5 of Paper II.. We see
larger variations only at large radii or low binding energies, where particles
are close to being stripped. Overall, the mean shift in relative energy from
orbit to orbit is fairly constant as a function of particle radius or initial
energy, and also from particle to particle, even after 5 orbits. This shows
why the energy truncation model produces a good approximation to the rate and
radial dependence of the actual mass loss, particularly in the slow mass-loss
case. The largest variations in $|\Delta\mathcal{E}|$ are close to the edge of
the system, where particles are close to being stripped.

The weak and declining dependence of $|\Delta\mathcal{E}|$ on particle radius
may seem counter-intuitive, since the particles at large radii are expected to
get larger velocity kicks in the satellite frame. Overall, the change in
relative energy does not scale with radius as expected for at least two
reasons. First, particles at small radii tend to be more bound, and therefore
can experience more acceleration and still remain bound, while particles at
large radii that receive large kicks will be lost from the system. We do not
track this unbound material in our model. This mass loss reduces the average
energy change in the material that remains bound, and may also explain why the
energy change decreases in magnitude quite rapidly close to the boundary of
the system. The second reason is that the main contribution to the relative
energy change is actually from the potential term, which decreases more or
less uniformly over the whole inner part of the satellite, as mass is stripped
primarily from the outside in. This term will deviate from a constant only in
the outer part of the satellite where some but not all of the mass is lost, so
here too this contributes to the rapid change in $|\Delta\mathcal{E}|$ at
large radii.

Figure 9:  The average change in relative energy of the satellite self-bound
particles ($\Delta\mathcal{E}$). Note that negative changes in energies
correspond to particles becoming less bound to the satellite. Specific
energies are calculated for both the Slow Simulation (top) and Fast Simulation
(bottom) after one, three and five orbits. Shaded regions show the $1-\sigma$
standard deviation in each bin. We also show the mass distribution of
particles as a function of initial radius (left) and initial energy (right).

We can also consider which particles are removed in phase space. In Fig. 10 we
show the average change in relative energy as a function of initial particle
position. Close to the phase-space boundary of the system (at a constant
relative energy), the change in relative energy is constant. The decrease in
relative energy in the centre of the system is due to the overall mass of the
system decreasing, causing particles to become less bound. This figure
demonstrates that $\Delta\mathcal{E}_{P}$ for the self-bound particles mainly
depends on their initial energy, especially for the Slow Simulation.

Figure 10:  The average particle change in relative energy,
$\overline{\Delta\mathcal{E}}$, as a function of the initial location in phase
space for the self-bound remnant. We plot contours of constant energy in an
NFW profile with dashed lines. The change of energy is constant along the edge
of the system, which roughly corresponds to a constant energy, indicating that
change in relative energy is primarily a function of particles initial energy.
The decrease in $\mathcal{E}$ in the center of the system (initial
$(r,v)\approx(0,0)$) is due to the decreased self-bound mass of the satellite.

In summary, there are several reasons why the energy-truncation model provides
a good estimate of the rate and radial dependence of mass loss, or
equivalently, why mass loss remains ordered by the original relative energy to
a significant degree. While simple arguments predict that tidal heating should
depend on particle radius and possibly also on energy, this variation is
reduced by selective mass loss at large radii and/or low energies, and the net
variation in the remaining bound mass is much smaller. Furthermore, the total
change in relative energy is actually dominated by the potential energy
change, which is relatively constant throughout the inner part of the
satellite.

Overall, the physical picture for mass-loss in our framework is that as the
satellite moves through the background potential, the tidal fields heat its
particle orbits and deform its self-potential (see Appendix B). The lowest
$\mathcal{E}$ particles rapidly escape, since they are no longer bound to the
satellite. Once this mass loss has occurred, however, the potential well of
the system decreases by a roughly constant amount, and the energies of all the
particles are consequently shifted while remaining ordered. Our model
approximates this process by assuming the initial distribution function is
lowered by a constant amount on each orbit. As discussed in Section 9, future
work will focus on improvements to this energy-truncation model, for instance
by adding scatter to the change in relative energy.

##  9 Discussion

In Paper I, we described an approach to modelling tidal mass loss based on
truncating and lowering the subhalo DF by a specified tidal energy. In Paper
II we used this approach to develop a model for mass loss with zero free
parameters, and showed that it provided a good description of the evolution of
subhaloes with NFW density profiles. In this work, we have demonstrated that
energy truncation and our mass loss model can be applied to a wide range of
other density profiles, with similar accuracy. We found the model does an
excellent job of describing the density, mass and velocity profile evolution
in all the tested cases. Additionally, our model naturally captures the effect
that cuspy profiles conserve central density, while cored profiles have a
large decrease in central density, as described in Peñarrubia et al. (2010).

Idealized simulations of individual subhaloes (such as the ones used in this
paper) require an assumption about the initial profile; this is typically
taken to be NFW, despite evidence that dark matter haloes are probably better
described by Einasto profiles (e.g. Navarro et al., 2010; Klypin et al.,
2016). Additionally, a recent study by Brown et al. (2020) indicates that the
universal profile might be due to the narrow range of initial conditions used
in simulations, and the diversity of dark matter profiles should be much
larger than previously thought. Regardless, the density profile of dark matter
haloes in cosmological $N$-body simulations is never resolved below the scale
of the softening length. To establish the true behaviour of haloes in their
innermost regions, we will likely need a theoretical model for the origin of
the universal density profile.

As discussed in Paper II, mass loss predictions are sensitive to the central
density of the system. Artificial disruption is likely a huge problem in
cosmological simulations on subhalo scales (van den Bosch et al., 2018; Errani
& Peñarrubia, 2019; Errani & Navarro, 2021; Green et al., 2021) and is thought
to be caused by artificial constant-density cores on the scale of the
resolution limit of the simulation that drive enhanced tidal disruption.
Controlled simulations suggest that centerally-divergent profiles can never be
fully disrupted by a smooth potential (e.g. Kazantzidis et al., 2004; Errani &
Peñarrubia, 2019). Similarly, Amorisco (2021) recently argued that the
conditions of (1) a centrally divergent density profile $\rho\sim r^{-1}$ and
(2) an isotropic phase space distribution are sufficient to ensure that
subhaloes can never be disrupted. These finding are in agreement with the
energy-truncation model, which predicts that systems will naturally always
retain a bound remnant as long as $\bar{\rho}(r)\to\infty$ as $r\to 0$.

These issues are particularly interesting in the context of microhaloes. These
structures form early in the universe, and their size is determined by the
free streaming scale of the dark matter particles. Assuming dark matter
particles have a mass of 100 GeV, the mass of the smallest microhaloes is
approximately an Earth mass (e.g. Diamanti et al., 2015). These early haloes
probably have cuspy profiles, with an inner slope of $\sim-1.5$ (e.g. Ogiya et
al., 2016). Most evidence suggests that the steep central cusps of microhaloes
are not conserved, and mergers are thought to drive profiles towards the
universal density profile, with an inner slope of $-1$. For instance, Ishiyama
(2014) show that the cusp slope gradually becomes shallower with increasing
halo mass. More recently, by using cosmological zoom-in simulations, Wang et
al. (2019) were able to simulate Earth-mass haloes at redshift zero, and find
that they look NFW (with an inner slope of $-1$).

Generally, it is thought that major mergers are responsible for this
flattening (e.g. Ogiya et al., 2016; Angulo et al., 2017; Delos et al., 2019).
Contrary to this picture, however, isolated simulations of major mergers show
that the inner densities of haloes are typically conserved (e.g. Kazantzidis
et al., 2006; Drakos et al., 2019a, b). A possible solution to this
discrepancy is that isolated binary merger simulations are overly simplistic,
and do not take into account the complex gravitational interactions to which
haloes are subject. However, it is also possible that high central densities
are not conserved in cosmological simulations due to numerical artefacts. The
central densities of these microhaloes are especially important to dark matter
annihilation calculations. For instance, Ishiyama (2014) found that
microhaloes had higher densities than expected from extrapolations of the low-
redshift concentration–mass–redshift relation, while Okoli et al. (2018)
showed that if these densities are conserved, it could increase estimates of
the dark matter annihilation boost factor by up to two orders of magnitude.

Overall, there is a growing body of evidence that the evolution of the inner
central density of haloes and their concentrations are not well understood.
This is particularly concerning since models of substructure evolution are key
ingredients in lensing predictions and dark matter annihilation constraints.
The consequences of this are examined in a companion paper, where we show that
while lensing constraints are not particularly sensitive to assumptions about
subhalo central density, dark matter annihilation signals depend very
sensitively on this quantity, and derived constraints on dark matter
properties are therefore still uncertain (Drakos et al., 2022, in prep).

The energy-truncation model used in this work is a promising approach to
modelling substructure evolution. As examined thoroughly in Paper II, the
energy truncation model includes a number of simplifying assumptions which
could be relaxed to make more accurate predictions. A straightforward
extension of the model would be allow some additional particle-to-particle
scatter in the energy change between orbits; equivalently, we could remove
particles from the DF using a more gradual cut-off in energy, rather than an
abrupt truncation. Additionally, we only consider mass loss in discrete steps
once per orbit, at the time of pericentric passage. Continuous models for
tidal mass loss are often created by dividing the orbital period into discrete
steps and assuming a fraction of the mass outside the tidal radius is lost at
each of these steps, according to a characteristic time scale for mass loss
(e.g. Taylor & Babul, 2001). We found in Paper II that a continuous mass-loss
model of this kind requires an additional free parameter, calibrated to
simulations, and thus we decided to use the discrete orbit-averaged model
instead. However, this is clearly a simplification, and we aim to improve the
mass-loss model in future work.

While our findings support the assertion that our energy-truncation model is
universal, we have only considered spherical, isotropic subhaloes, evolving in
an NFW main potential. Examining more complicated ICs (including anisotropic
and multicomponent systems) will be the focus of future work. Another
interesting extension would be to include a disk potential to the host halo,
as it has been shown that including an embedded central disk potential to
dark-matter-only simulations increases the rate of subhalo disruption (e.g.
Garrison-Kimmel et al., 2017). For example, recent work by Webb & Bovy (2020)
considers more realistic external tidal fields, including baryonic bulge and
disk components, and finds that the baryonic components lead to more subhalo
disruption and lower mass subhaloes. Since the energy-truncation model works
well for a wide range of collisionless systems, it may also be useful for
predicting the evolution of multiple-component systems. For instance, tidal
stripping is thought to be a likely mechanism for creating observed dark-
matter-deficient ultra-diffuse galaxies (e.g. Ogiya et al., 2021b; Ogiya et
al., 2021a). In future work, we will examine the evolution of these galaxies
in energy space.

Overall, a growing body of evidence suggests that the central density and
concentration of dark matter haloes are not well understood. Physically-based
models for subhalo evolution are critical for correctly predicting density
profile evolution at small radii, as these regions cannot be resolved in
simulations. The theoretically-based energy-truncation model offers a
convenient tool to study the evolution of collisionless systems as they are
tidally stripped, allowing us to work towards a complete understanding of how
subhaloes evolve, and to ultimately place more robust and accurate constraints
on the properties of dark matter.

## Acknowledgements

The authors thank the anonymous referee for useful comments. NED acknowledges
support from NSERC Canada, through a postdoctoral fellowship. JET acknowledges
financial support from NSERC Canada, through a Discovery Grant. The
simulations for this work were carried out on computing clusters provided by
Compute Ontario (https://computeontario.ca/) and Compute Canada
(www.computecanada.ca).

_Software_ : Gadget-2 (Springel, 2005), Icicle (Drakos et al., 2017), numpy
(Harris et al., 2020), matplotlib (Hunter, 2007), scipy (Virtanen et al.,
2020).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Amorisco (2021) Amorisco N. C., 2021, arXiv e-prints, p. arXiv:2111.01148
  * Angulo et al. (2017) Angulo R. E., Hahn O., Ludlow A. D., Bonoli S., 2017, MNRAS, 471, 4687
  * Baltz et al. (2009) Baltz E. A., Marshall P., Oguri M., 2009, J. Cosmology Astropart. Phys., 2009, 015
  * Binney & Tremaine (1987) Binney J., Tremaine S., 1987, Galactic dynamics. Princeton University Press
  * Bonaca et al. (2020) Bonaca A., et al., 2020, ApJ, 892, L37
  * Boylan-Kolchin & Ma (2007) Boylan-Kolchin M., Ma C.-P., 2007, MNRAS, 374, 1227
  * Brown et al. (2020) Brown S. T., McCarthy I. G., Diemer B., Font A. S., Stafford S. G., Pfiefer S., 2020, MNRAS, 495, 4994
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Carlberg (2020) Carlberg R. G., 2020, ApJ, 889, 107
  * Choi et al. (2009) Choi J.-H., Weinberg M. D., Katz N., 2009, MNRAS, 400, 1247
  * Delos (2019) Delos M. S., 2019, Phys. Rev. D, 100, 063505
  * Delos et al. (2019) Delos M. S., Bruff M., Erickcek A. L., 2019, Phys. Rev. D, 100, 023523
  * Diamanti et al. (2015) Diamanti R., Catalan M. E. C., Ando S., 2015, Phys. Rev. D, 92, 065029
  * Drakos et al. (2017) Drakos N. E., Taylor J. E., Benson A. J., 2017, MNRAS, 468, 2345
  * Drakos et al. (2019a) Drakos N. E., Taylor J. E., Berrouet A., Robotham A. S. G., Power C., 2019a, MNRAS, 487, 993
  * Drakos et al. (2019b) Drakos N. E., Taylor J. E., Berrouet A., Robotham A. S. G., Power C., 2019b, MNRAS, 487, 1008
  * Drakos et al. (2020) Drakos N. E., Taylor J. E., Benson A. J., 2020, MNRAS, 494, 378
  * Drakos et al. (2022) Drakos N. E., Taylor J. E., Benson A. J., 2022, Do assumptions about the central density of halos effect lensing and dark matter annihilation calculations?, _In prep._
  * Einasto (1965) Einasto J., 1965, Trudy Inst. Astrofiz. Alma-Ata, 5, 87
  * Errani & Navarro (2021) Errani R., Navarro J. F., 2021, MNRAS, 505, 18
  * Errani & Peñarrubia (2019) Errani R., Peñarrubia J., 2019, MNRAS, p. 2998
  * Errani et al. (2021) Errani R., Navarro J. F., Ibata R., Peñarrubia J., 2021, arXiv e-prints, p. arXiv:2111.05866
  * Gao et al. (2008) Gao L., Navarro J. F., Cole S., Frenk C. S., White S. D. M., Springel V., Jenkins A., Neto A. F., 2008, MNRAS, 387, 536
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel S., et al., 2017, MNRAS, 471, 1709
  * Gnedin & Ostriker (1999) Gnedin O. Y., Ostriker J. P., 1999, ApJ, 513, 626
  * Green & van den Bosch (2019) Green S. B., van den Bosch F. C., 2019, MNRAS, 490, 2091
  * Green et al. (2021) Green S. B., van den Bosch F. C., Jiang F., 2021, MNRAS, 503, 4075
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hayashi et al. (2003) Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T., 2003, ApJ, 584, 541
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ishiyama (2014) Ishiyama T., 2014, ApJ, 788, 27
  * Kampakoglou & Benson (2007) Kampakoglou M., Benson A. J., 2007, MNRAS, 374, 775
  * Kazantzidis et al. (2004) Kazantzidis S., Magorrian J., Moore B., 2004, ApJ, 601, 37
  * Kazantzidis et al. (2006) Kazantzidis S., Zentner A. R., Kravtsov A. V., 2006, ApJ, 641, 647
  * King (1962) King I., 1962, AJ, 67, 471
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Limousin et al. (2005) Limousin M., Kneib J.-P., Natarajan P., 2005, MNRAS, 356, 309
  * Miller et al. (2020) Miller T. B., van den Bosch F. C., Green S. B., Ogiya G., 2020, MNRAS, 495, 4496
  * Mo et al. (2010) Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolution
  * Navarro et al. (2004) Navarro J. F., et al., 2004, MNRAS, 349, 1039
  * Navarro et al. (2010) Navarro J. F., et al., 2010, MNRAS, 402, 21
  * Ogiya et al. (2016) Ogiya G., Nagai D., Ishiyama T., 2016, MNRAS, 461, 3385
  * Ogiya et al. (2019) Ogiya G., van den Bosch F. C., Hahn O., Green S. B., Miller T. B., Burkert A., 2019, MNRAS, 485, 189
  * Ogiya et al. (2021a) Ogiya G., van den Bosch F. C., Burkert A., 2021a, arXiv e-prints, p. arXiv:2111.12104
  * Ogiya et al. (2021b) Ogiya G., Taylor J. E., Hudson M. J., 2021b, MNRAS, 503, 1233
  * Okoli et al. (2018) Okoli C., Taylor J. E., Afshordi N., 2018, J. Cosmology Astropart. Phys., 8, 019
  * Peñarrubia et al. (2008a) Peñarrubia J., McConnachie A. W., Navarro J. F., 2008a, ApJ, 672, 904
  * Peñarrubia et al. (2008b) Peñarrubia J., Navarro J. F., McConnachie A. W., 2008b, ApJ, 673, 226
  * Peñarrubia et al. (2009) Peñarrubia J., Navarro J. F., McConnachie A. W., Martin N. F., 2009, ApJ, 698, 222
  * Peñarrubia et al. (2010) Peñarrubia J., Benson A. J., Walker M. G., Gilmore G., McConnachie A. W., Mayer L., 2010, MNRAS, 406, 1290
  * Sereno et al. (2016) Sereno M., Fedeli C., Moscardini L., 2016, J. Cosmology Astropart. Phys., 2016, 042
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Stref et al. (2019) Stref M., Lacroix T., Lavalle J., 2019, Galaxies, 7, 65
  * Stücker et al. (2021) Stücker J., Angulo R. E., Busch P., 2021, arXiv e-prints, p. arXiv:2107.13008
  * Taylor & Babul (2001) Taylor J. E., Babul A., 2001, ApJ, 559, 716
  * Tormen et al. (1997) Tormen G., Bouchet F. R., White S. D. M., 1997, MNRAS, 286, 865
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2019) Wang J., Bose S., Frenk C. S., Gao L., Jenkins A., Springel V., White S. D. M., 2019, arXiv e-prints, p. arXiv:1911.09720
  * Webb & Bovy (2020) Webb J. J., Bovy J., 2020, MNRAS, 499, 116
  * Widrow & Dubinski (2005) Widrow L. M., Dubinski J., 2005, ApJ, 631, 838
  * van Kampen (2000) van Kampen E., 2000, preprint,  (arXiv:astro-ph/0002027)
  * van den Bosch & Ogiya (2018) van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066
  * van den Bosch et al. (2018) van den Bosch F. C., Ogiya G., Hahn O., Burkert A., 2018, MNRAS, 474, 3043

##  Appendix A Accuracy of profile predictions

Fig. 11 shows the residuals in the density and mass profiles of the energy-
truncation model compared to the simulation results. Residuals are calculated
as the difference between the model predictions and simulations, divided by
the model predictions. We note that there is very little mass within $\sim
r_{\rm unit}$—particularly for the EinLow and KingLow simulations—which makes
the relative errors in the mass profiles untrustworthy at low radii.
Nonetheless, these residual plots offer important insights into the accuracy
of the energy-truncation model.

The Slow Simulation (top), generally captures the model to within $\sim 10$
per cent, except when $r$ approaches the truncation radius. An exception is
the KingLow profile, whose central density is very sensitive to the total mass
of the system. In general, for the Slow Simulation, the energy-truncation
model over-predicts central density and mass of the stripped satellite. The
Fast Simulations (bottom), shows a larger disagreement with the energy-
truncation model (but the central density is still fairly well predicted, to
within $\sim 20$ per cent or better). Since this simulation was chosen as a
case where the energy-truncation model does not work as well (see Papers I and
II), this is expected. In general, large deviations between the model and
simulations can be seen at radii larger than $\sim 2\,r_{\rm unit}$. At these
larger radii, the energy-truncation model predicts too much mass for the Slow
Simulation, and too little mass for the Fast Simulation. This result is
consistent with the mass-loss curves shown in Fig. 5.

(a)  (b)

Figure 11:  Residuals in the density profiles and mass profiles of the tidally
stripped systems, when comparing the energy-truncation model to the Slow (top)
and Fast (bottom) Simulations, as shown in Fig. 11. The shaded regions
correspond to 10 per cent (dark grey) and 20 per cent (light grey) agreement.
In general, the model does a fairly good job at predicting the central density
(within $\sim 2\,r_{\rm unit}$). At larger radii, there are significant
deviations between the energy-truncation model and the self-bound remnant from
simulations. The exact origin of these deviations are unclear, but are likely
tied to the detailed mass-loss timescales of each particle, and how the bound
particles of the satellite are defined.

This paper has mainly focussed on predicting the central density of
satellites, which the energy-truncation model does fairly well. The material
at large subhalo radii are more dynamically complicated, and are sensitive to
how we define bound particles. Future work will focus on understanding the
detailed time scales of mass loss for all particles, and whether some of the
bound mass is in the process of leaving the system.

##  Appendix B The deforming bowl picture of tidal mass loss

Stücker et al. (2021) recently proposed that there is a natural energy at
which particles will be lost, calculated from the “boosted potential". They
advocate for a model for tidal stripping in which mass loss is explained by a
lowering of the escape energy caused by the tidal field termed the “deforming
bowl” picture. This naturally leads to the concept of a “truncation energy".

In tidal stripping analyses, we often only consider the self potential of the
satellite. However, the satellite exists within the large scale potential of
the host it is orbiting, as illustrated in Fig. 12. This figure shows an
(infinitely extended) NFW potential of mass $M_{\rm sat}$ located at
$R_{P}=50\,r_{\rm unit}$ of the host halo ($M_{H}=300M_{\rm sat}$).

Figure 12:  The global potential around an NFW satellite located at
$R=50\,r_{\rm unit}$ of its host halo. The potential from the host halo
strongly effects the local potential of the satellite.

To account for the affect of the host potential on the satellite, Stücker et
al. (2021) define the boosted potential as:

$\Phi_{\rm boost}(\mathbf{x})=\Phi(\mathbf{x})+\mathbf{x}\cdot\mathbf{a_{0}}$
(23)

where $\Phi$ is the total (host plus satellite) potential and $\mathbf{a_{0}}$
is an additional apparent acceleration. Thus, in our example, the boosted
potential subtracts the local large scale gradient of the host halo from the
total potential.

Under the assumption of spherical symmetry for both the host and satellite
halo, and assuming that the satellite is located at pericentre, Stücker et al.
(2021) define the boosted potential of the satellite system as:

$\displaystyle\Phi_{\rm boost}(R)=$ $\displaystyle\Phi_{H}(R)+\Phi_{\rm
sat}(r)$ (24)
$\displaystyle-(R-R_{p})\dfrac{GM_{H}(<R_{p})}{R_{p}^{2}}-\Phi_{H}(R_{p})\,\,\,,$

where $R=R_{p}+r$. The final term has been subtracted to give the self-
potential of the satellite. Solving for the saddle points of Equation (24)
gives:

$\bar{\rho}_{\rm
sat}(r)=\pm\left[\dfrac{R_{p}}{r}-\dfrac{R_{p}^{3}}{r(R_{p}\pm
r)^{2}}\dfrac{M_{H}(<R_{p}\pm
r)}{M_{H}(<R_{p})}\right]\bar{\rho}_{H}(R_{p})\,\,\,.$ (25)

Fig. 13 shows the boosted potential of an NFW satellite. For this example, the
saddle points are located at $r=-7.3\,r_{\rm unit}$ and $8.4\,r_{\rm unit}$.
In the deforming bowl picture, particles are lost as they spill over the top
of this bowl.

Figure 13: The potential of an NFW satellite (dashed grey line) and the
boosted potential (solid black line). In the “deforming bowl" picture,
particles outside the saddle points of the boosted potential are lost.

This picture of mass-loss can be easily integrated into our energy-truncation
model with the observation that Equation (25) is in the form of Equation (5),
with

$\eta_{\rm boost}=\pm\left(\dfrac{R_{p}}{r_{\rm
lim}}-\dfrac{R_{p}^{3}}{r(R_{p}\pm r_{\rm lim})^{2}}\dfrac{M_{H}(<R_{p}\pm
r_{\rm lim})}{M_{H}(<R_{p})}\right)\,\,\,.$ (26)

Other common definitions of $\eta$, as discussed in detail in Paper II,
include the well-known Roche ($\eta=2$) and Jacobi limits ($\eta=3$),

$\eta_{1}=2-\dfrac{{\rm d}\ln M}{{\rm d}\ln R}\,\,\,,$ (27)

which describes a satellite and host with extended bodies, and

$\eta_{2}=\dfrac{\omega^{2}}{\omega_{c}^{2}}-\dfrac{1}{\omega_{c}^{2}}\dfrac{{\rm
d}^{2}\Phi}{{\rm d}R^{2}}\,\,\,,$ (28)

which includes the centrifugal force.

In Paper II we found that in the energy-truncation framework, $\eta_{2}$
overestimates mass-loss rates and work the best for circular orbits, while
$\eta_{1}$ underestimates mass loss and work the best for radial orbits. The
value we adopt, $\eta_{\rm eff}$, is the orbital average of the instantaneous
value of $\eta_{2}$ (Equation (7)). In practice, we found $\eta_{\rm eff}$
(which lies between $\eta_{1}$ and $\eta_{2}$) worked best to describe mass-
loss.

There are important differences between Equation (26) and the alternate $\eta$
definitions listed above. First, there are two $\eta_{\rm boost}$ values,
reflecting the fact that the tidal field is stronger on the side of the
satellite closer to the host halo. We will take $\eta_{\rm boost}$ to be the
negative root, as this corresponds to stronger tidal forces. Secondly, since
the value of $r_{\rm lim}$ depends on the satellite model, $\eta_{\rm boost}$
does not only depend on the host halo and orbit, but also on the details of
the satellite profile.

In Fig. 14 we compare $\eta_{\rm boost}$ to $\eta_{\rm eff}$ as well as
commonly used values in the literature. We consider both the Slow and Fast
Simulation, and use an infinitely extended NFW profile to calculate the
$\eta_{\rm boost}$ values. We find that the boosted potential calculation
predicts an $\eta$ value between $\eta_{1}$ and $\eta_{\rm eff}$. Since we
interpret $\eta_{\rm eff}$ as being the “true" value, this suggests that
$\eta_{\rm boost}$ underestimates mass-loss in our framework.

Figure 14:  Values of $\eta$ used in the mass-loss prescription, Equation (6).
From left to right (1) $\eta_{\rm eff}$ is the value we use in Equation (7)
shown to match well with simulations (Paper II) so can be considered close to
the "true" value (2) $\eta_{1}$, Equation (27), which works well for radial
orbits (3) $\eta_{2}$, Equation (28), which works well for circular orbits and
(4) $\eta_{\rm boost}$, Equation (26), which is calculated from the boosted
potential. For comparison, we also show the Roche ($\eta=2$) and Jacobi
($\eta=3$) limits. $\eta_{\rm boost}$ is lower than the “true" value,
$\eta_{\rm eff}$, and therefore under-predicts mass loss in the energy-
truncation framework.

Although $\eta_{\rm boost}$ (and effectively the truncation energy) calculated
from this formulation of the boosted potential underestimates the best-fit
value (which is expected to be roughly equivalent to $\eta_{\rm eff}$), the
boosted potential framework seems to be a natural extension to the energy-
based tidal stripping model we advocate. A potential explanation for why
$\eta_{\rm boost}$ is too low is because we neglected the centrifugal term;
i.e. we could subtract a term $\frac{1}{2}|\omega\times\mathbf{R}|^{2}$ from
Equation (24) to account for the corotating potential. We leave more detailed
explorations of the boosted potential to future work.


# Semicomplete Arithmetic Sequences, Division of Hypercubes, and the Pell
Constant

Z. Hoelscher

###### Abstract -

In this paper we produce a few continuations of our previous work on
partitions into fractions. Specifically, we study strictly increasing integer
sequences $\\{n_{j}\\}$ such that there are partitions for all integers less
than the floor of $G$, where
$G=\frac{n_{1}}{j}+\frac{n_{2}}{j}+\cdots+\frac{n_{j-1}}{j}+\frac{n_{j}}{j}$,
and all summands are distinct terms drawn from
$\frac{n_{1}}{j}+\frac{n_{2}}{j}+\cdots+\frac{n_{j-1}}{j}+\frac{n_{j}}{j}$. We
call such sequences “semicomplete”. We find that there are only three
semicomplete arithmetic sequences. We also study sequences that give the
maximum number of pieces that an $M$ dimensional hypercube can be cut into
using $N-1$ hyperplanes. We find that these are semicomplete in one, two,
three, and four dimensions. As an aside, we use one of our generating
functions to produce what appears to be a new identity for the Pell constant,
a number which is closely connected to the density of solutions to the
negative Pell equation.

Keywords : Integer partitions, $q$-series, Pell equation, Pell constant

Mathematics Subject Classification (2020) : 05A17, 11B65

##  1 Introduction

Integer partitions have a long history of study, with mathematicians from
Euler to Ramanujan producing interesting results. Past mathematicians have
also studied complete sequences. These are defined in the work of Brown [1] as
sequences $\\{f_{i}\\}^{\infty}_{i=1}$ such that all terms are positive and
every $n\in\mathbb{N}$ can be written as shown below.

$n=\sum^{\infty}_{i=1}\alpha_{i}f_{i}\hskip 14.45377pt\alpha_{i}\in\\{0,1\\}$
(1.0.1)

In this paper we build on our previous work, where we studied partitions of
integers into fractions with constant denominators and distinct numerators
drawn from sets of even or odd integers [2]. We noted that for any positive
integer $j$, one can write $j^{2}=\sum_{n=1}^{j}(2n-1)$. After dividing both
sides of the equation by $j$, one can represent $j$ as a sum of fractions.

$j=\frac{1}{j}+\frac{3}{j}+\frac{5}{j}+\cdots+\frac{2j-1}{j}$ (1.0.2)

We then proved that when $j>2$, any $k\in\mathbb{N},k<j$ can be written as the
sum of some combination of distinct terms from the series that sums to produce
$j$. We noted that this is not always possible when the numerators are even
integers, prompting us to consider what other sequences allow for similar
behavior. Such a generalization of the problem is considered in this paper.
Suppose we have a series of the form given below:

$G=\frac{n_{1}}{j}+\frac{n_{2}}{j}+\cdots+\frac{n_{j-1}}{j}+\frac{n_{j}}{j}$
(1.0.3)

We look for solutions for natural numbers $k<\lfloor G\rfloor$, where
$k=\frac{a_{1}}{j}+\frac{a_{2}}{j}+\cdots+\frac{a_{m-1}}{j}+\frac{a_{m}}{j}$,
and $a_{1},\ldots,a_{m}$ are distinct coefficients drawn from
$\\{n_{1},n_{2},\ldots,n_{j}\\}$. We also require that $k,j\in\mathbb{N}$,
$j>2$. We say $\\{n_{j}\\}$ is semicomplete if it enables such solutions for
all $k$, where the terms of the sequence strictly increase, and all terms are
positive. We refer to such sequences as semicomplete because a semicomplete
sequence is not necessarily complete. For example, we prove that the cake
numbers are both complete and semicomplete, yet the odd integers are only
semicomplete. This is because the odd integers enable sums for all $kj$ under
the restrictions above, though one cannot find integers such as 2 as the sum
of distinct odd summands. This makes semicomplete sequences perhaps a bit less
obvious, as one does not have to be able to find partitions for every natural
number. Such questions lead naturally to the following theorem.

###### Theorem 1.0.1

The only semicomplete arithmetic sequences are $\\{1,3,5,7,9\ldots\\}$,
$\\{2,3,4,5,6\ldots\\}$, and $\\{1,2,3,4,5\ldots\\}$.

We prove this theorem in Section 2.

In our previous paper, we conjectured that the cake numbers and lazy caterer’s
sequence are both semicomplete [2]. We prove this now, and generalize these
sequences to higher dimensions. Let $C_{M}^{N}$ denote the maximum number of
pieces that an $M$-dimensional hypercube can be cut into using $N-1$
hyperplanes. Let $\\{C_{M}^{N}\\}$ be the sequence of such terms for a given
value of $M$ where $N$ varies. (Observe that when $M=2$ we have the lazy
caterer’s sequence, and when $M=3$, we have the cake numbers.) Note that we
use $N-1$ for the number of cuts so that $N=1$ gives the first term of the
sequence, $N=5$ gives the fifth term, et cetera, rather than $N=0$ for the
first term and $N=4$ for the fifth.

###### Theorem 1.0.2

For values of $M$ equal to one, two, three, or four, $\\{C_{M}^{N}\\}$ is
semicomplete.

We provide a proof for this in Section 3. We also give an interesting
conjecture on these sequences.

###### Conjecture 1.0.3

$\\{C_{M}^{N}\\}$ is semicomplete for any value of $M\geq 1$,
$M\in\mathbb{N}$.

If true, this would imply that an infinite number of semicomplete sequences
exist, hence one cannot list them all. To solve the full semicompleteness
problem, one would then have to find the general necessary and sufficient
conditions for any type of sequence to be semicomplete.

Perhaps a more interesting problem would be to search for sequences that, like
the odd integers, are semicomplete but not complete. We know of one example,
$\\{2,3,4,5,6.\dots\\}$, though it would be interesting to find more.

We note that if one does not require $k$ to be an integer, one can write the
generating function for the case of odd numerators as shown below.

$\prod_{n=1}^{j}(1+q^{2n-1})=\sum_{kj=0}^{j^{2}}f_{O_{j}}(k,j)q^{kj}$ (1.0.4)

As a somewhat interesting aside, we give what appears to be a new $q$-series
identity for the Pell constant, $\mathcal{P}_{Pell}$. This constant is closely
connected to the product $\prod_{n=1}^{j}(1+q^{2n-1})$ as well as the density
of solutions to the negative Pell equation $x^{2}-Dy^{2}=-1$ [3]. It is
currently unknown whether the Pell constant is transcendental, though it has
been proven to be irrational [3].

###### Theorem 1.0.4

$F(q)=\sum_{m=0}^{j+1}\sum_{n=0}^{j+1}\frac{(-1)^{\frac{3n}{2}}q^{\frac{1}{2}(n^{2}-2n+2)}(-1)^{\frac{m}{2}}q^{\frac{1}{2}(m^{2}-2m)}}{q+1}{j+1\brack
n}_{q}{j+1\brack m}_{q}$ (1.0.5)

$\mathcal{P}_{Pell}=1-\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{j\rightarrow\infty}\;$}}F(-\frac{1}{2})\approx
0.58057$

We provide a proof for this in Section 4. We note that ${j+1\brack n}_{q}$ is
a Gaussian binomial coefficient, which is the $q$-analog of a binomial
coefficient. Gaussian binomial coefficients are defined as shown below [4],
where $m\geq n$.

${m\brack
n}_{q}=\frac{(1-q)(1-q^{2})\cdots(1-q^{m})}{(1-q)(1-q^{2})\cdots(1-q^{n})(1-q)(1-q^{2})\cdots(1-q^{m-n})}$
(1.0.6)

One should note that when $n>m$, ${m\brack n}_{q}$ is defined as zero. We can
then restate our conjecture from our previous paper [2] in terms of $F(q)$,
where we note that $F(q)$ is a polynomial in $q$ when fully simplified.

###### Conjecture 1.0.5

Let $G(q)$ contain only the terms from $F(q)$ where the exponent of $q$ is
$kj$, $k,j\in\mathbb{N},k<j,j>2$. The coefficients of $G(q)$ are always either
unimodal or bimodal.

##  2 Proof for Theorem 1.0.1

As we require that we must have solutions for all $k<\lfloor G\rfloor$ where
$j>2$, $k,j\in\mathbb{N}$, we can prove that a sequence does not work by
finding a value of $k$ that has no solution for some value of $j$. While
somewhat tedious, this process works well to eliminate all sequences other
than the three semicomplete arithmetic sequences. We can then readily show
these three are semicomplete.

Proof.  When $j=3$, we have a series of the form

$G=\frac{a}{3}+\frac{a+b}{3}+\frac{a+2b}{3}$ (2.0.1)

where $a,b\in\mathbb{N}$. Note that $a$ is the first term in the sequence of
numerators and $b$ is the difference between consecutive terms in that
sequence.

Let $a=1$:

$G=\frac{1}{3}+\frac{1+b}{3}+\frac{1+2b}{3}$ (2.0.2)

If $b=0$, we have the sequence $\\{1,1,1,1,\ldots\\}$ as numerators. In this
case $G$ is always exactly one, hence there are no natural numbers $k<\lfloor
G\rfloor$. We thus see that $\\{1,1,1,1,\ldots\\}$ is not semicomplete. If
$b=1$, we have the consecutive integers as numerators. It is easy to see that
any integer less than the sum of the first $n$ consecutive integers can be
found as the sum of some combination of terms from the first $n$ consecutive
integers, hence this sequence is semicomplete.

If $b=2$ our numerators are the odd integers, which we have previously proven
to be semicomplete [2]. If $b>2$ all terms in the series but $\frac{1}{3}$ are
greater than one, hence it is impossible to find a combination of terms that
sum to $k=1$. This then proves that any sequence with $a=1$ and $b>2$ is not
semicomplete.

We now consider sequences where $a>1$. If $b=0$ and $a=2$ there is no solution
for $k=1$, as $1$ is not an integral multiple of $\frac{2}{3}$. If $a=3$ and
$b=0$ there are solutions for all $k$ when $j=3$, but not when $j=4$, hence
$\\{3,3,3,3,\ldots\\}$ is not semicomplete. If $a>3$ and $b=0$ all terms in
the series

$G=\frac{a}{3}+\frac{a+b}{3}+\frac{a+2b}{3}$ (2.0.3)

are greater than one, hence there is no combination for $k=1$. All sequences
$\\{a,a,a,a,\ldots\\}$, where $a>3$ are thus incomplete.

If $a=2$ and $b=1$ we have the sequence $\\{2,3,4,5,\ldots\\}$. One can see
that any integer except $1$ and $(2+3+4+5+\cdots+(j+1))-1$ can be found as the
sum of some combination of terms from $\\{2,3,4,5,\ldots\\}$. We know that
$\frac{1}{j}$ and $\frac{(2+3+4+5+\cdots+(j+1))-1}{j}$ cannot be integers when
$j>2$, hence every integer $k<\lfloor G\rfloor$ can be found as the sum of
some combination of terms. We know $\frac{(2+3+4+5+\cdots+(j+1))-1}{j}$ cannot
be an integer through the argument shown below.

$G=\frac{(j+1)(j+2)}{2j}-\frac{1}{j}=\frac{j+3}{2}$ (2.0.4)
$\frac{(2+3+4+5+\cdots+(j+1))-1}{j}=\frac{j+3}{2}-\frac{1}{j}=\frac{j}{2}-\frac{1}{j}+\frac{3}{2}$
(2.0.5)

We know that $\frac{j}{2}-\frac{1}{j}+\frac{3}{2}$ cannot be an integer when
$j>2$, hence $\frac{(2+3+4+5+\cdots+(j+1))-1}{j}$ cannot be an integer when
$j>2$.

If $a=2$ and $b>1$ then there is no combination for $k=1$ when $j=3$, as every
term in the series but $\frac{2}{3}$ is greater than $1$. Such sequences are
thus incomplete. If $a=3$ and $b=1$ or $b=2$, there is no combination for
$k=2$ when $j=3$. If $a=3$ and $b=3$ there is no combination for $k=1$ when
$j=4$. If $a=3$ and $b>3$ there is no combination for $k=2$ when $j=3$. If
$a>3$ then all terms in the series are greater than $1$ when $j=3$, regardless
of the value of $b$, hence there is no combination for $k=1$. Such sequences
are thus incomplete. $\Box$

##  3 Proof for Theorem 1.0.2

Proof.  In one dimension: It is easy to see that in this case we are cutting a
line into pieces with points. This then results in
$\\{C_{1}^{N}\\}=\\{1,2,3,\ldots\\}$, which is the sequence of consecutive
integers. We know this sequence is semicomplete.

In two dimensions:

$\\{C_{2}^{N}\\}=\bigg{\\{}\frac{N^{2}-N+2}{2}\bigg{\\}}$ (3.0.1)

The sum of the first $t$ terms can be found through induction.

$\sum_{N=1}^{t}\frac{N^{2}-N+2}{2}=\frac{1}{6}t(t^{2}+5)$ (3.0.2)

We see that when $t$ is sufficiently large, the $(t+1)^{TH}$ term is less the
sum of the first $t$ terms. We have

$\frac{(t+1)^{2}-(t+1)+2}{2}<\frac{1}{6}t(t^{2}+5)$ (3.0.3)

where $3<t$.

Suppose one takes a set of the first $t$ terms and one sees that any integer
less than the sum of those $t$ terms can be found as the sum of some
combination of terms from that set. Now take an integer $I$ that satisfies the
following inequality:

$\sum_{N=1}^{t}\frac{N^{2}-N+2}{2}<I<\sum_{N=1}^{t+1}\frac{N^{2}-N+2}{2}$
(3.0.4)

We know that the last term is less than the sum of all previous terms, where
$t$ is sufficiently large. If every integer less than the sum of the previous
terms can be found as the sum of some combination of the previous terms, one
can find a partition for any integer $I$ by summing the $(t+1)^{TH}$ term
$C_{2}^{t+1}$ and the partition for $I-C_{2}^{t+1}$, as
$I-C_{2}^{t+1}<\sum_{N=1}^{t}\frac{N^{2}-N+2}{2}$. By induction, one then sees
that any integer less than the sum of $\\{C_{2}^{N}\\}$ can be found as the
sum of some combination of the terms from this sequence, if this can be done
for the first $t$ terms, where $t\in\mathbb{N}$, $t>3$. (If you can do this
for the first $t$ terms, you can do it for the first $t+1$ terms, and thus the
first $t+2$ terms, et cetera.) This method closely follows that used by Brown
to show whether a sequence is complete [1]. Here we take the first four terms
of the sequence, then examine the integers that can be found as sums of the
terms.

$\\{C_{2}^{N}\\}=\\{1,2,4,7\\}$ (3.0.5) $1+2+4+7=14$ (3.0.6)

We see we have $1=1$, $2=2$, $3=1+2$, $4=4$, $5=1+4$, $6=2+4$, $7=7$, $8=1+7$,
$9=7+2$, $10=1+2+7$, $11=4+7$, $12=1+4+7$, and $13=2+4+7$. This then confirms
that the lazy caterer’s sequence is complete, and hence also semicomplete.
(Note that one can individually check cases to see that it is semicomplete
when using the first 3, 2, or 1 terms. We just needed to apply this technique
to a set of more than 3 terms.) To illustrate that this holds upon adding
another term, we examine the first five terms of the sequence.

$\\{C_{2}^{N}\\}=\\{1,2,4,7,11\\}$ (3.0.7) $1+2+4+7+11=25$ (3.0.8)

We already know we have partitions for $1,2,\ldots,14$. Now we need to find
these for $15,16\ldots,25$. As an example, we take the case of $I=20$. We see
that $20-11=9=2+7$, hence $20=11+2+7$.

In higher dimensions: There is a known recurrence relation that describes such
sequences.

$C_{M}^{N}=C_{M}^{N-1}+C_{M-1}^{N-1}$ (3.0.9) $C_{M}^{1}=1$

For a discussion of this recurrence, see [5].

In three dimensions:

$C_{3}^{N}=\bigg{\\{}\frac{N^{3}-3N^{2}+8N}{6}\bigg{\\}}$ (3.0.10)
$\sum_{N=1}^{t}\frac{N^{3}-3N^{2}+8N}{6}=\frac{t(t+1)(t^{2}-3t+14)}{24}$
(3.0.11)

We then have

$\frac{(t+1)^{3}-3(t+1)^{2}+8(t+1)}{6}<\frac{t(t+1)(t^{2}-3t+14)}{24}$
(3.0.12)

where $4<t$.

In four dimensions:

$C_{4}^{N}=\bigg{\\{}\frac{N^{4}-6N^{3}+23N^{2}-18N+24}{24}\bigg{\\}}$
(3.0.13)
$\sum_{N=1}^{t}\frac{N^{4}-6N^{3}+23N^{2}-18N+24}{24}=\frac{t(t^{4}-5t^{3}+25t^{2}+5t+94)}{120}$
(3.0.14)

We then have

$\frac{(t+1)^{4}-6(t+1)^{3}+23(t+1)^{2}-18(t+1)+24}{24}<\frac{t(t^{4}-5t^{3}+25t^{2}+5t+94)}{120}$
(3.0.15)

where $5<t$. One can then apply the same process used in two dimensions to
prove semicompleteness in three and four dimensions. For three dimensions one
must manually check semicompleteness for $t\leq 4$, and for four dimensions,
one must manually check semicompleteness for $t\leq 5$. $\Box$

##  4 Proof for Theorem 1.0.4

Proof.

$F(q)=\sum_{m=0}^{j+1}\sum_{n=0}^{j+1}\frac{(-1)^{\frac{3n}{2}}q^{\frac{1}{2}(n^{2}-2n+2)}(-1)^{\frac{m}{2}}q^{\frac{1}{2}(m^{2}-2m)}}{q+1}{j+1\brack
n}_{q}{j+1\brack m}_{q}$ (4.0.1)

We can split this double sum into the product of two sums.

$F(q)=\biggl{(}\sum_{m=0}^{j+1}\frac{(-1)^{\frac{3m}{2}}q^{\frac{1}{2}(m^{2}-2m+2)}}{q+1}{j+1\brack
m}_{q}\biggr{)}\biggl{(}\sum_{m=0}^{j+1}(-1)^{\frac{m}{2}}q^{\frac{1}{2}(m^{2}-2m)}{j+1\brack
m}_{q}\biggr{)}$ (4.0.2)

By noting the definition for a binomial coefficient and manipulating the
series, we can arrive at a more convenient form.

$\binom{m}{2}=\frac{m!}{2!(m-2)!}=\frac{1}{2}m(m-1)$ (4.0.3)
$F(q)=\biggl{(}\frac{q}{q+1}\biggr{)}\biggl{(}\sum_{m=0}^{j+1}(-1)^{m}\biggl{(}-\frac{1}{q}\biggr{)}^{\frac{m}{2}}q^{\binom{m}{2}}{j+1\brack
m}_{q}\biggr{)}\biggl{(}\sum_{m=0}^{j+1}\biggl{(}-\frac{1}{q}\biggr{)}^{\frac{m}{2}}q^{\binom{m}{2}}{j+1\brack
m}_{q}\biggr{)}$ (4.0.4)

We have the following identity from the work of Koekoek and Swarttouw [6]. It
follows from the $q$-binomial theorem.

$(a;q)_{n}=\sum_{k=0}^{n}(-a)^{k}q^{\binom{k}{2}}{n\brack k}_{q}$ (4.0.5)

This identity then enables us to rewrite $F(q)$.

$F(q)=\frac{q}{q+1}\biggl{(}\sqrt{-\frac{1}{q}};q\biggr{)}_{j+1}\biggl{(}-\sqrt{-\frac{1}{q}};q\biggr{)}_{j+1}$
(4.0.6)

The identity given below is known [6].

$(a^{2};q^{2})_{n}=(a;q)_{n}(-a;q)_{n}$ (4.0.7) $\therefore\hskip
5.69046ptF(q)=\frac{q}{q+1}\biggl{(}\sqrt{-\frac{1}{q}};q\biggr{)}_{j+1}\biggl{(}-\sqrt{-\frac{1}{q}};q\biggr{)}_{j+1}=\frac{q}{q+1}\biggl{(}-\frac{1}{q};q^{2}\biggr{)}_{j+1}$
(4.0.8)

The $q$-analog of the Pochhammer symbol can be defined as the following
product for positive $k$.

$(a;q)_{n}=\prod_{k=0}^{n-1}(1-aq^{k})$ (4.0.9)
$F(q)=\frac{q(-\frac{1}{q};q^{2})_{j+1}}{q+1}=\frac{q}{q+1}(1+q^{-1})(1+q)(1+q^{3})(1+q^{5})\cdots=\prod_{n=1}^{j}(1+q^{2n-1})$
(4.0.10)

This product is a generating function for partitions of $kj$ into distinct odd
integers drawn from $\\{1,3,5,\ldots,2j-1\\}$.

$\prod_{n=1}^{j}(1+q^{2n-1})=\sum_{kj=0}^{j^{2}}f_{O_{j}}(k,j)q^{kj}$ (4.0.11)

The Pell constant $\mathcal{P}_{Pell}$ can be written as shown below [3].

$\mathcal{P}_{Pell}=1-\prod_{k=0}^{\infty}\biggl{(}1-\frac{1}{2^{2k+1}}\biggr{)}$
(4.0.12)
$\lim_{j\to\infty}F\biggl{(}-\frac{1}{2}\biggr{)}=\lim_{j\to\infty}\biggl{(}\prod_{n=1}^{j}\biggl{(}1+\biggl{(}-\frac{1}{2}\biggr{)}^{2n-1}\biggr{)}\biggr{)}=\prod_{n=0}^{\infty}\biggl{(}1-\frac{1}{2^{2n+1}}\biggl{)}=1-\mathcal{P}_{Pell}$
(4.0.13)

Hence
$\mathcal{P}_{Pell}=1-\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{j\to\infty}\;$}}(F(-\frac{1}{2}))$.
$\Box$

## References

  * [1] J. L. Brown Jr., A note on complete sequences of integers, Amer. Math. Monthly, 68 (1961), 557–560, available online at the URL: https://www.jstor.org/stable/2311150?seq=1
  * [2] Z. Hoelscher and E. Palsson, Counting restricted partitions of integers into fractions: symmetry and modes of the generating function and a connection to $\omega(t)$, PUMP J. Undergrad. Res., 3 (2020), 277–307, available online at the URL: https://journals.calstate.edu/pump/article/view/2428
  * [3] P. Stevenhagen, The number of real quadratic fields having units of negative norm, Experimental Mathematics, 2 (1993), 121–136, available online at the URL: https://projecteuclid.org/euclid.em/1048516217
  * [4] K. O’Hara, Unimodality of Gaussian coefficients: A constructive proof, _J. Combin. Theory Ser. A_ , 53 (1990), 29-52, available online at the URL: https://www.sciencedirect.com/science/article/pii/009731659090018R
  * [5] N.J.A. Sloane, Sequence A000125, The On-Line Encyclopedia of Integer Sequences, available online at the URL: https://oeis.org
  * [6] R. Koekoek and R. F. Swarttouw, The Askey-Scheme of Hypergeometric Orthogonal Polynomials and its $q$-Analogue, Technische Universiteit Delft, Faculty of Technical Mathematics and Informatics Report, (1998).