A core-halo pattern of entropy creation in gravitational collapse

TL;DR

This paper introduces a kinetic theory model revealing a core-halo entropy pattern during gravitational collapse, showing entropy destruction in the core and creation in the halo, which may help identify structure formation.

Contribution

It presents a novel kinetic theory approach demonstrating core-halo entropy patterns in gravitational collapse without prior assumptions, aiding structure formation analysis.

Findings

Entropy destruction in the core and creation in the halo during collapse

Core-halo pattern emerges without prior assumptions

Proposes a new scheme for identifying structure formation

Abstract

This paper presents a kinetic theory model of gravitational collapse due to a small perturbation. Solving the relevant equations yields a pattern of entropy destruction in a spherical core around the perturbation, and entropy creation in a surrounding halo. This indicates collisional "de-relaxation" in the core, and collisional relaxation in the halo. Core-halo patterns are ubiquitous in the astrophysics of gravitational collapse, and are found here without any of the prior assumptions of such a pattern usually made in analytical models. Motivated by this analysis, the paper outlines a possible scheme for identifying structure formation in a set of observations or a simulation. This scheme involves a choice of coarse-graining scale appropriate to the structure under consideration, and might aid exploration of hierarchical structure formation, supplementing the usual density-based methods for highlighting astrophysical and cosmological structure at various scales.

Tables

Table 1. Table 1: Errors in the approximation η ​ ( k ) = k J ​ σ , 𝜂 𝑘 subscript 𝑘 J 𝜎 \eta(k)=k_{\text{J}}\sigma, as calculated in Wren ( 2018 ) for a range of values of k / k J . 𝑘 subscript 𝑘 J k/k_{\text{J}}. To two significant figures, the error is as predicted by the − 3 ​ σ ​ k 2 / ( 2 ​ k J ) 3 𝜎 superscript 𝑘 2 2 subscript 𝑘 J -3\sigma k^{2}/(2k_{\text{J}}) correction from Eq. ( 56 ). The results make clear that, for any k ≪ k J , much-less-than 𝑘 subscript 𝑘 J k\ll k_{\text{J}}, the approximation η ​ ( k ) ​ t = k J ​ σ ​ t 𝜂 𝑘 𝑡 subscript 𝑘 J 𝜎 𝑡 \eta(k)t=k_{\text{J}}\sigma t is very good except for extremely large time-scales t ​ k J ​ σ . 𝑡 subscript 𝑘 J 𝜎 tk_{\text{J}}\sigma.

$k/k_{\text{J}}$	${10}^{-1}$	${10}^{-2}$	${10}^{-3}$	${10}^{-4}$
Approximation error $\left[k_{\text{J}}\sigma -\eta \left(k\right)\right]/k_{\text{J}}\sigma$	$1.5\text{×}{10}^{-2}$	$1.5\text{×}{10}^{-4}$	$1.5\text{×}{10}^{-6}$	$1.5\text{×}{10}^{-8}$
Time-scale $t{k}_{\text{J}}\sigma$ before significant approximation errors	$6.7\text{×}{10}^1$	$6.7\text{×}{10}^3$	$6.7\text{×}{10}^5$	$6.7\text{×}{10}^7$

Equations

$$f(1,2)=f(1)f(2)+\frac{1}{N}g(1,2)+O\left(\frac{1}{N^2}\right)$$

$$f(1,2)\approx f(1)f(2)+\frac{1}{N}g(1,2).$$

$$\text{a}(1,2)=-Gm(N-1)\frac{({\text{x}}_1-{\text{x}}_2)}{|{\text{x}}_1-{\text{x}}_2{|}^3}\approx -GmN\frac{({\text{x}}_1-{\text{x}}_2)}{|{\text{x}}_1-{\text{x}}_2{|}^3},$$

$$\begin{gathered} \frac{\partial f(1)}{\partial t}+{\text{v}}_1\cdot \frac{\partial f(1)}{\partial {\text{x}}_1}+\int \text{a}(1,2)f(2)\text{d}(2)\cdot \frac{\partial f(1)}{\partial {\text{v}}_1} \\ =-\frac{1}{N}\frac{\partial }{\partial {\text{v}}_1}\int \text{a}(1,2)g(1,2)\text{d}(2), \end{gathered}$$

$$\frac{\partial g(1,2)}{\partial t}+\Bigg{\{}\textbf{{v}}_{1}\boldsymbol{\cdot}\frac{\partial g(1,2)}{\partial\textbf{{x}}_{1}}+\int\!\textbf{{a}}(1,3)f(3)\,\text{d}(3)\boldsymbol{\cdot}\frac{\partial g(1,2)}{\partial\textbf{{v}}_{1}}\\ +\frac{\partial f(1)}{\partial\textbf{{v}}_{1}}\boldsymbol{\cdot}\int\textbf{{a}}(1,3)\,g(3,2)\,\text{d}(3)\\ +\left(\textbf{{a}}(1,2)-\int\!\textbf{{a}}(1,3)f(3)\,\text{d}(3)\right)\boldsymbol{\cdot}\frac{\partial f(1)}{\partial\textbf{{v}}_{1}}f(2)+\ (1)\leftrightarrow(2)\Bigg{\}}=0,$$

$$\frac{\partial f(1)}{\partial t}+{\text{v}}_1\cdot \frac{\partial f(1)}{\partial {\text{x}}_1}+\int \text{a}(1,2)f(2)\text{d}(2)\cdot \frac{\partial f(1)}{\partial {\text{v}}_1}=0.$$

$$M({\text{v}}_j)=\frac{\exp \left[-\text{nicefrac}{v}_j^{2}2{\sigma }^2\right]}{(2\pi {\sigma }^2{)}^{\text{nicefrac}32}}.$$

$$k_{\text{J}}^2\equiv \frac{4\pi Gmn}{{\sigma }^2}=\frac{4\pi GmN}{{\sigma }^{2}V}.$$

$$f(1)=(1-ϵ)f_0(1)+ϵf_1(1),$$

$$g(1,2)=(1-ϵ)g_0(1,2)+ϵg_1(1,2).$$

$$f_0({\text{x}}_1,{\text{v}}_1)=\{\begin{array}{ll}\frac{\exp \left[-\text{nicefrac}{\text{v}}_1^{2}2{\sigma }^2\right]}{V(2\pi {\sigma }^2{)}^{\text{nicefrac}32}}& \text{if }|{\text{x}}_1|<R\\ 0& \text{if }|{\text{x}}_1|\ge R\end{array}.$$

$$\int \text{a}(1,2)g_0(1,2)\text{d}(2)=GmN\int \frac{\mathbf{r}}{r^3}{g}_0(r,{\text{v}}_1,{\text{v}}_2){\text{d}}^3\mathbf{r}{\text{d}}^3{\text{v}}_2=0,$$

$$f_{1,\text{init}}(1)={\delta }^{(3)}({\text{x}}_1)M({\text{v}}_1),$$

$$\frac{\partial f_1(1)}{\partial t}+{\text{v}}_1\cdot \frac{\partial f_1(1)}{\partial {\text{x}}_1}+\int \text{a}(1,2)f_1(2)\text{d}(2)\cdot \frac{\partial f_0(1)}{\partial {\text{v}}_1}=0;$$

$$\frac{\partial g_{0}(1,2)}{\partial t}+\Bigg{\{}\textbf{{v}}_{1}\boldsymbol{\cdot}\frac{\partial g_{0}(1,2)}{\partial\textbf{{x}}_{1}}+\frac{\partial f_{0}(1)}{\partial\textbf{{v}}_{1}}\boldsymbol{\cdot}\int\textbf{{a}}(1,3)\,g_{0}(3,2)\,\text{d}(3)\\ +\textbf{{a}}(1,2)\boldsymbol{\cdot}\frac{\partial f_{0}(1)}{\partial\textbf{{v}}_{1}}f_{0}(2)+\ (1)\leftrightarrow(2)\Bigg{\}}=0\,;$$

$$\frac{\partial g_{1}(1,2)}{\partial t}+\Bigg{\{}\textbf{{v}}_{1}\boldsymbol{\cdot}\frac{\partial g_{1}(1,2)}{\partial\textbf{{x}}_{1}}+\int\!\textbf{{a}}(1,3)f_{1}(3)\,\text{d}(3)\boldsymbol{\cdot}\frac{\partial g_{0}(1,2)}{\partial\textbf{{v}}_{1}}\\ +\frac{\partial f_{0}(1)}{\partial\textbf{{v}}_{1}}\boldsymbol{\cdot}\int\textbf{{a}}(1,3)\,g_{1}(3,2)\,\text{d}(3)+\frac{\partial f_{1}(1)}{\partial\textbf{{v}}_{1}}\boldsymbol{\cdot}\int\textbf{{a}}(1,3)\,g_{0}(3,2)\,\text{d}(3)\\ +\textbf{{a}}(1,2)\boldsymbol{\cdot}\frac{\partial f_{0}(1)}{\partial\textbf{{v}}_{1}}f_{1}(2)+\textbf{{a}}(1,2)\boldsymbol{\cdot}\frac{\partial f_{1}(1)}{\partial\textbf{{v}}_{1}}f_{0}(2)\\ -\int\!\textbf{{a}}(1,3)f_{1}(3)\,\text{d}(3)\boldsymbol{\cdot}\frac{\partial f_{0}(1)}{\partial\textbf{{v}}_{1}}f_{0}(2)+\ (1)\leftrightarrow(2)\Bigg{\}}=0.$$

$$S\equiv -N\int f(1)\ln \left[f(1)\right]\text{d}(1).$$

$$\begin{gathered} \frac{\mathrm{d}S}{\mathrm{d}t}=-N\int {\left(\frac{\partial f(1)}{\partial t}\right)}_{\text{coll}}\ln \left[f(1)\right]\text{d}(1) \\ =\int \ln \left[f(1)\right]\frac{\partial }{\partial {\text{v}}_1}\cdot \int \text{a}(1,2)g(1,2)\text{d}(1,2). \end{gathered}$$

$${\left(\frac{\partial S_{{\text{x}}_1}}{\partial t}\right)}_{\text{flow}}=-N\frac{\partial }{\partial {\text{x}}_1}\cdot \int {\text{v}}_{1}f(1)\ln \left[f(1)\right]{\text{d}}^3{\text{v}}_1,$$

$$\begin{gathered} {\left(\frac{\partial S_{{\text{x}}_1}}{\partial t}\right)}_{\text{creation}}=\int \ln \left[f(1)\right]\frac{\partial }{\partial {\text{v}}_1}\cdot \int \text{a}(1,2)g(1,2)\text{d}(2){\text{d}}^3{\text{v}}_1 \\ =-N\int \ln \left[f(1)\right]{\left(\frac{\partial f(1)}{\partial t}\right)}_{\text{coll}}{\text{d}}^3{\text{v}}_1. \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}S}{\mathrm{d}t}=\int \ln \left[f_0(1)\right]\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g(1,2)\text{d}(1,2) \\ +\int \ln \left[1+\frac{ϵf_1(1)}{f_0(1)}+O\left({ϵ}^2\right)\right]\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)ϵg_1(1,2)\text{d}(1,2) \\ =-\frac{1}{2{\sigma }^2}\int v_1^2\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g(1,2)\text{d}(1,2) \\ +{ϵ}^2\int \frac{f_1(1)}{f_0(1)}\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g_1(1,2)\text{d}(1,2)+O({ϵ}^3), \end{gathered}$$

$$\begin{gathered} 0=\frac{\mathrm{d}}{\mathrm{d}t}\left\{\int \left[-\frac{G{m}^2\left(N-1\right)}{2}\int \frac{f(3)}{|{\text{x}}_1-{\text{x}}_3|}\text{d}(3)+\frac{1}{2}m{v}_1^2\right]f(1)\text{d}(1)\right\} \\ =\frac{G{m}^2\left(N-1\right)}{2}\int \left[\frac{{\text{v}}_3\cdot \frac{\partial f(3)}{\partial {\text{x}}_3}f(1)}{|{\text{x}}_1-{\text{x}}_3|}+\frac{f(3){\text{v}}_1\cdot \frac{\partial f(1)}{\partial {\text{x}}_1}}{|{\text{x}}_1-{\text{x}}_3|}\right]\text{d}(1,3) \\ -\frac{1}{2}m\int v_1^{2}f(2)\text{a}(1,2)\cdot \frac{\partial f(1)}{\partial {\text{v}}_1}\text{d}(1,2) \\ -\frac{m}{2N}\int v_1^2\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g(1,2)\text{d}(1,2) \\ =-\frac{m}{2N}\int v_1^2\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g(1,2)\text{d}(1,2), \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}S}{\mathrm{d}t}={ϵ}^2\int \frac{f_1(1)}{f_0(1)}\frac{\partial }{\partial {\text{v}}_1}\cdot \int \text{a}(1,2)g_1(1,2)\text{d}(1,2)+O({ϵ}^3) \\ \equiv -{ϵ}^{2}N\int \left[\frac{f_1(1)}{f_0(1)}\right]{\left(\frac{\partial f_1(1)}{\partial t}\right)}_{\text{coll}}\text{d}(1)+O({ϵ}^3), \end{gathered}$$

$$\frac{\mathrm{d}{S}_{\text{acg}}}{\mathrm{d}t}\sim {ϵ}^2{\int }_{\mathcal{K}}\frac{f_{1,\text{a}}(1)}{f_0(1)}\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g_{1,\text{a}}(1,2)\text{d}(1,2),$$

$${\bar{f}}_1(1)\equiv {\bar{f}}_1({\text{k}}_1,{\text{v}}_1)\equiv \int f_1({\text{x}}_1,{\text{v}}_1){\mathrm{e}}^{-\mathrm{i}{\text{k}}_1\cdot {\text{x}}_1}{\text{d}}^3{\text{x}}_1,$$

$${\widetilde{\stackrel{ˉ}{f}}}_1(1)\equiv {\widetilde{\stackrel{ˉ}{f}}}_1({\text{k}}_1,{\text{v}}_1,\omega )\equiv {\int }_0^{\infty }{\bar{f}}_1({\text{k}}_1,{\text{v}}_1,t){\mathrm{e}}^{\mathrm{i}\omega t}\text{d}t,$$

$$\begin{gathered} \int {\text{d}}^3{\text{x}}_1\left[\int \text{d}(2)\text{a}(1,2)g_{1,\text{a}}(1,2)\right]{\mathrm{e}}^{-\mathrm{i}{\text{k}}_1\cdot {\text{x}}_1} \\ =-\mathrm{i}\int \frac{{\text{d}}^3{\text{k}}_2}{(2\pi {)}^3}\frac{{\text{k}}_2{\overline{\stackrel{ˉ}{g}}}_{1,\text{a}}({\text{k}}_1-{\text{k}}_2,{\text{v}}_1,{\text{k}}_2,{\text{v}}_2)}{k_2^2}, \end{gathered}$$

$$\begin{gathered} {ϵ}^2{\int }_{\mathcal{K}}\frac{f_{1,\text{a}}(1)}{f_0(1)}\frac{\partial }{\partial {\text{v}}_1}\cdot \text{a}(1,2)g_{1,\text{a}}(1,2)\text{d}(1,2) \\ =-4\pi GmN\mathrm{i}{ϵ}^2{\int }_{\mathcal{K}}{\text{d}}^3{\text{v}}_1{\text{d}}^3{\text{v}}_2\frac{{\text{d}}^3{\text{k}}_1}{(2\pi {)}^3}\frac{{\text{d}}^3{\text{k}}_2}{(2\pi {)}^3} \\ \frac{\partial \left({\bar{f}}_{1,\text{a}}({\text{k}}_1,{\text{v}}_1)/f_0({\text{v}}_1)\right)}{\partial {\text{v}}_1}\cdot \frac{{\text{k}}_2{\overline{\stackrel{ˉ}{g}}}_{1,\text{a}}(-{\text{k}}_1-{\text{k}}_2,{\text{v}}_1,{\text{k}}_2,{\text{v}}_2)}{k_2^2}, \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{S}_{\text{acg}}}{\mathrm{d}t}\equiv -4\pi GmN\mathrm{i}{ϵ}^2{\int }_{\mathcal{K}(\beta k_{\text{J}})}{\text{d}}^3{\text{v}}_1\frac{{\text{d}}^3{\text{k}}_1}{(2\pi {)}^3}\frac{{\text{d}}^3{\text{k}}_2}{(2\pi {)}^3} \\ \frac{\partial \left({\bar{f}}_{1,\text{a}}({\text{k}}_1,{\text{v}}_1)/f_0({\text{v}}_1)\right)}{\partial {\text{v}}_1}\cdot \frac{{\text{k}}_2{\gamma }_{\text{a}}(-{\text{k}}_{+},{\text{v}}_1,{\text{k}}_2)}{k_2^2}, \end{gathered}$$

$${\gamma }_{\text{a}}({\text{k}}_1,{\text{v}}_1,{\text{k}}_2)\equiv \int {\overline{\stackrel{ˉ}{g}}}_{1,\text{a}}({\text{k}}_1,{\text{v}}_1,{\text{k}}_2,{\text{v}}_2){\text{d}}^3{\text{v}}_2.$$

Full text

††pubyear: 2018††pagerange: 1–29

1 Introduction

An early landmark in the study of kinetic theory entropy and gravitational collapse was the consideration by Antonov (1962)111There is an English translation by Antonov in Goodman & Hut (1985). of the thermodynamics of a model in which self-gravitating particles are confined within a (reflecting) sphere, in particular examining the “Antonov instability” associated with the absence of a global state of maximum kinetic theory entropy. Following this, Lynden-Bell & Wood (1968) looked at the link between thermodynamics, entropy and the formation of a core-halo pattern, showing numerically that, in their model of finite volume, the inner part of the system loses energy to the outer part, but that the kinetic energy – the temperature – of the inner part increases as its particles fall into the potential energy well. This negative specific heat capacity can drive a continuing “gravothermal catastrophe” which increases the flow of energy from the higher-temperature core to the lower-temperature halo. This can be illustrated (see also, for example, Binney & Tremaine, 2008, p572) as an expression of the virial theorem, albeit making an artificial division of the system into core and halo as an assumption for the argument, rather than its conclusion. This paper describes a kinetic theory model and analysis which avoids that artificial division, and sees a core-halo pattern emerge naturally in terms of a quantity we will construct, the asymptotic course-grained entropy creation rate, which indicates the effects of two-body collisions, and the rate of the system’s collisional relaxation. We then consider the potential physical implications of this result in terms of astrophysical and cosmological structure, in particular for supplementing the identification of structure based on patterns in density.

There are now a variety of approaches to considering the kinetic theory of gravitational collapse in cosmology and astrophysics. Some examples of useful sources include: Vereshchagin & Aksenov (2017) for a cosmological context; Binney & Tremaine (2008), the standard text on galactic dynamics; Merritt (2013) on galactic nuclei; and Heggie & Hut (2003) on star clusters.

Collisionless dynamics neglects the specific interactions between specific particles, focusing only on the evolution of their distribution under the “Vlasov equation” in the “mean-field” created by the average effect of all particles. In our current context, it is worth noting Lynden-Bell (1967)’s characterisation of “violent relaxation” through a coarse-graining of the collisionless Vlasov equation. Violent relaxation sees a coarse-grained Boltzmann entropy increasing not through entropy creation but via so-called “phase-mixing”. The fine-grained Boltzmann entropy remains unchanged, as for any evolution of the Vlasov equation. The evolution of entropy-like “H-functions” during violent relaxation was explored in Tremaine et al. (1986). Violent relaxation was more recently considered in, for example, Chavanis et al., 1996 and Dehnen, 2005.

Collisional dynamics is usually approached by adding a “collisional” term to the Vlasov equation to give, for example, the Fokker-Planck equation (Cohn, 1980, and as focused on in Binney & Tremaine, 2008), the Klimontovich equation (see, for example, Campa et al., 2009), the Smoluchowski equation (see Chavanis et al., 2002 on the emergence of a core-halo pattern in that context), Lenard-Balescu-type equations (see, for example, Chavanis, 2012), a generalized Landau equation (as in, for example, Chavanis, 2013) or a hierarchy of “BBGKY” equations linking $n$-particle distribution functions for $n=1,2,3...$ (see, for example, Gilbert, 1968, and other references mentioned below). The BBGKY hierarchy is often truncated to provide a closed set of equations. A form of truncated BBGKY hierarchy is used in this paper. It is also possible to consider collisional dynamics from the point of view of applied mathematics, as reviewed in Villani (2002).

Statistical mechanics also provides approaches to modelling distributions of self-gravitating particles – see, for example, Padmanabhan (1990) or Campa et al. (2014). Chavanis (2006) reviews the statistical mechanics of self-gravitating systems, illustrating their core-halo structure. One statistical mechanics approach is to use the one-dimensional Hamiltonian Mean Field (HMF) as a toy model to explore systems with long-range interactions. Exploration of the HMF, and its analogies with astrophysics, in Staniscia et al. (2009) illustrates the thermodynamics that can be associated with self-gravitating systems and, as set out in Levin et al. (2014), the appearance of a core-halo form in the HMF, starting from a particularly simple (so-called “water bag”) type of initial condition.

As described in, for example, Binney & Tremaine (2008) or Mo et al. (2010), there is also a widely-used approach of modelling astrophysical self-gravitating systems through fairly ad hoc, but useful, density distributions. These frequently have a core-halo form. Katz (1980) considers gravitational instabilities in connection with a selection of such models.

As mentioned, the approach used in this paper is based on the well-known BBGKY hierarchy. The BBGKY hierarchy is named from the initials of its pioneers: Bogolioubov (1946), Born & Green (1946), Kirkwood (1946), and Yvon (1935). It is reviewed in, for example, Balescu (1997), Huang (1987), and in an astrophysical context in Gilbert (1968). A notable astrophysical application was made in Weinberg (1993), which highlighted the importance of large-scale fluctuations of growing amplitude in the relaxation of a self-gravitating system. A broadly similar approach to ours, differing considerably in context and detail, is found in Heyvaerts (2010), which examines the collisional evolution and entropy of otherwise stable self-gravitating systems.

In this paper, we explore how much, and where, entropy is created and destroyed during gravitational collapse. The underlying notion of entropy is the Boltzmann entropy of the one-particle distribution function in kinetic theory phase space. To facilitate the derivation of Boltzmann entropy and the mapping of core-halo patterns, this paper uses physical space and physical velocity co-ordinates, and the Fourier transforms of the space co-ordinates. An alternative approach of so-called angle-action (or action-angle) variables is often used in the context of less homogeneous systems – see for example, Heggie & Hut (2003), Binney & Tremaine (2008), Heyvaerts (2010), and Chavanis (2012, 2013).

Section 2 describes the kinetic theory model used in this paper. The approach is to use the first two equations of the BBGKY hierarchy, under an assumption that the number of particles is sufficiently large that account need not be taken of the effect of collisions on one-particle distribution functions, while the collisional term is still useful for calculating the rate of creation of entropy. Our model system consists of a homogeneous Maxwellian distribution with a small central, nearly point-like, perturbation. Care is taken to define the underlying distribution so it is held equilibrium by external forces – an approach equivalent to the well-known Jeans swindle.

Section 3 constructs the rate of entropy creation, and then extracts the “asymptotic coarse-grained” part of that rate. It is “asymptotic” in the sense that it uses the term in the entropy creation rate with the strongest exponential time dependence, and so will asymptotically over time become the dominant part of entropy creation (provided the perturbation is small enough that perturbation theory is still valid when that dominance begins). “Coarse-grained” means that only small wave-number parts are retained – these represent the fastest-growing asymptotic parts, and also make calculations relatively tractable.

Section 4 then addresses the relevant evolution equations: recalling the well-known (Landau, 1946) solution for the first order perturbation of the distribution function, and then dealing with the zeroth and first order perturbations of the correlation function. Section 5 then calculates the creation rate for the asymptotic coarse-grained entropy in our system, both over all space, and for its distribution in space, before noting some simple variants of the main model and further avenues for exploration. Section 6 discusses physical implications, proposing a use for our approach in identifying structure formation, before a brief conclusion in Section 7.

Many detailed considerations and calculations are set out in the appendices to the paper. Appendix A explains a technicality in the definition of our initial point-like perturbation. Appendix B looks at the plasma dispersion function, which plays a key role in the formula for the distribution function associated with the perturbation. It gives an asymptotic formula for the “Landau zeros" of the plasma dispersion function, motivating a proof that, as often assumed, the residue sum rule for the inverse Laplace transform applies for the one-particle distribution function. Appendix C shows that the coarse-grained asymptotic number density and entropy density both have the same shape, strongly peaked near the initial central perturbation’s location: so the entropy density’s pattern is not a useful supplement to the number density’s. Appendix D calculates the correlation function associated with the underlying homogeneous distribution function, whilst Appendix E derives the equation for the correlation function’s perturbation. Appendices F and G provide two different approaches for finding an integral involving that correlation perturbation, which is needed to obtain the asymptotic coarse-grained entropy rate. Appendix F employs a propagator method, whilst Appendix G stays closer to Landau (1946)’s technique used to derive the perturbed distribution function. Appendix H sets out some detailed work needed to complete the entropy creation rate calculations. Further details of many calculations, and associated numerical integrations, are set out in a Mathematica notebook at Wren (2018).

2 The BBGKY hierarchy and the modelled system

2.1 Distributions and the BBGKY equations

This subsection establishes notation and recalls well-known equations and terminology. Following Gilbert (1968), we let $f(j)$ be the one-particle distribution function (DF), the probability density that a given particle is at the phase space point $(\textbf{{x}}_{j},\textbf{{v}}_{j}).$ This implies that $\int\!f(j)\,d(j)=1,$ where $d(j)$ is short-hand for $\text{d}^{3}\textbf{{x}}_{j}\text{d}^{3}\textbf{{v}}_{j},$ and, unless otherwise indicated, in this paper integrals will always be over the whole range of the relevant variable(s) of integration. Similarly, let $f(1,2)$ be be the two-particle distribution function, the probability that two given distinct particles are respectively at the phase space points $(\textbf{{x}}_{1},\textbf{{v}}_{1})$ and $(\textbf{{x}}_{2},\textbf{{v}}_{2}).$

Let $N$ be the total number of particles, which we will assume is very large. We define the (two-particle) correlation function, $g(1,2)\equiv(N-1)\left[f(1,2)-f(1)f(2)\right]$ and note that

[TABLE]

so that in practice we will assume

[TABLE]

Let $\textbf{{a}}(1,2)$ be $(N-1)$ times the acceleration of particle $1$ due to particle $2,$ in other words the acceleration if all the other particles were at position $2$ in phase space:

[TABLE]

which $G$ is Newton’s gravitational constant and $m$ is the small mass of each particle (we assume that each particle has the same mass).

For self-gravitating particles, from Gilbert (1968) we then have the pair of equations, truncated to leading order in ${\nicefrac{{1}}{{N}}},$

[TABLE]

and, also truncating to disregard three-particle correlations,

[TABLE]

where, for brevity, we used the abbreviation $+\ (1)\leftrightarrow(2)$ to indicate that we need to add terms which repeat all the terms in the braces, but with the variables $(1)=(\textbf{{x}}_{1},\textbf{{v}}_{1})$ and $(2)=(\textbf{{x}}_{2},\textbf{{v}}_{2})$ swapped over. Eqs. (4) and (5) represent the first two equations of the BBGKY hierarchy. The form of the acceleration, from Eq. (3), implies that Eqs. (4) and (5) suffer from short-range “ultraviolet” divergences, but these are not relevant when we focus on only large-scale, or coarse-grained, effects.

If the collisional term on the right-hand side of Eq. (4) is omitted, it becomes the Vlasov equation,

[TABLE]

The collisional assumption for our model is that $N$ is large enough, as for many physical systems, so the effect of the correlation function $g(1,2)$ on the one-particle DF $f(1)$ is minimal, even integrated over the whole of the time period we shall consider. In contrast, the one-particle DF will still drive the evolution of the correlation function. We will only take account of the correlation function’s effect on the one-particle DF when formulating our definition of entropy, or more precisely our definition of entropy creation. In that context, we will see in Section 3 that the collisional term will play a key role because, as is well known, Boltzmann entropy, as for any functional of the DF222Referred to as a Casimir functional., is invariant under the Vlasov equation.

For notational convenience we will often write expressions such as $x_{j}=|\textbf{{x}}_{j}|,v_{j}=|\textbf{{v}}_{j}|$ and write the Maxwellian velocity distribution as

[TABLE]

We can also identify some key parameters for our system, supplementing the Maxwellian velocity parameter $\sigma.$ Suppose there is a volume $V$ associated with our system (this volume will be specified in the next subsection). We write $n\equiv{\nicefrac{{N}}{{V}}}$ for the corresponding average number density. Recalling that we have assumed all particles have the same mass, we write that mass as $m.$ Our system has a characteristic length scale $k_{\text{J}}^{-1},$ where $k_{\text{J}}$ is the Jeans wave-number,

[TABLE]

Our system also has a characteristic time scale, its dynamical time, given by $(k_{\text{J}}\sigma)^{-1}.$

2.2 The perturbed model

Our model consists of an underlying DF and a perturbation. This subsection sets out our model, and recalls the BBGKY equations corresponding to the underlying and perturbation DFs. In light of our normalisation convention that DFs are probability densities, in particular integrating to unity over all phase space, we set

[TABLE]

where $\epsilon\ll 1$ is the perturbation parameter. Similarly we will also write the correlation function as

[TABLE]

We now define our underlying one-particle DF $f_{0}$. As in, for example, Binney & Tremaine (2008), we assume that the velocity dependence of our underlying DF is Maxwellian and that it is also homogeneous across a large spherical $V={\nicefrac{{4\pi R^{3}}}{{3}}}$ beyond which it vanishes,

[TABLE]

We will regard $R$ and $V$ as so large that, for many practical purposes, we are taking the limit $R\to\infty.$

The large, but finite, volume $V$ is needed to give us a non-zero probability density for the location of a particle at a given point. If we interpreted this as meaning that there is no mass beyond the volume $V,$ then the underlying distribution would itself not be in equilibrium, but would be undergoing gravitational collapse. Instead, we assume that it is held in equilibrium by external accelerations. These accelerations are set to be such that the acceleration integral terms in the BBGKY equations, Eqs. (4)-(6) are not limited to the volume $V$ but extend over all space. We are, in effect, assuming that there is mass beyond the volume $V,$ but that its response to the perturbation will be ignored. As noted in Binney & Tremaine (2008, p. 403), this kind of approach is a version of the well-known Jeans swindle (Jeans, 1902). In particular, this allows us to assume that $f_{0}$ is time-invariant under an evolution governed by the Vlasov equation, Eq. (6), because it enables us to disregard that equation’s acceleration integral.

We assume, and later confirm, that there is a “translation-invariant” $g_{0}$ which depends on position only through the distance $r\equiv|\mathbf{r}|\equiv|\textbf{{x}}_{2}-\textbf{{x}}_{1}|.$ Such a choice implies that

[TABLE]

where the last equality follows from the anti-symmetry of the integrand in $\mathbf{r}.$ This means that such a $g_{0}$ makes our $f_{0}$ not only time-invariant under the Vlasov equation, Eq. (6), but also collisionally time-invariant under evolution via the full one-particle BBGKY equation, Eq. (4). We also assume that $g_{0}$ is time invariant, in keeping with our aim that unperturbed functions represent an equilibrium state. Section 4.2 checks that there is indeed a translation- and time-invariant choice of $g_{0},$ consistent with Eq. (5).

We now define the first order perturbation $f_{1}$ of the one-particle distribution function. We assume the perturbation consists of particles of the same mass $m$ as the underlying distribution. The perturbation is defined by the initial value $f_{1,\text{init}}$ at $t=0$ of the perturbation $f_{1}.$ We shall use the idealised expression

[TABLE]

representing a sharp perturbation concentrated entirely at the origin $\textbf{{x}}_{1}=\textbf{{0}},$ with a Maxwellian velocity distribution which has the same parameter $\sigma$ as the underlying DF.333In Subsection 5.3 and Appendix H.4, we consider the case of choosing a different Maxwellian parameter for the perturbation. The formulation of $f_{1,\text{init}}$ via a spatial Dirac delta function is an approximation. The delta function takes an infinite value at the origin and so, in principle, is not compatible with perturbation theory. This technicality is dealt with in Appendix A. We will also take the initial perturbation to be uncorrelated, that is $g_{1}=0$ at time $t=0.$

We now, as is standard, write the BBGKY equations in terms of the underlying and perturbation functions, $f_{0},f_{1},g_{0}$ and $g_{1}.$ Using the Jeans swindle’s implications for acceleration integrals of $f_{0}$ and $g_{0},$ noted in and before Eq. (12), we have: the first order444In this context “order” refers to perturbation order in $\epsilon.$ Vlasov equation,

[TABLE]

the zeroth order correlation equation

[TABLE]

and the first order correlation equation

[TABLE]

Note that in the last two equations although $g_{0}$ is invariant under translations of both variables, it is not homogeneous in a single variable, implying that we cannot drop terms involving integrals of $g_{0}$ of forms like $\int\!\!\textbf{{a}}(2,3)\,g_{0}(1,3)\,\text{d}(3).$

In summary, our assumptions are as follows. That the total number of particles $N$ is large enough for the collisional assumption, described after Eq. (6), to hold. That the perturbation parameter $\epsilon,$ see Eq. (9), is small enough for first order perturbation theory to work. That the volume $V,$ of Eq. (11) is large enough for the effects of the exterior beyond $V$ to be ignored, at least in a large region around the centre of the initial perturbation. All these assumptions can be presumed to have a finite, but conceivably very long, lifetime from the introduction of the initial perturbation.

It is well known (see, for example, Binney & Tremaine, 2008, and Eq. (37) below) that the solution to the first order Vlasov equation, Eq. (14), is dominated by components of $f_{1}$ with small wave-numbers, which grow exponentially with time. We will call these the asymptotically-dominant, or simply asymptotic, parts of the solution.

3 Asymptotic coarse-grained entropy creation

3.1 A Boltzmann entropy rate formula

We now derive a formula for the rate of creation of the standard Boltzmann entropy in our model. This will motivate a more tractable definition of asymptotic coarse-grained entropy creation in the following subsection.

For a one-particle distribution function, and $N$ particles in total, the standard definition of entropy is given by the Boltzmann entropy

[TABLE]

It is well known that the overall entropy remains constant if $f$ is governed by the Vlasov equation, Eq. (6), and the creation of total entropy over time comes entirely from Eq. (4)’s collisional term, with

[TABLE]

For a given physical point $\textbf{{x}}_{1},$ similar arguments show that the rate of flow of entropy into $\textbf{{x}}_{1}$ is given by

[TABLE]

which for our perturbation, $f=(1-\epsilon)f_{0}+\epsilon f_{1},$ will be of order $\epsilon,$ and that the entropy creation rate is given, from Eq. (4), by

[TABLE]

Superficially, the sign of the entropy creation rate depends on whether $\ln[f(1)]$ is positive or negative, that is whether $f(1)$ is greater or less than $1.$ However, this is misleading: the velocity derivative in the collisional term implies that we can replace $\ln[f(1)]$ on the right-hand side of Eq. (20) with $\ln[f(1)/C]$ for any positive constant $C$ which is independent of velocity. This indicates that the entropy creation rate measures the tendency of the collisional term to push the DFs $f(\textbf{{x}}_{1},\textbf{{v}}_{1})$ at fixed $\textbf{{x}}_{1}$ towards a constant value independent of $\textbf{{v}}_{1}$ (the logarithm implying this measurement is relative to $f(\textbf{{x}}_{1},\textbf{{v}}_{1})$’s size). In summary, the entropy creation rate probes the tendency of collisions to flatten the velocity distribution – in other words, the rate of collisional relaxation at $\textbf{{x}}_{1},$ with a positive entropy rate indicating increasing collisional relaxation.

For now, we consider only the total creation of entropy over all space $V.$ From Eq. (18), we find

[TABLE]

where the $g_{0}$ term in the second line was eliminated using Eq. (12), making Eq. (21) exact up to, and including, order $\epsilon^{2}.$

We can see that the first summand in the last expression of Eq. (21) is zero, as follows. By conservation of energy, the total energy change for Eq. (4) must be zero, as shown in Irving & Kirkwood (1950) (and see also Martys, 1999). Making sure we count the gravitational potential from each particle interaction only once, we then have, after division by $N,$

[TABLE]

where we used the Vlasov equation, Eq. (4), to substitute for ${\nicefrac{{\partial f(1)}}{{\partial t}}}$ and ${\nicefrac{{\partial f(3)}}{{\partial t}}},$ omitting the terms from that equation which vanish as they give total derivatives, and then used substitution and partial integration to reach the final result. We can interpret Eq. (22) as reflecting the well-known result (see, for example, Binney & Tremaine, 2008, pp557-558) that two-particle collisional interactions do not result in the particles becoming bound. This implies that all purely collisional interactions begin and end with the particles relatively far apart with effectively zero mutual potential energy and hence the collision itself does not alter the two particles’ total kinetic energy.

Applying Eq. (22) to Eq. (21), we find,

[TABLE]

where we have defined the first order collisional term $\left({\nicefrac{{\partial f_{1}(1)}}{{\partial t}}}\right)_{\text{coll}},$ which arises from perturbation expansion of Eq. (4)’s right-hand side. Note that there is no order $\epsilon$ term in this equation for the rate of entropy creation, so under-densities (negative $\epsilon$) have the same entropy creation as equal and opposite over-densities (positive $\epsilon$).

As for Eq. (20), we can interpret Eq. (23) in terms of collisional relaxation. The integrand in Eq. (23)’s final expression consists of a weighting (the factor in square brackets) multiplied by the first order collisional term. This means that Eq. (23)’s entropy creation rate measures the (weighted average) tendency of collisions to suppress the perturbation (for positive rates) or enhance it (for negative rates). For example, if, at a phase space point $(1),$ we have $f_{1}(1)$ positive, then a positive entropy creation rate at that point implies the collisional term is negative, tending the eliminate the perturbation. The weighting $\left[{\nicefrac{{f_{1}(1)}}{{f_{0}(1)}}}\right]$ used in the integral’s averaging of such point rates over all phase space recognises changes affecting the perturbation $f_{1}$ relative to the size of the underlying distribution $f_{0}$ at the velocity concerned. It is worth recalling that in our model, the ${\nicefrac{{1}}{{N}}}$ factor in Eq. (4)’s collisional term implies that its tendency to enhance or suppress the perturbation is very slight – our model assumed very large $N$ and hence a very long relaxation time.

3.2 Defining the asymptotic coarse-grained entropy

As indicated at the end of Section 2, it is well known the perturbing DF $f_{1}$ is dominated by exponentially-growing asymptotic components, associated with small wave-numbers. It seems plausible, and we shall confirm in Subsections 4.2 and 4.3 below, that $g_{1}$ also has this behaviour. Motivated by Eq. (23), we introduce a definition of the rate of creation of asymptotic coarse-grained entropy along the following lines

[TABLE]

where $\mathcal{K}$ indicates that we do the integral over some region involving only small wave-numbers and the subscript “a” indicates that we only take the asymptotically-dominant term of each of $f_{1}$ and $g_{1}.$ We will indeed assume that $f_{1,\text{a}}$ and $g_{1,\text{a}}$ are also defined so that their Fourier transforms are non-zero only for small wave-numbers compatible with the region $\mathcal{K}.$

We accordingly now introduce Fourier transforms to isolate the small wave-number components. We also define our Laplace transform convention which we will need subsequently. In writing transforms, a single bar on a function indicates a Fourier transform with respect to one space variable, for example

[TABLE]

and a double bar, such as in $\bar{\bar{g}}_{1}(1,2),$ indicates a Fourier transform with respect to two space variables using the same convention. A tilde then further indicates a Laplace transform, as in

[TABLE]

or, with the same Laplace convention, $\accentset{\cong}{g}_{1}(1,2).$ These Fourier and Laplace conventions are as in Binney & Tremaine (2008), and it is worth noting that the Fourier and Laplace exponentials have differing signs. For concision, we will generally suppress both the time variable $t$ and its Laplace conjugate $\omega$ in our functions.

To define the region of interest, $\mathcal{K},$ we convert the expression of Eq. (24) into an integral over the Fourier space of wave-numbers. To do this, we use a standard approach to express the acceleration via the Fourier transform of Poisson’s equation for a point mass. This gives the Fourier transform with respect to $\textbf{{x}}_{1},$

[TABLE]

where $k_{2}=|\textbf{{k}}_{2}|$ is the wave-number associated with $\textbf{{k}}_{2}$ (and similarly below for other $k_{j}$). We can now see that Eq. (24) gives

[TABLE]

using the convolution theorem and integrating by parts.

We now specify more precisely the notion of “small” wave-numbers. Recall from Eq. (8) that we have a fundamental dynamical length scale in our model, corresponding to the Jeans wave-number, $k_{\text{J}}.$ For $\beta\ll 1,$ let $\mathcal{K}({\beta k_{\text{J}}})$ be the region of Fourier-transformed two-particle phase space which includes points with arbitrary velocities $\textbf{{v}}_{1}$ and $\textbf{{v}}_{2}$ and includes only small wave-numbers $\textbf{{k}}_{1},$ $\textbf{{k}}_{2}$ and $\textbf{{k}}_{+}\equiv\textbf{{k}}_{1}+\textbf{{k}}_{2},$ with $0<k_{1},\,k_{2},\,k_{+}<\beta k_{\text{J}}\ll k_{\text{J}}.$ We coarse-grain by setting our functions $f_{1,\text{a}}$ and $g_{1,\text{a}}$ to be zero for any argument with wave-number greater than or equal to $\beta k_{\text{J}}.$ Motivated by Eq. (28)’s last expression, we finally make our definition of the asymptotic coarse-grained entropy creation rate as

[TABLE]

where, as before, the subscript “a” indicates that for each of $\bar{f}_{1}$ and $\bar{\bar{g}}_{1}$ we keep only the part with the asymptotically-dominant growth, which will come from the poles of the Laplace transforms $\tilde{\bar{f}}_{1}$ and $\accentset{\cong}{g}_{1}$ with the most positive imaginary parts; we also wrote $\textbf{{k}}_{+}\equiv\textbf{{k}}_{1}+\textbf{{k}}_{2};$ and, for brevity, we defined

[TABLE]

In order to be consistent with our assumption that the system is within a finite volume $V$ of radius $R,$ we should require $\beta$ to be chosen so that $R^{-1}\ll k_{\text{J}}\beta\ll k_{\text{J}}.$ In Appendix H, we note after Eq. (120) that the contribution of very small $k_{1},k_{2},k_{+}\sim R^{-1}$ does not materially affect the entropy calculation for $\beta$ satisfying $R^{-1}\ll k_{\text{J}}\beta,$ so we do not need to account for the infrared cut-off at such small wave-numbers, and can take [math] as the lower limit for the wave-number integrals.

Note that because $S_{\text{acg}}$ in Eq. (29) is a measure of entropy creation, although it introduces a form of coarse-graining by including only small wave-numbers, it excludes phase-mixing which occurs in some other definitions of coarse-grained entropy (see, for example, the original reference in Lynden-Bell, 1967). As the name suggests, phase-mixing arises from particles with different velocities mixing closely together – an effect associated with entropy flow rather than entropy creation.

4 Solving the evolution equations

4.1 The Vlasov perturbation equation

The standard approach for deriving the first order perturbation $f_{1}$ from the Vlasov perturbation equation, Eq. (14), originated in Landau (1946). Its application to self-gravitating particles is reviewed in, for example, Binney & Tremaine (2008). This approach takes both Fourier and Laplace transforms of the distributions and the associated equations.

With the convention set out in Subsection 3.2 above, the Fourier-transformed initial perturbation is

[TABLE]

Fourier- and Laplace-transforming the Vlasov perturbation equation, Eq. (14), then gives

[TABLE]

where Fourier transforming the acceleration term was handled via its relationship with the Poisson equation for a point mass, along the lines described above Eq. (27). Following the standard method of Landau (1946), integrating Eq. (32) with respect to $\textbf{{v}}_{1}$ gives the first order perturbation for $\textbf{{k}}_{1}\neq 0$ as

[TABLE]

We used definitions, for $\omega$ with $\operatorname{Im}\omega>0,$

[TABLE]

and

[TABLE]

where the so-called plasma dispersion function $Z$ is defined by

[TABLE]

These three functions $Y,P$ and $Z$ can each be extended to $\operatorname{Im}\omega\leq 0$ by analytic continuation. Appendix B recalls and explores relevant properties of these functions.

As discussed in Binney & Tremaine (2008), the key property is that, for a given $0<k_{1}<k_{\text{J}},$ there is exactly one value of $\omega$ with $\operatorname{Im}\omega>0,$ satisfying the dispersion relation $k_{1}^{2}=k_{\text{J}}^{2}P(k_{1},\omega).$ For each such $0<k_{1}<k_{\text{J}}$ and $\omega,$ we define the positive real number $\eta(k_{1})$ by $\omega\equiv\mathrm{i}\eta(k_{1}),$ which is shown by the solid black line in Figure 1. For other values of $k_{1},$ there is no such $\omega$ with $\operatorname{Im}\omega>0.$

We now look at $\bar{f}_{1},$ which is the inverse Laplace transform of $\tilde{\bar{f}}_{1}.$ To do this inverse transform, we use the well-known residue formula recalled in Eq. (72). Appendix B.4 demonstrates that this formula does indeed work for $\tilde{\bar{f}}_{1}$ (this is often assumed without proof).

The residue formula is the sum of terms each with exponential time-dependency with rate $\omega,$ for each solution $\omega$ of the dispersion relation. This implies that, for a given wave-number $k_{1},$ the asymptotically-dominant part of $\bar{f}_{1}(\textbf{{k}}_{1},\textbf{{v}}_{1},t)$ – the term with the fastest growth in the limit of large times – is given by the value of $\omega$ with the largest positive imaginary value. We saw just above that, for $0<k_{1}<k_{\text{J}},$ this value is purely imaginary, and we then have

[TABLE]

for this asymptotically fastest-growing part of $\bar{f}_{1}(\textbf{{k}}_{1},\textbf{{v}}_{1},t).$ In Appendix C, we see that the asymptotic coarse-grained number density (by volume) has a strong positive central peak, with a much weaker, oscillating tail, as shown in Figure 3 (Right). The entropy density pattern’s provides no information additional to the pattern of the number density, because they are proportional. This helps motivate us to find out whether the entropy creation rate gives a different density pattern: we shall find this is indeed the case, and that the new pattern has a clear core-halo configuration.

4.2 The zeroth order correlation equation

We now solve the zeroth order correlation equation, Eq. (15), to calculate the Fourier transform, $\bar{\bar{g}}_{0},$ of the correlation function associated with the time-invariant equilibrium distribution function $f_{0}.$ In Subsection 2.2, we chose $g_{0}$ to be both time and translation invariant, implying that we can write

[TABLE]

and in turn this implies

[TABLE]

where $\bar{G}_{0}$ is the Fourier transform with respect to $\mathbf{r},$ of $G_{0},$ $\textbf{{k}}_{-}={\nicefrac{{(\textbf{{k}}_{2}-\textbf{{k}}_{1})}}{{2}}},$ and we continue to write $\textbf{{k}}_{+}=\textbf{{k}}_{1}+\textbf{{k}}_{2}.$

We substitute Eq. (38) into the zeroth order correlation equation, Eq. (15), and use time and translation invariance to get

[TABLE]

In Appendix D, we Fourier transform Eq. (40) and show that $\bar{G}_{0}(\textbf{{k}}_{-},\textbf{{v}}_{1},\textbf{{v}}_{2})=\bar{q}(\textbf{{k}}_{-})\operatorname{\mathcal{M}}(\textbf{{v}}_{1})\operatorname{\mathcal{M}}(\textbf{{v}}_{2}),$ where, neglecting a term which is irrelevant for $k_{1}<k_{\text{J}},$ from Eq. (81) we have

[TABLE]

As mentioned in Appendix D, this expression for $\bar{q}$ was derived in Kandrup (1983).

4.3 The first order correlation function

We now discuss the first order correlation equation, Eq. (16), for $g_{1},$ or rather for its Fourier and Laplace transform, $\accentset{\cong}{g}_{1}.$ Appendix E works through the Fourier and Laplace transforms of the first order correlation equation, finding that

[TABLE]

Given the zeroth and first order distribution functions, Eq. (42) is an integral equation for $\tilde{\bar{g}}_{1}$ in $\textbf{{v}}_{1}$ and $\textbf{{v}}_{2},$ with $\omega,\textbf{{k}}_{1}$ and $\textbf{{k}}_{2}$ as parameters.

General integral equations can only be solved numerically, a procedure which would be cumbersome given that our equation, Eq. (42), has multiple parameters. However, it is well known that integral equations such as this can be solved by use of a propagator approach, which, in essence, derives a Green’s function for the equation. This kind of approach is used in, for example, Ichimaru (1973) in the context of plasmas and Heyvaerts (2010) employing angle-action variables in the context of self-gravitating particles. The propagator approach to solving Eq. (42) is set out in Appendix F. In principle, this could give us results (which might involve new special functions of origin similar to $P$ and $Y$) for any $k_{1},k_{2}$ and $\omega.$ However, as indicated by Eq. (29), we focus on the case where $0<k_{1},k_{2}\ll k_{\text{J}}.$ This gives us tractable analytical solutions.

Another method of addressing Eq. (42) is based on the approach of Landau (1946) which was used to derive Eq. (33) for $\tilde{\bar{f}}_{1}.$ Along those lines, we can rearrange Eq. (42) and integrate it with respect to velocities. This approach is followed through in Appendix G for $0<k_{1},k_{2}\ll k_{\text{J}}.$ It is more complicated than the derivation of $\tilde{\bar{f}}_{1},$ requiring calculation of a number of integrals through solving a set of seven simultaneous equations.

Using either the propagator method of Appendix F, or the Landau method of Appendix G, we get a power series in $k_{j}$ for $\gamma_{\text{a}}$ as set out in Eq. (106). This provides the key function needed to calculate the rate of asymptotic coarse-grained entropy creation.

5 The rate of asymptotic coarse-grained entropy creation

5.1 The total entropy creation

We are now in a position to calculate the rate of asymptotic coarse-grained entropy creation over all space, using the results of Section 4 in the definition set out in Eq. (29). We start by looking at the factor in Eq. (29) that comes from the velocity derivative related to the one-particle distribution function. We find,

[TABLE]

Appendix B.1 reviews a relevant asymptotic series for small wave-numbers, and this provides a good approximation for the derivative, which is set out in Eq. (119).

Calculations in Appendix H.1 evaluate the asymptotic coarse-grained entropy creation rate of Eq. (29), using the formula for $\gamma_{\text{a}}$ from Eq. (106), and Eq. (119)’s series approximation. To leading order in $\epsilon$ and $k_{j},$ we find that we have

[TABLE]

We wrote $N_{1}\equiv\epsilon N$ for the number of particles associated with the perturbation, and recalled that $n={\nicefrac{{N}}{{V}}}$ is the average number density of the system. While $\epsilon$ is a measure of the size of the perturbation relative to the total number of particles, $N_{1}$ is a more absolute measure of that size. We also wrote $B\equiv{\nicefrac{{4\pi}}{{3\left(k_{\text{J}}\beta\right)^{3}}}}$ for the volume of a sphere associated with the coarse-graining scale $k_{\text{J}}\beta.$

The total (net) asymptotic coarse-grained entropy creation from the initial time [math] to some later time $t>0$ is then clearly

[TABLE]

Note that our asymptotic coarse-grained entropy decreases with time. This is consistent with the second law of thermodynamics because our system is not isolated – as set out following Eq. (11), it is subject to Jeans swindle forces from outside the system (beyond the radius $R$). These forces ensure that the acceleration integrals of the perturbation evolution equations, Eqs. (14)-(16), are not limited to the volume $V$ but may be taken to extend over all space. If, as a thought experiment, we imagine these forces as being created by some actual physical machine, its generation of these forces must produce entropy which at least offsets the negative entropy creation within the system. The negative entropy creation can also be viewed as indicating that, without Jeans swindle forces, it would be entropically-favourable for the system to contract under gravity.

5.2 The distribution of entropy creation in space

We now look at the distribution of asymptotic coarse-grained entropy creation in space. Clearly our system is spherically symmetric, so we look primarily at the shell density of the entropy creation rate, measuring the entropy creation rate on a thin sphere of a given radius, $r,$ centred on the initial perturbation. Referring back to Section 3, we see that the position of the entropy creation is indicated by the variable $\textbf{{x}}_{1},$ and that to extract that entropy creation’s spatial distribution we need to introduce a factor of $\delta(x_{1}-r)$ into the integral of Eq. (23). On Fourier-transforming, this corresponds to inserting a convolution with $4\pi r^{2}\operatorname{sinc}(k_{1}r)=4\pi r\sin(k_{1}r)/k_{1}$ into the integral of Eq. (29). Writing $\textbf{{k}}_{12}=\textbf{{k}}_{1}+\textbf{{k}}_{2}$ and $\textbf{{k}}_{01}=\textbf{{k}}_{0}+\textbf{{k}}_{1},$ we therefore have,

[TABLE]

where $S_{\text{acg}}^{\circ}(r)$ is the asymptotic coarse-grained entropy within a sphere of radius $r,$ and we defined $\mathcal{K}^{\prime}(\beta k_{\text{J}})$ to be the region $0<k_{0},k_{1},k_{2}<k_{\text{J}}\beta.$ Note that $k_{0}$ only appears in the integrand’s first and second factors because the associated convolution links together these two factors, whilst a convolution via $k_{1}$ links together the second and third.

Calculations in Appendix H.2 show that, at leading order,

[TABLE]

where $\hat{S}_{\text{acg}}^{\circ}$ might be termed the leading order entropy-creation pattern function and is shown in Figure 2 (Left), from calculations in Wren (2018). We have defined $\hat{S}_{\text{acg}}^{\circ}$ so that it corresponds to the numerical factor of $-0.0116$ in Eq. (44).

To compare the entropy creation density by shell of Eq. (47) with the total net entropy creation of Eq. (44), we should integrate Eq. (47) over $r$ to obtain the entropy creation in a region, for example the total entropy creation in either the core or the halo. As detailed in Appendix H.2, comparing Eq. (44) with the $r$-integrated Eq. (47), we see the size of either the core or the halo is around $\beta^{-2}\gg 1$ times larger than the total net entropy creation.

From Figure 2 (Left), we can see that there is a core, where entropy is destroyed, which spans from $r=0$ to r\approx{\nicefrac{{3.5}}{{k_{\text{J}}\beta}}}, and a halo, where entropy is created, which spans from r\approx{\nicefrac{{3.5}}{{k_{\text{J}}\beta}}} to, somewhere around, roughly, r\approx{\nicefrac{{7}}{{k_{\text{J}}\beta}}} to {\nicefrac{{8}}{{k_{\text{J}}\beta}}}. It is plausible that, beyond that radius, volumes of decreasing and increasing entropy alternate, but these are highly suppressed compared with the core and the halo.

Note that the notions of core and halo are scale dependent with respect to $\beta.$ A shell at radius $r,$ which is outside the core and halo at scale $\beta_{1},$ will be inside the halo for some smaller $\beta_{2}<\beta_{1},$ and inside the core for some yet smaller $\beta_{3}<\beta_{2}.$ The core-halo pattern is therefore more subtle than a simple core-halo model, with the boundaries of those two regions depending on the scale factor $\beta.$

Recall that the core and halo only indicate local destruction and creation of entropy – they take no account of entropy flowing from one place to another with the motion of particles. The core (resp. halo) is where particles’ associated asymptotic coarse-grained entropy tends to be destroyed (resp. created) by collisions. As an aside, recall that, given the long-range nature of gravity, these collisions are not necessarily close-range, but may be with distant particles outside the core (resp. halo). Referring to the paragraphs after Eq. (23), we also see that the core (resp. halo) is where collisions tend gently to enhance (resp. suppress) the perturbation. corresponding to collisional de-relaxation (resp. relaxation).

As confirmed in Appendix H.2, the absolute value of the total entropy creation in either the core or halo are of very similar size. This is to be expected, as we know from Subsection 5.1 that they must essentially cancel out. Appendix H.2 also notes that this absolute value is very close to being $\beta^{-2}$ times the size of the total net entropy creation. Since $\beta^{-2}\gg 1,$ the core-halo pattern is therefore much more pronounced that the total net entropy creation.

It might be asked if similar calculations to those in this and the previous subsection can be performed for different definitions of asymptotic coarse-grained entropy creation, in particular with different coarse-graining from that set out in the paragraph following Eq. (28). That set three constraints, which may be summarised as $0<k_{1},k_{2},k_{+}<k_{\text{J}}\beta.$ In Appendices H.1 and H.2, we consider the implications of relaxing any one of the three upper constraints. We get results for entropy creation totalled over all space, but the integrations needed for the distribution over space are no longer so tractable: their integrands no longer contain just the small wave-numbers needed for our analytical calculations. Moreover, Appendix H.1 notes that those alternatives match the time dependence of our perturbations less well, again suggesting a preference for the approach to constraints which we have adopted.

We could alternatively use a coarse-graining which matches the asymptotic time behaviour of our model yet more closely – albeit one which in a physical context we might be unlikely to select if we did not have an analytical expression for that behaviour. As set out in Appendix H.1, this is to choose $\mathcal{K}_{2}=\left\{(1,2):k_{1}^{2}+k_{2}^{2}+k_{+}^{2}<2k_{J}^{2}\beta^{2}\right\},$ which has the property that it captures all wave-numbers for which the exponential increase with time is estimated (to second order in $k_{j}$) as being faster than a given rate. The resulting entropy creation rate is as in Eq. (44), but with $-0.0116$ replaced by $$-0.0125$.$

As described in Appendix H.2, for $\mathcal{K}_{2}$ the leading order space distribution of the entropy creation rate satisfies the formula of Eq. (47), but with the entropy pattern function shown in Figure 2 replaced by that of Figure 4. We again see a core of entropy destruction surrounded by a halo of entropy creation, here with also small outer shells of entropy destruction and creation (compare with Figure 2 where there are perhaps such shells but of smaller amplitude than the estimated integration errors). The $\mathcal{K}_{2}$ coarse-graining, and another “taxicab” coarse-graining outlined at the end of Appendix H.2, suggest that a core-halo pattern of entropy destruction and creation may be generic for a broad class of coarse-grainings which are symmetrical with respect to $k_{1},k_{2}$ and $k_{+}.$

5.3 Further avenues for exploration, including of alternative systems

In Appendix H.4, we look at entropy creation when the Maxwellian parameter associated with the initial perturbation differs from that of the underlying perturbation. Excepting very large parameters, which our approximation methods cannot address, the leading order entropy creation rate is as in Eq. (47) and Figure 2, with no dependence on the initial perturbation’s Maxwellian parameter. This applies in particular for a perturbation with all its particles initially stationary.

It is also possible to vary the assumption $g_{1}(1,2,t=0)=0$ that the initial perturbation is uncorrelated. As discussed in Appendix H.4, choosing the same correlation as for the underlying distribution, $g_{1}(1,2,t=0)=g_{0}(1,2,t=0),$ makes no difference to our results.

We note some further possible avenues for investigation, looking at alternative systems to the one investigated here. It is an open question as to how analytically-tractable they might be, or whether approaches which are more numerical than used here might be needed.

This paper examined a single, almost point-like, initial perturbation. This might be replaced by a “dipole” pair of such perturbations, or a ”multipole” arrangement of many point-like perturbations.

More varied alternatives could also be considered. They might avoid having to impose a Jeans swindle restriction of the system to a finite volume, and the associated total net destruction of entropy.

One system that could be investigated would be the evolution, possibly with a small central perturbation, of a spherically-symmetric, but non-homogeneous, distribution of self-gravitating particles. This would roughly model a collapsing halo of dark matter, gas or stars. An introduction to such distributions can be found in Binney & Tremaine (2008, section 4.3).

Another alternative is to seek to model the evolution under gravity of a system comprising a small central perturbation in an underlying razor-thin disc of particles, the underlying disc being in equilibrium and also rotating. This might approximate a galactic disc, for example. A discussion of equilibrium distributions for razor-thin discs can be found in Binney & Tremaine (2008, section 4.5) or in Kalnajs (1976).

The assumption of equilibrium for those two distributions might avoid the requirement for a fixed exterior to keep the initial underlying distribution artificially in equilibrium. This would then perhaps show the dominance of entropy creation over its destruction within the whole system.

A different avenue of investigation would be to model the expanding universe, with a small perturbation being introduced. The universe’s expansion might then avoid the need for defining an artificial limit to the volume considered, with, perhaps, the expansion’s Hubble horizon providing a more natural finite limit. A start might be made with the formulation of kinetic theory equations set out for a Newtonian expanding universe in Kandrup (1983). However, this would not capture general relativity’s limitation of gravitational effects to a sphere of influence travelling out from the source at the speed of light, an effect which might define a sharper finite volume limit than the Hubble horizon. A general relativistic approach would therefore be preferable, see, for example, Vereshchagin & Aksenov (2017) or Andréasson (2011). The post-Newtonian approximation might be useful – see, for example, Poisson & Will (2014) for a general introduction, or Agón et al. (2011) and Ramos-Caro et al. (2012) for consideration of kinetic theory and the post-Newtonian approximation.

Near the initial central perturbation even a general relativistic system could, presumably, look very like the Newtonian system examined here if the central perturbation were not too dense. We might see asymptotic coarse-grained entropy destruction locally dominating over time, and, perhaps, being offset by distant entropy creation associated with the initial central perturbation’s expanding sphere of influence.

There is also the possibility of the central mass being a black hole. The accreting particles’ kinetic theory entropy would be distinct from the Bekenstein-Hawking entropy (Bekenstein, 1973; Hawking, 1976) of the black hole itself, and there might or might not be any noteworthy relationship between these entropies. A starting point for such a study might be provided by the very recent consideration of the kinetic theory of collisionless gas accreting on to a Schwarzschild black hole in Rioseco & Sarbach (2017).

6 Physical discussion

It is well known that the Universe has a multi-scale hierarchical structure, in which core-halo patterns are ubiquitous. The identification of observed or simulated astrophysical structure typically involves considering features of especially high or low densities, in physical space, or phase space. There is no unambiguous definition of structure in this context, which can result in different methods giving different results – for example, see Onions et al. (2012) on identifying sub-haloes near the centre of dark matter haloes, Behroozi et al. (2015) on major halo mergers, and Libeskind et al. (2018) on classifying elements of the cosmic web. This suggests that complementary methods for identifying structure, or structure formation, may be helpful.

We outline a possible scheme for doing this using kinetic-theory entropy creation. In dealing with observations or simulations, we will need to first construct a phase space distribution function (DF) from the data. Practically speaking, this will need involve some smoothing in phase space. We then also need to choose a scale of interest in position space for coarse-graining in order to help identify the structures we are aiming to explore. The coarse-graining scale must be at least as large, and might be much larger, than the scale for the practically-necessitated smoothing.

Given now our DF $f(1)$, smoothed for practical reasons, and then coarse-grained to our scale of interest, there are then two alternative approaches:

One approach is to rely solely on the DF, starting with the total entropy change at a point and subtracting the entropy flow of Eq. (19), giving

[TABLE] 2. 2.

The approach closer to that used for our model in Section 5 is to construct from the observations or simulation a two-point correlation function $g(1,2),$ coarse-grain to a scale of interest, and then apply Eq. (20),

[TABLE]

A comparison of methods for estimating correlations is given in Kerscher et al. (2000), albeit for a correlation function that depends only on the distance between two points, instead of depending, as here, on the position and velocity within phase space. The estimation will, of course, involve an element of practically-necessitated smoothing.

By analogy with our model – see the remark after Eq. (23) – it is possible that this scheme would characterise under-densities, for example cosmic voids or their centres, as undergoing structure formation. It is possible that, in some circumstances, transitional regions between over-dense structure and these under-densities might be the regions of maximum entropy creation.

For the model we constructed in Sections 2-5, there was no question of phase-mixing, as explained at the end of Section 3. However, once we apply our scheme to observations or a simulation, the practical requirement to smooth the particle distribution implies phase-mixing is a possibility, and so even truly collisionless microscopic processes might give rise to macroscopic entropy creation.

It is beyond the scope of this paper to identify if the above scheme would work in practice. This might depend upon the details of the particular simulation or set of observations, the approach to smoothing, and the choice of approach and scale for coarse-graining.

An obstacle to robust results – how challenging an obstacle is to be determined – could well arise from the need to identify differences between possibly relatively similar quantities. For approach 1, this requirement is explicit in the form of Eq. (48). For approach 2, it arises because of the requirement to calculate correlation functions. Feasibility of our scheme is therefore, perhaps most crucially, dependent on the size of simulation or observation errors relative to the underlying entropy creation and destruction.

The scheme is feasible and useful if, at whatever scale is focused upon, we can robustly detect entropy creation and/or destruction. It is also feasible, but presumably of less use, if we can robustly exclude entropy creation and/or destruction.

Our scheme’s coarse-graining might help draw out properties of different scales. For example, as is well known, dark matter is essentially collisionless, whereas dark matter haloes are mutually collisional. At scales relevant to dark matter haloes, would we see larger entropy creation? Varying the coarse-graining scale might also help elucidate interaction between various levels of hierarchical structure. If we wanted to include gas in our system, as the dominant baryonic matter component, we might also add terms to encompass hydrodynamic entropy creation in the gas.

7 Conclusion

As mentioned in the introduction, a core-halo model was described in Binney & Tremaine (2008, p572), making an artificial distinction between the core and the halo during gravitational collapse. A similar argument (Binney & Tremaine, 2008, pp377-378) draws out more clearly that the creation of entropy takes place predominantly within the halo, not the core. If we further assume that the core’s structure scales with its radius, then its phase space volume varies like $r^{\nicefrac{{3}}{{2}}}$ and it is easily seen that entropy is in fact destroyed within the (shrinking) core.

In the current paper, we have constructed an analytical kinetic theory perturbation model for the beginning of gravitational collapse. We introduced an asymptotic coarse-grained entropy, which in our model is associated with the system’s fastest-growing modes, and indicates the rate of their collisional relaxation.

Overall for our model, which is not an isolated system, we see net entropy destruction. However, this is a higher order, more suppressed, effect compared with a pattern of entropy destruction and creation. Entropy destruction occurs in a “core” around the central perturbation, with equal and opposite entropy creation in a “halo” extending for a finite radius beyond that core, as shown in Figure 2. The physical scale for the core-halo pattern depends on the coarse-graining parameter chosen: the coarser the graining, the bigger the physical scale.

In the core, collisions enhance the perturbation in a process of collisional “de-relaxation.” Conversely, in the halo, collisional relaxation suppresses the perturbation. In our linear perturbation model, the effect of such collisional evolution on the perturbation (and hence the entropy creation) is well defined, but small in size compared with its collisionless evolution (and the entropy flow).

A core-halo pattern of gravitational collapse, well known from simulations and observations, is generally set “by hand” in analytical models. As far as the author has been able to determine, this is the first time an analytical kinetic theory model has produced a core-halo pattern.

This motivates a scheme for measuring structure formation in observations or simulations, via patterns of entropy creation and destruction, as set out in Section 6. The feasibility of this scheme in the contexts of various observations or simulations is the key unanswered question arising from this paper. Because the main difficulty is likely to arise from the size of observation or simulation errors relative to entropy creation and destruction, feasibility is likely to improve along with the precision of observations and simulations.

Acknowledgement

The author is very grateful to the anonymous referee for his or her extremely helpful suggestions, including on the paper’s structure and style, and for prompting inclusion of more discussion of physical implications and of alternative coarse-grainings.

Appendix A Approximating the initial delta function perturbation by a Gaussian

As noted after Eq. (13), the formulation of $f_{1,\text{init}}$ via a delta function is strictly speaking not compatible with perturbation theory. This is dealt with by regarding the Dirac delta function as an approximation of a Gaussian in $\textbf{{x}}_{1},$

[TABLE]

for some relatively small width $w>0$ and, having fixed $w,$ then taking $\epsilon$ to be small enough to ensure perturbation theory works. Note that, with our conventions, the Fourier transform of the Dirac delta function of Eq. (13) is $1,$ while the Fourier transform of the Maxwellian of Eq. (50) which it approximates is $\mathrm{e}^{-{\nicefrac{{k_{1}^{2}\,w^{2}}}{{2}}}}\approx 1-{\nicefrac{{k_{1}^{2}\,w^{2}}}{{2}}}.$ So, for small $k_{1}w,$ the Fourier transform is effectively $1.$ We shall be most interested in wave-numbers $k_{1}<k_{\text{J}}\beta,$ for some $\beta\ll 1.$ If we are given $\beta,$ then we need only insist that $w\ll k_{\text{J}}^{-1}\beta^{-1},$ to get $k_{1}w\ll 1.$

From the growth with time of $\bar{f}_{1,\text{a}}$ (and hence $f_{1,\text{a}}$) set out in Eq. (37), and from Eqs. (11) and (50), we can see that for first order perturbation theory to be valid we need $\epsilon$ and/or $t$ to be of small enough order such that

[TABLE]

When considering distributions of entropy creation in space, as in Subsection 5.2, in order to maintain our Gaussian approximation we must require that radius $r$ considered satisfies $w\ll r\ll R,$ recalling that $R$ is the radius of the large volume $V.$ To also satisfy Eq. (51) again requires small enough $\epsilon$ and/or $t.$

Appendix B The plasma dispersion function

This appendix explores the properties of the plasma dispersion function defined by

[TABLE]

for $\operatorname{Im}z>0,$ and by analytic continuation for $\operatorname{Im}z\leq 0.$ It is easy to see that this relates to the definition of $Z(k,\omega)$ in Eq. (36) via $Z(k,\omega)\equiv Z({\nicefrac{{\omega}}{{\sqrt{2}k\sigma}}})\,.$ The key source for this appendix is Fried & Conte (1961), and see also Binney & Tremaine (2008, app. C.3). The plasma dispersion function can alternatively be defined by

[TABLE]

where $\operatorname{erf}$ is the usual error function, and $\operatorname{erfc}$ the usual complementary error function.

We saw definitions of $Y(k,\omega)$ and $P(k,\omega)$ in Eqs. (34) and (35), which immediately give related definitions of $Y(z)$ and $P(z)$ along the same lines as for $Z(k,\omega)$ and $Z(z).$ The notation $Z$ and $P$ is fairly standard and is used in Fried & Conte (1961), but the function $Y$ in Fried & Conte (1961) is different from ours. A related function, the Fadeeva function $w(z)={\nicefrac{{Z(z)}}{{\mathrm{i}\sqrt{\pi}}}}$ is often considered, and its properties are discussed in, for example, DLMF (2014, Ch. 7), which is the online companion to Olver et al. (2010).

B.1 Asymptotic series for small wave-numbers

We will be particularly interested in the properties of $Z(k,\omega)$ and the associated functions for small $k,$ that is the properties of $Z(z)$ for large $z.$ As $z\to\infty,$ we have an asymptotic series for $Z,$

[TABLE]

which can be derived from Eq. (52), and analytic continuation for $\operatorname{Im}z\leq 0,$ using standard results for the moments of the Gaussian. The approximation excludes the “tails” of the integral as $|s|\to\infty,$ to ensure that $|{\nicefrac{{s}}{{z}}}|<1$ for the series expansion. For large $|z|,$ the error is made small by the integral’s exponential function.

As mentioned, the series Eq. (54) is not convergent, but asymptotic. As a power series, it does not converge for any finite $z,$ because the ratio of a term over its predecessor is ${\left(j-{\nicefrac{{1}}{{2}}}\right)}/{z^{2}},$ which goes to infinity for any finite $z.$ The utility of the series arises because taking the first few terms of the series can give a very good approximation for large but finite $z.$ Heuristically, the first term which is not used in the approximation provides an estimate of the error. Asymptotic series are discussed in depth in, for example, Bender & Orszag (1999). From Eq. (54), we can also find similar asymptotic series for $Z(k,\omega),Y(k,\omega),$ and $P(k,\omega),$ for small $k>0.$

B.2 The plasma dispersion relation and its zeros

The plasma dispersion relation, or here simply dispersion relation, is the equation

[TABLE]

and, given $k\geq 0,$ we call a solution $\omega$ a dispersion zero. From the series in Eq. (54), for small $k>0$ and $\operatorname{Im}\omega>0,$ it can be found, order by order, that the (unique, Binney & Tremaine, 2008) dispersion zero with positive imaginary part is given by

[TABLE]

which is checked in Wren (2018).

The numerically-calculated values of $\eta$ used in Figure 1 are obtained in Wren (2018) by calculating $Z(z)$ using Eq. (53) for $-1<\operatorname{Im}z<1,$ and using a continued fraction method from Fried & Conte (1961) for $\operatorname{Im}z\geq 1,$ continuing the fraction using $20$ terms. For $\operatorname{Im}z\leq-1,$ the value of $Z(z)$ follows from the continued fraction method and use of the result (Binney & Tremaine, 2008, eq. C.25) that, for real $x$ and $y,$

[TABLE]

where $Z^{*}$ denotes the complex conjugate of $Z.$ (There is a typo omitting that conjugation of $Z$ in the corresponding formula in Fried & Conte, 1961.)

Binney & Tremaine (2008) discusses the dispersion zeros. As mentioned, given $k,$ there is at most one dispersion zero with positive imaginary part; this zero only occurs for $0<k<k_{\text{J}}.$ The only real dispersion zeros occur for $k=0$ and $k=k_{\text{J}}$ (and have value [math]), while, for any $k>0,$ there are an infinite number of Landau dispersion zeros with $\operatorname{Im}z<0$ . Let $\omega=|\omega|\left[\cos(\theta)-\mathrm{i}\sin(\theta)\right]$ be a Landau zero, where $0<\theta<\pi,$ and then using Eq. (54) and the dispersion relation Eq. (55), we have

[TABLE]

where we have assumed $k$ is fixed and so the final order term is written in terms of $\omega.$

For $0\leq\theta<{\nicefrac{{\pi}}{{4}}}$ or ${\nicefrac{{3\pi}}{{4}}}<\theta\leq\pi$ we have $-\cos(2\theta)<0,$ and so, when $|\omega|\gg k\sigma$ the value of the exponential in Eq. (58) is close to zero; on the other hand, for ${\nicefrac{{\pi}}{{4}}}<\theta<{\nicefrac{{3\pi}}{{4}}},$ when $|\omega|\gg k\sigma$ the value of the exponential in Eq. (58) is very large.555A line, such as $\theta={\nicefrac{{\pi}}{{4}}}$ or $\theta={\nicefrac{{3\pi}}{{4}}},$ with this kind of behaviour is often described as a Stokes line. Therefore, when $|\omega|\gg k\sigma,$ we must have $\theta$ close to either ${\nicefrac{{\pi}}{{4}}}$ or ${\nicefrac{{3\pi}}{{4}}}.$ We concentrate on $\theta\approx{\nicefrac{{\pi}}{{4}}},$ the other choice being very similar.

B.3 Landau zeros of relatively large size

In this subsection, we look at the case when we have Landau dispersion zeros of large size, that is where $|\omega|\gg k_{\text{J}}\sigma.$ This part of Appendix B motivates Appendix B.4’s key step in verifying that we can apply the residue formula, mentioned before Eq. (37), to the one-particle distribution function.

We assume $|\omega|\gg\sigma k,$ and treat $k>0$ (which is not assumed to be smaller than $k_{\text{J}}$) as fixed. Write $\phi={\nicefrac{{\pi}}{{4}}}-\theta,$ giving us $\omega=|\omega|\mathrm{e}^{\mathrm{i}\left(\phi-{\nicefrac{{\pi}}{{4}}}\right)},$ and then, from Eq. (58), we have

[TABLE]

Taking the phases, Eq. (59) implies

[TABLE]

where, in Appendix B only, $n$ is an integer, and we approximated $\cos(2\phi)=1+\operatorname{O}(\phi^{2}).$ Taking the absolute value of Eq. (59), we then have

[TABLE]

although the left-hand side of this equation is small, the factor of ${\nicefrac{{|\omega|^{2}}}{{\sigma^{2}k^{2}}}}$ in the exponent on the right-hand side requires that for the equality to hold we must have

[TABLE]

to avoid the exponential on the right-hand side of Eq. (61) being very large or very small. From Eq. (61), we then have that

[TABLE]

The results of Eqs. (62) and (63) are similar in form to expressions in DLMF (2014, eq. 7.13.4) for the zeros of the $\operatorname{erfc}$ operator, which occurs in Eq. (53). Numerical calculations in Wren (2018) confirm the accuracy of these approximations.

Together with the similar results for $\theta\sim{\nicefrac{{3\pi}}{{4}}},$ Eq. (63) demonstrates the well-known fact that the Landau zeros with negative imaginary parts and large absolute values have real parts of very slightly larger size than their imaginary parts, and that the difference in size between the real and imaginary sizes grows increasingly small as the absolute value of the Landau zero grows.

B.4 The residue approach for the inverse Laplace transform of the one-particle distribution function

We shall now confirm that the (Fourier- and) Laplace-transformed one-particle distribution function $\tilde{\bar{f}}_{1}$ of Eq. (33) can be inverse Laplace-transformed using the well-known formula

[TABLE]

where $p$ ranges over all poles $\omega_{p}$ of $h.$ The applicability of this well-known formula to $\tilde{\bar{f}}_{1}$ is assumed in similar contexts in, for example, Ichimaru (1973) and Binney & Tremaine (2008), but is not demonstrated in those texts. We treat $k>0$ as fixed, and, as in the previous subsection, do not need to assume that it is less than $k_{\text{J}}.$

Recall that the inverse Laplace transform of a function $h(\omega)$ is given by

[TABLE]

for $c>0$ such that $h(t)\mathrm{e}^{-ct}\to 0$ as $t\to 0.$ The standard argument used to justify Eq. (64) applies to an analytic function with poles, $h(\omega),$ that tends to zero as $|\omega|\to\infty$. Jordan’s Lemma then implies that the contour integral around a large semi-circle dropped from the straight line from $-X+\mathrm{i}c$ to $X+\mathrm{i}c$ also tends to zero, and therefore

[TABLE]

where $\mathcal{C}_{X}$ is the closed contour formed by the straight line from $-X+\mathrm{i}c$ to $X+\mathrm{i}c$ followed by the semi-circle dropped below. The condition for $c$ implies that any pole of $h$ lies within any contour $\mathcal{C}_{X}$ for large enough $X,$ and Eq. (64) then follows from the residue theorem.

Returning to the inverse Laplace transform of $\tilde{\bar{f}}_{1},$ recall that Appendix B.3 showed the Landau zeros extend out to infinity. For fixed wave-number and velocity, the first term of Eq. (33) clearly tends to zero as $|\omega|\to\infty,$ so we focus entirely on the second term, which, because of the Landau zeros, is not even bounded as $\omega\to\infty$. This lack of boundedness means we cannot directly apply the standard argument recalled above for $h.$ The residue formula Eq. (64) will none the less hold if we can construct a sequence of dropped-below semi-circles growing in radius to infinity, with (the second term of) $\tilde{\bar{f}}_{1}(\textbf{{k}},\textbf{{v}},\omega)$ tending to zero on that sequence of semi-circles. This follows by using the standard argument, but in Eq. (66) confining our attention to that sequence of semi-circles. As in the previous subsection, we will concentrate on dealing with Landau zeros having positive real parts – very similar steps deal with the conjugate Landau zeros with negative real parts, and we do not need to set those out explicitly. In the notation of Appendix B.4, we are assuming that $-{\nicefrac{{\pi}}{{4}}}\leq\phi\leq{\nicefrac{{3\pi}}{{4}}}.$

We have the second term of $\tilde{\bar{f}}_{1}(\textbf{{k}},\textbf{{v}},\omega)$ as

[TABLE]

For $\operatorname{Im}\omega\geq 0,$ that is ${\nicefrac{{\pi}}{{4}}}\leq\phi\leq{\nicefrac{{3\pi}}{{4}}}\,,$ we see from Eq. (54) that $Y(k,\omega)$ tends to zero as $\omega\to\infty,$ and $\omega\,Y(k,\omega)\to-1,$ so, Eq. (67)’s middle expression shows that the second term of $\tilde{\bar{f}}_{1}(\textbf{{k}},\textbf{{v}},\omega)$ tends to zero.

Consider next the case $-{\nicefrac{{\pi}}{{4}}}\leq\phi<0.$ From Eq. (54), we see that, for large $\omega,$ we then have $Y(k,\omega)$ dominated by the first, $\tau,$ term of Eq. (54), which is like $\mathrm{e}^{-\omega^{2}/(2k^{2}\sigma^{2})},$ with $\omega^{2}$ having negative real part. So $Y(k,\omega)$ tends exponentially to infinity as $\omega\to\infty,$ and therefore Eq. (67)’s final expression shows that the second term of $\tilde{\bar{f}}_{1}(\textbf{{k}},\textbf{{v}},\omega)$ again tends to zero as $\omega\to\infty.$

The remaining case is for $0\leq\phi\leq{\nicefrac{{\pi}}{{4}}}.$ For this case, we now construct a sequence of semi-circles on which, as $\omega\to\infty,$ the dispersion expression $k^{2}-k_{\text{J}}^{2}P(k,\omega)$ is bounded below. In constructing this sequence, we will choose the semi-circles to pass between the Landau zeros. We have from Eq. (58),

[TABLE]

We will call that final expression’s first summand, $\left({\nicefrac{{\sqrt{2\pi}k_{\text{J}}^{2}|\omega|}}{{\sigma k}}}\right)\,\exp\left[\cdots\right],$ the exponential summand. We can see that if $\sin(2\phi)>({2\sigma^{2}k^{2}}/{|\omega|^{2}})\ln({2\sqrt{2\pi}k_{\text{J}}^{2}|\omega|}/{k^{3}\sigma}),$ then the exponential summand’s absolute value is less than $k^{2}/2\,,$ implying that $|k^{2}-k_{\text{J}}^{2}P(k,\omega)|>k^{2}/2.$ We choose the radius of our semi-circles, $|\omega|,$ to be sufficiently large that the condition on $\sin(2\phi)$ only fails for $\phi$ small enough that $\cos(2\phi)\approx 1$ is a very close approximation. Note that, for such $\phi,$ we have $\phi\leq\operatorname{O}(\ln[|\omega|]\,|\omega|^{-2}).$

To deal with the small $\phi$ case, we now consider the phase from Eq. (68). Set our large $|\omega|=2\sigma k\sqrt{\pi\left(n-{\nicefrac{{1}}{{8}}}+{\nicefrac{{1}}{{2}}}\right)}=2\sigma k\sqrt{\pi\left(n+{\nicefrac{{3}}{{8}}}\right)}\,,$ where $n$ is a large positive integer. This puts each of our semi-circles roughly midway between two successive Landau zeros. Using $\cos(2\phi)=1+\operatorname{O}(\phi^{2}),$ we now have that the exponential summand’s phase is

[TABLE]

which, since $\phi$ is small, and $\ln[n]^{2}\,n^{-1}\to 0$ as $n\to\infty\,,$ implies the exponential summand’s real part is positive, and hence $|k^{2}-k_{\text{J}}^{2}P(k,\omega)|>k^{2}>k^{2}/2.$ We have therefore shown that, on semi-circles of radius $|\omega|=2\sigma k\sqrt{\pi\left(n+{\nicefrac{{3}}{{8}}}\right)},$ for sufficiently large $n,$ we have $|k^{2}-k_{\text{J}}^{2}P(k,\omega)|>k^{2}/2.$ The other factors of the second term of $\tilde{\bar{f}}_{1}(\textbf{{k}},\textbf{{v}},\omega)$ taken together tend to zero as $|\omega|\to\infty\,,$ so we have shown that the second term tends to zero on our sequence of semi-circles for $0\leq\phi<{\nicefrac{{\pi}}{{4}}}.$ We have now shown this in turn for all relevant cases – which were ${\nicefrac{{\pi}}{{4}}}\leq\phi\leq{\nicefrac{{3\pi}}{{4}}}\,,$ $-{\nicefrac{{\pi}}{{4}}}\leq\phi<0,$ and $0\leq\phi<{\nicefrac{{\pi}}{{4}}}$ – giving us the condition we noted in the paragraph before Eq. (67), and allowing us to apply the residue formula of Eq. (64) to $\tilde{\bar{f}}_{1}(\textbf{{k}},\textbf{{v}},\omega).$

We can calculate the residue associated with the dispersion relation in Eq. (33) using the well-known relation that, if $h(z)$ is any holomorphic function with a simple zero at $z=z_{0},$ meaning that $h^{\prime}(z_{0})\neq 0,$ then

[TABLE]

Assuming, $\omega\neq 0$ is a dispersion zero, we find

[TABLE]

where we used $Z^{\prime}(z)=-2-2zZ(z)$ from Binney & Tremaine (2008, eq. C.26), and the dispersion relation itself.

We now check that the dispersion poles are all simple for $k>0.$ The right-hand side of Eq. (71) can only be zero if $\omega^{2}$ is real. From the properties of the dispersion relation discussed in the paragraph after Eq. (57), this is only possible for $0\leq k\leq k_{\text{J}},$ and when $\omega$ is purely imaginary. As shown by the blue dotted line in Figure 1, the right-hand side of Eq. (71) only vanishes for $k=0.$ For $k\neq 0,$ we therefore have only simple poles, which means we can use Eq. (70) and also enables us to simplify Eq. (64) to

[TABLE]

where the sum ranges over all the poles of $\tilde{\bar{f}}_{1}.$ Asymptotically over time the fastest growing part of $\bar{f}_{1}$ is therefore that associated with the pole $\omega_{p}=\mathrm{i}\eta.$666It is possible for other terms to give initially faster-growing parts. For example, for $a,b>0,$ suppose $\omega=a-\mathrm{i}b$ is a Landau zero, and hence so is $-a-\mathrm{i}b.$ We can also approximate $a\approx b$ if $|\omega|$ is large. Pairing terms, we get time-dependence like $\mathrm{i}\,\mathrm{e}^{-at}\sin\left(at\right).$ It is easy to see that, for $at\ll 1,$ this quantity is fast growing, while for $t\gg{\nicefrac{{1}}{{a}}},$ it will be highly suppressed. This, with Eq. (71), gives us the result for $\bar{f}_{1}$ quoted in Eq. (37).

Appendix C Asymptotic coarse-grained number density and entropy density of the first order perturbation function

Integrating Eq. (37) with respect to $\textbf{{v}}_{1},$ we get the (Fourier-transformed) asymptotic number density function,

[TABLE]

where we used the dispersion relation, and the series approximation from Eq. (56). We now make our coarse-graining explicit and treat $\bar{n}_{1,\text{a}}$ as vanishing unless $0<k<k_{\text{J}}\beta,$ where, as usual, $\beta\ll 1.$ If we then inverse Fourier transform, we get

[TABLE]

This is shown in Figure 3 (Right). The pattern of the asymptotic coarse-grained number density (by volume) is a strong positive peak, with a much weaker tail. As described in the caption, Figure 3 (Left) shows the number density by shell.

We now consider the entropy density in space,

[TABLE]

and find that the Fourier-transformed asymptotically-dominant entropy density in space, to leading order in $\epsilon,$ is

[TABLE]

where the integration can be done by hand, and is also calculated using Mathematica in Wren (2018).

Because the the Fourier-transformed asymptotic entropy density in Eq. (76) and number density in Eq. (73) both have constant leading order in $k_{1}$, they are proportional, implying that the corresponding asymptotic coarse-grained quantities in position space are also proportional to leading order in $k_{\text{J}}\beta.$ At leading order therefore, the asymptotic coarse-grained entropy density pattern does not add any new information.

Appendix D Calculations for the zeroth order correlation function

In this appendix, we solve Eq. (40), which was derived from the zeroth order correlation equation. Fourier transforming Eq. (40) gives

[TABLE]

This can be solved by the ansatz $\bar{G}_{0}(\textbf{{k}}_{-},\textbf{{v}}_{1},\textbf{{v}}_{2})=\bar{q}(\textbf{{k}}_{-})\operatorname{\mathcal{M}}(\textbf{{v}}_{1})\operatorname{\mathcal{M}}(\textbf{{v}}_{2})\,,$ which gives us

[TABLE]

implying that our solution is of the form

[TABLE]

where $\lambda(\textbf{{k}}_{-})$ is an arbitrary function, and $C$ is an arbitrary constant which allows for the possibility that $\textbf{{k}}_{-}=\textbf{{0}}$ in Eq. (78). Viewing Eq. (78) as a differential equation in $\mathbf{r}=\textbf{{x}}_{2}-\textbf{{x}}_{1},$ the constant $C$ arises because Eq. (78) represents a total derivative, while the first two summands of Eq. (79)’s right-hand side represent, respectively, a particular integral and a complementary function.

Note that, since $f_{0}(1,2)=f_{0}(1)f_{0}(2)+({\nicefrac{{1}}{{N}}})g_{0}(1,2)$ to first order in ${\nicefrac{{1}}{{N}}},$ integrating over all $(2)$ gives us

[TABLE]

and we therefore must have $\int\!q(\mathbf{r})\,\text{d}^{3}\mathbf{r}=0.$ The first and second terms of $\bar{q}(\textbf{{k}}_{-})$ have well-defined values at $\textbf{{k}}_{-}=\textbf{{0}}.$ The integral over all space of their inverse Fourier transforms can therefore be evaluated by setting $\textbf{{k}}_{-}=\textbf{{0}}.$ The result is $-{\nicefrac{{1}}{{V}}}$. We cannot evaluate the third, $C,$ term at $\textbf{{k}}_{-}=\textbf{{0}},$ but we can see that its inverse Fourier transform is ${\nicefrac{{C}}{{(2\pi)^{3}}}}.$ Therefore to have $q$ integrating to zero over all space, or more precisely over the volume $V,$ we must have that $C={\nicefrac{{(2\pi)^{3}}}{{V^{2}}}}$ and so

[TABLE]

As the focus in this paper will be on small $k\ll k_{\text{J}},$ the function $\lambda(\textbf{{k}}_{-})$ will not affect our conclusions, so we disregard it for the remainder of this paper, and have the result set out in Eq. (41).

As an aside, we can find the inverse Fourier transform, $q,$ of $\bar{q},$ directly from inverse Fourier transforming Eq. (78) to get an equation which can be straightforwardly solved, along lines set out in Kandrup (1983), to get

[TABLE]

Appendix E The first order correlation equation

In this appendix, we Fourier transform Eq. (16), the first order correlation equation, with respect to $\textbf{{x}}_{1}$ and $\textbf{{x}}_{2}$ which gives us

[TABLE]

where, in the context of Fourier transforms, we read $(j)$ as $(\textbf{{k}}_{j},\textbf{{v}}_{j}).$ Using Eq. (39) and the $\bar{q}$ ansatz noted before Eq. (41), we have

[TABLE]

Recalling, as set out after Eq. (13), that our initial perturbation is uncorrelated, $g_{1}(1,2,t=0)=0,$ Laplace transforming Eq. (84) then gives Eq. (42).

Appendix F The first order correlation equation using a propagator

This appendix solves the BBGKY equation for the first order correlation function, or more precisely the related function $\gamma_{\text{a}}$ from Eq. (30), using a propagator approach, which, in essence, derives a Green’s function for the equation. This well-known approach is used in, for example, Ichimaru (1973) for plasmas and, in an “angle-action” approach different to ours, in Heyvaerts (2010) for self-gravitating particles.

F.1 Propagators

Note that the (Fourier-transformed) equation for the correlation function, given in Eq. (84), is of the form

[TABLE]

for $t>0.$ For clarity below, we have made the time argument of $\bar{\bar{g}}_{1}$ explicit, and we have written

[TABLE]

and similarly for $H_{2},$ and we have also collected all the other terms, none of which involve $\bar{\bar{g}}_{1},$ as a “driving” term $D_{1}$ on the right-hand side of Eq. (85),

[TABLE]

for $t>0,$ and zero otherwise. Note that $\textbf{{k}}_{1}$ and $\textbf{{k}}_{2}$ simply parametrise Eqs. (85) and (86), whilst, in contrast, there is differentiation and integration with respect to the velocity variables, which play the active role in the equations. For brevity, where needed we will write variables $1^{\prime}\equiv(\textbf{{k}}_{1},\textbf{{v}}_{1}^{\prime})$ and $2^{\prime}\equiv(\textbf{{k}}_{2},\textbf{{v}}_{2}^{\prime}).$

To solve Eq. (85) with a propagator, we adopt the standard approach of writing

[TABLE]

with, for $j=1,2,$

[TABLE]

and the initial condition $\mathcal{G}(j,\textbf{{v}}_{j}^{\prime},0)={\delta^{(3)}}(\textbf{{v}}_{j}-\textbf{{v}}_{j}^{\prime}).$ This equation has exactly the same form as the Fourier-transformed Vlasov equation for $\bar{f}_{1}(\textbf{{k}}_{j},\textbf{{v}}_{j},t)\,,$ the only difference being the initial condition. We now solve Eq. (89), using the same approach that took us to Eq. (32) and then Eq. (33) to find

[TABLE]

The technical discussion of Appendix B.4, shows that, as for $\tilde{\bar{f}}_{1},$ we can apply Eq. (72), to get, for $0<k_{j}<k_{\text{J}},$ that the asymptotically-dominant term of $\mathcal{G}$ is

[TABLE]

where, just as for $\bar{f}_{1,\text{a}},$ we used Eqs. (70) and (71) to find the square bracketed factor.

F.2 Time dependence of distribution and correlation functions

We now want to look at the dominant time dependence of distribution functions and correlation functions. Throughout this appendix, we assume all functions of a time argument vanish when that time is less than $0,$ and consider only their behaviour for times greater than or equal to $0.$ We write $y\preccurlyeq z$ to mean that the proportionate growth of $y$ with time is no faster than the proportionate growth of $z$ with time (for times greater than or equal to [math]). For example, $t\preccurlyeq\mathrm{e}^{t}$ or $\mathrm{e}^{t}\preccurlyeq\mathrm{e}^{2t},$ but also $\mathrm{e}^{t}\preccurlyeq\mathrm{e}^{t},$ and, because we ignore constants of proportionality, $\mathrm{e}^{t}\preccurlyeq 10\,\mathrm{e}^{2t},$ and even $10\,\mathrm{e}^{t}\preccurlyeq\mathrm{e}^{t}.$ Note also that $${10}^{100}$t\preccurlyeq\mathrm{e}^{t},$ because we are comparing rates of growth, not absolute size. Factors without a time dependency can be omitted.

We show that $\bar{\bar{g}}_{1}(1,2,t)\preccurlyeq\mathrm{e}^{2k_{\text{J}}\sigma\,t}.$ From the discussion of the dispersion relation after Eq. (36), and the detailed discussion in Appendix B.2, we see that $\mathcal{G}(\textbf{{k}}_{j},\textbf{{v}}_{j},\textbf{{v}}_{j}^{\prime},t)\preccurlyeq\mathrm{e}^{\eta(k_{j})\,t},$ where, as before, $\eta(k_{j})$ is the unique positive imaginary part of a dispersion zero for wave-number $0<k_{j}<k_{\text{J}}.$ For a fixed $k_{j},$ from Eq. (37) we have $\bar{f}_{1}(\textbf{{k}}_{j},\textbf{{v}}_{j},t)\preccurlyeq\mathrm{e}^{\eta(k_{j})\,t}.$ From the Fourier-transformed first order correlation equation Eq. (84), we therefore have

[TABLE]

the last relation following because $\eta(k_{j})\leq k_{\text{J}}\sigma.$

Implicitly assuming still that all functions of time vanish for negative times, we have

[TABLE]

So, from Eqs. (88) and (93), we have, for wave-numbers $k_{j}$ in the range $0<k_{j}<k_{\text{J}},$

[TABLE]

where in the second line we omitted all the factors without a time dependency.

We therefore have that the relative growth with time $t>0$ of the first order correlation function, $\bar{\bar{g}}_{1}(1,2,t)$, is no faster than $\mathrm{e}^{2k_{\text{J}}\sigma\,t}.$ If $\eta(k_{1})>{\nicefrac{{k_{\text{J}}\sigma}}{{2}}}$ and $\eta(k_{2})>{\nicefrac{{k_{\text{J}}\sigma}}{{2}}},$ which will be the case for small wave-numbers, then we can also see from Eq. (94) that $\bar{\bar{g}}_{1,\text{a}}(1,2,t),$ for a given $k_{1}$ and $k_{2}$ will grow like $\mathrm{e}^{\left[\eta(k_{1})+\eta(k_{2})\right]\,t}.$

F.3 Addressing the correlation equation using propagators

From Eq. (88), and the discussion immediately following Eq. (94), we have that, for $\eta(k_{1})>{\nicefrac{{k_{\text{J}}\sigma}}{{2}}}$ and $\eta(k_{2})>{\nicefrac{{k_{\text{J}}\sigma}}{{2}}},$ the asymptotically-dominant part of $\bar{\bar{g}}_{1}$ is given by

[TABLE]

where the asterisk denotes a time convolution. For brevity, we now write $E=\eta(k_{1})+\eta(k_{2}).$ Note that, for $\omega$ with $\operatorname{Im}\omega>E,$ the Laplace transform of $\mathrm{e}^{Et}$ is

[TABLE]

This means that the convolution at the end of Eq. (95) is an inverse Laplace transform

[TABLE]

Clearly ${\nicefrac{{\mathrm{i}}}{{\left(\omega-\mathrm{i}E\right)}}}$ is bounded as $\omega\to\infty.$ From the discussion of Appendix B.4, given $\textbf{{k}}_{1}$ and $\textbf{{k}}_{2},$ then for $\textbf{{k}}_{j}=\textbf{{k}}_{1},\textbf{{k}}_{2}$ or $\textbf{{k}}_{+},$ $\tilde{\bar{f}}_{1}(\textbf{{k}}_{j},\textbf{{v}},\omega)$ is bounded on a sequence of semi-circles. Since $\tilde{D}_{1}$ is the sum of a finite number of terms each depending on $\tilde{\bar{f}}_{1},$ we can apply the residue formula for the inverse Laplace transform, Eq. (72), to $\tilde{D}_{1}$ and therefore to ${\nicefrac{{\mathrm{i}\,\tilde{D}_{1}(\textbf{{k}}_{1},\textbf{{v}}_{1}^{\prime},\textbf{{k}}_{2},\textbf{{v}}_{2}^{\prime},\mathrm{i}E)}}{{\left(\omega-\mathrm{i}E\right)}}}$ overall. Our assumption that $\eta(k_{1})>{\nicefrac{{k_{\text{J}}\sigma}}{{2}}}$ and $\eta(k_{2})>{\nicefrac{{k_{\text{J}}\sigma}}{{2}}},$ implies that $E>k_{\text{J}}\sigma,$ and so the asymptotically-dominant term comes from the simple pole at $\omega=\mathrm{i}\,E,$ where the residue is $\mathrm{i}.$ So, from the residue formula, we have that the asymptotically-dominant term of the convolution is

[TABLE]

Note that the time-dependency here comes entirely from the ${\nicefrac{{\mathrm{i}}}{{\left(\omega-\mathrm{i}E\right)}}}$ factor on the right-hand side of Eq. (97) – this means that the driving term enters into Eq. (98)’s right-hand side solely through its Laplace transform $\tilde{D}_{1}(\omega),$ and the question of whether to restrict attention to the asymptotically-dominant part of $D_{1}(t)$ does not arise.

Referring back to Eq. (95), we now have

[TABLE]

We now need to evaluate the velocity integrals in Eq. (99). To do this we use the Laplace transform of from Eq. (87) and drop the delta function terms, which are of order $V^{-1},$ while the remaining terms are order $V^{0}.$

Using the expression for $\tilde{\bar{f}}_{1}$ in Eq. (33), we now have

[TABLE]

where we handled the terms corresponding to velocity derivatives of $\tilde{\bar{f}}_{1}$ via integration by parts. To help in doing the velocity integrals, write

[TABLE]

We also have a useful relation

[TABLE]

Doing the velocity integrations from Eq. (100), we get

[TABLE]

We can also recall an expression for $\bar{q}$ from Eq. (41).

From Eq. (95), we have

[TABLE]

Using the dispersion relation, this implies that

[TABLE]

The series expansion of $I_{D_{1}}$ for small $k_{j}$ is obtained from Eq. (103) via Mathematica computer algebra in Wren (2018), making use of the dispersion relation, and $Z^{\prime}(z)=-2-2zZ(z)$ from Binney & Tremaine (2008, eq. C.26). Evaluating Eq. (105) in Wren (2018), we then find that, for small $k_{j},$ we have

[TABLE]

where on the right-hand side we replaced $\textbf{{k}}_{2}$ by $\textbf{{k}}_{+}-\textbf{{k}}_{1},$ which puts the result in the form most helpful for use in Appendix H.

Appendix G The Landau approach to the first order correlation equation

To provide an alternative to, and to check, the propagator method set out in Appendix F for finding $\gamma_{\text{a}},$ Eq. (42), this appendix takes an approach based on that used in Landau (1946), and Subsection 4.1, to derive Eq. (33) for $\tilde{\bar{f}}_{1}.$ Along those lines, we solve Eq. (42) by rearranging the equation and integrating it with respect to velocities. It is lengthier than the derivation of $\tilde{\bar{f}}_{1},$ requiring calculation of a number of integrals through solving a set of seven simultaneous equations.

We integrate Eq. (42) with respect to both $\textbf{{v}}_{1}$ and $\textbf{{v}}_{2}$ to get

[TABLE]

where we have written

[TABLE]

where $\tilde{D}_{1}(1,2)$ is the Laplace transform of the driving term as set out in Eq. (87).

Now divide Eq. (42) by $1-\frac{\textbf{{k}}_{1}\boldsymbol{\cdot}\textbf{{v}}_{1}+\textbf{{k}}_{2}\boldsymbol{\cdot}\textbf{{v}}_{2}}{\omega}$ and integrate with respect to $\textbf{{v}}_{2},$ to get

[TABLE]

The leading order of the integral involving the driving term in $k_{j}$ is $-1,$ which implies this is also the leading order of $\tilde{\gamma}(1,\textbf{{k}}_{2}).$ As Eq. (109) suggests, consider $\tilde{\gamma}(1,\textbf{{k}}_{2})$ as a power series in $\textbf{{k}}_{1}\boldsymbol{\cdot}\textbf{{v}}_{1},$ times $\operatorname{\mathcal{M}}(\textbf{{v}}_{1}).$ From this power series, and the property that integration of powers of $\textbf{{k}}_{1}\boldsymbol{\cdot}\textbf{{v}}_{1}$ times the Maxwellian vanishes for odd powers, it can be seen that $\tilde{\Gamma}^{(0)}$ and $\tilde{\Gamma}^{(1)}$ are both of order zero in $k_{j},$ while $\tilde{\Gamma}^{(2)}$ and $\tilde{\Gamma}^{(3)}$ are both of order two in $k_{j}$ and $\tilde{\Gamma}^{(4)}$ and $\tilde{\Gamma}^{(5)}$ are both of order four.

Write

[TABLE]

Multiplying Eq. (109) by $\textbf{{k}}_{1}\boldsymbol{\cdot}\textbf{{v}}_{1}$ and integrating with respect to $\textbf{{v}}_{1},$ we have

[TABLE]

Now multiplying Eq. (109) by $(\textbf{{k}}_{1}\boldsymbol{\cdot}\textbf{{v}}_{1})^{2}$ and integrating with respect to $\textbf{{v}}_{1},$ we get

[TABLE]

Multiplying Eq. (109) by $(\textbf{{k}}_{1}\boldsymbol{\cdot}\textbf{{v}}_{1})^{3}$ and integrating with respect to $\textbf{{v}}_{1},$ we get

[TABLE]

We now solve simultaneously Eq. (107), Eqs. (111)-(113), and Eqs. (111)-(113) with $\textbf{{k}}_{1}$ and $\textbf{{k}}_{2}$ swapped, making seven equations in total. We write

[TABLE]

where

[TABLE]

[TABLE]

noting that the first entry in $\boldsymbol{\nu}$ is $\mu^{(0)}$ rather than $\nu^{(0)},$ and M is then defined via terms on the right-hand side of Eqs. (107) and (111)-(113).

Write $\textbf{{L}}=\omega\,\textbf{{1}}-\textbf{{M}},$ where 1 is the $7\times 7$ identity matrix, and we then have

[TABLE]

giving us, in particular, $\tilde{\Gamma}^{(s)}(\textbf{{k}}_{2},\textbf{{k}}_{1})$ for $s=0,..3.$ In Wren (2018), Mathematica is used to make find $\boldsymbol{\Gamma}$ from Eq. (117), and to substitute its components back into Eq. (109). The residue at $\omega=\mathrm{i}\,\eta(k_{1})+\mathrm{i}\,\eta(k_{2}),$ is calculated, noting that the driving function $\tilde{D}_{1}(1,2)$ has residue zero at that value of $\omega.$ We now have a closed equation for $\gamma_{\text{a}}(1,\textbf{{k}}_{2},t),$ explicitly to second order in $k_{j},$ which, after some manipulation in Wren (2018), gives the same result as obtained in Appendix F’s Eq. (106).

Appendix H Entropy creation rate calculations

The first two parts of this appendix cover detailed calculations needed for, respectively, Subsection 5.1 and Subsection 5.2 in the main text. The third part checks consistency between total and local entropy creation calculations. The fourth part considers variants to the paper’s main model, as mentioned in Subsection 5.3.

H.1 Calculating the total entropy creation rate

Recall Eq. (29) for the asymptotic coarse-grained entropy creation rate. An expression for $\gamma_{\text{a}}(-\textbf{{k}}_{+},\textbf{{v}}_{1},\textbf{{k}}_{2})$ can be obtained from Eq. (106) as

[TABLE]

by mapping $(\textbf{{k}}_{1},\textbf{{k}}_{2})\mapsto(-\textbf{{k}}_{+},\textbf{{k}}_{2}),$ which also induces the mapping $\textbf{{k}}_{+}\mapsto-\textbf{{k}}_{1}.$ From Eqs. (43) and (56), the factor involving a velocity derivative is

[TABLE]

The $\textbf{{v}}_{1}$ integral in Eq. (29) is evaluated in Wren (2018), effectively using identities for a multivariate normal distribution, giving

[TABLE]

The terms in the braces on the first line of Eq. (120)’s final expression are of order $-2$ in $k_{j},$ while the other terms in the braces are of order $0.$ Taking account of the $\text{d}^{3}\textbf{{k}}_{1}\,\text{d}^{3}\textbf{{k}}_{2}$ factors, the overall integral is therefore of leading order $4$ in $k_{j}.$ In principle, our approach of not accounting for entropy exterior to the volume $V$ means that we should set a lower limit for $k_{j}$ in the integrals, as well as the upper limit implied by the region $\mathcal{K}({\beta k_{\text{J}}}).$ However, because the integral is order four in $k_{j}$ and $R^{-1}\ll k_{\text{J}}\beta,$ we can safely omit this lower limit with negligible effect on the final answer.

The order $-2$ integrand gives

[TABLE]

where the final result of [math] follows by swapping the variables $\textbf{{k}}_{1}$ and $\textbf{{k}}_{2}$ in the integral, noting that this preserves $\mathcal{K}({\beta k_{\text{J}}}).$

From Eq. (120), we can therefore write

[TABLE]

where we took advantage of the integrand being homogeneously of order [math] in $k_{j},$ and $\textbf{{k}}_{j}$ was scaled by $k_{\text{J}}\beta$ to become dimensionless. The variable $\theta$ is the angle $\textbf{{k}}_{2}$ makes with $\textbf{{k}}_{1}$: the factor associated with $k_{2}^{2}$ is then $2\pi$ rather than $4\pi.$ Note that we can write the pre-factors before the integral in Eq. (122) as

[TABLE]

where $N_{1}$ and $B$ were defined after Eq. (44).

Assuming that $t$ is sufficiently small that we are willing to make the approximation $\left[\eta(k_{1})+\eta(k_{+})+\eta(k_{2})\right]t\approx 3k_{\text{J}}\sigma\,t,$ the exponential can be factored out from the integral. The resulting integral is evaluated numerically in Wren (2018), getting a result of $-0.0115541$, with an estimated error of $$9.94634\text{\times}{10}^{-6}$.$ The overall result for ${\nicefrac{{\mathrm{d}\,S_{\text{acg}}}}{{\mathrm{d}t}}}$ is shown in Eq. (44).

Note that the approximation $\left[\eta(k_{1})+\eta(k_{+})+\eta(k_{2})\right]t\approx 3k_{\text{J}}\sigma\,t,$ no longer holding at very late times is not a fundamental difficulty – it would be straightforward to calculate a time-dependent entropy creation formula which explicitly accounts for accurate values of $\eta(k_{j})\,t$ in the numerical integration. It is also straightforward to estimate the maximum error from the approximation $\eta(k_{j})\,t=k_{\text{J}}\sigma\,t,$ using the values of $\eta(k_{\text{J}}\beta)$ as calculated in Wren (2018) for drawing Figure 1. Table 1 confirms that, for values of $k\ll k_{\text{J}},$ the approximation $\eta(k)\,t=k_{\text{J}}\sigma\,t$ is very accurate until extremely late times $t\gg{\nicefrac{{1}}{{k_{\text{J}}\sigma}}},$ by which time, following Eq. (51), $\epsilon$ will need to have been very small indeed for the perturbative regime to remain valid.

As mentioned at the end of Subsection 5.2, we can try to alter the definition of asymptotic coarse-grained entropy, by relaxing the constraints that define the region $\mathcal{K}({\beta k_{\text{J}}}).$ That region sets three constraints, which may be summarised as $0<k_{1},k_{2},k_{+}<k_{\text{J}}\beta.$ In Wren (2018) are calculations analogous to those in this part of Appendix H, but for the cases where only two of the three upper constraints have effect. Because, for example, $k_{+}\leq k_{1}+k_{2},$ all three wave-numbers will still be small, allowing our analytical approximations. The results are as for Eq. (44), but with $-0.0116$ replaced by $-1.54\text{\times}{10}^{-3}$ if we only have the constraints $0<k_{1},k_{2}<k_{\text{J}}\beta,$ or by -0.0240$\beta^{-2}$ if we only have the constraints $0<k_{1},k_{+}<k_{\text{J}}\beta,$ or by 0.0239$\beta^{-2}$ (that is, positive entropy creation) if we only have $0<k_{2},k_{+}<k_{\text{J}}\beta.$ (Estimated errors in numerical integration are around 1\text{\times}{10}^{-5}$,$1\text{\times}{10}^{-4}$\beta^{-2},$ and 1\text{\times}{10}^{-5}$\beta^{-2},$ respectively.) The factors of $\beta^{-2}$ arise when the order $-2$ integrand of Eq. (121) no longer vanishes on integration, because the integration limits are no longer symmetric in $k_{1}$ and $k_{2},$ corresponding to having differing coarse-graining approaches for those two wave-numbers.

The wave-number $k_{1}$ corresponds to the entropy creation’s physical location. The choice above which leads to positive total net entropy creation, of directly constraining only $0<k_{2},k_{+}<k_{\text{J}}\beta,$ leads to only an indirect constraint on $k_{1}$ from $\textbf{{k}}_{1}=\textbf{{k}}_{+}-\textbf{{k}}_{2}$ and is therefore physically a particularly contrived form of coarse-graining with respect to $x_{1}.$

We prefer applying all three upper constraints, as better reflecting coarse-graining that might be done for practical observational or simulation reasons, or to study a particular scale. The symmetry in treatment of $k_{1},k_{2}$ and $k_{+}$ also better reflects asymptotic behaviour over time, as seen, for example, in the exponential factor of Eq. (120). Other choices of region $\mathcal{K},$ including the approaches above which directly constrain only two of the three wave-numbers, are therefore inappropriate for defining asymptotic entropy creation. We also note another significant advantage of our choice of $\mathcal{K}$ at the end of Appendix H.2.

Another approach chooses a coarse-graining which matches the asymptotic behaviour yet move closely: for small $k_{j}$ it captures essentially all wave-numbers which asymptote faster than a given rate. The asymptotic rate is given by the $\eta(k_{1})+\eta(k_{2})+\eta(k_{+})$ factor in Eq. (120)’s exponential. From Eq. (56) this is approximated to second order in $k_{j}$ by $3k_{\text{J}}\sigma-{3\sigma(k_{1}^{2}+k_{2}^{2}+k_{+}^{2})}/{2k_{\text{J}}},$ and, as shown in Table 1, this approximation is very good for all $k_{j}\leq 0.1k_{\text{J}}.$ We choose our coarse-graining to be defined by $0<k_{1}^{2}+k_{2}^{2}+k_{+}^{2}<2k_{\text{J}}^{2}\beta^{2},$ with the choice of $2$ as a factor for the constraint being somewhat arbitrary, and being made as it gives a similar overall scale to our original constraint $0<k_{1},k_{2},k_{+}<k_{\text{J}}\beta\,,$ for example capturing a similar volume of $k_{j}$ space (Wren, 2018). The part of the total net entropy associated with the order $-2$ integrand of Eq. (121) once again vanishes, because the constraint is symmetrical in $\textbf{{k}}_{1}$ and $\textbf{{k}}_{2}.$ As calculated in Wren (2018), this coarse-graining $0<k_{1}^{2}+k_{2}^{2}+k_{+}^{2}<2k_{\text{J}}^{2}\beta^{2}$ produces a total net entropy creation rate as in Eq. (44), but with $-0.0116$ replaced by -0.0125$,$ the latter being subject to an estimated integration error of around 1\text{\times}{10}^{-5}$.$

H.2 Calculating the distribution of the entropy creation rate over space

The distribution of the entropy creation rate over space is given by Eq. (46). We now evaluate that equation. We calculated the velocity integral in Eq. (120), expanding explicitly to order [math] in $k_{j}.$ Following that calculation, although there was also an order $-2$ term, Eq. (121) showed that it vanished on integration by $\textbf{{k}}_{1}$ and $\textbf{{k}}_{2}.$ However, this need not now be the case for Eq. (46), because of the new role of $\textbf{{k}}_{0}$ and the $\operatorname{sinc}$ factor disrupting the symmetry which ensured this term vanished in the previous subsection. We therefore look at the product of Eq. (118) and Eq. (119), in the latter substituting $\textbf{{k}}_{1}\mapsto\textbf{{k}}_{01},$ to get get the integrand corresponding to that of Eq. (120). Doing the velocity integral only to order $-2$ in $k_{j},$ we find

[TABLE]

We chose to pull $(8\pi^{4})^{-1}$ out of the integral in order to get our first factor of the same form as the corresponding term in Eq. (44), choosing the sign to make positive (resp. negative) values of our integral correspond to entropy creation (resp. destruction). We also kept one factor of $\beta$ inside the integral, to ensure the integral is a dimensionless function $\hat{S}_{\text{acg}}^{\circ}$ of the dimensionless quantity $k_{\text{J}}\beta r,$ and of order zero in $\beta$ (after the integration). The numerically-calculated function $\hat{S}_{\text{acg}}^{\circ}$ will be called the (leading order) entropy-creation pattern function, and is evaluated in Wren (2018) for varying values of $k_{\text{J}}\beta r.$ The overall result is shown in Eq. (47) and Figure 2. If we integrate ${\nicefrac{{\partial^{2}S_{\text{acg}}^{\circ}(r)}}{{\partial t\,\partial r}}}$ with respect to $r$ to get entropy creation in a generic region, we get a result of order $\beta^{4},$ compared with order $\beta^{6}$ for the total net entropy creation over all space. In other words, the total net entropy creation is suppressed by a factor of $\beta^{2}\ll 1$ compared with local entropy creation in a generic region.

Calculations in Wren (2018) confirm that the absolute values of $\hat{S}_{\text{acg}}^{\circ}$ integrated over either the core or halo are of very similar size, 0.0110$,$ with the core and halo values agreeing to around $0.01\%.$ This similarity is to be expected, as we know from Eq. ([121](#A8.E121)) that the total net entropy creation over both must essentially vanish. Note that $0.0110$ is very close in absolute value to -0.0116$,$ which is the corresponding factor for the total net entropy creation over all space, reinforcing the point made in the previous paragraph that the core-halo pattern is much more prominent than the total net entropy creation.

We now look again at the alternative definitions of asymptotic coarse-grained entropy creation which were considered at the end of Appendix H.1. Comparing Eqs. (29) and (46), the three constraints which now apply to our basic definition are $0<k_{01},k_{2},k_{12}<k_{\text{J}}\beta.$ However, if we now relax any one of the upper constraints, then at least one of $k_{01},k_{2},k_{12}$ will be unconstrained above. For example, if we relax the constraint that $k_{01}<k_{\text{J}}\beta,$ then there is nothing to prevent $k_{01}$ taking any large value, in this case because $0<k_{2},k_{12}<k_{\text{J}}\beta$ provides no constraint on $k_{0}.$ This means that we cannot apply the approach of this paper, which is based on analytically-tractable series for $\gamma_{\text{a}}(-\textbf{{k}}_{12},\textbf{{v}}_{1},\textbf{{k}}_{2})$ and $\bar{f}_{1,\text{a}}(\textbf{{k}}_{01},\textbf{{v}}_{1}),$ derived from expansion in their $\textbf{{k}}_{j}$ arguments for small wave-number. This is the further reason, referred to at the end of Appendix H.1, for preferring our definition of asymptotic coarse-grained entropy creation.

The coarse-graining $0<k_{1}^{2}+k_{2}^{2}+k_{+}^{2}<2k_{\text{J}}^{2}\beta^{2}$ (which corresponds to $0<k_{01}^{2}+k_{2}^{2}+k_{12}^{2}<2k_{\text{J}}^{2}\beta^{2}$ in Eq. (124)’s spatial distribution equation) also gives us a tractable analytical treatment for small wave-numbers. It yields a spatial distribution as in Eq. (47), but with the resulting entropy pattern function, calculated in Wren (2018), being $S_{\text{acg,}2}^{\circ}$ of Figure 4. Another “taxicab” coarse-graining, which constrains $k_{1}+k_{2}+k_{+},$ is additionally explored in Wren (2018). This again gives negative total net entropy creation, at next-to-leading order, and a leading order pattern of an entropy-destroying core, and an entropy-creating halo. Obscured by estimated integration error, there appear to be smaller-amplitude outer shells of entropy destruction and creation beyond the halo.

H.3 Checking consistency between the total entropy creation and its distribution

We can also follow a similar route to that of Appendix H.2 to calculate the next-to-leading order term for the distribution of entropy creation. This follows from taking the order [math] terms of Eq. (120), instead of the order $-2$ terms we considered above. It can be found that, at this next-to-leading order, we have

[TABLE]

where $\hat{S}_{\text{acg}}^{\circ\text{ntl}}$ is a dimensionless next-to-leading order entropy-creation pattern function, which is shown in Figure 5 as numerically calculated in Wren (2018). Note that, when integrated over a generic region, ${\nicefrac{{\partial^{2}S_{\text{acg}}^{\circ\text{ntl}}(r)}}{{\partial t\,\partial r}}}$ is of the same order in $\beta$ as the total net entropy creation over all space, and suppressed by a factor of $\beta^{2}$ compared with its leading order equivalent.

From Eq. (44), integrating Eq. (125)’s shell density over all radii should give us

[TABLE]

in order to ensure consistency between those equations. The size of the error bars in Figure 5 (Left) might give us some pause about using our numerical calculations to do the integration in Eq. (126). None the less, as calculated in Wren (2018), approximating the integral by summing and appropriately scaling the values shown in Figure 5 (Left) gives a result of $-0.0115,$ surprisingly close to $-0.0116.$ By doing this for a range of radii including only the central sphere of destruction and the surrounding shell of creation we also get a result of $-0.0115,$ suggesting that there is essentially no net entropy destruction outside the sphere and that innermost shell.

H.4 Modifications to the main model

A simple variant of our model is to allow the initial perturbation $f_{1,\text{init}}$ to have a different Maxwellian parameter $\sigma_{1}$ from the parameter $\sigma$ of the underlying perturbation $f_{0}.$ Following the approach of Subsection 4.1, we find, compare with Eq. (33),

[TABLE]

where the $\operatorname{\mathcal{M}}_{1}$ and $Y_{1}$ (only) are defined using $\sigma_{1}$ rather than $\sigma.$ Note that the residue with largest positive imaginary part still comes from the dispersion relation $k_{1}^{2}-k_{\text{J}}^{2}P(k_{1},\omega)=0,$ dependent on $\sigma.$ Hence, we have, compare with Eq. (37),

[TABLE]

with $\sigma_{1}$ only entering in through $Y_{1}(k_{1},\mathrm{i}\,\eta(k_{1})).$ By considering the series of Eq. (54) for $z={\nicefrac{{\mathrm{i}\,\eta(k_{1})}}{{\sqrt{2}k_{1}\sigma_{1}}}}\approx{\nicefrac{{\mathrm{i}\,k_{\text{J}}\sigma}}{{\sqrt{2}k_{1}\sigma_{1}}}},$ we can see that if $\beta\sigma_{1}\ll\sigma,$ then that $z$ will be large, and we can safely apply this paper’s small $k$ approximation approach.777If this $\beta\sigma_{1}\ll\sigma$ condition fails, then $\sigma_{1}/\sigma$ must be very large, and our perturbation is very different from a point-like perturbation as usually understood – the perturbing particles’ typical velocity is such that it more resembles an explosion from a point. Calculations in Wren (2018), following the approaches of Appendices F or G, and then Appendices H.1 and H.2, show that, at leading order in $\beta,$ we get a core-halo pattern exactly as in Eq. (47), with the same pattern function shown in Figure 2. In particular, there is no $\sigma_{1}$ dependency at leading order.

As for the main, $\sigma_{1}=\sigma,$ model, the total net entropy creation is at a higher order in $\beta$ than the core-halo pattern. Calculations in Wren (2018) give the total net entropy creation as in Eq. (44), but with the replacement

[TABLE]

The estimated integration errors for each of the two numerical coefficients on the right-hand side are around 1\text{\times}{10}^{-5}$.$ As we would expect, if we set $\sigma_{1}=\sigma$ in the right-hand side of Eq. ([129](#A8.E129)), we recover its left-hand side, -0.0116$.$ Note that we always get negative total net entropy creation. The requirement that $\beta\sigma_{1}\ll\sigma$ implies that the size of the total net entropy creation remains much less than the size of the entropy destruction in the core (or the size of its creation in the halo).

Taking the limit $\sigma_{1}\to 0$ represents an initial perturbation with all its particles initially stationary.888It is also straightforward to derive the equivalent of Eq. (128) for ${\bar{f}_{1,\text{init}}}(\textbf{{k}}_{1},\textbf{{v}}_{1})={\delta^{(3)}}(\textbf{{v}}_{1}),$ and then see, from the power series expansion of $Y_{1}(k_{1},\omega)$ along the lines of Eq. (54), that this does indeed correspond to the limit $\sigma_{1}\to 0.$ Since local leading order entropy creation is independent of $\sigma_{1},$ we still have exactly the same core-halo equation and pattern as for our main model. The next-to-leading order entropy creation is then given by

[TABLE]

where $S_{\text{acg}}^{\circ\text{\,stat\,ntl}}$ is a dimensionless next-to-leading order entropy-creation pattern function for the initially stationary perturbation, which is shown in Figure 6 as numerically calculated in Wren (2018).

We can also consider varying the perturbation’s initial correlation function $g_{1}(1,2,t=0),$ which was assumed to vanish. Suppose the initial perturbation correlation is altered to be the same as for the underlying DF, so $g_{1}(1,2,t=0)=g_{0}(1,2,t=0)$, the latter being set out in, and just before, Eq. (41). The rule for Laplace transforming time derivatives implies this adds a new term, a constant times $g_{0}(1,2,t=0),$ to the Laplace transform $\tilde{D}_{1},$ as found in Eq. (100). From Eq. (100), the dispersion relation, and Eq. (41), it can be seen that the resulting additional term in $I_{D_{1}}$ is of leading order $4$ in $k_{j}.$ Via Eq. (105), this gives rise to an order $3$ term in $\gamma_{\text{a}},$ and does not affect the terms explicitly set out in Eq. (118). Therefore this new choice of initial perturbation correlation produces the same leading and next-to-leading order entropy creation as our main model with its uncorrelated initial perturbation, and does not affect the results in this paper.