A pure non-neutral plasma under an external harmonic field: equilibrium thermodynamics and chaos

Experimentally, a non-neutral plasma can be effectively confined within external magnetic and electric fields of a Penning trap, whose Hamiltonian can be written as follows [9],

$$H(\mathbf{r},\mathbf{p})=\sum_{i}\frac{1}{2m}\left[\mathbf{p}_{i}-q\mathbf{A}(\mathbf{r}_{i})\right]^{2}+q\phi_{T}(\mathbf{r}_{i})+\sum_{i<j}\frac{q^{2}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}.$$

(1)

However, this type of confinement does not avoid the evaporation of constituents during the system dynamical evolution [38] neither the Hamiltonian (1) obeys the standard form

$$H(q,p)=\sum_{i,j}\frac{1}{2}a^{ij}\left(q\right)p_{i}p_{j}+V(q),$$

(2)

required for the application of Riemannian approach of Hamiltonian chaos [41]. Although this geometric framework will not be considered in this work, we think convenient to address a more simple Hamiltonian model:

$$H(\mathbf{r},\mathbf{p})=\sum_{i}\frac{1}{2m}\mathbf{p}^{2}_{i}+\frac{1}{2}m\omega^{2}\mathbf{r}^{2}_{i}+\sum_{i<j}\frac{q^{2}}{\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|},$$

(3)

which avoids all previous difficulties, and it still enables us to study the chaoticity of a non-neutral plasma due to the long-range character of Coulombian forces. Hamiltonian (3) represents a gas with $N$ identical non-relativistic point particles of mass $m$ and charge $q$, interacting through Coulomb’s force under an external harmonic field with frequency $\omega$. This particular situation guarantees the fully confinement of the non-neutral plasma and avoids the incidence of evaporation. Consequently, the thermo-statistical description can be performed assuming a rigorous thermodynamic equilibrium.

The equilibrium one-body distribution function (DF) for the present model system is given by Maxwell-Boltzmann profile [42]:

$$f(\mathbf{r},\mathbf{p})=A\exp\left[-\beta\varepsilon(\mathbf{r},\mathbf{p})\right].$$

(4)

Here, $\varepsilon(\mathbf{r},\mathbf{p})$ denotes the mechanical energy of individual particles:

$$\varepsilon(\mathbf{r},\mathbf{p})=\frac{1}{2m}\mathbf{p}^{2}+q\varphi(\mathbf{r})+\frac{1}{2}m\omega^{2}\mathbf{r}^{2},$$

(5)

$\beta=1/kT$ is the inverse temperature parameter, $k$ the Boltzmann constant, and $\varphi(\mathbf{r})$, the electrostatic potential. Particles distribution $n(\mathbf{r})$ can be calculated as follows:

	$\displaystyle n(\mathbf{r})$	$\displaystyle=\int f(\mathbf{r},\mathbf{p})\frac{d^{3}\mathbf{p}}{(2\pi\hbar)^{3}}$
		$\displaystyle=\left(\frac{2m\pi}{\hbar^{2}\beta}\right)^{\frac{3}{2}}A\exp\left\{-\beta\left[\frac{1}{2}m\omega^{2}\mathbf{r}^{2}+q\varphi(\mathbf{r})\right]\right\},$

which obeys the normalization condition:

$$N[f]=\int f(\mathbf{r},\mathbf{p})\frac{d^{3}\mathbf{r}d^{3}\mathbf{p}}{(2\pi\hbar)^{3}}=\int n(\mathbf{r})d^{3}\mathbf{r}=N.$$

(7)

The total energy $U$ is given by:

$$U=\frac{3N}{2\beta}+\int\frac{1}{2}m\omega^{2}\mathbf{r}^{2}n(\mathbf{r})d^{3}\mathbf{r}+\frac{1}{2}\int q\varphi(\mathbf{r})n(\mathbf{r})d^{3}\mathbf{r},$$

(8)

where the first term represents the kinetic energy:

$$\int\frac{\mathbf{p}^{2}}{2m}f(\mathbf{r},\mathbf{p})\frac{d^{3}\mathbf{r}d^{3}\mathbf{p}}{(2\pi\hbar)^{3}}=\int\frac{3}{2\beta}n(\mathbf{r})d^{3}\mathbf{r}=\frac{3}{2}\frac{N}{\beta}.$$

(9)

The electrostatic potential $\varphi(\mathbf{r})$ is related to the particles density $n(\mathbf{r})$ throughout Poisson equation, $\Delta\varphi(\mathbf{r})=-4\pi qn(\mathbf{r})$. Introducing the effective field $u(\mathbf{r})$:

$$u(\mathbf{r})=\frac{1}{2}m\omega^{2}\mathbf{r}^{2}+q\varphi(\mathbf{r}),$$

(10)

and the dimensionless potential $\Phi(r)$:

$$\Phi(\mathbf{\mathbf{r}})=\beta\left[u(0)-u(\mathbf{r})\right],$$

(11)

the particles density can be rewritten as:

$$n(\mathbf{r})=n_{0}\exp\left[\Phi(\mathbf{r})\right],$$

(12)

where $n_{0}$ is the central density:

$$n_{0}=\left(\frac{2m\pi}{\hbar^{2}\beta}\right)^{\frac{3}{2}}A\exp\left[-\beta u(0)\right].$$

(13)

Additionally, it is convenient to introduce the characteristic radius constant $r_{c}$:

$$\frac{1}{r_{c}^{2}}=q^{2}n_{0}\beta,$$

(14)

the dimensionless radius variable $\xi=r/r_{c}$, and the auxiliary constant $\lambda$:

$$\lambda=\frac{3m\beta\omega^{2}}{4\pi}r_{c}^{2}=\frac{3m\omega^{2}}{4\pi q^{2}n_{0}}.$$

(15)

The above definitions enable us to rephrase Poisson equation with spherical symmetry as follows:

$$\frac{1}{\xi^{2}}\frac{d}{d\xi}\left[\xi^{2}\frac{d}{d\xi}\Phi(\xi)\right]=4\pi\left[e^{\Phi(\xi)}-\lambda\right],$$

(16)

which should be solved by numerical integration using the following boundary conditions at the origin:

$$\Phi(0)=0,\>\frac{d}{d\xi}\Phi(0)=0.$$

(17)

Additionally, it is necessary to take into account the asymptotic behavior of the dimensionless potential for large distances. For spherical solutions, the electrostatic potential can be expressed as follows:

$$\frac{d}{dr}\varphi(r)=-\frac{qN(r)}{r^{2}},$$

(18)

where $N(r)$ is the number of particles enclosed inside a sphere of radius $r$. Considering the first derivative of the dimensionless potential:

$$\frac{d\Phi(\xi)}{d\xi}=\beta r_{c}\left[m\omega^{2}r+q\frac{d\varphi(r)}{dr}\right],$$

(19)

and introducing the dimensionless parameter $\eta$:

$$\eta=\frac{q^{2}N\beta}{r_{c}},$$

(20)

one obtains the following asymptotic expression for $\eta$:

$$\eta=\lim_{\xi\rightarrow+\infty}\left[\frac{4\pi}{3}\lambda\xi^{3}+\xi^{2}\frac{d}{d\xi}\Phi(\xi)\right].$$

(21)

At zero temperature, the pure non-neutral plasma in a magnetic trap adopts the called Brillouin steady state [7], a profile of uniform density where electrostatic repulsion forces are exactly compensated with electric and magnetic fields of the Penning trap. An analogous situation also appears in the zero temperature limit for the present model. Considering the vanishing of the resulting force when $T\rightarrow 0$:

$$\mathbf{F}(\mathbf{r})=q\mathbf{E}(\mathbf{r})-m\omega^{2}\,\mathbf{r}=0,$$

(22)

one can apply the divergence to obtain the Brillouin density:

$$n_{B}=\frac{3m\omega^{2}}{4\pi q^{2}}$$

(23)

associated with the presence of the external harmonic field with frequency $\omega$. This profile of constant density is not infinitely extended because of the number of particles $N$ is bound. The right expression is given by:

$$n(r)=\left\{\begin{array}[]{cc}n_{B},&\mbox{ if }r<R_{B},\\ 0,&\mbox{ if }r\geq R_{B},\\ \end{array}\right.$$

(24)

where $R_{B}$ is the Brillouin radius:

$$n_{B}\frac{4\pi}{3}R^{3}_{B}=N,$$

(25)

which defines the linear size of the system at zero temperature.

For the sake of convenience, we shall refer numerical results of this study by using a set of characteristic units associated with Brillouin state. The linear distances and particles densities will be referred into units of Brillouin radius $R_{B}$ and density $n_{B}$. Using these units, the particles density (12) can be rewritten as:

$$\bar{n}(\bar{r})=\frac{1}{\lambda}\exp\left[\Phi(\bar{r})\right],$$

(26)

where $\bar{n}(r)=n(r)/n_{B}$ and $\bar{r}=r/R_{B}=\xi r_{c}/R_{B}$. It was considered here the relation:

$$\frac{n_{0}}{n_{B}}=\frac{1}{\lambda},$$

(27)

which is derived from definition (15). According to Eq.(26), the Brillouin steady state appears in the limit $\lambda\rightarrow 1$, where $\bar{n}(\bar{r})=1$ and $\Phi(\bar{r})=0$. Additionally, one can introduce the characteristic units:

$$U_{B}=\frac{q^{2}N^{2}}{R_{B}}\mbox{ and }T_{B}=\frac{q^{2}N}{R_{B}k}$$

(28)

for the energy and the temperature, respectively. The inverse temperature $\beta$ and the characteristic radius $r_{c}$ defined in Eq.(14) can be rewritten into Brillouin units as:

$$\bar{\beta}=\frac{T_{B}}{T}=\left(\frac{4\pi}{3}\lambda\eta^{2}\right)^{\frac{1}{3}}\mbox{ and }\bar{r}_{c}=\frac{r_{c}}{R_{B}}=\left(\frac{4\pi}{3}\frac{\lambda}{\eta}\right)^{\frac{1}{3}}.$$

(29)

Let us now obtain the working expressions for the total energy (8) in Brillouin units. The kinetic energy contribution is expressed as:

$$\bar{K}=\frac{K}{U_{B}}=\frac{3}{2}\frac{1}{\bar{\beta}}.$$

(30)

The total potential energy associated with the external harmonic field:

$$\frac{W}{U_{B}}=\frac{1}{2}\frac{\bar{r}_{c}^{3}}{\bar{\beta}}\mathbb{Q},$$

(31)

where $\mathbb{Q}$ denotes the auxiliary integral:

$$\mathbb{Q}=\int^{+\infty}_{0}\xi^{2}\exp[{\Phi(\xi)}]4\pi\xi^{2}d\xi.$$

(32)

Finally, let us obtain the expression for the total electrostatic energy. Rephrasing this contribution in terms of dimensionless potential $\Phi(\mathbf{r})$:

$$V=\frac{1}{2}Nq\varphi(0)-\frac{1}{2}W-\int\frac{1}{2\beta}\Phi(\mathbf{r})n(\mathbf{r})d^{3}\mathbf{r}$$

(33)

and rewriting it into Brillouin energy units $U_{B}$, one obtains:

$$\bar{V}=\frac{V}{U_{B}}=\frac{3}{8\pi}\frac{\bar{r}_{c}^{2}}{\lambda}\mathbb{C}-\frac{1}{4}\frac{\bar{r}_{c}^{3}}{\bar{\beta}}\mathbb{Q}-\frac{3}{8\pi}\frac{\bar{r}_{c}^{3}}{\lambda\bar{\beta}}\mathbb{P},$$

(34)

where $\mathbb{C}$ and $\mathbb{P}$ are the following auxiliary integrals:

	$\displaystyle\mathbb{C}$	$\displaystyle=$	$\displaystyle\int^{+\infty}_{0}\exp[{\Phi(\xi)}]4\pi\xi d\xi,$		(35)
	$\displaystyle\mathbb{P}$	$\displaystyle=$	$\displaystyle\int^{+\infty}_{0}\Phi(\xi)\exp[{\Phi(\xi)}]4\pi\xi^{2}d\xi.$		(36)

The sum of all these contributions is:

$$\bar{U}=\frac{U}{U_{B}}=\frac{3}{2}\frac{1}{\bar{\beta}}+\frac{1}{4}\frac{\bar{r}_{c}^{3}}{\bar{\beta}}\mathbb{Q}+\frac{3}{8\pi}\frac{\bar{r}_{c}^{2}}{\lambda}\left(\mathbb{C}-\frac{\bar{r}_{c}}{\bar{\beta}}\mathbb{P}\right).$$

(37)

Hamiltonian equations for the model (3) are given by:

$$\dot{\mathbf{r}}_{i}=\frac{1}{m}\mathbf{p}_{i},\;\dot{\mathbf{p}}_{i}=-m\omega^{2}\mathbf{r}_{i}+\sum_{j\neq i}\frac{q^{2}}{\left|\mathbf{r}_{ij}\right|^{3}}\mathbf{r}_{ij},$$

(38)

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ is the separation vector oriented from $j$-th to $i$-th particle. Numerical integration of this conservative system is better performed using some sympletic algorithm [43, 44]. An efficient and precise sympletic algorithm was proposed by Casetti [45]:

	$\displaystyle\tilde{\mathbf{r}}_{i}=\mathbf{r}_{i}(t),\>\tilde{\mathbf{p}}_{i}=\mathbf{p}_{i}(t)-\frac{1}{2}\Delta t\frac{\partial}{\partial\tilde{\mathbf{r}}_{i}}V(\tilde{\mathbf{r}}),$		(39)
	$\displaystyle\mathbf{r}_{i}(t+\Delta t)=\tilde{\mathbf{r}}_{i}+\Delta t\frac{1}{m}\tilde{\mathbf{p}}_{i},\>\mathbf{p}_{i}(t+\Delta t)=\tilde{\mathbf{p}}_{i}-\frac{1}{2}\Delta t\frac{\partial}{\partial{\mathbf{r}}_{i}}V\left[\mathbf{r}(t+\Delta t)\right],$
	$\displaystyle\hat{\mathbf{p}}_{i}=\mathbf{p}_{i}(t+\Delta t),\>\hat{\mathbf{r}}_{i}=\mathbf{r}_{i}(t+\Delta t)+\frac{1}{2m}\Delta t\hat{\mathbf{p}}_{i},$
	$\displaystyle\mathbf{p}_{i}(t+2\Delta t)=\hat{\mathbf{p}}_{i}-\Delta t\frac{\partial}{\partial\hat{\mathbf{r}}_{i}}V\left[\hat{\mathbf{r}}\right],\>\mathbf{r}_{i}(t+2\Delta t)=\hat{\mathbf{r}}_{i}+\frac{1}{2m}\Delta t\mathbf{p}_{i}(t+2\Delta t),$

which was employed by Pettini and co-workers in the astrophysical situation [41, 46]. This same algorithm will be considered in this work.

The instability of trajectories can be studied using the tangent dynamics, that is, the linearization of separation dynamics of two trajectories that are infinitely close:

$$\dot{\boldsymbol{\xi}}_{i}=\frac{1}{m}\boldsymbol{\eta}_{i},\>\dot{\boldsymbol{\eta}}_{i}=-m\omega^{2}\boldsymbol{\xi}_{i}+\sum_{j\neq i}\frac{q^{2}}{\left|\mathbf{r}_{ij}\right|^{3}}\left\{\boldsymbol{\xi}_{ij}-3\left(\mathbf{n}_{ij}\cdot\boldsymbol{\xi}_{ij}\right)\mathbf{n}_{ij}\right\},$$

(40)

where $\mathbf{n}_{ij}$ and $\boldsymbol{\xi}_{ij}$ are given by:

$$\mathbf{n}_{ij}=\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|}\mbox{ and }\boldsymbol{\xi}_{ij}=\boldsymbol{\xi}_{i}-\boldsymbol{\xi}_{j}.$$

(41)

Numerical integration of the above equations is coupled to Hamilton equations (38). For the numerical integration of the tangent dynamics, we shall employ a second-order Euler scheme:

	$\displaystyle\boldsymbol{\xi}_{i}\left(t+\Delta t\right)$	$\displaystyle=$	$\displaystyle\boldsymbol{\xi}_{i}\left(t\right)+\frac{1}{m}\boldsymbol{\eta}_{i}\left(t\right)\Delta t+\frac{1}{2m}\dot{\boldsymbol{\eta}}_{i}\left(t\right)\Delta t^{2}+\frac{1}{6m}\ddot{\boldsymbol{\eta}}_{i}\left(t\right)\Delta t^{3},$
	$\displaystyle\boldsymbol{\eta}_{i}\left(t+\Delta t\right)$	$\displaystyle=$	$\displaystyle\boldsymbol{\eta}_{i}\left(t\right)+\dot{\boldsymbol{\eta}}_{i}\left(t\right)\Delta t+\frac{1}{2m}\ddot{\boldsymbol{\eta}}_{i}\left(t\right)\Delta t^{2},$		(42)

which requires the second time derivative $\ddot{\boldsymbol{\eta}}_{i}$ of speed separation $\boldsymbol{\eta}_{i}$:

	$\displaystyle\ddot{\boldsymbol{\eta}}_{i}$	$\displaystyle=-{\omega}^{2}{\boldsymbol{\eta}}_{i}+\sum_{j\neq i}\frac{q^{2}}{\left|\mathbf{r}_{ij}\right|^{3}}\left[-3(\mathbf{n}_{ij}\cdot\boldsymbol{\kappa}_{ij})\left\{\boldsymbol{\xi}_{ij}-3\left(\mathbf{n}_{ij}\cdot\boldsymbol{\xi}_{ij}\right)\mathbf{n}_{ij}\right\}+\dot{\boldsymbol{\xi}}_{ij}\right.$			(43)
		$\displaystyle\qquad\qquad\left.-3\left(\mathbf{n}_{ij}\cdot\dot{\boldsymbol{\xi}}_{ij}\right)\mathbf{n}_{ij}-3\left(\dot{\mathbf{n}}_{ij}\cdot\boldsymbol{\xi}_{ij}\right)\mathbf{n}_{ij}-3\left(\mathbf{n}_{ij}\cdot\boldsymbol{\xi}_{ij}\right)\dot{\mathbf{n}}_{ij}\right].$

For the sake of briefly, we have considered here the following notations:

$$\boldsymbol{\kappa}_{ij}=\frac{\mathbf{v}_{ij}}{|\mathbf{r}_{ij}|},\>\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}\mbox{ and }\dot{\mathbf{n}}_{ij}=\left\{\boldsymbol{\kappa}_{ij}-(\mathbf{n}_{ij}\cdot\boldsymbol{\kappa}_{ij})\mathbf{n}_{ij}\right\}.$$

(44)

Finally, the calculation of Lyapunov exponent is performed considering the limit:

$$\lambda=\lim_{t\rightarrow+\infty}\frac{1}{t}\ln\frac{\left\|\Delta x(t)\right\|}{\left\|\Delta x(0)\right\|},$$

(45)

where the norm of deviation of trajectories $\left\|\Delta x(t)\right\|$ is defined as follows:

$$\left\|\Delta x(t)\right\|^{2}=\sum_{i}\frac{1}{2m}\boldsymbol{\eta}^{2}_{i}(t)+\frac{1}{2}m\omega^{2}\boldsymbol{\xi}^{2}_{i}(t).$$

(46)

For a computational viewpoint, the calculation the limit (45) can be rephrased as the average of a discrete time series with period $\tau$. Denoting the $n$-th time instant $t_{n}=n\tau$ and the corresponding exponential dispersion $\hat{\lambda}_{n}$ as:

$$\hat{\lambda}_{n}=\frac{1}{\left(t_{n}-t_{n-1}\right)}\ln\frac{\left\|\Delta x(t_{n})\right\|}{\left\|\Delta x(t_{n-1})\right\|},$$

(47)

the momentum of $m$-th order $\lambda^{(m)}_{n}$ of the exponential dispersion $\hat{\lambda}_{n}$ can be expressed as follows:

$$\lambda^{(m)}_{n}=\lambda^{(m)}_{n-1}+\frac{1}{n}\left[\left(\hat{\lambda}_{n}\right)^{m}-\lambda^{(m)}_{n-1}\right].$$

(48)

Accordingly, the infinite time limit (45) is now replaced by its discrete version as:

$$\lambda=\lim_{n\rightarrow+\infty}\lambda^{(1)}_{n}.$$

(49)

Of course, the Lyapunov exponent $\lambda$ can be estimated by the $n$-th value $\lambda^{(1)}_{n}$. Its statistical error $\delta\lambda_{n}$ can be expressed in terms of the square dispersion $\sigma^{2}_{n}=\lambda^{(2)}_{n}-\left(\lambda^{(1)}_{n}\right)^{2}$ as follows:

$$\delta\lambda_{n}=\sigma_{n}/\sqrt{n\tau/\tau_{d}},$$

(50)

with $\tau_{d}$ being the decorrelation time (the necessary time interval that the dynamics produces a statistical independent microscopic configuration). Taking into consideration the physical meaning of Lyapunov exponent $\lambda$, one could estimate the decorrelation time $\tau_{d}$ using the chaotization time, $\tau_{d}\sim\tau_{ch}\sim 1/\lambda$. This consideration enable us to estimate the relative error as:

$$\frac{\delta\lambda_{n}}{\lambda^{(1)}_{n}}=\frac{\sigma_{n}}{\lambda^{(1)}_{n}}\frac{1}{\sqrt{n\lambda^{(1)}_{n}\tau}}.$$

(51)

For a general physical observable $O$, e.g.: the total kinetic energy $K$, the calculation of its temporal expectation can be implemented using discrete time series as follows:

$$O_{n}=O_{n-1}+\frac{1}{n}\left(\hat{O}_{n}-O_{n-1}\right),$$

(52)

where $\hat{O}_{n}=O\left(t_{n}\right)$ is its corresponding value at the time instant $t_{n}$.

The short-range interacting systems obey the extensive thermodynamic limit $N\rightarrow\infty$:

$$N\rightarrow\infty\Rightarrow\frac{S}{N}=\emph{const},\,\frac{U}{N}=\emph{const},\,\frac{V}{N}=\emph{const},$$

(53)

which guarantees the intensive character of some relevant observables and thermodynamic parameters like temperature and particles density. For systems with long-range interactions, the additivity and separability of extensive systems are non necessarily applicable. This type of arguments should be considered with care for each particular situation. For the concrete case of astrophysical system composed of non-relativistic point particles that interact among them through Newtonian gravitation, it has been claimed the relevance of the following thermodynamic limit [47]:

$$N\rightarrow\infty\Rightarrow\frac{S}{N}=\emph{const},\,\frac{U}{N^{\frac{7}{3}}}=\emph{const},\,VN=\emph{const},$$

(54)

which accounts for a non-extensive behavior. One should expect a similar thermodynamic limit for the case of a pure non-neutral plasma considered in this study. In order to check this possibility, let us consider the arguments already employed in the astrophysical case. Considering the characteristic linear dimension of the present situation as the Brillouin radius $r_{c}\sim R_{B}$, as well as the characteristic momentum $p_{c}$ derived from the Brillouin energy $U_{B}=q^{2}N^{2}/R_{B}$:

$$p_{c}=\sqrt{\frac{2mU_{B}}{N}},$$

(55)

these quantities can be employed to estimate the volume of phase space $\Omega$ as:

$$\Omega\sim\frac{1}{N!}\left(\frac{p_{c}r_{c}}{\hbar}\right)^{3N},$$

(56)

which leads to the following estimation for the entropy $S=k\ln\Omega$:

$$S\sim\frac{Nk}{2}\ln\left(\frac{8m^{3}q^{6}NR_{B}^{3}e^{2}}{\hbar^{6}}\right).$$

(57)

According to the expression (57), the entropy per particle $S/N$ remains finite when $N\rightarrow\infty$ after assuming the non-extensive scaling laws:

$$U/N^{7/3}=const\mbox{ and }RN^{1/3}=const.$$

(58)

This limit evidences the non-extensive character of the pure non-neutral plasma described by the Hamiltonian (1), which is essentially the same behavior associated with astrophysical case. Considering the scaling behavior of the particles density $n$, one can obtain the scaling behavior of the frequency $\omega$ of the external harmonic field

$$n/N^{2}=const\Rightarrow\omega/N=const,$$

(59)

where the expression of Brillouin density (23) was taken into consideration.

Since the number of particles of the model system under consideration is finite, the particles density (12) must be a bound function. For large distances, the particles density should be a monotonous decreasing function. Analysing the Poisson problem (16), the behaviour of the density profile, or even the one associated with the dimensionless potential $\Phi(\xi)$, is driven by the sign of the function $(\exp\left[\Phi(\xi)\right]-\lambda)$. If this function is positive definite, the dimensionless potential will be a monotonous increasing function, and hence, this possibility cannot describe a density profile with a finite number of particles. Therefore, this function has to be negative definite in order to describe a monotonous decreasing function for the dimensionless $\Phi(\xi)$. The maximum will always take place at the origin $\xi=0$, so that, the integration parameter $\lambda\geq 1$. The lower bound $\lambda=1$ corresponds to the Brillouin steady state, that is, the zero temperature limit. For the sake of convenience, let us introduce the auxiliary parameter $\delta=\ln(\lambda)=\ln\left(n_{B}/n_{0}\right)$, which is a measure of the central density $n_{0}$ in units of the Brillouin density $n_{B}$.

We show in figure 1 the temperature dependencies (in units of the Brillouin temperature $T_{B}$) of two relevant observables, the auxiliary parameter $\delta$ and the specific heat (in units of the Boltzmann constant $k$). According to these results, the auxiliary parameter $\lambda$ exhibits a wide range of values to describe the physics associated with the present situation. By itself, the mathematical behavior of the parameter $\delta$ evidences difficulties of numerical integration of Poisson problem (16) in the low temperature region. On the other hand, temperature dependence of the specific heat (per particles) evidences the existence of two asymptotic limits: the zero temperature limit where appears the Brillouin steady state, and the high temperature limit, where the system behaves as a ideal gas of harmonic oscillators.

In the low temperature limit, we observe that the specific heat approaches to the asymptotic value $C/Nk=3/2$, which corresponds to the heat capacity of an ideal gas at constant volume $V_{B}=4\pi R^{3}_{B}/3$. Such a capacity is explained by the contribution of kinetic degrees of freedom only. The repulsive electrostatic forces among the charged particles is fully compensated with the linear forces of the external harmonic field. The contribution of spatial degrees of freedom to the heat capacity is vanishing at zero temperature. In the limit of high temperatures, the density of the system decreases so much that electrostatic potential energy is almost negligible, so that, the system behaves here as an ideal gas of harmonic oscillators. According to the equipartition theorem, the specific heat capacity asymptotically approaches the value $C/Nk=3$. During the transition from the low to high temperature limits, one observes an effective unfrozen of the oscillatory degrees of freedom. It is worthy to comment that the present behavior cannot be associated to the occurrence of a phase transition. Actually, it is rather analogous to the unfrozen of oscillatory and vibrational degrees of freedom in gases and solids due to quantum effects [42].

Let us discuss now some other relevant thermodynamic observables. We have highlighted in figure 1 a series of states (the red points labelled with Latin letters $a-k$). The same ones were uniformly located between the two asymptotic values of the specific heat capacity. For these points, we have obtained the radial distribution profiles $n\left(r\right)$ shown in figure 2. The step function (24) associated to the Brillouin state is established when the temperature $T$ approaches to zero. For nonzero temperatures, the radial distribution profiles turn smooth functions (without discontinuities), which decreases monotonically with the growth of radial coordinate $r$. The associated central density $n_{0}$ of these profiles decreases with the growth of temperature (this behavior is also evidenced by the temperature dependence of the quantity $\delta=\ln(n_{B}/n_{0})$ shown in figure 1). At the limit of hight temperatures $T\gg T_{B}$, the radial profiles turn the Gaussian profile:

$$n_{G}\left(r\right)=\left(\frac{m\omega^{2}}{4\pi kT}\right)^{3/2}N\exp\left[-m\omega^{2}r^{2}/2kT\right]$$

(60)

associated to an ideal gas of harmonic oscillators. Accordingly, the central density $n_{0}$ adopts in this asymptotic limit the analytical form $n_{0}=n_{B}\left(T_{B}/T\right)^{3/2}/6\sqrt{\pi}$. We show in figure 3 the dependencies of three contributions of the total energy $U=K+V_{ext}+V_{elect}$ (the total kinetic energy $K$, the potential energy of the harmonic field $V_{ext}$ and the total electrostatic energy $V_{elect}$) versus the total energy $U$. For low energies, the kinetic energy drops to zero, while both potential energy contributions exhibit finite values due to the system adopts the Brillouin profile (24). For low energies, the electrostatic potential energy approaches to zero and there is an asymptotic equipartition between the kinetic energy and the potential energy of the external harmonic field.

The non-neutral plasma model (3) does not exhibit microcanonical phase transitions neither negative heat capacities reported in other long-range interacting systems [48, 49, 50]. For a better understanding, let us perform a brief comparison of this non-neutral plasma model with its gravitational counterpart:

$$H(\mathbf{r},\mathbf{p})=\sum_{i}\frac{1}{2m}\mathbf{p}^{2}_{i}+\frac{1}{2}m\omega^{2}\mathbf{r}^{2}_{i}-\sum_{i<j}\frac{Gm^{2}}{\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|}.$$

(61)

For the sake of briefly, we shall avoid to enter into mathematical details concerning the thermo-statistical analysis of this second system (the mathematical treatment is quite similar). The gravitational equivalents of Brillouin density and radius are given by:

$$n_{B}=\frac{3\omega^{2}}{4\pi Gm}\mbox{ and }R_{B}=\left(GM/\omega^{2}\right)^{1/3}.$$

(62)

The gravitational radius $R_{B}$ represents the radius of the region where gravitation of this system dominates its microscopic dynamics. The gravitational radius $R_{B}$ is considered to defined the associated units for energy and temperature:

$$U_{B}=\frac{GM^{2}}{R_{B}}\mbox{ and }T_{B}=\frac{GMm}{kR_{B}}.$$

(63)

We show in figure 4 the temperature versus energy dependence (the called microcanonical caloric curve) of both models using their respective characteristic units. As naturally expected, these models exhibit the same asymptotic limit for large energies (or large temperatures), where they recover the thermodynamic behavior of the ideal gas of harmonic oscillators. However, they disagree in their respective behaviors at low energies.

The repulsive character of Coulomb forces among the charge particles of the non-neutral plasma enables the existence of Brillouin steady state at zero temperature. On the contrary, the attractive character of gravitation suppress the possibility that the model (61) reaches the zero temperature point. This system cannot exhibit temperatures below the critical value $T_{\beta}=0.43T_{B}$, precisely, because of the corresponding caloric curve exhibits a minimum at the point $U_{\beta}=0.215U_{B}$. This models also exhibit a minimum energy at the value $U_{\alpha}=0.016U_{B}$. It is worth noticing that the heat capacity $C=dU/dT$ diverges at the critical point $U_{\beta}$. For energies within the interval $U_{\alpha}<U<U_{\beta}$, the system temperature $T$ grows when the energy $U$ is decreased, thus evidencing the existence of a branch with negative heat capacities. At the critical point $U_{\alpha}$, the (negative) heat capacity approaches to zero, and the system undergoes the called gravothermal collapse [51]. According to the common understanding of this collective phenomenon (a microcanonical phase transition), the internal pressures of the system are unable to balance its own gravitational field, so that, the system undergoes a collapse that leads the formation of structures with very dense cores. The description of these post-collapse configurations, however, requires additional physical considerations beyond the classical statistical mechanics description of the model (61), e.g.: consideration of quantum effects, the effective linear size of constituents, etc.

Despite of the obvious differences, the non-neutral plasma (3) and its gravitational counterpart (61) share many analogies. Both the Coulomb and Newtonian forces do not exhibit characteristic lengths, so that, the stability of their thermodynamic description crucially depends on the existence of external factors, e.g.: the dynamic influence of the potential harmonic field. In particular, the thermodynamic description of both model systems depend on the frequency constant $\omega$ of the harmonic field (this constant parameter, in particular, determines the Brillouin units $R_{B}$, $n_{B}$, $U_{B}$, $T_{B}$ that enter in their low energy thermodynamic behaviors). Without the external influence of the harmonic field, the constituents will escape from both systems, thus ruling out the occurrence of a rigorous thermodynamic equilibrium. Finally, both long-range interacting systems belong to the same class of non-extensive systems, since their observables follow in the thermodynamic limit $N\rightarrow+\infty$ the same scaling laws (58) and (59)¹¹1Readers may object that derivation of the scaling laws (58) and (59) crucially depends on whether one demands or not the extensive character of the entropy. Actually, the same scaling behavior appears in the framework of the Thomas-Fermy theory for the asymptotic behaviors of atomic energy $U(Z)$ and density $\rho^{Z}(r)$ for an atom of charge $Z$ sufficiently large [52, 53, 54], $U\left(Z\right)=C_{TF}Z^{7/3}+O\left(Z\right)$ and $\rho_{TF}^{Z}\left(r\right)=Z^{2}\varrho\left(Z^{1/3}r\right)$. In general, this type of scaling behavior corresponds to a three-dimensional system of non-relativistic point particles that interact among them throughout $1/r$ potential without mattering about their quantum or classical behavior..

In the framework of non-neutral plasmas, the collisions among the constituting particles is a mechanism that contributes both the chaoticity of microscopic dynamics as well as the system relaxation (at macroscopic level) towards equilibrium. If collision events constitute the main source behind the chaoticity of the present model (1), the associated Lyapunov exponent $\lambda$ should be comparable to the frequency of collisions $\nu$ (collisions per unit of time). For the sake of convenience, let us perform a qualitative description for the present situation, which will be considered later to analyze the results obtained from the dynamical simulations. An elementary estimation is obtained in terms of the effective cross section $\sigma$ and the gas density $n$ as follows, $\nu\sim N\sigma nv_{rel}$, with $v_{rel}$ being the relative velocity among particles. The cross section $\sigma$ among charged particles can be estimated as the Rutherford scattering cross section, $\sigma\sim q^{4}/m^{2}\left\langle v^{2}\right\rangle^{2}$, with $\left\langle v^{2}\right\rangle$ being the square dispersion of velocity. Denoting the mean square velocity as $v=\sqrt{\left\langle v^{2}\right\rangle}$, the temperature dependence of both the mean square velocity $v$ the relative velocity among particles is $v_{rel}\sim\sqrt{2}v\propto\sqrt{kT/m}$. Taking into account the asymptotic dependence of the central density $n_{0}$ for low and large temperatures ($n_{0}\sim n_{B}$ and $n_{0}\propto n_{B}\left(T_{B}/T\right)^{3/2}$, respectively), one obtains the estimation

$$\nu\propto\left\{\begin{array}[]{cc}\omega(T_{B}/T)^{3/2},&\mbox{ if }T/T_{B}\ll 1,\\ \omega(T_{B}/T)^{3},&\mbox{ if }T/T_{B}\gg 1,\\ \end{array}\right.$$

(64)

with $\omega$ being the frequency parameter of the external harmonic field. The previous result predict that the frequency of collisions $\nu$ decreases with the growth of temperature, but it exhibits different behaviors in the low and high temperature limits. Moreover, the characteristic timescale between collision events $\tau_{col}=1/\nu$ exhibits the same order of characteristic timescale $\tau_{dyn}=1/\omega$ of microscopic dynamics. Accordingly, the system as a whole will be affected by the incidence of collisions in the timescale $\tau_{c.rlx}\propto N\tau_{col}$, which represents the characteristic timescale of the collisional relaxation.

One can suppose that the model system (1) obeys ergodicity and mixing properties that are necessary to perform a statistical description in terms of microcanonical ensemble. Besides, one can also admit that the system will arrive at an thermodynamic equilibrium by starting from an arbitrary initial condition in a finite relaxation time $\tau_{relax}<+\infty$. Since Lyapunov exponent (45) requires infinite time limit $t\rightarrow+\infty$, this indicator will characterize the equilibrium state. In our dynamical simulations, we have started the dynamical evolution from equilibrium situation associated with Maxwell-Boltzmann distribution (4). Proceeding thus, we shall avoid the system initial evolution towards equilibrium and speed up all calculations.

During random generation of initial conditions using distribution (4), the statistical fluctuations are large for values of $N$ relatively small, overall, for low values of temperatures. In these cases, the resulting initial configuration can significantly differ from the equilibrium configuration of interest. For this reason, we have considered a thermalization procedure to force the system to reach a given equilibrium state with temperature $T$. Along this process, the kinetic energy of each particles is periodically re-scaled in a way that the expectation value of the dynamical temperature $T_{dyn}=2K/3N$ acquires the value $T$ of interest, $\left\langle T_{dyn}\right\rangle=T$. A particular illustration of this procedure for the system with size $N=64$ is shown in figure 5, which was developed to force the equilibrium state with $T=10^{-5}T_{B}$.

Our first task is to contrast the predictions of the thermo-statistical analysis and the results obtained from the dynamical simulations. We show in figure 6 a direct comparison among the temporal expectation values

$$\left\langle O\right\rangle_{dyn}=\lim_{T\rightarrow+\infty}\frac{1}{T}\int^{T}_{0}O\left[\mathbf{q}(t),\mathbf{p}(t)\right]dt$$

(65)

of the three contributions of the total energy $U$ (points) and their corresponding statistical expectation values (solid lines):

$$\left\langle O\right\rangle_{stat}=\int O\left(\mathbf{q},\mathbf{p}\right)dP\left(\mathbf{q},\mathbf{p}\right)$$

(66)

already shown in figure 3. The agreement of these predictions is excellent taking into consideration that our dynamical simulations were restricted to a number of particles $N=2048$, while the statistical estimations were performed by invoking thermodynamic limit $N\rightarrow+\infty$. The present results evidence (1) the applicability of the thermodynamic limit approximation $N\rightarrow+\infty$ to obtain statistical expectation values for finite systems with moderate number of particles $N\sim 10^{3}$ and (2) the ergodic character of the non-neutral plasma model (3) considered in the present study, namely, the equalization between the statistical and dynamical expectation values of a given observable $O$:

$$\left\langle O\right\rangle_{dyn}=\left\langle O\right\rangle_{stat}.$$

(67)

Let us now analyse the relationship between the chaoticity of the microscopic dynamics and its thermodynamic behavior. Specifically, we have performed extensive simulations to obtain the dependence of Lyapunov exponent $\lambda$ on the system temperature $T$. Results for a number of particles $N=2048$ are shown in figure 7. In our simulations, the time variable $t$ was referred to into units of $\tau_{c}=1/\omega$, and hence, the Lyapunov exponent $\lambda$ (the maximum one) is expressed to in units of constant frequency $\omega$ of the external harmonic field. For comparative purposes, we have also included in this figure the dependence of heat capacity $C$ versus temperature $T$ already shown in figure 1. According to these results, the Lyapunov exponent decreases (and probably drops to zero) when the system approaches its two asymptotic limits of low and high temperatures, while it exhibits a local maximum during the transition around the temperature value $T\simeq 0.1T_{B}$. The observed temperature dependence of the Lyapunov exponent differs from the one considered by the rate of collisions (64). Of course, one cannot expect a direct identification between these quantities, but the growth of Lyapunov exponent in the low temperature limit cannot be explained in terms of the particles collisions. The mechanism of collisions turns more effective for low temperatures, which is in contradiction with the reduction of the system chaoticity observed when the temperature decreases.

For a better understanding of the system chaoticity, we have studied the dependence of the Lyapunov exponent $\lambda$ on the system size $N$. Results of extensive simulations considering the ranges $N=64-2048$ and $T/T_{B}=\left[10^{-5}-10\right]$ are shown in figure 8. According to these results, the Lyapunov exponent exhibits the same qualitative dependence on the temperature $T$, but the overall values of this quantity grow with the system size $N$. A simple analysis evidences that the Lyapunov exponent $\lambda$ (in units of $\omega$) grows with $N$ following the power-law $N^{1/4}$. This scaling law is evidenced by the data-collapse after re-scaling $\lambda/\omega$ by the factor $N^{1/4}$ (the bottom panel of figure 8).

The present results evidence that the characteristic chaotization timescale of this non-neutral plasma model is $\tau_{ch}=1/\lambda\propto 1/\omega N^{1/4}=\tau_{dyn}/N^{1/4}$, which is considerably smaller than the characteristic timescale of its microscopic dynamics $\tau_{dyn}$, $\tau_{ch}\ll\tau_{dyn}$. The observed chaoticity cannot be explained in terms of particles collisions because of their characteristic timescales considerably differ between them. The size dependence of the Lyapunov exponent should be explained by some type of collective influence of the system as a whole. Taking into consideration precedent studies in the context of astrophysical model [46, 40], the observed chaoticity should be explained in terms of the phenomenon of parametric resonance. The verification of this hypothesis requires the application of Riemannian approach of Hamiltonian chaos [41], which is beyond the scope of the present study.

The chaoticity of a Hamiltonian system is understood as a consequence of non-linearity of its microscopic dynamics. Apparently, the incidence of non-linear effects in this concrete situation reaches its maximum during the transition between the two asymptotic limits of the system thermodynamic behavior. According to the heat capacity $C$ versus temperature dependence, during the transit between Brillouin limit towards the ideal gas of harmonic oscillators limit, it takes place the unfreezing of the oscillatory degrees of freedom. Therefore, the non-linear effects of microscopic dynamics that explain the chaoticity shown in figures 7 and 8 must be in someway associated to this process of unfreezing of the oscillatory degrees of freedom.