Dynamical origin of non-thermal states in galactic filaments

TL;DR

This study uses numerical simulations of self-gravitating particles to demonstrate that galactic filaments naturally evolve into non-thermal steady states with density profiles and temperature inversions similar to observations, indicating a dynamical origin.

Contribution

It introduces a comprehensive numerical analysis showing that non-thermal states in galactic filaments can arise dynamically from gravitational collapse and perturbations, linking simulations to observed features.

Findings

Non-thermal steady states match observed filament profiles.

Temperature inversion occurs in cold collapse scenarios.

Non-thermal states may be common in systems with long-range interactions.

Abstract

Observations strongly suggest that filaments in galactic molecular clouds are in a non-thermal state. As a simple model of a filament we study a two-dimensional system of self-gravitating point particles by means of numerical simulations of the dynamics, with various methods: direct $N$-body integration of the equations of motion, particle-in-cell simulations and a recently developed numerical scheme that includes multiparticle collisions in a particle-in-cell approach. Studying the collapse of Gaussian overdensities we find that after the damping of virial oscillations the system settles in a non-thermal steady state whose radial density profile is similar to the observed ones, thus suggesting a dynamical origin of the non-thermal states observed in real filaments. Moreover, for sufficiently cold collapses the density profiles are anticorrelated with the kinetic temperature, i.e., exhibit temperature inversion, again a feature that has been found in some observations of filaments. The same happens in the state reached after a strong perturbation of an initially isothermal cylinder. Finally, we discuss our results in the light of recent findings in other contexts (including non-astrophysical ones) and argue that the same kind of non-thermal states may be observed in any physical system with long-range interactions.

Tables

Table 1. Table 1: Values of the α 𝛼 \alpha exponent obtained fitting Eq. ( 4 ) to the final radial density profiles for cold collapses of initially Gaussian profiles with r 0 = 1 subscript 𝑟 0 1 r_{0}=1 .

$2K/\left|W\right|$ (t = 0)	$\alpha$ (MPC+PIC)	$\alpha$ (PIC)	$\alpha$ ($N$-body)
0.0	3.0	3.0	1.9
0.1	3.5	4.2	2.8
0.3	4.0	7.0	4.0
0.5	5.2	8.7	4.8
0.7	6.1	8.5	7.3

Equations

$$\{\begin{array}{l}{\partial }_t\rho +\nabla \cdot (\rho \mathbf{u})=0,\\ {\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\nabla \Phi -\frac{1}{\rho }\nabla P,\\ \Delta \Phi =4\pi G\rho ,\end{array}$$

$$P=\kappa {\rho }^{\gamma },$$

$$\nabla \cdot ({\rho }^{-1}\nabla {\rho }^{\gamma })+\frac{4\pi G}{\kappa }\rho =0,$$

$$\rho (r)=\frac{{\rho }_c{r}_c^{\alpha }}{(r_c^2+r^2{)}^{\alpha /2}},$$

$$M_{\ell }(r)=2\pi {\int }_0^r\rho (r^{\prime })r^{\prime }\mathrm{d}{r}^{\prime }=\frac{2\pi {\rho }_c{r}_c^{\alpha }}{(2-\alpha )}\left[{\left(r^2+r_c^2\right)}^{1-\alpha /2}-r_c^{2-\alpha }\right].$$

$$r_c=\sqrt{2/\pi G\beta {\rho }_c},$$

$$\mathcal{H}=\sum _{i=1}^N\frac{{\left|{\mathbf{p}}_i\right|}^2}{2m}+G{m}^2\sum _{i=1}^N\sum _{j\ne i}^N\log \frac{\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|}{r_s}\text{leavevmode},$$

$${\ddot{\mathbf{r}}}_i=-Gm\sum _{j\ne i}^N\frac{{\mathbf{r}}_i-{\mathbf{r}}_j}{{\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|}^2}\text{leavevmode}.$$

$$t_{\ast }\equiv \sqrt{\frac{2r_{\ast }^2}{G{M}_{\ell }}},$$

$$v_{\ast }\equiv r_{\ast }/t_{\ast }.$$

$${\ddot{\mathbf{r}}}_i=-Gm\sum _{j\ne i}^N\frac{{\mathbf{r}}_i-{\mathbf{r}}_j}{{ϵ}^2+{\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|}^2}\text{leavevmode};$$

$${\rho }_{i,j}={\tilde{\rho }}_{i,j}+\frac{1}{\Delta x\Delta y}\sum _{k=1}^{n_{i,j}}{m}_k,$$

$${\hat{\Phi }}_{i,j}={\hat{\rho }}_{i,j}{\hat{\mathcal{G}}}_{i,j},$$

$$\begin{array}{r}\frac{\partial \Phi }{\partial x}({\mathbf{r}}_k)=\frac{{\Phi }_{i+1,j}-{\Phi }_{i-1,j}}{2\Delta x}+\frac{{\Phi }_{i+1,j}+{\Phi }_{i-1,j}-2{\Phi }_{i,j}}{\Delta x^2}\delta x_k\\ +\frac{{\Phi }_{i+1,j+1}-{\Phi }_{i-1,j+1}+{\Phi }_{i-1,j-1}-{\Phi }_{i+1,j-1}}{4\Delta x\Delta y}\delta y_k.\end{array}$$

$${\mathbf{u}}_{ij}=\frac{1}{n_{ij}}\sum _{k=1}^{n_{ij}}{\mathbf{v}}_k.$$

$${\mathbf{v}}_k^{\prime }={\mathbf{u}}_{ij}+{\mathcal{R}}_{ij}\cdot \delta {\mathbf{v}}_k,$$

$${\mathbf{P}}_{ij}=\sum _{k=1}^{n_{ij}}m{\mathbf{v}}_k=\sum _{k=1}^{n_{ij}}m{\mathbf{v}}_k^{\prime }$$

$$K_{ij}=\frac{1}{2}\sum _{k=1}^{n_{ij}}m{\mathbf{v}}_k^2=\frac{1}{2}\sum _{k=1}^{n_{ij}}m{\mathbf{v}}_k^{\prime 2}$$

$$\sin {\vartheta }_{ij}=-\frac{2a_{ij}{b}_{ij}}{a_{ij}^2+b_{ij}^2};\cos {\vartheta }_{ij}=\frac{a_{ij}^2-b_{ij}^2}{a_{ij}^2+b_{ij}^2}.$$

$$a_{ij}=\sum _{k=1}^{n_{ij}}{\mathbf{r}}_k\wedge ({\mathbf{v}}_k-{\mathbf{u}}_{ij});b_{ij}=\sum _{k=1}^{n_{ij}}{\mathbf{r}}_k\cdot ({\mathbf{v}}_k-{\mathbf{u}}_{ij}).$$

$${\mathcal{P}}_{ij}=1-\exp \left[-(v_{ij}\Delta t{n}_{ij}{d}_{ij}/\Delta x\Delta y{)}^2\right].$$

$$\rho (r)={\rho }_c\exp \left(-\frac{r^2}{2r_0^2}\right),$$

$$K=\sum _{i=1}^N\frac{m{\left|{\mathbf{v}}_i\right|}^2}{2}$$

$$W=\sum _{i=1}^N{\mathbf{r}}_i\cdot {\mathbf{F}}_i$$

$${\mathbf{F}}_i=-G{m}^2\sum _{j\ne i}^N\frac{{\mathbf{r}}_i-{\mathbf{r}}_j}{{\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|}^2}$$

$$W=-\frac{G{m}^2}{2}N(N-1).$$

$$\xi =\frac{K_r}{K_{\phi }},$$

$$\beta (r)=1-\frac{{\sigma }_{\phi }^2(r)}{{\sigma }_r^2(r)}\text{leavevmode},$$

$$\rho (r)=\frac{r_c^4{\rho }_c}{(r_c^2+r^2{)}^2},$$

$$M_{\ell }(r)=\frac{\pi {\rho }_c{r}^2{r}_c^2}{r^2+r_c^2}=\frac{M_{\ell }{r}^2}{r^2+r_c^2},$$

Full text

††pagerange: LABEL:firstpage–LABEL:lastpage††pubyear: 0000

1 Introduction

Many astrophysical systems are found in states out of thermal equilibrium. This is not surprising especially when gravity is the dominant interaction, given the overwhelmingly large relaxation times111Strictly speaking, “true” thermodynamic equilibrium is impossible in isolated and finite self-gravitating systems in three dimensions, due to the non-confining nature of the potential; in this case the relaxation time denotes the time needed to effectively erase the memory of the initial condition. However, in the following we shall be mainly concerned with systems that are effectively two-dimensional, where the gravitational potential becomes confining and a thermal equilibrium state is possible. that grow with the number of constituents of a self-gravitating system (see e.g. Binney & Tremaine 2008). Recent observations, especially those obtained in submillimetre band by the Herschel Space Observatory, revealed that filaments in molecular clouds are not an exception to this rule, i.e., exhibit density profiles that are not consistent with thermal equilibrium (Arzoumanian et al. 2011; see also the discussion in Toci & Galli 2015a). Filaments are elongated structures present in galactic molecular clouds with densities ranging from $10^{4}$ to $10^{6}$ particles per ${\rm cm^{3}}$ that extend for tens of parsecs in length and that typically host star forming cores (see e.g. Klessen 2011; Myers 2015). Molecular filaments are expected to be short-lived ($\approx 10^{7}$ yr) and, at least in their initial stages, gravitationally supported structures (Gómez & Vázquez-Semadeni 2014, see also Burkert & Hartmann 2004) embedded in galactic magnetic fields (Gehman et al. 1996; Crutcher 2012; Li et al. 2013; Toci & Galli 2015b, and references therein). Another intriguing feature of the non-thermal states of filaments is that the radial density $\rho(r)$ and kinetic temperature $T(r)$ profiles are often anticorrelated, i.e., the kinetic temperature of the gas is higher at larger radii where the density is smaller; this emerges indirectly, from the fact that density profiles are well fitted by polytropic solutions yielding $T$ versus $\rho$ power laws with negative exponents (Toci & Galli 2015a), and also from some direct observations of velocity dispersion (Goodman et al. 1998; Pineda et al. 2010). This feature, dubbed “temperature inversion”, is often found also in astrophysical plasmas encompassing a broad range of sizes and temperatures, from, for example, cool cores of galaxy clusters (Wise et al. 2004) to the hot ($5\times 10^{6}$ K) plasma torus surrounding Jupiter’s moon Io (Meyer-Vernet et al. 1993, 1995; Saur et al. 2004). Stellar coronae also show temperature profiles increasing with radius (from $10^{3}$ K at the photosphere up to $10^{6}$ K for our Sun) where instead the density rapidly falls off with radius (see e.g.  Aschwanden 2005; Golub & Pasachoff 2009, and references therein). Observations have been reported of elliptical galaxies where the velocity dispersion increases with the distance from the nucleus (Loubser et al., 2008; Veale et al., 2017). Transient anticorrelated density and temperature profiles have been observed also in one-dimensional self-gravitating systems (Worrakitpoonpon, 2011). Typical non-thermal features that might bear a relation with temperature inversion are core–halo structures, that are ubiquitous in long-range-interacting systems (Levin et al., 2014), having been observed in charged particle beams (Nghiem et al., 2014; Davidson & Qin, 2001; Holloway et al., 2008), in one-dimensional toy models of gravity (Cuperman et al., 1969; Teles et al., 2011; Schulz et al., 2013), and in two- and three-dimensional models of self-gravitating systems (Ruffini et al., 2015; Danilov, 1988; Drukier et al., 1999). Recently it has been suggested that non-thermal states with temperature inversion may be typical of any long-range-interacting many-body system (Casetti & Gupta 2014; Teles et al. 2015, 2016), being long-time quasi-stationary states of the dynamics arising after a violent relaxation process, and that could be observed also in condensed matter systems (Gupta & Casetti 2016).

The observation of non-thermal states in filaments raised questions about the physical processes that may support such nonequilibrium states, and many ideas have been put forward, typically involving the rôle of local turbulence and of magnetic fields (see e.g. Toci & Galli 2015a, b and references therein). In this paper we argue that ingredients other than gravitational interactions may not be necessary at all, as such states may simply be the end states of the dynamics of the self-gravitating particles composing the filament, after a cold collapse as well as after a strong perturbation of an initially thermal state. To this end, we present a dynamical study of a very simple model of a filament, where both dissipative effects and magnetic fields are completely neglected, assuming our system is infinitely extended with perfect cylindrical symmetry, so that its dynamics reduces to that of a two-dimensional system of $N$ self-gravitating particles. We study the dynamics of the collapse of Gaussian overdensities, showing that after the damping of the virial oscillations the system settles in a quasi-stationary state whose density profile is well described by the empiric law used to model the radial density profile of observed filaments in Arzoumanian et al. (2011) and discussed by Toci & Galli (2015a) showing its relation to polytropic solutions of the fluid equations for the gas in the filament. For sufficiently cold collapses, the end states we find do exhibit temperature inversion. The same kind of non-thermal states with temperature inversion is observed after a strong perturbation of an initially thermal state.

The paper is organized as follows. In Section 2 we first recall previously studied fluid models, then introduce our particle-based model and our numerical approach. In Section 3 we study the collapse of cylindrical Gaussian overdensities, with different initial values of the virial ratio and of the spread of the Gaussian, and in Section 4 we study the dynamics after a strong perturbation of an initial state in thermal equilibrium. In Section 5 we discuss our results and their connection with those found for other long-range-interacting many-body systems, also in non-astrophysical contexts, and finally in Section 6 we draw our conclusions.

2 Model and numerical techniques

2.1 Fluid models

Neglecting the presence of galactic magnetic fields, in the fluid picture the dynamical evolution of self-gravitating filaments is given in terms of mass density $\rho$, pressure $P$, gravitational potential $\Phi$ and velocity field $\mathbf{u}$ by

[TABLE]

where $G$ is the gravitational constant. Within the assumption of cylindrical symmetry, it is natural to express the operators in Eqs. (1) in cylindrical coordinates $(r,z,\varphi)$ and to set $\partial_{z}\ =\partial_{\varphi}=0$. Stationary solutions of Eqs.  (1) (i.e., with $\partial_{t}\rho=\partial_{t}\mathbf{u}=0$ and222As noted by Toci & Galli (2015a), stationary solutions with a nonzero velocity field are possible and might be relevant when modeling accretion phenomena or the so-called “varicose” instabilities. $(\mathbf{u}\cdot\nabla)\mathbf{u}=0$, and thus $\nabla\Phi=-\rho^{-1}\nabla P$) are typically obtained by imposing an equation of state relating pressure to density and solving for $\rho$ the differential equation obtained by combining the second and the third equations after having expressed $P$ as function of $\rho$. One of the most common forms used to model filamentary clouds is the polytropic equation of state (Rouleau & Bastien 1990; Arcoragi et al. 1991; Shadmehri 2005; Breysse et al. 2014; Lou 2015)

[TABLE]

where the constant $\kappa$ is connected to the system’s entropy, and the so-called polytropic exponent $\gamma$ is related to the polytropic index $n$ by $\gamma=1+1/n$. In this case, the density profile supporting the polytropic equation of state is obtained by solving

[TABLE]

that after exploiting the cylindrical symmetry and introducing dimensionless coordinates reduces to the standard cylindrical Lane-Emden equation (Toci & Galli 2015a). In polytropic stationary models, the density profile is related to the temperature profile as $\rho\propto T^{n}$.

A good approximation of the solution of Eq. (3) for infinite cylinders is given by the softened power-law radial density profile (Toci & Galli 2015a; Federrath 2016; Chira et al. 2016; Xiong et al. 2017)

[TABLE]

associated to the radial mass per unit length profile

[TABLE]

In the equations above $\rho_{c}$ and $r_{c}$ are the (finite) core density and core radius, respectively. Note that the total mass per unit length at height $z$, $M_{\ell}$, obtained extending to $+\infty$ the integral in Eq. (5), is finite for $\alpha>2$. Note also that if $\alpha=4$ in Eq. (4) we would obtain the density profile for the isothermal cylinder (Stodólkiewicz 1963; Ostriker 1964a, b). In this case the core radius $r_{c}$ becomes

[TABLE]

where $\beta=m_{0}/k_{B}T$ fixes the width of the equilibrium Maxwellian velocity distribution at temperature $T$ with mean molecular mass $m_{0}$, and $k_{B}$ is the Boltzmann constant. The parameter $\beta$ can be also expressed as function of the total mass333Note that $M_{\ell}$ is a finite quantity at variance with the diverging mass of an isothermal sphere. The total mass of the isothermal filament, however, is indeed infinite, since the integral of $M_{\ell}$ over $z$ diverges. per unit length $M_{\ell}$ as $\beta=2/GM_{\ell}$.

Although early observations with limited spatial extent yielded radial density profiles consistent with an isothermal cylinder, the density profiles from recent observations with high dynamic range like those by the Herschel Space Observatory, especially of filaments located in nearby molecular clouds, are well described at large radii by considerably softer power laws than the isothermal case, i.e., $\alpha\approx 2$ or even smaller in Eq. (4), indicating that these filaments are in a non-thermal state (Toci & Galli 2015a).

2.2 Particle-based model

Our aim is to investigate how non-thermal stationary states may emerge from the dynamics. Hence, instead of studying the fluid equations (1) and their stationary states, we use a particle-based model and study its dynamics by integrating its equations of motion. We still assume our system is infinitely extended and uniform along its longitudinal ($z$) axis, so that its dynamics is that of a system of infinite straight massive wires, whose mutual interaction is described by a two-dimensional (logarithmic) gravitational potential444We note that a similar model is used also in the context of unbunched charged particle beams in linear accelerators (see e.g. Riabko et al. 1995; Wangler et al. 1998; Levin et al. 2008, and references therein). (see e.g. Katz & Lynden-Bell 1978; Aly 1994; Teles et al. 2010). Hence, the dynamics in a plane transverse to the $z$ direction is governed by the Hamiltonian of a two-dimensional system of $N$ self-gravitating particles,

[TABLE]

where $\mathbf{r}_{i}$ and $\mathbf{p}_{i}$ are, respectively, the position and momentum of the $i$-th particle in the plane, $r_{s}$ is an arbitrary scale length needed to make the argument of the logarithm dimensionless, and the mass $m$ of each particle is related to the total mass per unit length $M_{\ell}$ of the filament by $m=M_{\ell}/N$. The equations of motion to be integrated are thus

[TABLE]

Various aspects of the dynamics of such a system have been studied by Chavanis (2006), Teles et al. (2010), Gabrielli et al. (2010), Marcos (2013), Levin et al. (2014), Marcos et al. (2017), Silvestre & Rocha Filho (2016). For this model it is convenient to introduce dimensionless quantities in terms of scale radius $r_{*}$, dynamical time $t_{*}$ and scale velocity $v_{*}$. We define the dynamical timescale of the system as

[TABLE]

where $r_{*}=r_{50}(0)$ (i.e. the Lagrangian radius containing half of $M_{\ell}$ at $t=0$), and the velocity scale as

[TABLE]

2.3 Numerical approach

2.3.1 Direct $N$-body simulations

First of all, we performed direct $N$-body simulations where the acceleration on each particle was evaluated by direct sum. In this case, the divergence of the inter-particle potential and force for vanishing interparticle separation was regularized by introducing a softening length555The results we are going to present in the following do not depend on the choice of the softening length, at least as long as $10^{-5}\leq\epsilon/r_{*}\leq 10^{-2}$. For a theoretical study of the effect of a softening of power-law interactions on the existence and lifetime of nonequilibrium quasi-stationary states see Gabrielli et al. (2010) and Marcos et al. (2017). $\epsilon=10^{-3}r_{*}$, so that Eq. (8) actually reads

[TABLE]

the above equations were then integrated with a third-order symplectic algorithm (see e.g. Kinoshita et al. 1991; Casetti 1995) with fixed timestep $\Delta t=t_{*}/100$ ensuring energy conservation up to 14 digits in double precision.

In order to effectively explore different initial conditions, we limited ourselves to rather small $N$’s, that is, $N=3\times 10^{4}$. Our aim is to understand whether the purely gravitational dynamics of the system defined by the Hamiltonian (7) is able to produce non-thermal states similar to those observed in real filaments. Since the two-body relaxation time, after which “collisional” effects would eventually destroy any non-thermal steady state, is expected to grow666This happens in any long-range-interacting system; see e.g. Campa et al. (2014) for a general discussion and Marcos (2013), Teles et al. (2010), Gabrielli et al. (2010) and Marcos et al. (2017) for a discussion more focused on gravity and general power-law interactions. with $N$, $N$ should be large enough to allow the formation and survival of such states. It has been shown by Marcos (2013) that the two-body relaxation time in a two-dimensional system of self-gravitating particles is of the order of $Nt_{*}$, that is, a factor of $\ln N$ longer than in a three-dimensional system with the same $N$, while the time needed to form a non-thermal steady state should not depend on $N$ (and we shall see in Sec. 3 that such a timescale is not larger than $100t_{*}$). Hence, we expect that a number of particles of the order of $10^{4}$ is already sufficient to describe the physics we are interested in. This notwithstanding, it would be interesting to check whether our results do depend on $N$. In principle one should run $N$-body simulations with much larger $N$. Being such a task computationally demanding, to increase the total number of particles of three orders of magnitudes, up to $N=1.5\times 10^{6}$, we resorted to a standard two-dimensional particle-in-cell (PIC) method (to be described in Sec. 2.3.2) and to a novel method capable of effectively taking close encounters (collisions) into account in a PIC framework (Sec. 2.3.3). These methods do not give a description of the dynamics of the model as accurate as the $N$-body method does, in that they do not properly model the physics at small length scales, but should still yield reasonable results, the more so because it has been shown by Marcos (2013) that in a two-dimensional self-gravitating system the slow relaxation processes that would eventually destroy the nonequilibrium stationary state are not dominated by the small-scale physics. Moreover, comparing the results obtained by evolving the same initial conditions with different numerical methods may help us to gain some deeper insights in the results themselves (see Sec. 5).

2.3.2 Particle-in-cell simulations

In a PIC approach, the system is represented by $N$ macro particles with mass $m_{k}$ moving in the plane. At every time step a fixed Cartesian grid of $N_{c}=N_{x}\times N_{y}$ cells is superimposed to the system, so that the mesh based mass density is obtained as

[TABLE]

where $\Delta x$ and $\Delta y$ are the cell’s side lengths, $n_{i,j}$ is the number of particles enclosed by the cell, and $\tilde{\rho}$ is a static background density. In our implementation the particles’ contribution to the density in a cell is evaluated with the simple nearest-grid-point method (NGP, Hockney & Eastwood 1981). In order to compute the resulting gravitational potential, the code solves the Fourier-space based Poisson equation

[TABLE]

where the hats over the quantities denote the Fourier transforms, and $\mathcal{G}$ is the Green’s function of the Laplace operator. Once the mesh based potential $\Phi_{i,j}$ is recovered by antitransforming $\hat{\Phi}_{i,j}$, the gravitational acceleration $\mathbf{a}(\mathbf{r}_{k})=-\nabla\Phi(\mathbf{r}_{k})$ acting on a particle sitting at position $\mathbf{r}_{k}$ inside the cell of indices $i,j$, is obtained via the gradient interpolation (see Fellhauer et al. 2000)

[TABLE]

In the equation above (and in the analogous one for the $y$ component of $\nabla\Phi$), $\delta x_{k}$ and $\delta y_{k}$ are the components of $\mathbf{r}_{k}-\mathbf{c}_{i,j}$, where $\mathbf{c}_{i,j}$ is the position of the cell’s center. As in the $N$-body case, the equations of motion of the particles are integrated with a third order symplectic integrator scheme with fixed timestep $\Delta t=t_{*}/100$ ensuring energy conservation up to 14 digits in double precision. For the simulations reported here we have used $N=1.5\times 10^{6}$ particles on a Cartesian grid with $N_{c}=128\times 128$ equal square cells.

2.3.3 Inclusion of collisions

One of the limits of a PIC approach is that it underestimates the contribution of close encounters (collisions) between particles, that are included in the computationally heavier $N$-body simulations, at least up to a distance scale of order $\epsilon$. To effectively overcome this limitation, following a suggestion of Bufferand et al. (2010); Bufferand et al. (2013) and Di Cintio et al. (2015, 2017) in the context of numerical modelling of plasma transport, we treat collisional processes (i.e.  particle scatterings) with the stochastic multi-particle collision scheme (hereafter MPC), originally introduced by Malevanets & Kapral (1999) in the field of mesoscopic fluid dynamics (see also Malevanets & Kapral 2004; Kapral 2008; Gompper et al. 2009). The MPC methods simulate inter-particle collisions at the cell level by performing a stochastic rotation of their velocity vectors $\mathbf{v}_{k}$, constrained by the conservation of the total momentum and kinetic energy of the cell. In a two dimensional implementation, at the beginning of each timestep $\delta t$, in the cell of indexes $i,j$ first of all the center of mass velocity is evaluated as

[TABLE]

Then an angle $\vartheta_{ij}$ is sampled from a uniform distribution in $(0,2\pi)$. The collision itself is simulated by rotating with probability one-half the relative velocities $\delta\mathbf{v}_{k}=\mathbf{v}_{k}-\mathbf{u}_{ij}$, as

[TABLE]

where $\mathcal{R}_{ij}$ is the 2D rotation matrix of an angle $\vartheta_{ij}$. Such a rotation guarantees the conservation of the total momentum

[TABLE]

and kinetic energy

[TABLE]

in the cell. By doing so, however, the total $z$-component $L_{z}$ of the angular momentum of the cell is not conserved. To overcome this inconvenience, following Ryder (2005) in our implementation of the MPC rule, we introduce the following constraints on $\vartheta_{ij}$

[TABLE]

The latter hold if the coefficients $a_{ij}$ and $b_{ij}$ are given as functions of particles’ positions and velocities in the cell by

[TABLE]

Since the effective collisionality may in principle depend on the local state of the system (i.e., particles’ velocities and density), the collision move is accepted only if a random probability $\mathcal{P}^{*}_{ij}$ sampled from a uniform distribution in $(0,1)$ is smaller than the cell dependent collision probability

[TABLE]

In the expression above, $v_{ij}$ and $d_{ij}$ are the average velocity and average interparticle distance in the cell, respectively. Note that the collision probability has no adjustable parameters. It must be noted that the “collisions” taken into account by this PIC+MPC scheme are not the “true” ones dictated by the gravitational interactions, but only a stochastic effective model of the latter.

3 Cold collapses

The first scenario we investigate for the dynamical generation of non-thermal states in our simple model of a filament is that of cold collapse. As initial condition for the positions we assume a cylindrically symmetric Gaussian overdensity with initial radial density profile given by

[TABLE]

where the central density $\rho_{c}$ is taken to be ten times larger than the homogeneous background $\rho$. The initial particle velocities are extracted from a position independent Maxwellian such that the initial virial ratio of the system is $(2K/|W|)_{t=0}<1$, where

[TABLE]

is the total kinetic energy of the system, and

[TABLE]

is the virial function, where

[TABLE]

is the force exerted on the $i$-th particle. We note that in the two-dimensional case $W$ does not depend on time,

[TABLE]

In the following we report on the results obtained with the Gaussian initial profiles (22), but those with the truncated initial distribution look almost the same.

As expected, at small times the system undergoes a phase of violent relaxation with strong oscillations of the virial ratio, stronger when the initial value of $2K/|W|$ is smaller (Fig. 1). Virial oscillations damp out on a timescale of some tens of dynamical times $t_{*}$, and at times $t>100t_{*}$ the system appears to be in an essentially steady state. In Fig. 1 only the results obtained with the MPC+PIC simulation protocol are reported but those obtained with simple PIC without collisions and $N$-body simulations are very similar.

Monitoring the kinetic temperature as a function of the position (i.e., the locally averaged kinetic energy) during the time evolution for a completely cold (i.e., with vanishing initial virial ratio) collapse, we observe that it wildly changes exhibiting highly irregular patterns during the violent relaxation phase, but then settles in an azimuthally symmetric pattern where temperature inversion is clearly evident, as shown in Fig. 2.

3.1 Density and temperature profiles

The density and kinetic temperature (azimuthally averaged) radial profiles in the steady states reached after a cold collapse of a Gaussian initial profile with $r_{0}=1$ are shown in Fig. 3, for all the three types of simulations we performed, PIC+MPC, PIC without collisions, and direct summation $N$-body. The qualitative properties of the density profiles are very similar: all exhibit a flat, essentially isothermal core and decay as a power law at larger radii. The power law is steeper for initially warmer and softer for initially colder systems, that also exhibit temperature inversion, while initially warmer systems do not show this feature and rather their kinetic temperature decreases with radius. Temperature inversion does depend on the simulation method: it is stronger in the $N$-body and especially in the PIC+MPC cases, while it is much less evident in the PIC without collisions simulations.

The properties of the final steady state also depend on the width of the initial Gaussian profile, as shown in Fig. 4 for MPC+PIC and simple PIC simulations: as a general rule, regardless of the presence of collisions, initially more concentrated systems produce softer decays at large radii and stronger temperature inversions in the final state.

Remarkably, as shown in Fig. 5, the density profiles in the final steady state are quite well described by the empiric relation777This may suggest that these states are indeed essentially polytropic, but an analysis of the velocity anisotropies (see Sec. 3.2) shows that this may be an oversimplified description. (4). The exponent $\alpha$ depends on the initial virial ratio and encompasses a broad range of values that includes the thermal value $\alpha=4$, being smaller for initially colder systems and varying from $\alpha=1.9$ for a completely cold collapse in an $N$-body simulations to $\alpha=8.7$ for a collapse with initial virial ratio equal to 0.5 in a PIC simulation.

The values of the fitted $\alpha$ exponent are also reported in Table 1 in the various cases of cold collapses of Gaussian initial states with $r_{0}=1$.

3.2 Velocity distribution anisotropy

Further insight into the properties of the violent relaxation process at the beginning of the collapse as well as into the features of the non-thermal state reached after the damping of the virial oscillations can be obtained by studying the anisotropies of the velocity distribution. Were the system in a thermal equilibrium state, the velocity distribution would be isotropic, but also non-thermal states can be isotropic, e.g., when the single-particle distribution function only depends on the total energy. To this end, in our direct $N$-body simulations we adopted standard diagnostic tools to reveal possible anisotropies. First, we measured the time evolution of the global anisotropy parameter (see e.g. Trenti & Bertin 2006 and references therein)

[TABLE]

where $K_{r}$ and $K_{\varphi}$ are the radial and tangential kinetic energies, respectively. In a system where the orbits are mostly radial (resp. tangential), $\xi$ is large (resp. small); an isotropic velocity distribution yields $\xi=1$. After the system had reached its steady state, we measured the azimuthally averaged radial profile of the anisotropy parameter

[TABLE]

where $\sigma_{r}$ and $\sigma_{\varphi}$ are the radial and tangential velocity dispersions, respectively. The parameter888Not to be confused with the $\beta$ entering Eq. 6, that is related to the temperature of the isothermal cylinder. $\beta$ measures the pressure anisotropy as a function of the radius $r$, as well as the kinematical structure of the particles’ orbits (Binney & Tremaine, 2008): in an isotropic state $\beta=0$, while positive (resp. negative) values of $\beta$ denote a radial (resp. tangential) anisotropy. Purely radial orbits yield $\beta=1$, while perfectly circular orbits would give $\beta=-\infty$. The results are shown in Fig. 6.

The time evolution of $\xi$ shows that the global velocity anisotropy is rather large at the beginning of the collapse (colder collapses are more anisotropic) and decreases with the damping of the virial oscillations; in warmer collapes the global amount of anisotropy is negligible at large times, while it remains quite relevant for initially colder states. The $\beta(r)$ profiles in the non-thermal steady state are all monotonously increasing with radius (with the only exception of the warmer collapse) and show that in the outer regions anisotropy is always present and low-angular momentum orbits dominate (again with the exception of the warmer collapse where the anisotropy becomes tangential at large radii). The fact that initially warmer states are also more isotropic, both in their time evolution and in their end states, is expected due to the larger amount of random motions with respect to the initially colder cases. The rather large amount of radial anisotropy found in the initially colder cases, especially in the outer regions, may be due to the fact that in cold collapses the velocities are initially extremely anisotropic, because in the initial stages of virial oscillations the motion is predominantly radial when starting from cylindrically symmetric initial conditions (this is witnessed by the fact that $\xi$ is very large at small times)999We note that in two dimensions we may expect larger anisotropies than in three dimensions, due to the smaller number of available tangential channels for the velocity to relax to an isotropic distribution.. Anisotropy is thus generated at the beginning of the virial oscillations and decreases during the relaxation to the non-thermal end state, but only in warmer collapses this process efficiently reduces the initial anisotropy to a small amount101010In three-dimensional collapses it has been shown that starting from highly symmetric initial conditions results in more anisotropic end states than in the case of less symmetric initial states (Trenti & Bertin, 2006; van Albada, 1982).. Since the initially colder collapses are also those where temperature inversion is more pronounced, we observe a correlation between temperature inversion and velocity anisotropy in our simulations; however, our results alone are not sufficient to suggest that a correlation between temperature inversion and anisotropy should hold in general in a low-dimensional self-gravitating system.

As mentioned in Sec. 2.1, density profiles given by Eq. (4) are very close to those obtained in a fluid picture assuming a polytropic equation of state. However, the latter states are isotropic, hence our results on velocity anisotropy suggest that an interpretation of the non-thermal states of filaments in terms of purely polytropic states may be an oversimplified picture, although it correctly describes the radial density profile.

4 Perturbation of a thermal state

The second dynamical scenario we consider is the evolution of a strong impulsive radial perturbation of an initially isothermal profile, i.e., given by the Ostriker-Stodólkiewicz self-gravitating cylinder–Eq. (4) with $\alpha=4$. Physically, this scenario may model, albeit in a very simplified way, the evolution of an already thermally relaxed filament that is suddenly perturbed by a shock wave. Moreover, such scenario is interesting because it bridges our results to other recent results obtained for models of long-range-interacting systems in non-astrophysical contexts (Casetti & Gupta 2014; Teles et al. 2015, 2016; Gupta & Casetti 2016), where perturbations of initially thermal states were considered.

The density profile of the isothermal cylinder is given by setting $\alpha=4$ in Eq. (4), thus obtaining

[TABLE]

so that the associated radial mass per unit length profile becomes

[TABLE]

where the core radius $r_{c}$ is linked to $T$ through Eq. (6) and $r_{*}=r_{c}$, see also Ostriker (1964a). Initial particle positions are obtained by sampling Eq. (30) for the radial coordinate $r$ and then assigning an angular coordinate $\varphi$ form a uniform distribution in the interval $(0,2\pi)$. After this, the velocity components $v_{r}$ and $v_{\varphi}$ are extracted consistently from the position independent isotropic Maxwellian at temperature $T$

[TABLE]

and finally converted in Cartesian coordinates along with their positions components. In order to simulate a radial compression, the radial component of particles’ velocities $v_{r}$ is instantaneously increased of a quantity $\delta v_{r}\approx 2\sqrt{k_{B}T}$. Alternatively, in some simulations particles’ positions have been displaced of a quantity $\delta_{r}\approx r_{*}/80$ along their radial coordinate; the results are very similar to those obtained with the velocity displacement, and in the following we report only on the latter.

The time evolution of the virial ratio for different strengths of the perturbation is shown in Fig. 7. As in the cold collapses, the system undergoes virial oscillations, that damp out on a timescale of tens of dynamical times, and then settles in a steady state. Kinetic temperature and density profiles of such a state are shown111111Preliminary results limited to the $N$-body case were already presented by Teles et al. (2015). in Fig. 8. It is apparent that the kinetic temperature and density radial profiles are anticorrelated: temperature inversion is even stronger than in the completely cold collapse case. As far as the relation between temperature and density is concerned, the steady states reached after the perturbation of a thermal cylinder appear as polytropic, at least when density is not too small: indeed, the $\rho$ versus $T$ relation is reasonably fitted by a power law (data not shown), and in general the quality of the fit improves for systems undergoing stronger perturbations. However, also in this case there are velocity anisotropies (data not shown), so that these states can not be considered as “true” polytropic states.

5 Discussion

As shown in Sec. 3.1, the steady states reached by our simple model of a filament after a cold collapse are definitely non-thermal, unless the initial condition is fine-tuned, and qualitatively very similar to the observed states, as their density profiles are in good agreement with Eq. (4) with values of the exponent $\alpha$ that depend on the initial condition, being in general smaller for colder collapses. Another indication of the fact that the steady states are non-thermal is the observed velocity anisotropy (Sec. 3.2). Moreover, colder collapses exhibit appreciable temperature inversion. The fact that such a simple model of a filament, neglecting all the dissipative effects, rotation, radiation, and magnetic fields, and retaining only self-gravity, exhibits steady states that are so similar to the observed ones suggests that the dynamics alone is sufficient to produce and support the non-thermal states observed in real filaments, and that other mechanisms are not strictly necessary. Moreover, the fact that all the three numerical methods ($N$-body, PIC and PIC+MPC) yield qualitatively similar results is an indication of the robustness of the results themselves. Quantitatively, the values of the exponent $\alpha$ governing the decay of the radial density profile at large radii obtained for colder collapses with $N$-body simulations seem compatible with those observed in real filaments (Toci & Galli 2015a), while all the other cases yield values that are larger than the observed ones (see Table 1). Still, the values of $\alpha$ depend not only on the initial conditions, but also on the choice of the simulation protocol. First of all, we note that the values of $\alpha$ closer to the observed ones are obtained by means of $N$-body simulations. This is an important result because the $N$-body method is the most correct one among the three considered, since it captures the physics at all length scales (larger than the softening length). Then, as witnessed by the values reported in Table 1, we observe that the general trend, i.e., colder collapses corresponding to final profiles decaying with softer power laws, is the same for all the three simulation techniques, but $N$-body and PIC simulations have the same spread in the values of $\alpha$ while the $\alpha$’s obtained with MPC+PIC simulations are closer to each other; moreover, as far as absolute values are concerned, $N$-body simulations yield the smallest $\alpha$’s, at least for colder collapses, while the values found in PIC simulations are the largest. This may suggest that the small-scale physics, including close encounters between particles, that is somewhat suppressed in PIC simulations and (only approximately, via a stochastic model) restored in MPC+PIC ones, is important for the processes that produce “hot” particles able to create a more diffuse halo at large radii. Consistently, temperature inversion is stronger in states obtained with MPC+PIC and $N$-body simulations, while is suppressed in PIC simulations. On the other hand, the fact that the values of $\alpha$ obtained with MPC+PIC simulations span a smaller interval and are closer to the thermal value may be a consequence of that including random collisions allows the system to get closer to thermal than in the other cases.

As recalled in the Introduction, self-gravitating systems do not reach thermal equilibrium (that in the case of two-dimensional systems exists also for a finite system, at variance with the three-dimensional case) when starting from nonequilibrium initial conditions because the time need to reach a thermal state is exceedingly long and grows with the number of particles $N$. The steady states we observe after the violent relaxation phase are rather non-thermal, quasi-stationary states that would become exactly stationary in the limit $N\to\infty$. Theoretically, such states are identified with stable stationary solutions of the Vlasov (or collisionless Boltzmann) equation, that describes the evolution of the single-particle distribution function for such systems for time scales smaller than the relaxation time (that becomes infinite in the $N\to\infty$ limit). In the last decades, it has been recognized that the presence of such quasi-stationary, non thermal states as well as of long relaxation times growing with $N$ and preventing an effective relaxation to thermodynamic equilibrium are common to all many-body systems with long-range interactions and not restricted at all to self-gravitating systems (Campa et al. 2014; see also Di Cintio et al. 2013; Di Cintio & Ciotti 2011). It has been recently suggested that non-thermal states with temperature inversion are a typical outcome when a generic many-body long-range-interacting system relaxes to a quasi-stationary states after undergoing collective oscillations (Casetti & Gupta 2014; Teles et al. 2015, 2016; Gupta & Casetti 2016). A possibly “universal” mechanism has been proposed to explain the ubiquity of temperature inversion by Teles et al. (2015). The mechanism can be summarized as follows: during the initial violent relaxation phase the interaction of the particles with the collective oscillations, that will eventually lead to the Landau-like damping of the oscillations themselves, may produce suprathermal tails in the velocity distribution function even when starting from a thermal state. In an inhomogeneous system, this may trigger a “velocity filtration” process (Scudder 1992a, b) that effectively broadens the velocity distribution function in the regions where the system is less dense, because only sufficiently fast particles can escape the potential well produced by the central concentration.

The systems studied by Teles et al. (2015) and Gupta & Casetti (2016) were initially in a thermal state, and then brought out of equilibrium by a sudden perturbation or a quench of an external parameter. The results we have shown in Sec. 4 confirm the general picture also for two-dimensional gravitational systems. The fact that also cold collapses, where the starting configuration is not a thermal one, may exhibit temperature inversion is a significant addition to this picture: it indeed suggests that collective oscillations and wave-particle interactions are the crucial ingredient to produce a final state with temperature inversion, although not a sufficient one. Our results with MPC+PIC simulations then suggest that collisions may even enhance the temperature inversion.

6 Concluding remarks

The main result of the work we have presented here is that a very simple model of a filament, that completely neglects dissipative effects, radiation, magnetic fields, and rotation, retaining only gravity and effectively mapping the dynamics of a filament onto that of a two-dimensional system of self-gravitating particles, is able to produce non-thermal steady states similar to those observed in real filaments, as end states emerging after the damping of virial oscillations in cold collapses. The radial density profiles of such states are well described by Eq. (4) and for sufficiently cold collapses $N$-body simulations yield values of the exponent $\alpha$ governing the decay of the density profile at large radii in the range of the observed ones.

Clearly, our results do not imply that mechanisms other than dissipationless dynamics alone are not at work and that they are not important in the real case; all the effects and processes we neglected in our very simple models may be important and a strategy to understand their rôle and relevance may be the progressive addition of other ingredients to our simple model, to see whether they may bring the non-thermal final states closer or not to the observed ones. The first step is naturally the inclusion of magnetic fields, that couple to the small fraction of ionized molecules in the filament and whose effect is thus transferred to the whole system by the interactions between the ionized particles and all the others. Work is in progress in this direction. In the fluid picture, polytropic solutions for infinite cylinders with magnetic fields were analyzed by Toci & Galli (2015b), while the effect of rotation and of an imposed temperature gradient were considered by Recchi et al. (2013, 2014).

Finally, the non-thermal states observed as end products of the dynamics of our simple model exhibit temperature inversion, that seems a common feature of a wide class of non-thermal states of any long-range-interacting system, hinting at a possible and intriguing relation between non-thermal states of galactic filaments and properties of nonequilibrium states that could be observed in other systems, even in condensed matter ones.

Acknowledgments

We thank Guido Ciraolo, Daniele Galli, and Tarcísio N. Teles for insightful discussions and important comments at an early stage of this work. We also thank an anonymous referee for very useful comments that helped us to improve the presentation of our work.