text:
# Inelastic co-tunneling through an excited state of a quantum dot

## I Introduction

Transport spectroscopy provides information on the excited states of quantum dots.  Outside the Coulomb blockaded regions (“diamonds”), peaks in the differential conductance can be observed, which linearly shift with the gate voltage.  Such a peak indicates that an excited state has entered the transport voltage window. Recently, De Franceschi et al. observed new features in the spectroscopy data of a vertical semiconductor quantum dot coupled relatively strong to the leads.  Inside several Coulomb diamonds with an even number of electrons they observed a step in the differential conductance with a peak at the step edge for a quantum dot with six electrons. In contrast to the peaks outside the diamond, this feature depends very weakly on the gate-voltage. This was interpreted as being due to inelastic co-tunneling through an excited state of the dot. It was pointed out that the sharpness of the step is determined either the temperature or by the broadening of the excited level, both of which are small energy scales compared to the widths of the first-order tunneling peaks. This provides an new spectroscopic tool for studying excited levels in quantum dots. Similar features were observed in the transport spectroscopy of carbon nanotubes obtained by Nygård et. al..  Here also, the gate voltage independent features at finite bias voltage were clearly seen in consecutive Coulomb diamonds with an even number of electrons. The observed peaks at the step edges were attributed to Kondo-physics. In both experiments a Kondo effect was observed at zero bias voltage in a magnetic field. The occurence of Kondo peaks in the differential conductance, although of different origin, indicates a relatively strong coupling to the leads.

Motivated by these experiments, we have calculated the co-tunneling current through an excited level of a quantum dot using a master equation approach. We investigate the possibility that the observed features arise from changes in the nonequilibrium occupations of the states in the dot. Such changes are caused by voltage-dependence of co-tunneling processes. We find that a peak structure in the differential conductance can develop at the onset of inelastic co-tunneling. We consider the possible regimes defined by the energy of the excited level relative to the addition energy. Two possible cases are distinguished:  (i) for low excitation energy, the excited state can decay only by inelastic co-tunneling due to Coulomb blockade; (ii) for higher excitation energy the excited state can decay by sequential tunneling. We take the spin degrees of freedom into account and derive master equations for the nonequilibrium occupations of the relevant states, which are averaged over the spin projections. The case (i) has recently been considered by M. Eto  without taking the spin degrees of freedom into account.

The plan of the paper is as follows. In Sec. II we introduce a model for the transport spectroscopy of a quantum dot. The low bias elastic co-tunneling current is discussed in Sec. III and in Sec. IV we consider the different regions of the Coulomb diamond at higher bias where inelastic co-tunneling dominates the transport. Our conclusions are presented in Sec. V.

## II Transport spectroscopy

We consider a quantum dot connected to two leads $`L`$ and $`R`$ and one gate electrode $`G`$. The dot is modeled by the Hamiltonian $`H=H_{\text{sp}}+H_{\text{el}}+H_{\text{ex}}`$. The discrete single particle states are labeled by $`d=1,2,\mathrm{}`$ and are spin degenerate, $`s=,`$:

$`H_{\text{sp}}={\displaystyle \underset{ds}{}}\epsilon _dn_{ds},`$

where $`n_{ds}=a_{ds}^{}a_{ds}`$ is the number operator of level $`ds`$. We describe the electron-electron interactions in the classical electron liquid picture (Coulomb blockade model) . The electrostatic energy depends on the total number of electrons $`N=_{ds}n_{ds}`$ in the dot and on the voltages on the electrodes, which couple to the island through the capacitances $`C_L,C_R,C_G`$ (total capacitance $`C=_{k=L,R,G}C_k`$):

$`H_{\text{el}}={\displaystyle \frac{1}{2}}N\left(N1\right)U+\left({\displaystyle \frac{1}{2}}U{\displaystyle \underset{k=L,R,G}{}}{\displaystyle \frac{C_k}{C}}V_k\right)N.`$

Here $`U=1/C`$ is the charging energy (we use units $`e=k=\mathrm{}=1`$). We also include an exchange term that lowers the energy of states with parallel spins:

$`H_{\text{ex}}=J{\displaystyle \underset{dd^{}s}{}}n_{ds}n_{d^{}s}`$

The leads are modeled as noninteracting quasi-particle reservoirs with spin degenerate states. For the left lead with states $`l`$ we write $`H_L=_{ls}\epsilon _ln_{ls}`$ where $`n_{ls}=a_{ls}^{}a_{ls}`$ and $`H_R`$ is defined similarly for the right lead (states labeled by $`r`$). The leads are assumed to be in equilibrium at chemical potential $`\mu _L=\mu _R+V`$ and $`\mu _R`$, respectively. They are coupled to the dot by tunnel junctions. The tunneling Hamiltonian of the left junction is

$`H_{TL}={\displaystyle \underset{ds}{}}\left(H_{TLd}^{}+H_{TLd}\right),H_{TLd}={\displaystyle \underset{ls}{}}t_{ld}a_{ls}^{}a_{ds},`$

and $`H_{TR}`$ is defined similarly for the right junction ($`L,lR,r`$). The temperature $`T`$ is assumed to be small: $`T\mathrm{\Gamma },V,\delta ,`$ where $`\mathrm{\Gamma }`$ is the typical broadening of an energy level in the dot due to the coupling to the leads and $`\delta `$ is the relevant level spacing introduced below.

We consider here the case where the number of electrons on the dot $`N`$ is even at zero bias voltage. Let us first consider the conditions for sequential tunneling transport through the dot. For transitions between the $`N`$ and $`N\pm 1`$ electron ground states only two orbitals need to be considered. Let $`|0`$ denote the $`N2`$ electron ground state and $`d`$ and $`d^{}`$ label the $`N/2^{\text{th}}`$ and $`N/2+1^{\text{st}}`$ orbital with level spacing $`\delta =\epsilon _{N/2+1}\epsilon _{N/2}`$. The singlet ground state then reads $`|N_{0,0}=a_d^{}a_d^{}|0`$. Here the subscripts in $`|N_{S,S_z}`$ indicate the spin $`S`$ on the dot and its projection $`S_z`$ in some arbitrary $`z`$-direction. The doublet ground states with $`N\pm 1`$ electrons are denoted by

$`|\left(N1\right)_{1/2,+1/2}`$ $`=`$ $`a_d^{}|0,|\left(N+1\right)_{1/2,+1/2}=a_d^{}^{}a_d^{}a_d^{}|0,`$ (2)
$`|\left(N1\right)_{1/2,1/2}`$ $`=`$ $`a_d^{}|0,|\left(N+1\right)_{1/2,1/2}=a_d^{}^{}a_d^{}a_d^{}|0,`$ (3)
The differences between the ground state energies $`E_N`$ are

$`E_{NN1}`$ $``$ $`E_NE_{N1}=\epsilon _{N/2}+\left(N{\displaystyle \frac{1}{2}}\right)U{\displaystyle \underset{k=LRG}{}}{\displaystyle \frac{C_k}{C}}V_k,`$
$`E_{N+1N}`$ $``$ $`E_{N+1}E_N=\epsilon _{N/2+1}+\left(N+{\displaystyle \frac{1}{2}}\right)U{\displaystyle \underset{k=LRG}{}}{\displaystyle \frac{C_k}{C}}V_k.`$
In order to have sequential tunneling transport it must be possible to inject an electron from the left lead and emit one to the right lead in a subsequent event:

$`\mu _L`$ $``$ $`E_{NN1}\mu _R`$ (6)
$`\mu _L`$ $``$ $`E_{N+1N}\mu _R`$ (7)
These transport conditions translate into a stability region in the plane of bias and gate voltages, which is shown in Fig. 1. Let us express the addition energies as a function of the position ($`V,V_G`$) in this diagram. In the limit of zero bias voltage (i.e. $`V=0`$ in Fig. 1), sequential tunneling transport is possible only at two isolated points: $`V_G=V_G^{NN1}`$ \[Eq. (6) is satisfied\] and $`V_G=V_G^{N+1N}`$ \[Eq. (7) is satisfied\]. At these degeneracy points the number of electrons on the dot can change by one. At the intermediate gate voltages $`V_G^{NN1}<V_G<V_G^{N+1N}`$ the number of electrons on the dot is stable and equal to $`N`$. Sweeping the gate voltage one thus observes the Coulomb blockade oscillations with peak spacing

$`V_G^{N+1N}V_G^{NN1}=\left(C/C_G\right)\left(\delta +U\right).`$

This involves the level spacing $`\delta `$ because the number of electrons on the dot is even. For a finite, fixed bias voltage each of these degeneracy points develops into a gate-voltage transport window, cf. Fig. 1. Equivalently, the transport conditions (6) and (7) require that the bias voltage exceeds a gate-voltage dependent threshold to allow the number of electrons on the dot to fluctuate. The threshold lines are (depicted in Fig. 1)

$`V_L^{NN1}`$ $``$ $`{\displaystyle \frac{C_G}{CC_L}}\left(V_GV_G^{NN1}\right),`$ (9)
$`V_R^{NN1}`$ $``$ $`+{\displaystyle \frac{C_G}{C_L}}\left(V_GV_G^{NN1}\right),`$ (10)
$`V_L^{N+1N}`$ $``$ $`+{\displaystyle \frac{C_G}{CC_L}}\left(V_G^{N+1N}V_G\right),`$ (11)
$`V_R^{N+1N}`$ $``$ $`{\displaystyle \frac{C_G}{C_L}}\left(V_G^{N+1N}V_G\right).`$ (12)
The region where the $`N`$ electron ground state is stable, the $`N`$th Coulomb diamond, is thus delimited by $`V_L^{NN1}\left(V_G\right),V_R^{N+1N}\left(V_G\right)<V<V_R^{NN1}\left(V_G\right),V_L^{N+1N}\left(V_G\right)`$ for gate voltages $`V_G^{NN1}<V_G<V_G^{N+1N}`$. The maximal bias voltage that can be applied without lifting the Coulomb blockade is

$$V_{\text{max}}^N=\frac{C_G}{C}\left(V_G^{N+1N}V_G^{NN1}\right)=\delta +U.$$
(13)
This is achieved at the intersection of the threshold lines (10) and (11) (the top of the diamond) at gate voltage $`V_{G\text{ max}}^N=V_G^{NN1}+(C_L/C)V_G^{N+1N}V_G^{NN1})`$. Similarly, $`V=V_{\text{max}}^N`$ at the intersection of lines (9) and (12) at $`V_G=V_G^{N+1N}\left(C_L/C\right)(V_G^{N+1N}V_G^{NN1})`$. The addition energies relative to the electro-chemical potentials can now be expressed in terms of distance of the bias voltage to the (gate-voltage dependent) thresholds value

$`E_{NN1}`$ $`=`$ $`\mu _R{\displaystyle \frac{C_L}{C}}\left(V_R^{NN1}\left(V_G\right)V\right)<\mu _R,`$ (15)
$`E_{N+1N}`$ $`=`$ $`\mu _L+\left(1{\displaystyle \frac{C_L}{C}}\right)\left(V_L^{N+1N}\left(V_G\right)V\right)>\mu _L.`$ (16)
Let us now consider the transport conditions for the $`N`$ electron excited state. We assume that the lowest lying excited state is a spin triplet ,

$`|N_{1,+1}^{}`$ $`=`$ $`a_d^{}^{}a_d^{}|0,`$
$`|N_{1,0}^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(a_d^{}^{}a_d^{}+a_d^{}^{}a_d^{}\right)|0,`$
$`|N_{1,1}^{}`$ $`=`$ $`a_d^{}^{}a_d^{}|0.`$
with energy $`E_N^{}=\mathrm{\Delta }+E_N`$. The prime indicates the excited state and the excitation energy $`\mathrm{\Delta }=\delta J`$ is just the spacing between the two active levels $`d`$ and $`d^{}`$ reduced by the exchange energy. We concentrate on the change in the current due to this excited state. The effect of excited states with higher energy (such as the spin singlet at $`E_N+\delta `$) is disregarded. In order to have transport by sequential tunneling, the transitions between the $`N+1`$ and $`N1`$ electron ground states and the $`N`$ electron excited state must be possible:

$`\mu _L`$ $``$ $`E_{N^{}N1}\mu _R,`$ (18)
$`\mu _L`$ $``$ $`E_{N+1N^{}}\mu _R.`$ (19)
where

$$E_{N^{}N1}=E_{NN1}+\mathrm{\Delta },E_{N+1N^{}}=E_{N+1N}\mathrm{\Delta }$$
(20)
The regions in the stability diagram where the conditions (18) and (7) are satisfied, are delimited by another four thresholds lines for the bias voltage:

$`V_\alpha ^{N^{}N1}\left(V_G\right)`$ $`=`$ $`V_\alpha ^{NN1}\left(V_G+{\displaystyle \frac{C}{C_G}}\mathrm{\Delta }\right),`$ (22)
$`V_\alpha ^{N^{}N+1}\left(V_G\right)`$ $`=`$ $`V_\alpha ^{NN+1}\left(V_G{\displaystyle \frac{C}{C_G}}\mathrm{\Delta }\right),`$ (23)
respectively, where $`\alpha =L,R`$. These are simply the regions in Fig. 1 where (18) and (7) are satisfied, shifted in gate voltage by $`\left(C/C_G\right)\mathrm{\Delta }`$ towards the center of the Coulomb diamond. Outside the Coulomb diamond the sequential tunneling transitions to and from the $`N`$ electron excited state gives rise to an increase in the current. By measuring the differential conductance a spectroscopy of excited states can be performed.  Within the Coulomb diamond the $`N`$ electron ground state is stable with respect to sequential tunneling. If the tunnel coupling to the leads is weak, the transport is suppressed within the Coulomb diamond. As the tunnel coupling is enhanced transport becomes possible as a result of higher order tunneling processes. At low bias voltage, $`V<\mathrm{\Delta }`$, only elastic co-tunneling is allowed by energy conservation. Charge is transported through the dot without exciting it. At higher bias voltage, $`V\mathrm{\Delta }`$, inelastic co-tunneling processes excite the dot and give rise to an additional current. Since these processes transfer energy between the dot and the leads, the co-tunneling is called inelastic.  The dot switches between the $`N`$ electron ground and excited state and in certain cases the $`N\pm 1`$ electron ground states can be occupied, too. These two regimes are considered separately in the next two sections. We point out that the model introduced here can be applied to a limited range of bias and gate voltages in an experimental situation where the capacitances and single particle energies provide a sufficiently accurate parametrization of the many-body spectrum of the dot.

## III Elastic co-tunneling

At low bias voltage $`V<\mathrm{\Delta }`$, a tunneling process of second order in the transition amplitude cannot populate the excited state. However, in this order of tunneling, transport is possible via elastic co-tunneling without exciting the dot. This is the dominant contribution to the current at sufficiently low temperatures and voltages. We do not consider here the integer-spin Kondo physics on a small energy scale $`T_K`$, which has recently been observed:  we assume $`T,VT_K`$. It was pointed out in Refs.  that elastic co-tunneling is the only process limiting the operation of a quantum dot as a single electron transistor when the level quantization $`\delta `$ is comparable with the charging energy $`U`$.

In a co-tunneling process two tunneling events occur in a very short time interval. For instance, an electron is first injected into the dot from the left lead and then emitted to the right lead. On this small time scale the two tunneling events are coherent. The energy of the total system (dot and leads) is sufficiently uncertain for the process not to be completely blocked. However, the transition amplitude for this process is inversely proportional to the energy difference between the initial and intermediate state of the total system.  The amplitude for transferring an electron in the opposite sequence (first emit an electron, then inject one) via another intermediate state must also be added if the final state of the system is the same. From the total second order amplitude (or matrix element) one obtains the tunneling rate by statistically averaging over all initial and final states of the leads.  Thus a current can pass through the dot although the $`N1`$ and $`N+1`$ electron states are only virtually accessible.

Let us describe the calculation of the tunneling rate for this simple case in order to skip the details of similar calculations below. Consider the initial state of the system where two electrons are in level $`l`$ of the left leads and the dot in the ground state, denoted by $`|N_{0,0}a_l^{}a_l^{}|`$ where $`|`$ denotes the state of the other electrons in the noninteracting leads. Applying the tunneling operators $`H_{TLd}^{}`$ and $`H_{RLd}`$ in the two possible sequences, one obtains the final state $`|N_{0,0}\widehat{F}|`$ where the operator $`\widehat{F}`$ is given in Table I. The amplitude for this process is written as $`\sqrt{2}M_{NN}`$ where the factor $`\sqrt{2}`$ is due to the spin degeneracy of the levels in the left lead. The explicit expression for $`M_{NN}`$ is given in Table II. Introducing the tunneling rate

$$\mathrm{\Gamma }_{NN}=2\pi \underset{l,r}{}\left|M_{NN}\right|^2f\left(\epsilon _l\mu _L\right)\left[1f\left(\epsilon _r\mu _R\right)\right]\delta \left(\epsilon _l\epsilon _r\right),$$
(24)
we write for the elastic co-tunneling current

$$I=2\mathrm{\Gamma }_{NN}.$$
(25)
The sum over all possible states of the leads with the Fermi distribution $`f\left(\epsilon \right)=1/\left(1+e^{\epsilon /T}\right)`$ can be converted into an integration and gives

$`\mathrm{\Gamma }_{NN}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}[\mathrm{\Gamma }_{dd}^L\mathrm{\Gamma }_{dd}^R\tau (E_{NN1},V)+\mathrm{\Gamma }_{d^{}d^{}}^L\mathrm{\Gamma }_{d^{}d^{}}^R\tau (E_{N+1N},V)`$ (27)
$`+\mathrm{\Gamma }_{dd^{}}^L\mathrm{\Gamma }_{d^{}d}^R\tau ^{}(E_{NN1},E_{N+1N},V)]`$
Here we introduced the rates $`\mathrm{\Gamma }_{dd}^L=2\pi _l\left|t_{ld}\right|^2\delta \left(\epsilon \epsilon _l\right)`$ and $`\mathrm{\Gamma }_{d^{}d^{}}^L`$ ($`dd^{}`$) for sequential tunneling into level $`d`$ and $`d^{}`$ and we assumed that these depend only weakly on the energy $`\epsilon `$. The same applies to the quantity $`\mathrm{\Gamma }_{dd^{}}^L=2\pi _lt_{ld}^{}t_{ld^{}}\delta \left(\epsilon \epsilon _l\right)`$ in the last term on the right-hand side of Eq. (27), which takes into account the interference between the two possible sequences for the tunneling events (“paths” in energy space); $`\mathrm{\Gamma }_{dd^{}}^R,\mathrm{\Gamma }_{dd}^R`$,$`\mathrm{\Gamma }_{d^{}d^{}}^R`$ are defined similarly ($`L,lR,r`$). All these quantities are treated as parameters here. In the expression for the rate (27) the explicit ($`V`$) and implicit dependence \[through the addition energies; cf. Eqs. (13)\] on the voltages is incorporated in the functions (explicitly evaluated in the limit $`T0`$)

$`\tau (E,V)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\epsilon {\displaystyle \frac{f\left(\epsilon \mu _L\right)\left[1f\left(\epsilon \mu _R\right)\right]}{\left(\epsilon E\right)^2}}`$ (29)
$`=`$ $`{\displaystyle \frac{V}{\left(\mu _LE\right)\left(\mu _RE\right)}}\mathrm{\Theta }\left(V\right),`$ (30)
$`\tau ^{}(E,E^{},V)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\epsilon {\displaystyle \frac{f\left(\epsilon \mu _L\right)\left[1f\left(\epsilon \mu _R\right)\right]}{\left|\epsilon E\right|\left|\epsilon E^{}\right|}}`$ (31)
$`=`$ $`{\displaystyle \frac{1}{\left|E^{}E\right|}}\left(\left|\mathrm{ln}{\displaystyle \frac{\mu _LE}{\mu _RE}}\right|+\left|\mathrm{ln}{\displaystyle \frac{\mu _LE^{}}{\mu _RE^{}}}\right|\right)\mathrm{\Theta }\left(V\right),`$ (32)
Here $`\mathrm{\Theta }`$ is the unit step function and $`lim_{E^{}E}\tau ^{}(E,E^{},V)=\tau (E,V)`$. The voltage dependence will be important below. The rate (27) diverges near the edges of the Coulomb diamond where sequential tunneling becomes dominant \[cf. condition (II)\]. Here the second order perturbation theory for the amplitudes breaks down. However, well within the Coulomb diamond we have $`\mathrm{\Gamma }_{NN}\mathrm{\Gamma }_{dd^{}}^{L,R},\mathrm{\Gamma }_{dd}^{L,R},\mathrm{\Gamma }_{d^{}d^{}}^{L,R}`$ and perturbation theory can be used. This is the regime of interest here.

## IV Inelastic co-tunneling

If the bias voltage can supply the necessary excitation energy, $`V\mathrm{\Delta }`$, then the excited state can be populated by inelastic co-tunneling events.  In such a coherent process, an electron is transferred from the left to the right lead, thereby loosing an amount $`\mathrm{\Delta }`$ of its excess energy, which is supplied by the transport voltage; this energy is transferred to the dot, which is thereby excited. The excited state can relax to the $`N`$ electron ground state by two sorts of processes depending on $`\mathrm{\Delta }`$. If $`\mathrm{\Delta }`$ is sufficiently small then there is an energy barrier for injecting or emitting electrons \[i.e. conditions (18) and (19) are not satisfied\]. Thus the excited state is also Coulomb blockaded. The dot can relax by three different inelastic co-tunneling processes. An electron can be transferred from the left to the right lead again, but now gain an energy $`\mathrm{\Delta }`$, thereby relaxing the dot. Similarly, an electron injected from either lead can be excited above the Fermi level of that same lead. The dot relaxes without a net charge transfer. This is the situation in the region labeled (I) in Fig. 1. For the current the elastic co-tunneling through the excited state is also important. If $`\mathrm{\Delta }`$ is sufficiently large, however, the excited state can decay to the ground state by two sequential tunneling processes. There are three possibilities. Either the dot can relax to the $`N1`$ electron ground state and then return to the $`N`$ electron ground state \[region (II) in Fig. 1\]; or the dot can be excited to the $`N+1`$ electron ground state from which it returns to the $`N`$ electron ground state \[region (II’)\]; or both processes can occur \[region (III)\]. For positive bias voltage these regions are delimited by the bias voltage threshold lines \[Eqs. (II) and (20)\] as follows:

(I) $`\mathrm{\Delta }`$ $`<`$ $`V`$ $`<`$ $`V_L^{N^{}N1},V_L^{N+1N^{}}`$ $`\mathrm{\Delta }<\frac{1}{3}V_{\text{max}}^N`$  (II) $`\mathrm{\Delta },V_R^{N^{}N1}`$ $`<`$ $`V`$ $`<`$ $`V_R^{NN1},V_L^{N+1N^{}}`$ $`\mathrm{\Delta }<\frac{1}{2}V_{\text{max}}^N`$  (II’) $`\mathrm{\Delta },V_L^{N+1N^{}}`$ $`<`$ $`V`$ $`<`$ $`V_L^{N+1N},V_L^{N^{}N1}`$ $`\mathrm{\Delta }<\frac{1}{2}V_{\text{max}}^N`$  (III) $`\mathrm{\Delta },V_R^{N^{}N1},V_L^{N+1N^{}}`$ $`<`$ $`V`$ $`<`$ $`V_R^{NN1},V_L^{N+1N}`$ $`\mathrm{\Delta }<V_{\text{max}}^N`$

Depending on the excitation energy $`\mathrm{\Delta }=\delta J`$ relative to the maximum bias voltage $`V_{\text{max}}^N=\delta +U`$ some of these regions do not exist. The different limiting cases are depicted in Fig. 2. When the excitation energy is sufficiently small, $`\mathrm{\Delta }<V_{\text{max}}^N/3`$, all these regions exist. From this we find that the Coulomb and/or exchange interaction should large relative to the level spacing: $`\delta <\frac{1}{2}U+\frac{3}{2}J`$. Region (I) ceases to exist for larger values of $`\mathrm{\Delta }`$ (or $`\delta `$) while regions (II) and (II’) continue to exists as long as $`\mathrm{\Delta }<V_{\text{max}}^N/2`$, i.e., $`\delta <U+2J`$. Region (III) always exists since $`\mathrm{\Delta }<V_{\text{max}}^N`$ is satisfied for $`U+J>0`$. The situation is clearly more complicated than at low voltage $`V<\mathrm{\Delta }`$. The transport in each region needs to be considered separately. Below we first consider the simplest region (I), where only the ground and excited state participate in the transport. Then we extend this analysis to the other regions.

### A Region I - relaxation by co-tunneling

This region exists only when the excited state lies sufficiently close to the ground state. There are four states involved in the transport, of which three are degenerate. In order to calculate the current we need the nonequilibrium occupation of the states in the dot. We assume that at the bias voltage where inelastic co-tunneling processes become important, $`V\mathrm{\Delta }`$, the lifetime of the states in the dots is much smaller than the typical time scale for building up Kondo-type correlations between the leads and dots system. In terms of energies we thus assume $`T_K\mathrm{\Delta }<V\mathrm{\Gamma }`$ where $`T_K`$ is the energy scale for such correlations. In this limit we can use the following master equations  to describe the occupation of the ground state $`\rho _{N_{0,0}}`$ and of the three triplet excited states $`\rho _{N_{1,m}^{}},m=0,\pm 1`$:

$`_t\rho _{N_{0,0}}`$ $`=`$ $`{\displaystyle \underset{i=1,2,3}{}}\mathrm{\Gamma }_{NN^{}}^{\left(i\right)}{\displaystyle \underset{m=0,\pm 1}{}}\rho _{N_{1,m}^{}}3\mathrm{\Gamma }_{N^{}N}\rho _{N_{0,0}},`$ (34)
$`_t\rho _{N_{1,0}^{}}`$ $`=`$ $`\mathrm{\Gamma }_{N^{}N}\rho _{N_{0,0}}+{\displaystyle \frac{1}{2}}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}{\displaystyle \underset{m=\pm 1}{}}\rho _{N_{1,m}^{}}\left[\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}+{\displaystyle \underset{i=1,2,3}{}}\mathrm{\Gamma }_{NN^{}}^{\left(i\right)}\right]\rho _{N_{1,0}^{}},`$ (35)
$`_t\rho _{N_{1,m}^{}}`$ $`=`$ $`\mathrm{\Gamma }_{N^{}N}\rho _{N_{0,0}}+{\displaystyle \frac{1}{2}}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}\rho _{N_{1,0}^{}}\left[{\displaystyle \frac{1}{2}}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}+{\displaystyle \underset{i=1,2,3}{}}\mathrm{\Gamma }_{NN^{}}^{\left(i\right)}\right]\rho _{N_{1,m}^{}},`$ (37)
where $`m=\pm 1`$ and $`\rho _{N_{0,0}}+_{m=0,\pm 1}\rho _{N_{1,m}^{}}=1`$. The current reads

$$I=\left(2\mathrm{\Gamma }_{NN}+3\mathrm{\Gamma }_{N^{}N}\right)\rho _{N_{0,0}}+\left(\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}+\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}+\mathrm{\Gamma }_{NN^{}}^{\left(1\right)}\right)\underset{m=0,\pm 1}{}\rho _{N_{1,m}^{}}.$$
(39)
The transitions between the states are schematically indicated in Fig. 3. The rates are calculated using the Golden Rule as was done for Eq. (27). In Table I we have listed the final state of the leads, the matrix element and the transition rate for each transition. The explicit expressions for the matrix elements are given in Table II. Let us now calculate these rates and discuss the transitions that they describe starting with the inelastic co-tunneling processes.

Each excited state of the triplet is populated by inelastic co-tunneling transition from the ground state with the same rate $`\mathrm{\Gamma }_{N^{}N}`$. An electron with energy $`\epsilon _l`$ and spin $`s=\pm 1/2`$ is first injected into the dot from lead $`L`$, bringing the dot in the intermediate state $`|N+1_{1/2,s}`$, and is subsequently emitted into the lead $`R`$ at an energy $`\epsilon _r=\epsilon _l\mathrm{\Delta }`$. The other sequence is also possible, in which case the intermediate state of the dot is $`|N1_{1/2,s}`$. The sum of the amplitudes for these processes yields the matrix element $`M_{N^{}N}`$, listed in Table II. The inelastic co-tunneling rate for exciting the dot is then reads \[cf. Eqs. (29) and (31)\]

$`\mathrm{\Gamma }_{N^{}N}`$ $`=`$ $`2\pi {\displaystyle \underset{l,r}{}}\left|M_{N^{}N}\right|^2f\left(\epsilon _l\mu _L\right)\left[1f\left(\epsilon _r\mu _R\right)\right]\delta \left(\epsilon _l\epsilon _r\mathrm{\Delta }\right)`$ (41)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\mathrm{\Gamma }_{d^{}d^{}}^L\mathrm{\Gamma }_{dd}^R[\tau (E_{NN1},V\mathrm{\Delta })+\tau (E_{N+1N}\mathrm{\Delta },V\mathrm{\Delta })`$ (43)
$`+2\tau ^{}(E_{NN1},E_{N+1N}\mathrm{\Delta },V\mathrm{\Delta })]`$
Conversely, each state of the triplet can relax to the singlet ground state by inelastic co-tunneling. The rate for these processes must be calculated separately because they result in a different final state of the leads. The matrix element $`M_{NN^{}}^{\left(1\right)}`$ in Table II describes the transfer of an electron from the left to the right lead starting in any of the three triplet sublevels of the excited state. The energy of the emitted electron is increased by the excitation energy $`\epsilon _r=\epsilon _l+\mathrm{\Delta }`$. The inelastic co-tunneling rate for this relaxation process is \[cf. Eq. (20)\]

$`\mathrm{\Gamma }_{NN^{}}^{\left(1\right)}`$ $`=`$ $`2\pi {\displaystyle \underset{l,r}{}}\left|M_{NN^{}}^{\left(1\right)}\right|^2f\left(\epsilon _l\mu _L\right)\left[1f\left(\epsilon _r\mu _R\right)\right]\delta \left(\epsilon _l+\mathrm{\Delta }\epsilon _r\right)`$ (45)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\mathrm{\Gamma }_{d^{}d^{}}^L\mathrm{\Gamma }_{dd}^R[\tau (E_{N^{}N1},V+\mathrm{\Delta })+\tau (E_{N+1N^{}}+\mathrm{\Delta },V+\mathrm{\Delta })`$ (47)
$`+2\tau ^{}(E_{N^{}N1},E_{N+1N^{}}+\mathrm{\Delta },V+\mathrm{\Delta })].`$
This rate explicitly depends on the voltage because the voltage window (both $`\mu _L`$ and $`\mu _R`$ appear in the Fermi functions). The second process involves the injection of an electron at energy $`\epsilon _l`$ from the left lead and the emission of an electron to an empty state in the same lead at a higher energy $`\epsilon _l^{}=\epsilon _l+\mathrm{\Delta }.`$ The matrix element $`M_{NN^{}}^{\left(2\right)}`$ is thus found by replacing $`rl^{}`$, cf. Table II. The tunneling rate is similarly found by replacing $`r,Rl^{},L`$ in Eq. (45) and from the result (47) we obtain:

$`\mathrm{\Gamma }_{NN^{}}^{\left(2\right)}={\displaystyle \frac{\mathrm{\Gamma }_{dd}^L}{\mathrm{\Gamma }_{dd}^R}}\mathrm{\Gamma }_{NN^{}}^{\left(1\right)}\left(V=0\right).`$

The third process involves the injection and emission of an electron to and from the right lead: replacing $`LR`$ we find:

$`\mathrm{\Gamma }_{NN^{}}^{\left(3\right)}={\displaystyle \frac{\mathrm{\Gamma }_{d^{}d^{}}^R}{\mathrm{\Gamma }_{d^{}d^{}}^L}}\mathrm{\Gamma }_{NN^{}}^{\left(1\right)}\left(V=0\right).`$

Let us now consider the elastic co-tunneling processes in which electrons are transferred from the left to the right lead. These involve four spin-degenerate intermediate states of the dot. In addition to the four ground states (II), there are four other intermediate excited states

$`|\left(N1\right)_{1/2,+1/2}^{}`$ $`=`$ $`a_d^{}^{}|0,|\left(N+1\right)_{1/2,+1/2}^{}=a_d^{}^{}a_d^{}^{}a_d^{}|0,`$
$`|\left(N1\right)_{1/2,1/2}^{}`$ $`=`$ $`a_d^{}^{}|0,|\left(N+1\right)_{1/2,1/2}^{}=a_d^{}^{}a_d^{}^{}a_d^{}|0.`$
The addition energies for the transitions from the $`N`$ electron excited state to these states are

$`E_{N^{}N1^{}}`$ $`=`$ $`E_{NN1}+\mathrm{\Delta }\delta =E_{NN1}+J,`$
$`E_{N+1^{}N^{}}`$ $`=`$ $`E_{N+1N}+\delta \mathrm{\Delta }=E_{N+1N}J.`$
These energies lie closer to the Fermi energies than $`E_{NN1}`$ and $`E_{N+1N}`$. Therefore the contributions of these intermediate states to elastic co-tunneling amplitudes are larger than those in the amplitude $`M_{NN}`$ and they cannot be disregarded. Three types of elastic co-tunneling processes must be considered separately. Two of these processes do not involve a flip of the spin of the transmitted electron. The two different final states of the leads as indicated in Table I. The first type of process occurs with the same rate $`\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}`$ for all states of the triplet. The second type is only possible for the two states with nonzero spin projection and occurs with a rate $`\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}`$. The rates are calculated from the matrix elements $`M_{N^{}N^{}}^{(1,2)}`$ in Table II and give \[compare this with Eq. (27)\]

$`\mathrm{\Gamma }_{N^{}N^{}}^{(1,2)}`$ $`=`$ $`2\pi {\displaystyle \underset{l,r}{}}\left|M_{N^{}N^{}}^{(1,2)}\right|^2f\left(\epsilon _l\mu _L\right)\left[1f\left(\epsilon _r\mu _R\right)\right]\delta \left(\epsilon _l\epsilon _r\right)`$ (51)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\{\mathrm{\Gamma }_{dd}^L\mathrm{\Gamma }_{dd}^R[\tau (E_{N^{}N1},V)+\tau (E_{N+1N^{}},V)\pm 2\tau ^{}(E_{N^{}N1},E_{N+1N^{}},V)]`$ (55)
$`+\mathrm{\Gamma }_{d^{}d^{}}^L\mathrm{\Gamma }_{d^{}d^{}}^R\left[\tau (E_{N^{}N1^{}},V)+\tau (E_{N+1^{}N^{}},V)\pm 2\tau ^{}(E_{N^{}N1^{}},E_{N+1^{}N^{}},V)\right]`$
$`+2Re\mathrm{\Gamma }_{dd^{}}^L\mathrm{\Gamma }_{d^{}d}^R[\tau ^{}(E_{N^{}N1},E_{N^{}N1^{}},V)\pm \tau ^{}(E_{N^{}N1},E_{N+1^{}N^{}},V)`$
$`+\tau ^{}(E_{N+1N^{}},E_{N+1^{}N^{}},V)\pm \tau ^{}(E_{N^{}N1^{}},E_{N+1N^{}},V)]\}.`$
These processes do not induce transitions between the states \[they don not contribute to Eqs. (IV A)\], but they do contribute to the transfer of electrons \[the rate appears in Eq. (39)\]. The third type of elastic co-tunneling process, involves a flip of the spin of the transmitted electron and induces transitions between the triplet states states. The triplet state $`|N_{1,0}^{}`$ can decay to both $`|N_{1,1}^{}`$ and $`|N_{1,1}^{}`$ with a rate equal to $`\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}`$. The latter two states can return to the former state with the same transition rate.

Because there is no Zeeman energy or spin-orbit coupling in our model, the transition rates between the states of the triplet balance each other, cf. Fig. 3. Therefore we can introduce the total occupations of the excited state $`\rho _N^{}_{m=0,\pm 1}\rho _{N_{1,m}^{}}`$summed over the spin projection and omit the spin indices $`\rho _N\rho _{N_{0,0}}`$. By adding the equations for the triplet states we obtain a closed set of master equations for the total occupations:

$`_t\rho _N`$ $`=`$ $`{\displaystyle \underset{i=1,2,3}{}}\mathrm{\Gamma }_{NN^{}}^{\left(i\right)}\rho _N^{}3\mathrm{\Gamma }_{N^{}N}\rho _N,`$ (57)
$`_t\rho _N^{}`$ $`=`$ $`3\mathrm{\Gamma }_{N^{}N}\rho _N{\displaystyle \underset{i=1,2,3}{}}\mathrm{\Gamma }_{NN^{}}^{\left(i\right)}\rho _N^{},`$ (58)
with $`\rho _N+\rho _N^{}=1`$. The elastic co-tunneling transitions between the triplet states have precisely cancelled. However, these processes do contribute to the current (39). The fact that the triplet state with zero spin projection transfers twice as much charge per unit time ($`2\times \frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}`$) as the other states ($`\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}`$) is precisely balanced by the spin-flip processes that can only occur for states with nonzero spin projection ($`\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}`$). The elastic co-tunneling contribution of each state of the triplet is thus proportional to the same rate (39) and we can also express the current in total occupations

$$I=\left(2\mathrm{\Gamma }_{NN}+3\mathrm{\Gamma }_{N^{}N}\right)\rho _N+\left(\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}+\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}+\mathrm{\Gamma }_{NN^{}}^{\left(1\right)}\right)\rho _N^{}.$$
(59)
Note that the inelastic co-tunneling relaxation processes with rates $`\mathrm{\Gamma }_{NN^{}}^{\left(2.3\right)}`$ do not appear here since they do not involve charge transfer. In the stationary limit $`lim_t\mathrm{}_t\rho _N^{},_t\rho _N=0`$ equations (IV A) are readily solved. The current can be written as

$$I=\frac{\left(\frac{1}{2}\mathrm{\Gamma }_{N^{}N^{}}^{\left(1\right)}+\mathrm{\Gamma }_{N^{}N^{}}^{\left(2\right)}2\mathrm{\Gamma }_{NN}\right)+2\mathrm{\Gamma }_{NN^{}}^{\left(1\right)}+_{i=2,3}\mathrm{\Gamma }_{NN^{}}^{\left(i\right)}}{_{i=1,2,3}\mathrm{\Gamma }_{N^{}N}^{\left(i\right)}+3\mathrm{\Gamma }_{N^{}N}}3\mathrm{\Gamma }_{N^{}N}+2\mathrm{\Gamma }_{NN}.$$
(60)
For $`V<\mathrm{\Delta }`$ the excitation rate $`\mathrm{\Gamma }_{N^{}N}`$ vanishes and we recover the elastic co-tunneling current 25. In Fig. 4 we have plotted the current, differential conductance and the co-tunneling rates as a function of bias voltages for the case $`\mathrm{\Delta }<V_{\text{max}}^N/3`$. In this limit region (I) covers most of the Coulomb diamond for $`V\mathrm{\Delta }`$. Apart from a clear kink in the current and a step in the differential conductance, a sharp peak is present at the step edge. This is clearly due to a change of the population of the excited state.

### B Regions (II), (II’) and (III) - relaxation by sequential tunneling

The remain three regions where sequential tunneling transitions from the excited state to $`N1`$ and/or $`N+1`$ electron ground state are possible are now easily included. The elastic co-tunneling processes involving the excited state must now be disregarded ( and their contributions to the current). The rates for these processes diverge at the boundaries of regime (I) signalling that perturbation theory breaks down. Instead we write a full set of master equations similar to (IV A) without the transitions between the triplet states and now including the occupations of the spin sublevels of the $`N1`$ and/or $`N+1`$ electron ground states. Proceeding as before we arrive at the following closed set of master equations for the total occupations, including $`\rho _{N\pm 1}_{m=\pm 1/2}\rho _{\left(N\pm 1\right)_{1/2,m}}`$:

$`_t\rho _{N1}`$ $`=`$ $`\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^R\rho _N^{}\left(\alpha _{}\mathrm{\Gamma }_{dd}^R+\alpha _+\mathrm{\Gamma }_{d^{}d^{}}^L\right)\rho _{N1},`$ (62)
$`_t\rho _{N+1}`$ $`=`$ $`\alpha _+\mathrm{\Gamma }_{dd}^L\rho _N^{}\left(\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^L+\alpha _+\mathrm{\Gamma }_{d^{}d^{}}^R\right)\rho _{N+1},`$ (63)
$`_t\rho _N`$ $`=`$ $`+\alpha _{}\mathrm{\Gamma }_{dd}^L\rho _{N1}+\alpha _+\mathrm{\Gamma }_{d^{}d^{}}^R\rho _{N+1}3\mathrm{\Gamma }_{N^{}N}\rho _N,`$ (64)
$`_t\rho _N^{}`$ $`=`$ $`\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^L\rho _{N1}+\alpha _+\mathrm{\Gamma }_{dd}^R\rho _{N+1}+3\mathrm{\Gamma }_{N^{}N}\rho _N`$ (66)
$`\left(+\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^R\rho _N^{}+\alpha _+\mathrm{\Gamma }_{dd}^L\rho _N^{}\right)\rho _N^{},`$
The current reads

$$I=\left(2\mathrm{\Gamma }_{NN}+3\mathrm{\Gamma }_{N^{}N}\right)\rho _N+\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^R\rho _N^{}+\alpha _+\left(\mathrm{\Gamma }_{dd}^R+\mathrm{\Gamma }_{d^{}d^{}}^R\right)\rho _{N+1}.$$
(67)
These expression covers regions (II) \[$`\alpha _{}=1,\alpha _+=0`$\], (II’) \[$`\alpha _{}=0,\alpha _+=1`$\] and (III) \[$`\alpha _\pm =1`$\]. The solution for all cases can be written in the form

$$I=\frac{2K+K^{}2\kappa \mathrm{\Gamma }_{NN}+2\mathrm{\Gamma }_{N^{}N}^{\left(1\right)}}{K+3\kappa \mathrm{\Gamma }_{N^{}N}}3\mathrm{\Gamma }_{N^{}N}+2\mathrm{\Gamma }_{NN}.$$
(68)
The coefficients appearing here are constants that characterize the differences between the three regions:

$`K`$ $`=`$ $`\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^R{\displaystyle \frac{\mathrm{\Gamma }_{dd}^L}{\mathrm{\Gamma }_{dd}^L+\mathrm{\Gamma }_{d^{}d^{}}^L}}+\alpha _+\mathrm{\Gamma }_{d^{}d^{}}^L{\displaystyle \frac{\mathrm{\Gamma }_{dd}^R}{\mathrm{\Gamma }_{dd}^R+\mathrm{\Gamma }_{d^{}d^{}}^R}},`$
$`K^{}`$ $`=`$ $`\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^R{\displaystyle \frac{\mathrm{\Gamma }_{d^{}d^{}}^L}{\mathrm{\Gamma }_{dd}^L+\mathrm{\Gamma }_{d^{}d^{}}^L}}+\alpha _+\mathrm{\Gamma }_{d^{}d^{}}^L{\displaystyle \frac{\mathrm{\Gamma }_{d^{}d^{}}^R}{\mathrm{\Gamma }_{dd}^R+\mathrm{\Gamma }_{d^{}d^{}}^R}},`$
$`\kappa `$ $`=`$ $`\alpha _{}{\displaystyle \frac{\mathrm{\Gamma }_{d^{}d^{}}^R}{\mathrm{\Gamma }_{dd}^L+\mathrm{\Gamma }_{d^{}d^{}}^L}}+\alpha _+{\displaystyle \frac{\mathrm{\Gamma }_{d^{}d^{}}^L}{\mathrm{\Gamma }_{dd}^R+\mathrm{\Gamma }_{d^{}d^{}}^R}}+1.`$
where $`K+K^{}=\alpha _{}\mathrm{\Gamma }_{d^{}d^{}}^R+\alpha _+\mathrm{\Gamma }_{d^{}d^{}}^L`$. In Fig. 5 we have plotted the co-tunneling rates, the occupations of the states and the differential conductance as a function of the bias voltages for $`\mathrm{\Gamma }_{dd}^\alpha =\mathrm{\Gamma }_{d^{}d^{}}^\alpha =\mathrm{\Gamma },\alpha =L,R`$ for a typical case where $`V_{\text{max}}^N/2\mathrm{\Delta }`$. In this limit only region (III) exists. The qualitative result is the same as in region (I). The essential condition for the peak in the conductance is that the total rate for charge transport through the excited state is bigger than the rate $`\mathrm{\Gamma }_{N^{}N}`$ at which this state is populated.

## V Conclusions

We have analyzed the finite bias transport within a Coulomb blockaded region of a quantum dot with an even number of electrons. The differential conductance through the dot was found to exhibit a step within the Coulomb diamond at a bias voltage equal to the excitation energy. At the step edge an additional peak appears that is caused by a change in the nonequilibrium populations of the dot states due the voltage dependence of the co-tunneling rates. This may explain the recent observations in the transport spectroscopy of semi-conductor quantum dots by De Franceschi et. al and of carbon nanotubes by Nygård et. al. As pointed out in Ref.  the sharp feature in the differential conductance can be used to do a spectroscopy of excited states with increased resolution. We have considered two well-separated regimes where the excited state is depopulated by inelastic co-tunneling and sequential tunneling, respectively. It is expected that in the crossover regime the results do not change qualitatively since the inelastic co-tunneling processes that excite the dot form the bottleneck for the transport. The rate for these processes as calculated here is correct also in the crossover regime. The calculations presented here can be improved by including (additive) contributions of more virtual states of the dot (more than two orbitals) to the transition amplitudes, which determine the co-tunneling rates. Also, the situation where a second $`N`$-electron excited state lies the Coulomb diamond can be treated in the same way as single excited state. The occupations of all transport states can be calculated from an extended set of master equations. One must now also include co-tunneling transitions between the two excited states. The rates that enter in these equations have the same form as those calculated here. Furthermore, the peak at the onset of inelastic co-tunneling is broadened by thermal fluctuations and level broadening effects, as pointed out in Ref. . The effect of a electron finite temperature can be included in our calculations. Level broadening effects require a nonperturbative scheme, which is beyond the scope of this paper. It is clear, however, that near the onset of inelastic co-tunneling, the broadening of the peak, which is determined the decay rate of the excited state (either by sequential or inelastic co-tunneling), is much larger than the rate for exciting the dot.

## VI Acknowledgements

We would like to thank S. De Franceschi for bringing this problem to our attention and for useful discussions. One of us (M.R.W) would like to thank M. Eto and Y. Tokura for valuable discussions. This work is part of the research program of the “Stichting voor Fundamenteel Onderzoek der Materie” (FOM), which is financially supported by the ”Nederlandse Organisatie voor Wetenschappelijk Onderzoek” (NWO) and the NEDO project NTDP-98.


format: text

text:
# 1 Introduction

## 1 Introduction

The principles of equilibrium statistical mechanics are well established. But, thermal equilibrium is a special case, and little is known about the properties of systems not in equilibrium, for example about the relaxation toward the stationary state. Some interesting problems in non-equilibrium systems are non-equilibrium phase transitions described by phenomenological rate equations, and the way the system relaxes to its steady state. As mean-field techniques, generally, do not give correct results for low-dimensional systems, people are motivated to study exactly-solvable stochastic models in low dimensions. Moreover, solving one-dimensional systems should in principle be easier. Exact results for some models on a one-dimensional lattice have been obtained, for example in \[1–12\]. Different methods have been used to study these models, including analytical and asymptotic methods, mean field methods, and large-scale numerical methods. Systems with more than one species have also been studied \[13–27\]. Many of the arguments are based on simulation results. There are, however, some exact results as well. For most of the models studied, the interaction is between nearest neighbors.

In , a 10-parameter family of stochastic models with interactions between nearest neighbors has been studied. In these models, the $`k`$-point equal time correlation functions $`n_in_j\mathrm{}n_k`$ satisfy linear differential equations involving no higher-order correlations. We call these models autonomous, in the sense that the evolution equations of $`n`$-point functions are closed (contain only $`n`$\- or less-point functions). These linear equations for the average number of the particles $`n_i`$ have been solved. The same models were studied on lattices with boundaries in . It was shown that these models may exhibit dynamic and static phase transitions. The same idea has been generalized to multi-species models  in one dimension with two-site interactions. There, conditions have been obtained that the Hamiltonian should satisfy in order that the evolution equation for correlation functions be closed. The set of equations for average densities can be written in terms of four matrices. These matrices are not determined uniquely from the Hamiltonian: there is a kind of gauge transformation one can apply on them which of course, does not change the evolution equation. A formal solution for the average densities of different species was found. The large-time behaviour of the average densities of different species was also studied. The time evolution equations for more-point functions, generally contain not only these four matrices, but also elements of the Hamiltonian, and to obtain a closed form for their solution is generally not easy.

The Glauber dynamics was originally proposed to study the relaxation of the Ising model near equilibrium states. It was also shown that, there is a relation between the kinetic Ising model at zero temperature and the diffusion annihilation model in one dimension. There is an equivalence between domain walls in the Ising model and particles in the diffusion annihilation model. Kinetic generalizations of the Ising model, for example the Glauber model or the Kawasaki model, are phenomenological models and have been studied extensively \[30–41\].

In this paper, autonomous multispecies systems with more-than-two-neighbor interactions are studied. Conditions necessary and sufficient for closedness of the evolution equations of the $`n`$-point functions are obtained. As an example, we obtain explicitly the average number of the particles at each site for one species and three-site interactions. This is then generalized to the case where more than three site interact. As another example, a generalization of the Glauber model is presented. In this generalization, the processes are the same as those of the ordinary Glauber model, but the rates depend on three free parameters, rather than one free parameter in the ordinary Glauber model. Finally this model is further generalized to the case where the number of interacting sites is more than three and the number of states at each site is more than two.

## 2 Models leading to closed set of evolution equations

The models addressed are multi-species exclusion reaction-diffusion models. That is, each site is a vacancy or contains one particle. There are several kinds of particles, but at any time at most one kind can be present at each site. Throughout the paper, the system is assumed to be translationally invariant. Consider first a case where the interaction is between three neighboring sites. Then the Hamiltonian describing the system can be written as

$$=\underset{i=1}{\overset{L}{}}H_{i,i+1,i+2}.$$
(1)
The number of sites is $`L`$ and the number of possible states in a site is $`q`$ (one of these states, for example the $`q`$-th one, may be the vacancy); different states of each site are denoted by $`\alpha `$, $`\alpha =1,\mathrm{}q`$. Introducing $`n_i^\alpha `$ as the number operator of the particles of type $`\alpha `$ in the site $`i`$, we have

$$\underset{\alpha =1}{\overset{q}{}}n_i^\alpha =1.$$
(2)
The average number of the particles of the type $`\alpha `$ in the site $`i`$ at the time $`t`$ is

$$n_i^\alpha =S|n_i^\alpha |P(t),$$
(3)
where $`|P(t):=\mathrm{exp}(t)|P(0)`$ represents the state of the system at the time $`t`$,

$$S|=\underset{L}{\underset{}{s|\mathrm{}s|}},$$
(4)
and

$$s|:=\underset{q}{\underset{}{(\mathrm{1\; 1}\mathrm{}1)}}.$$
(5)
So, the time evolution of $`n_i^\alpha `$ is given by

$$\frac{\mathrm{d}}{\mathrm{d}t}n_i^\alpha =S|n_i^\alpha |P(t).$$
(6)
The only terms of the Hamiltonian $``$ entering the above equation are $`H_{i,i+1,i+1}`$, $`H_{i1,i,i+1}`$, and, $`H_{i1,i,i+1}`$. The result of acting any matrix $`Q`$ on the bra $`s|`$ is equal to that of acting the diagonal matrix $`\stackrel{~}{Q}`$ on the same bra, provided each diagonal element of the matrix $`\stackrel{~}{Q}`$ is equal to the sum of all elements of the corresponding column in the matrix $`Q`$. So, the actions of $`(11n^\alpha )H`$, $`(1n^\alpha 1)H`$, and $`(n^\alpha 11)H`$ on $`s|s|s|`$ are equal to the actions of three diagonal matrices on $`s|s|s|`$. We use the symbol $``$ to denote the equality of the action on $`s|s|s|`$. We have

$`(n^\alpha 11)H`$ $`{\displaystyle \underset{\beta \gamma \lambda }{}}{}_{}{}^{1}𝒜_{\beta \gamma \lambda }^{\alpha }n^\beta n^\gamma n^\lambda `$
$`(1n^\alpha 1)H`$ $`{\displaystyle \underset{\beta \gamma \lambda }{}}{}_{}{}^{2}𝒜_{\beta \gamma \lambda }^{\alpha }n^\beta n^\gamma n^\lambda `$
$`(11n^\alpha )H`$ $`{\displaystyle \underset{\beta \gamma \lambda }{}}{}_{}{}^{3}𝒜_{\beta \gamma \lambda }^{\alpha }n^\beta n^\gamma n^\lambda ,`$ (7)
where $`{}_{}{}^{i}𝒜_{\beta \gamma \lambda }^{\alpha }`$’s are defined as

$`{}_{}{}^{1}𝒜_{\beta \gamma \lambda }^{\alpha }`$ $`:=s_\tau s_\omega H_{\beta \gamma \lambda }^{\alpha \tau \omega }`$
$`{}_{}{}^{2}𝒜_{\beta \gamma \lambda }^{\alpha }`$ $`:=s_\tau s_\omega H_{\beta \gamma \lambda }^{\tau \alpha \omega }`$
$`{}_{}{}^{3}𝒜_{\beta \gamma \lambda }^{\alpha }`$ $`:=s_\tau s_\omega H_{\beta \gamma \lambda }^{\tau \omega \alpha }.`$ (8)
An implicit summation (from 1 to $`q`$) over the a same subscript and superscript is always assumed. From these, (6) takes the form

$$\dot{n}_i^\alpha ={}_{}{}^{1}𝒜_{\beta \gamma \lambda }^{\alpha }n_i^\beta n_{i+1}^\gamma n_{i+2}^\lambda +{}_{}{}^{2}𝒜_{\beta \gamma \lambda }^{\alpha }n_{i1}^\beta n_i^\gamma n_{i+1}^\lambda +{}_{}{}^{3}𝒜_{\beta \gamma \lambda }^{\alpha }n_{i2}^\beta n_{i1}^\gamma n_i^\lambda $$
(9)
Generally, the right-hand side of (9) contains one-, two-, and three-point functions. (Note that $`n^\alpha `$’s are not independent.) We want to obtain a condition that only one point functions appear in the right-hand side. To do this, we consider the expression

$$u=f_{\alpha \beta \gamma }n_i^\alpha n_j^\beta n_k^\gamma ,$$
(10)
and ask for the condition that the right-hand side is expressible in terms of linear combinations of $`n`$’s, provided

$$s_\alpha n_l^\alpha =1.$$
(11)
It is obvious that if

$$f_{\alpha \beta \gamma }={}_{1}{}^{}F_{\alpha }^{}+{}_{2}{}^{}F_{\beta }^{}+{}_{3}{}^{}F_{\gamma }^{},$$
(12)
then the right hand side of (10) is expressible in terms of linear combinations of $`n`$’s. To prove that this form for $`f`$ is necessary as well, we just count the number of independent variables in $`f`$’s satisfying the desired property. One can write $`n_l^q`$ in terms one 1 and other $`n_l^\alpha `$’s. Then it is seen that a general cubic form of $`n`$’s is expressible in terms of $`q^3`$ independent forms of $`n`$’s, each containing no more than three $`n`$’s. Of these, $`1+3(q1)`$ expressions (the monomials of degree zero and one of $`3(q1)`$ independent variables) are desirable. The coefficients of other monomials should be zero. So, from $`q^3`$ independent variables in $`f`$, there remains only $`3(q1)+1`$ independent variables in $`f`$’s satisfying the desired condition. It seems that the right-hand side of (12) contains more independent variables, $`3q`$. But we note that the transformation

$${}_{i}{}^{}F_{\alpha }^{}{}_{i}{}^{}F_{\alpha }^{}+{}_{i}{}^{}B,$$
(13)
does not change the right-hand side of (12), provided

$$\underset{i}{}{}_{i}{}^{}B=0.$$
(14)
This means that there are $`31`$ redundant variables in the expression (12). So (12) contains actually the correct number of independent variables, and hence is the most general form of $`f`$ with the desired property.

So, in order that (9) be expressible in terms of only one-point functions, one must have

$${}_{}{}^{i}𝒜_{\beta \gamma \lambda }^{\alpha }={}_{1}{}^{i}𝒜_{\beta }^{\alpha }+{}_{2}{}^{i}𝒜_{\gamma }^{\alpha }+{}_{3}{}^{i}𝒜_{\lambda }^{\alpha }.$$
(15)
It is noted that $`{}_{j}{}^{i}𝒜`$’s are not determined uniquely. Applying the gauge transformation

$${}_{j}{}^{i}𝒜_{\beta }^{\alpha }{}_{j}{}^{i}𝒜_{\beta }^{\alpha }+{}_{j}{}^{i}B_{}^{\alpha }={}_{j}{}^{i}𝒜_{\beta }^{\alpha }+{}_{j}{}^{i}B_{}^{\alpha }s_\beta ,\text{with}\underset{j}{}{}_{j}{}^{i}B_{}^{\alpha }=0$$
(16)
does not change the right-hand side of (15).

If (15) is satisfied, (6) takes the form

$`\dot{n}_i^\alpha =`$ $`({}_{1}{}^{1}𝒜_{\beta }^{\alpha }+{}_{2}{}^{2}𝒜_{\beta }^{\alpha }+{}_{3}{}^{3}𝒜_{\beta }^{\alpha })n_i^\beta `$
$`+({}_{2}{}^{1}𝒜_{\beta }^{\alpha }+{}_{3}{}^{2}𝒜_{\beta }^{\alpha })n_{i+1}^\beta +({}_{1}{}^{2}𝒜_{\beta }^{\alpha }+{}_{2}{}^{3}𝒜_{\beta }^{\alpha })n_{i1}^\beta `$
$`+{}_{3}{}^{1}𝒜_{\beta }^{\alpha }n_{i+2}^\beta +{}_{1}{}^{3}𝒜_{\beta }^{\alpha }n_{i2}^\beta .`$ (17)
(15) in fact guarantees that the time-evolution equations of $`n`$-point functions contain only $`n`$\- and less-point functions. In the simplest case, the one-species, each site is vacant or occupied by only one kind of particles. Then, the matrices $`{}_{j}{}^{i}𝒜`$ are two by two.

One can do the same arguments for the case where more than three neighboring sites interact. Suppose the number of interacting sites is $`N`$. One defines

$${}_{}{}^{i}𝒜_{\beta _1\mathrm{}\beta _N}^{\alpha _i}:=\left(\underset{li}{}s_{\alpha _l}\right)H_{\beta _1\mathrm{}\beta _N}^{\alpha _1\mathrm{}\alpha _N}.$$
(18)
To ensure that in the time-evolution equation of one-point functions only one-point functions appear, one must have

$${}_{}{}^{i}𝒜_{\beta _1\mathrm{}\beta _N}^{\alpha _i}=\underset{j}{}{}_{j}{}^{i}𝒜_{\beta _j}^{\alpha _i}.$$
(19)
Here too, the gauge transformation (16) does not change the right-hand side of (19), and hence the physics of the problem. It is to be noted that (19) is also sufficient for $`n`$-point functions evolution equations to contain no more than $`n`$-point functions.

### 2.1 Some special cases

We now consider some special cases.

#### 2.1.1 single species case

In this case the matrices $`{}_{j}{}^{i}𝒜`$ are two by two. The time-evolution equation for $`n_k`$ will then be

$$\dot{n}_i=\alpha n_i+\beta n_{i+1}+\beta ^{}n_{i1}+\gamma n_{i+2}+\gamma ^{}n_{i2}+\delta .$$
(20)
$`\delta `$ can be eliminated using the redefinition

$$x_i:=n_i\frac{\delta }{\beta +\beta ^{}+\gamma +\gamma ^{}\alpha }.$$
(21)
Then, introducing the generating function

$$G(z,t)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}x_iz^i,$$
(22)
one arrives at,

$$\dot{G}(z,t)=\left(\alpha +\beta z^1+\beta ^{}z+\gamma z^2+\gamma ^{}z^2\right)G(z,t),$$
(23)
the solution to which is

$$G(z,t)=\mathrm{exp}\left[t\left(\alpha +\beta z^1+\beta ^{}z+\gamma z^2+\gamma ^{}z^2\right)\right]G(z,0).$$
(24)
Using

$$\mathrm{exp}\left[\frac{u}{2}\left(z+z^1\right)\right]=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{I}_k(u)z^k,$$
(25)
one arrives at

$$x_k(t)=e^{\alpha t}\underset{j,l=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{I}_{kj2l}\left(2t\sqrt{\beta \beta ^{}}\right)\mathrm{I}_l\left(2t\sqrt{\gamma \gamma ^{}}\right)(\beta ^{}/\beta )^{(kj2l)/2}(\gamma ^{}/\gamma )^{l/2}x_j(0).$$
(26)
A similar procedure can be done for more-than-three-neighboring-sites interactions. The main difference will be the number of modified Bessel functions appearing in the expression.

#### 2.1.2 Generalizations of the Glauber model

Consider a two-state three-neighbor interaction of the form

$${}_{}{}^{1}𝒜_{\beta \gamma \lambda }^{\alpha }={}_{}{}^{3}𝒜_{\beta \gamma \lambda }^{\alpha }=0.$$
(27)
This means that, similar to the Glauber model, any site interacts only with its neighboring sites. The interactions are

$`AAA`$ $`A\mathrm{}A`$ $`\mu _1`$
$`\mathrm{}\mathrm{}\mathrm{}`$ $`\mathrm{}A\mathrm{}`$ $`\mu _2`$
$`A\mathrm{}A`$ $`AAA`$ $`\lambda _1`$
$`\mathrm{}A\mathrm{}`$ $`\mathrm{}\mathrm{}\mathrm{}`$ $`\lambda _2`$
$`AA\mathrm{}`$ $`A\mathrm{}\mathrm{}`$ $`\alpha _1`$
$`\mathrm{}\mathrm{}A`$ $`\mathrm{}AA`$ $`\alpha _2`$
$`A\mathrm{}\mathrm{}`$ $`AA\mathrm{}`$ $`\beta _1`$
$`\mathrm{}AA`$ $`\mathrm{}\mathrm{}A`$ $`\beta _2.`$ (28)
This is a generalization of the Glauber model. For the ordinary Glauber model

$`\mu _1`$ $`=\mu _2=1\mathrm{tanh}{\displaystyle \frac{J}{k_\mathrm{B}T}}`$
$`\lambda _1`$ $`=\lambda _2=1+\mathrm{tanh}{\displaystyle \frac{J}{k_\mathrm{B}T}}`$
$`\alpha _1`$ $`=\alpha _2=\beta _1=\beta _2=1`$ (29)
The criterion (19) for the closedness of the evolution-equation for one-point functions results in the following relations between the rates in the generalized Glauber model.

$`\mu _i\alpha _i`$ $`=\beta _j\lambda _j\text{for any }i,j`$
$`\alpha _1+\beta _1`$ $`=\alpha _2+\beta _2`$ (30)
So there are four independent variables in terms of them the above eight parameters can be expressed. One can write the expressions as

$`\mu _1`$ $`=ABCC^{}`$
$`\mu _2`$ $`=A+B`$
$`\lambda _1`$ $`=A+B+C+C^{}`$
$`\lambda _2`$ $`=AB`$
$`\alpha _1`$ $`=ABC^{}`$
$`\alpha _2`$ $`=A+B+C`$
$`\beta _1`$ $`=A+B+C^{}`$
$`\beta _2`$ $`=ABC.`$ (31)
Of course one of the parameters can be absorbed through a time-rescaling. The resulting evolution-equation for the average particle-number is

$$\dot{n}_i=A+B2An_i+Cn_{i+1}+C^{}n_{i1},$$
(32)
which is easy to solve.

This generalized Glauber model can further be generalized in two directions: when the number of the interacting sites is more than 3, and when the number of states of each site is more than 2. The first case means that the interaction is in a block of length $`N`$, resulting in the change of the state of a single specific site in that block. This rate of change depends on the states of this site and the states of the other $`N1`$ sites. Let’s label this specific site of the block by $`0`$. (Usually the length of the block is considered to be an odd integer ($`2k+1`$), and the evolving site is assumed to be the central one.) Denote the state of the site $`i`$ by $`\sigma _i`$, where $`\sigma _i`$ can take the values $`1`$ (particle) or $`0`$ (vacancy). Then the evolution equation for the average particle number is

$`\dot{n}_0=`$ $`{\displaystyle \underset{\stackrel{}{\sigma }}{}}R(1,\stackrel{}{\sigma })n_0{\displaystyle \underset{i0}{}}\left[1n_i+\sigma _i(2n_i1)\right]`$
$`+{\displaystyle \underset{\stackrel{}{\sigma }}{}}R(0,\stackrel{}{\sigma })(1n_0){\displaystyle \underset{i0}{}}\left[1n_i+\sigma _i(2n_i1)\right].`$ (33)
Here the state of other interacting sites is denoted by $`\stackrel{}{\sigma }`$, $`R(\sigma _0,\stackrel{}{\sigma })`$ is the rate of change of the state of site $`0`$, from $`\sigma _0`$ to $`1\sigma _0`$, when the states of the other interacting sites is $`\stackrel{}{\sigma }`$.

We are looking for those rates $`R(\sigma _0,\stackrel{}{\sigma })`$ that make the right-hand side of this evolution equation a linear combination of $`n_i`$’s. The claim is that the general form of these rates is

$$R(\sigma _0,\stackrel{}{\sigma })=A+(1)^{\sigma _0}\left(B+\underset{i0}{}C_i\sigma _i\right).$$
(34)
Inserting this ansatz into the evolution equation (2.1.2) we get

$$\dot{n}_0=B+A12n_0+\underset{i0}{}C_in_i.$$
(35)
So it is clear that the ansatz (34) leads into a closed set of evolution-equations for the average particle number. It remains to prove that this ansatz is the most general one satisfying this property. To see this, one considers (2.1.2). In the right-hand side of this equation there are $`2^N`$ terms (the expectation of monomials in terms of $`n_i`$’s). Of these, we wish that the coefficients of all be zero, except for the coefficients of the constant term, and linear terms. So, there are $`2^N(N+1)`$ equations to be satisfied for the rates (consisting themselves of $`2^N`$ unknowns). This shows that the rates satisfying the desired property contain $`N+1`$ independent variables, and it is clear that the ansatz (34) contains $`N+1`$ independent variables. So it is the most general solution.

Now consider the second generalization, when the number of possible states at each site is more than 2, say $`q`$, and each block consists of $`N`$ sites. The state of the site $`i`$ is denoted by $`\sigma _i`$, which can take $`q`$ values. That site the state of which is denoted by $`0`$. The rate of change of the state of he site $`0`$ from $`\beta `$ to $`\alpha `$ is denoted by $`R_{\beta ,\stackrel{}{\sigma }}^\alpha `$. The evolution equation for the average numbers is then

$`\dot{n}_0^\alpha =`$ $`{\displaystyle \underset{\stackrel{}{\sigma }}{}}{\displaystyle \underset{\beta \alpha }{}}R_{\beta ,\stackrel{}{\sigma }}^\alpha n_0^\beta {\displaystyle \underset{i0}{}}n_i^{\sigma _i}`$
$`{\displaystyle \underset{\stackrel{}{\sigma }}{}}{\displaystyle \underset{\beta \alpha }{}}R_{\alpha ,\stackrel{}{\sigma }}^\beta n_0^\alpha {\displaystyle \underset{i0}{}}n_i^{\sigma _i}.`$ (36)
Defining

$$R_{\alpha ,\stackrel{}{\sigma }}^\alpha :=\underset{\beta \alpha }{}R_{\alpha ,\stackrel{}{\sigma }}^\beta ,$$
(37)
the evolution equation reads

$$\dot{n}_0^\alpha =\underset{\stackrel{}{\sigma }}{}R_{\beta ,\stackrel{}{\sigma }}^\alpha n_0^\beta \underset{i0}{}n_i^{\sigma _i}.$$
(38)
It is easy to see that for the right-hand side of this equation be expressible in terms of one-point functions, one should have

$$R_{\beta ,\stackrel{}{\sigma }}^\alpha ={}_{0}{}^{}R_{\beta }^{\alpha }+\underset{i0}{}{}_{i}{}^{}R_{\sigma _i}^{\alpha }$$
(39)
Note that this is nothing but

$${}_{}{}^{0}𝒜_{\beta ,\stackrel{}{\sigma }}^{\alpha }=\underset{i}{}{}_{i}{}^{0}𝒜_{\sigma _i}^{\alpha },$$
(40)
obtained in (19). In this case, other $`{}_{}{}^{i}𝒜`$’s vanish. The evolution equation of the one-point functions is then

$$\dot{n}_0^\alpha =\underset{i}{}{}_{i}{}^{}R_{\beta }^{\alpha }n_i^\beta .$$
(41)
Previous discussions shows that this is the most general form of the rates for them the evolution-equations of the one-point functions contain only one-point functions. Here too, there is a gauge freedom in choosing $`{}_{i}{}^{}R`$’s, namely

$${}_{i}{}^{}R_{\beta }^{\alpha }{}_{i}{}^{}R_{\beta }^{\alpha }+{}_{i}{}^{}B_{}^{\alpha }s_\beta ,\text{with}\underset{i}{}{}_{i}{}^{}B_{}^{\alpha }=0.$$
(42)
This is the same gauge freedom encountered with earlier, (16).


text:
# kinetic equation model for self-gravitating systems

## 1  Introduction

There is a long tradition of applications the kinetic theory to the gas of self-gravitating particles -. The first regular approach of that kind was developed by Chandrasekhar , and appeared to be very useful in astrophysical applications. Indeed, the kinetic theory is certainly applicable for globular star clusters, which are rather well-isolated systems on the relevant time-scales, and contain sufficiently many (typically of order $`N10^5`$) stars . Therefore, the kinetic approach, which typically studies a closed, macroscopic system, has a definite range of validity here. It allows to identify correctly the collisional time-scale -, and predict the moderate and late stages of the evolution for globular clusters .

The general method for constructing the collisional Fokker-Planck equation for statistical systems with $`1/r`$ interaction was proposed in  (see also the text-book descriptions in ). This equation adopts a more general strategy for Boltzmann’s kinetic equation to the considered long-range interacting system. After suitably dealing with the peculiarities of the long-range interaction, one ends up with the collisional Fokker-Planck equation. The main difficulty with this equation is that it is technically rather involved and resistant to reliable analytical treatments. Thus, one has to change the level of description and move to hydrodynamical models of the stellar dynamics (see e.g.,  and refs. therein). Instead of the one-particle distribution function, which is the basic object of the kinetic theory, the hydrodynamical approach operates with a few coarse-grained hydrodynamic variables, but by its meaning it is applicable only to the dense central regions of the self-gravitating system, where collisions are certainly dominating.

In this contribution we present another approach to the kinetic theory of self-gravitating systems. This is an approach of kinetic models, where instead of going immediately to the hydrodynamical level of description, one is properly modelling the collisional part of the kinetic equation, trying to preserve its main qualitative characteristics, but to make the system analytically tractable. One of the most successful kinetic models was proposed by Bhatnager, Gross and Krook (BGK) , and since that time appeared to be very useful in the kinetic theory. The reason of its success is that the model offers a minimal way to incorporate the necessary conservation laws (those of particle number, momentum and energy) with a simple relaxation mechanism. For a self-gravitating statistical system this approach provides natural methods to consider simultaneously the collision-dominated (overdamped) regions of the system, and its almost collisionless outer domains (underdamping).

In the present contribution we will restrict ourselves to the presentation of the BGK-method as applied to self-gravitating systems (section II), and to straightforward analytical investigations which will recover within a single setup several known and also new results for the steady regime of behavior (section III). More detailed investigation of the situation, as well as the presentation of more elaborated BGK-type models, are planned for future.

## 2 Bhatnager-Gross-Krook kinetic equation

The statistical dynamics of $`N`$ point particles with unit mass interacting through Newton’s inverse-square law, will be considered in the classical kinetic approach. This basically involves two important assumptions: (i) The field acting on a given particle (test particle) can be represented as a mean, self-consistent field plus a contribution from two-particle collisions. (ii) As usual, a collision between a pair of particles is taken to be independent of the others, and local (namely, it leads to sudden changes in the corresponding momenta, whereas the coordinates can be considered as fixed).

The mean-field assumption is reasonable for a system having long-range forces, since they suppress fluctuations. On the other hand, the second assumption concerning the simple cumulative effect of independent collisions was subjected to a certain criticism (see e.g., ), because the long-range unshielded gravitational coupling involve many simultaneously interacting particles. Nevertheless, an impressive amount of observational and numerical material has been obtained, which advocates the use of the present approach al least on certain time-scales .

According to our assumptions the contribution from the direct interaction between the particles is added on the right-hand side of the Liouville equation for the one-particle distribution function $`f(\stackrel{}{r},\stackrel{}{v},t)`$ as the corresponding collision integral $`𝒞[f]`$:

$$\frac{\mathrm{d}f}{\mathrm{d}t}\frac{f}{t}+\stackrel{}{v}\frac{f}{\stackrel{}{r}}\frac{\varphi }{\stackrel{}{r}}\frac{f}{\stackrel{}{v}}=𝒞[f]$$
(1)
The mass of each particle is taken be unity, $`m=1`$, and the potential $`\varphi (\stackrel{}{r})`$ can be viewed as an external despite its self-consistent closure with the Poison equation:

$$\mathrm{\Delta }\varphi =4\pi Gd\stackrel{}{v}f(\stackrel{}{r},\stackrel{}{v},t)=4\pi Gn(\stackrel{}{r})$$
(2)
As usual in kinetic theory, the distribution function $`f(\stackrel{}{r},\stackrel{}{v},t)`$ is normalized on the total number of particles, namely

$`{\displaystyle d\stackrel{}{v}f(\stackrel{}{r},\stackrel{}{v},t)}=n(\stackrel{}{r})`$ (3)
is the density of particles, and

$`{\displaystyle d\stackrel{}{r}n(\stackrel{}{r})}=N`$ (4)
is the total number of particles.

Instead of deriving the collision integral of the kinetic equation from the underlying microscopic dynamics, Bhatnager, Gross and Krook proposed the following model for the collisional integral :

$$𝒞_{BGK}[f]=\nu (r)[f(\stackrel{}{r},\stackrel{}{v},t)f_0(\stackrel{}{r},\stackrel{}{v},t)],$$
(5)
where $`\nu (r)>0`$, and $`f_0`$ is the best local Maxwellian

$$f_0(\stackrel{}{r},\stackrel{}{v},t)=\frac{n(\stackrel{}{r},t)}{(2\pi T(\stackrel{}{r},t))^{3/2}}\mathrm{exp}\left(\frac{(\stackrel{}{v}\stackrel{}{u}(\stackrel{}{r},t))^2}{2T(\stackrel{}{r},t)}\right),$$
(6)
where $`n(\stackrel{}{r},t)`$, $`T(\stackrel{}{r},t)`$, $`\stackrel{}{u}(\stackrel{}{r},t)`$ are defined by the following conditions:

$$n(\stackrel{}{r},t)=d\stackrel{}{v}f(\stackrel{}{r},\stackrel{}{v},t),$$
(7)
$$n(\stackrel{}{r},t)\stackrel{}{u}(\stackrel{}{r},t)=d\stackrel{}{v}\stackrel{}{v}f(\stackrel{}{r},\stackrel{}{v},t),$$
(8)
$$3n(\stackrel{}{r},t)T(\stackrel{}{r},t)=d\stackrel{}{v}(\stackrel{}{v}\stackrel{}{u}(\stackrel{}{r},t))^2f(\stackrel{}{r},\stackrel{}{v},t).$$
(9)
Due to this conditions the collision integral (5) conserves probability, momentum and energy, and in the spatially homogeneous case $`T=T(\stackrel{}{r},t)`$ it enforces convergence of the distribution function $`f(\stackrel{}{r},\stackrel{}{v},t)`$ to the corresponding Gibbs distribution. One should notice that BGK collision integral (5) is still non-linear due to these self-consistency relations. However, this non-linearity is based on the first three moments of $`f`$ only, which means that BGK scheme will be friendly to moments equations and related procedures. Another attractive property of BGK model is that though being proposed ad hoc, it was later derived from the Boltzmann equation in both rarefied (low-density), and dense limits . In particular, relations with more formal approximation methods (e.g., those proposed by Chapman-Enskog or Grad) were established . Thus, there is a general expectation that this approximation will reasonably mimic the main properties of the original Boltzmann equation.

The general conditions of invariance and convergence do not specify $`\nu (r)`$, which has a meaning of the inverse characteristic relaxation time. It should not depend on $`\stackrel{}{v}`$, since else the conservation laws cannot be satisfied. Obviously, the concrete choice of $`\nu (\stackrel{}{r},t)`$ should be determined demanding that Eq. (5) is as close as possible to the original Boltzmann equation. A standard assumption is that $`\nu `$ depends on the coordinate $`\stackrel{}{r}`$ and time $`t`$ through the $`\stackrel{}{r}`$-dependence of the first three moments, $`\nu =\nu (n,\stackrel{}{u},T)`$, and a general procedure was proposed to derive its concrete form . We will not go into details of that procedure, but we will show that already rather simple qualitative considerations determine its rough form, and then we will directly take for it the expression for the inverse characteristic time, which was obtained from the full Boltzmann equation in Refs. . To determine $`\nu `$ one notices that in the rarefied case it is proportional to $`n`$, and then due to obvious dimensional reasons one has

$$\nu (\stackrel{}{r},t)\frac{n(\stackrel{}{r},t)}{T(\stackrel{}{r},t)^{3/2}}.$$
(10)
In this context one could wonder why only $`T`$ and not $`\varphi `$ should enter $`\nu `$. Since $`\nu `$ belongs to the collision integral, according to the common initial assumption the influence of the field is not taken into account there. Indeed, in the full Boltzmann equation binary collisions are considered in such a way that only microscopic two-particle interaction enters there, and not the mean-field. This is inevitable to find tractable schemes, and means that collisions occur on such time and length-scales, where the effect of the mean-field is not relevant.

More precisely, we notice the binary collision relaxation time obtained through the complete Boltzmann equation

$$t_r=\frac{0.065v^2^{3/2}}{nmG^2\mathrm{ln}(0.4N)}$$
(11)
with the right hand side calculated at typical radius $`R`$, and where we momentarily recovered the dimensional units. Thus, we just take $`1/\nu `$ equal to the local relaxation time, so

$`\nu (r)={\displaystyle \frac{1}{t_r(r)}}={\displaystyle \frac{nmG^2\mathrm{ln}(0.4N)}{0.065v^2^{3/2}}}{\displaystyle \frac{15.4\mathrm{ln}(0.4N)G^2m^{5/2}n(r)}{T(r)^{3/2}}}=\gamma {\displaystyle \frac{n(r)}{T(r)^{3/2}}},`$ (12)
$`\gamma =15.4\mathrm{ln}(0.4N)G^2m^{5/2}`$ (13)
is the inverse relaxation time at position $`r`$ . To estimate whether collisions are relevant compared with the motion in the mean-field, one has to compare $`t_{rel}`$ with the dynamical time-scale $`t_d`$. For a cluster with mass $`M=Nm`$ and typical radius $`R`$ the dynamical time is

$$t_d=\frac{R}{\overline{v}}=\frac{R^{3/2}}{\sqrt{GM}}$$
(14)
The situation $`t_r(r)t_d`$ means that collisions dominate (overdamping), and $`t_r(r)t_d`$ indicates underdamping. This qualitative criterion will be applied to determine the transition from the outer weakly-damped part of the system (halo) to the inner strongly-damped part (core).

### 2.1 H-theorem and convergence to equilibrium

Let us now briefly show that if the conditions for equilibrium are satisfied, then the BGK-collisional term (5) indeed induces convergence toward the Gibbs distribution. To this end we will consider the time-behavior of the Boltzmann-Gibbs entropy:

$$S(t)=d\stackrel{}{r}d\stackrel{}{v}f(\stackrel{}{r},\stackrel{}{v},t)\mathrm{ln}f(\stackrel{}{r},\stackrel{}{v},t),$$
(15)
$$\dot{S}=d\stackrel{}{r}d\stackrel{}{v}\nu (r)[ff_0]\mathrm{ln}f+d\stackrel{}{r}d\stackrel{}{v}\left[\stackrel{}{v}\frac{f}{\stackrel{}{r}}\frac{\varphi }{\stackrel{}{r}}\frac{f}{\stackrel{}{v}}\right]\mathrm{ln}f$$
(16)
Using the fact that $`d\stackrel{}{v}[ff_0]\mathrm{ln}f_0=0`$, which is conservation of probability, momentum and energy, one can rewrite the first term in r.h.s. of Eq. (16) as

$$d\stackrel{}{r}\nu (r)d\stackrel{}{v}\left[f\mathrm{ln}\frac{f}{f_0}+f_0\mathrm{ln}\frac{f_0}{f}\right]0,$$
(17)
since

$$d\stackrel{}{v}f\mathrm{ln}\frac{f}{f_0}0$$
(18)
for any distributions $`f`$, $`f_0`$. The second term term in r.h.s. of Eq. (16) can be presented as an integral over the corresponding surface in the coordinate space

$$d\stackrel{}{S}d\stackrel{}{v}\stackrel{}{v}f\mathrm{ln}f$$
(19)
If there are no global currents this integral is zero, and we obtain the H-theorem:

$$\dot{S}0$$
(20)
Monotonically increasing at finite times $`S`$ will saturate in the infinite-time limit: $`f=f_0`$. Substituting this expression in l.h.s of the kinetic equation we will get the Gibbs distribution

$$T=\mathrm{const}.,n(\stackrel{}{r})=\frac{1}{Z}e^{\beta \varphi (\stackrel{}{r})}$$
(21)
If there are global currents in the system, then the actual distribution is not known a priori, and the full kinetic equation has to be solved.

### 2.2 Kinetic equation in spherical coordinates

Since we will consider non-equilibrium distribution functions, which are in a sense close to be isotropic, it is reasonable to employ spherical coordinates:

$`x=r\mathrm{sin}\theta \mathrm{cos}\phi ,y=r\mathrm{sin}\theta \mathrm{sin}\phi ,z=r\mathrm{cos}\theta ,`$ (22)
In the local Cartesian base $`\{\stackrel{}{e}_r`$, $`\stackrel{}{e}_\theta `$, $`\stackrel{}{e}_\phi \}`$ (e.g., $`\stackrel{}{e}_r`$ is a unit vector, which points out from $`r`$ to $`r+\mathrm{d}r`$) velocity will have components

$$\stackrel{}{v}=v_r\stackrel{}{e}_r+v_\theta \stackrel{}{e}_\theta +v_\phi \stackrel{}{e}_\phi ,$$
(23)
where

$`v_x=v_r\mathrm{sin}\theta \mathrm{cos}\phi +v_\theta \mathrm{cos}\theta \mathrm{cos}\phi v_\phi \mathrm{sin}\phi ,`$ (24)
$`v_y=v_r\mathrm{sin}\theta \mathrm{sin}\phi +v_\theta \mathrm{cos}\theta \mathrm{sin}\phi +v_\phi \mathrm{cos}\phi ,`$ (25)
$`v_z=v_r\mathrm{cos}\theta v_\theta \mathrm{sin}\theta ,`$ (26)
In the spherical coordinates the kinetic equation will read

$`{\displaystyle \frac{f}{t}}`$ $`+v_r{\displaystyle \frac{f}{r}}+{\displaystyle \frac{v_\theta }{r}}{\displaystyle \frac{f}{\theta }}+{\displaystyle \frac{v_\phi }{r\mathrm{sin}\theta }}{\displaystyle \frac{f}{\phi }}+\left[{\displaystyle \frac{v_\theta ^2+v_\phi ^2}{r}}{\displaystyle \frac{\varphi }{r}}\right]{\displaystyle \frac{f}{v_r}}`$ (27)
$`+{\displaystyle \frac{\mathrm{cot}\theta }{r}}\left[v_\phi ^2{\displaystyle \frac{f}{v_\theta }}v_\phi v_\theta {\displaystyle \frac{f}{v_\phi }}\right]{\displaystyle \frac{v_r}{r}}\left[v_\theta {\displaystyle \frac{f}{v_\theta }}+v_\phi {\displaystyle \frac{f}{v_\phi }}\right]=𝒞_{BGK}[f]`$
We will consider non-equilibrium states with radial symmetry. We will assume that $`\varphi (r,\theta ,\phi )=\varphi (r)`$, and the distribution function $`f`$ is non-equilibrium exclusively through its radial properties:

$$f(r,\theta ,\phi ,v_r,v_\theta ,v_\phi )=f(r,v_r)\frac{1}{2\pi T(r)}\mathrm{exp}\left[\frac{v_\theta ^2+v_\phi ^2}{2T(r)}\right],$$
(28)
namely there is no dependence on the angular coordinates, and the angular velocities have relaxed to equilibrium. In fact, Eq. (28) is the minimal Ansatz, which allows to study non-equilibrium states within the simplest spherical geometry.

For the collision term one has

$$𝒞_{BGK}[f]=\nu (r)\left[f(r,v_r)f_0(r,v_r)\right]\times \frac{1}{2\pi T(r)}\mathrm{exp}\left[\frac{v_\theta ^2+v_\phi ^2}{2T(r)}\right],$$
(29)
$$f_0(r,v)=\frac{n(r)}{[2\pi T(r)]^{1/2}}\mathrm{exp}\left(\frac{(vu)^2}{2T(r)}\right),$$
(30)
where the mean radial velocity $`u`$ describes an average motion (current of particles) along the radial direction.

Let us now substitute Eq. (28) to Eq. (27), and integrate by $`v_\theta `$, $`v_\phi `$. The result will read

$$\frac{f}{t}+v_rf^{}\varphi ^{}\frac{f}{v_r}+\frac{2T(r)}{r}\left[f(r)\frac{v_r}{T(r)}+\frac{f}{v_r}\right]=\nu (r)\left[f(r,v_r)f_0(r,v_r)\right]$$
(31)
where $`f^{}=f(r,v_r)/r`$, $`\varphi ^{}=\mathrm{d}\varphi /\mathrm{d}r`$.

### 2.3 Moments

Eq. (31) is still complicated integro-differential equation. To subtract the relevant information from it one considers the subtracted and non-subtracted moments

$$M_l(r)=dv_r(v_ru)^lf(r,v_r),$$
(32)
$$_l(r)=dv_rv_r^lf(r,v_r),$$
(33)
These two types of moments correspond to different experimental situations, since $`M_l(r)`$ are measured in the comoving frame, whereas $`_l`$ at the rest. We will use them in parallel choosing the one or the another as appears most convenient.

The corresponding equations read:

$`\dot{M}_l`$ $`+`$ $`l\dot{u}M_{l1}+{\displaystyle \frac{1}{r^2}}[r^2(M_{l+1}+uM_l)]^{}+l\left(\varphi ^{}{\displaystyle \frac{2T}{r}}\right)M_{l1}+lu^{}[M_l+uM_{l1}]`$ (34)
$`=`$ $`\nu (r)[M_l\omega (l)nT^{l/2}(l1)!!],`$
where $`\omega (l)`$ is zero for $`l`$ odd, and equal to $`1`$ for $`l`$ even.

$`\dot{}_l`$ $`+`$ $`{\displaystyle \frac{1}{r^2}}(r^2_{l+1})^{}+l\left(\varphi ^{}{\displaystyle \frac{2T}{r}}\right)_{l1}`$ (35)
$`=`$ $`\nu (r)[_ln{\displaystyle \underset{p=0}{\overset{l}{}}}u^{lp}T^{p/2}\omega (p)(p1)!!],`$
Recall that, by construction, the low moments

$`M_0=n(r),M_1=0,`$ $`M_2=n(r)T(r),`$ (36)
$`_0=n(r),_1=n(r)u(r),`$ $`_2=n(r)[T+u^2],`$ (37)
are the same for $`f_0(r,v_r)`$ and $`f(r,v_r)`$. Since this fact expresses conservation of probability, momentum and energy, the equations (34, 35) with $`l=0,1,2`$ will be universal, namely they do not depend on the concrete form of the collision integral. Notice also that the third moments are connected as:

$$_3=M_3+nu^3+3unT.$$
(38)
Here $`_3`$ is energy current, which consists of heat current $`M_3`$, transfer of kinetic energy $`nuu^2`$, and work done by pressure $`3unT=3uP`$.

Let us write down the first members of Eq. (34)

$`l=0:`$ $`r^2\dot{n}(r)+(r^2nu)^{}=0,`$ (39)
$`l=1:`$ $`n\dot{u}+(nT)^{}+n\varphi ^{}+nuu^{}=0,`$ (40)
$`l=2:`$ $`n\dot{T}+{\displaystyle \frac{1}{r^2}}(r^2M_3)^{}+2u^{}nT+unT^{}=0,`$ (41)
$`l=3:`$ $`\dot{M}_3+3nT\dot{u}+{\displaystyle \frac{1}{r^2}}[r^2(M_4+uM_3)]^{}+3u^{}[M_3+nuT]`$ (42)
$`+3nT(\varphi ^{}{\displaystyle \frac{2T(r)}{r}})=\nu (r)M_3,`$
$`l=4:`$ $`\dot{M}_4+4M_3\dot{u}+{\displaystyle \frac{1}{r^2}}[r^2(M_5+uM_4)]^{}+4u^{}[M_4+uM_3]`$ (43)
$`+4M_3(\varphi ^{}{\displaystyle \frac{2T(r)}{r}})=\nu (r)(M_43nT^2)`$
Eq. (39) expresses conservation of the mass inside the sphere $`r`$. When $`u=0`$ Eq. (40) expresses hydrostatic equilibrium; when $`u0`$ it contains flow terms, that should be easy to understand. The equation for $`l=2`$ expresses the conservation of energy; for obtaining it we have also inserted Eq. (39).

Somewhat different, but of course, totally equivalent expressions will be obtained for the $``$-moments.

Eqs. (39-43) are inherently attached to the Poisson equation:

$$\frac{1}{r^2}(r^2\varphi ^{})^{}=4\pi Gn(r).$$
(44)
The density $`n(r)`$ should satisfy to the condition of integrability,

$`4\pi {\displaystyle _0^{\mathrm{}}}drn(r)r^2=N.`$ (45)
Notice that equations with $`l2`$ can be written in a form, which does not contain $`\varphi ^{}`$ explicitly. In that respect their combination can be viewed as equations of state. For example, using Eq. (40), the for result for $`l=3`$ reads

$$\dot{M}_3+\frac{1}{r^2}[r^2(M_4+uM_33nT^2)]^{}+3u^{}M_3+3nTT^{}=\nu (r)M_3.$$
(46)
## 3 Steady state.

As was already indicated in the introduction, we will consider the stationary (steady) situation. This case, where all quantities are time-independent, is worth to study for its own sake, since there are clear observational evidences that typical globular clusters do have such a phase in their evolution . On the other hand, it is a necessary step towards understanding of more general, time-dependent situation.

In the situation, where the known conditions for equilibrium are satisfied (e.g., absence of energy and/or probability currents), the self-gravitating system under study has the Gibbs distribution as its only steady state. This point was already illustrated by us when considering the H-theorem in section 2.1. In certain range of temperatures this distribution is physical (e.g., it is at least metastable), but it loses its stability for temperatures lower than a certain critical temperature . This is the notorious phenomenon of gravo-thermal instability (collapse) . However, there are situations, where one expects that the very reasons for the existence of the equilibrium can be invalid. The most typical situation of that kind appears with binary star formation and tightening processes in the central region of a star cluster . These processes are accompanied with a release of energy, so that the rest of the system (namely its outer with respect to a relatively small part, where the binary-formation process occurs) can be considered as being subjected to sources of energy put in the central part. The existence of energy currents from the central part leads to a non-equilibrium steady state. Additionally, they can stabilize the behavior of the self-gravitating system preventing its collapse, since the mechanism of the gravothermal instability is connected with spontaneous energy currents towards the center .

A somewhat more exotic example of a non-equilibrium steady state can be provided by the existence of a black hole in the central part of a star cluster . Then one has a stationary current of consumed stars, which is directed towards the center, as well as a related current of energy.

Our setup with the moments equations is especially suitable for describing the above non-equilibrium states, since there are explicit expressions for the energy current, $`_3`$ and the current of particles $`u`$.

### 3.1 Ideal hydrodynamics

This scheme is determined by a condition $`M_3=0`$ (no heat current). The nonequilibrium character is displayed only through the stationary current of particles $`u`$.

Eqs. (39-41) read

$`(r^2nu)^{}=0,(nT)^{}+n\varphi ^{}+nuu^{}=0,2u^{}nT+unT^{}=0.`$ (47)
Recall that the usual gibbsian equilibrium is a particular case of this situation, which is realized for $`u=0`$, and an additional assumption that $`T=\mathrm{const}`$. In this case Eqs. (47, 47) are trivial, and Eq. (47) leads to the Gibbs distribution.

Let us first assume that $`u`$ does not depend of $`r`$. Then Eq. (47) predicts that $`T=\mathrm{const}`$. With Eq. (47) one gets

$$n=\frac{c}{u}\frac{1}{r^2},c=\mathrm{const},$$
(48)
Notice that similar to the isothermal case this distribution is non-integrable at infinity, and therefore predicts infinite mass. Being combined with the Poisson equation (44), Eq. (47) offers an exact solution:

$$\varphi =2T\mathrm{ln}r+\varphi _0.$$
(49)
To study the case $`u\mathrm{const}`$ one first notice that Eq. (47) can be exactly integrated

$$T=\frac{c_1}{u^2},c_1=\mathrm{const}$$
(50)
Eqs. (48, 50) combined with (47-47) lead to

$`c_1ur^2(u^3r^2)^{}\varphi ^{}+u^{}u=0,u(r^2\varphi ^{})^{}=4\pi Gc.`$ (51)
To study self-similar solutions, one makes an Ansatz $`u=u_0r^\alpha `$, and gets from Eq. (51)

$$c_1u_0^2(3\alpha +2)r^{2\alpha 1}+\varphi ^{}+\alpha u_0^2r^{2\alpha 1}=0$$
(52)
Consulting with the Poisson equation one observes that for small distances the last term in Eq. (52) can be neglected, and as the result one has the following asymptotic solution at small distances $`\alpha =1/2`$ and

$`\varphi ^{}{\displaystyle \frac{7c_1}{2u_0^2}}r^2+{\displaystyle \frac{8\pi Gc}{u_0}}r^{3/2},T{\displaystyle \frac{c_1}{u_0^2}}r^1,n{\displaystyle \frac{c}{u_0}}r^{5/2}.`$ (53)
It is seen that the density is singular at the origin, but still integrable. The first term in $`\varphi ^{}`$ indicates the presence of the central mass $`7c_1/(2u_0^2G)`$. This implies the following choice: $`u_0<0`$, $`c<0`$, and then the obtained solution describes a stationary consumption regime of stars by the central mass. The energy current $`_3=nu^3+3nuT`$ is also directed towards the center.

### 3.2 Overdamped situation

Starting since this moment, we will into the consideration an energy current, described by $`M_3`$. We make a self-consistent assumption that there is a range of distances close to the inner part (i.e., the core), where collisions dominate. In this overdamped regime $`\eta =(\mathrm{dynamic}\mathrm{time})/(\mathrm{kinetic}\mathrm{time})`$ is large, and the large $`\nu (r)`$ in the r.h.s. of Eqs. (42, 43) imposes the Gaussian behavior for the corresponding moments:

$`M_3=𝒪(1/\eta ),M_43nT^2=𝒪(1/\eta )`$ (54)
Using the last result as $`M_43nT^2`$ one gets from Eq. (42)

$`3n(r)T(r)T^{}(r)=\nu (r)M_3(r)`$ (55)
as a constitutive equation for the overdamped case. Let us consider separately cases with $`u=0`$, and $`u0`$.

The case: $`u=0`$

Then Eq. (41) reduces to

$$M_3r^2=\nu _3=\mathrm{const}$$
(56)
If $`\nu _3>0`$ there is a solution

$$T=\left(\frac{2\gamma \nu _3}{21}(\frac{1}{r}+\frac{1}{r_0})\right)^{2/7}$$
(57)
with some integration constant $`r_0`$. The condition $`\nu _3>0`$ means a positive outward flow of energy. One thus gets for small $`r`$

$$TT_0r^{2/7},T_0=\left(\frac{2\gamma \nu _3}{21}\right)^{2/7}.$$
(58)
Taking this into account one easily gets from the Poisson equation and hydrostatic equilibrium

$$n(r)n_0r^{16/7},n_0=\frac{1}{4\pi G}\frac{90}{49}T_0,\frac{16}{7}=2.2857143,$$
(59)
$$\varphi ^{}(r)\varphi _0r^{9/7},\varphi _0=\frac{18}{7}T_0.$$
(60)
Notice that density is singular at the origin, but integrable. Thus, this solution can be intended to describe a stationary regime without a central mass (as follows from Eq. (60)), but with a stationary tightening of binaries at the very central region. In our setup this last process is reflected through a constant outward heat current.

The mass insider a sphere of radius $`r`$ scales as $`M(r)r^{5/7}`$. The inverse is $`rM^{7/5}`$. These are so-called Lagrange radii. If one chooses a set of equidistant $`M`$ values, e.g. $`M=i/20`$, $`i=1,\mathrm{},20`$, then the radii will be closer to each other near the center than in the outer parts of the cluster. This happen one since $`16/7>2`$. If, on the other hand, the density would be finite at the center, then one would have $`rM^{1/3}`$, implying large separations between Lagrange radii near $`M=0`$.

The case: $`u0`$

Here one can still use Eq. (55) for $`M_3`$, which now does not reduce to $`\mathrm{const}/r^2`$. Having substituted this equation to Eq. (41), one obtains

$`l=0:`$ $`(r^2un)^{}=0,l=1:(nT)^{}+n\varphi ^{}+u^{}un=0,`$ (61)
$`l=2:`$ $`{\displaystyle \frac{3}{\gamma }}{\displaystyle \frac{1}{r^2}}[r^2T^{5/2}T^{}]^{}+unT^{}+2u^{}nT=0,`$ (62)
For small distances the self-consistent solution of these equations reads

$`n(r)n_0r^{12/5},T(r)T_0r^{2/5},u(r)u_0r^{2/5},\varphi ^{}(r){\displaystyle \frac{14}{5}}T_0r^{7/5},`$ (63)
where $`6T_0^{5/2}=5\gamma u_0n_0`$. Notice that there is no central mass (as follows from Eq. (63)). Moreover, $`u_0`$ has to be positive and $`u(0)=0`$. Therefore, the obtained solution can be suitable to describe evaporation phenomena.

### 3.3 Underdamped situation

There is yet another extremal situation realized at sufficiently large distances from the center, where the behavior of the system is almost collisionless (underdamping).

The case with $`u=0`$

In the underdamping regime one expects that all moments with $`k3`$ will be equally important, therefore they have to behave in the nearly similar way. Looking for a self-consistent, long-distance solution of the hierarchy (39-43) one gets:

$`n(r){\displaystyle \frac{n_0}{r^{7/2}}},T(r){\displaystyle \frac{T_0}{r}},\nu (r)\gamma {\displaystyle \frac{n_0T_0^{3/2}}{r^2}},`$ (64)
$`M_3={\displaystyle \frac{\nu _3}{r^2}},M_4{\displaystyle \frac{\nu _4}{r^2}}+\gamma {\displaystyle \frac{\nu _3n_0T_0^{3/2}}{r^3}},`$ (65)
where $`\nu _4`$, $`n_0`$, $`T_0`$ are constant. At this stage of our asymptotic analysis they will remain unfixed. For $`\varphi ^{}`$ one has from the Poisson equation:

$$\varphi ^{}(r)=\frac{GM}{r^2}\frac{8\pi Gn_0}{r^{5/2}},$$
(66)
where the first term is the standard asymptotic limit for a finite-mass cluster with $`M`$ as the total mass. The second term follows from Eq. (64). The behavior of $`M_5`$ is qualitatively the same as of $`M_4`$ as seen from Eq. (43). Our expression for the density given by Eq. (64) coincides with that obtained in  for the halo of a globular cluster by means of qualitative physical considerations.

Having substituted Eqs. (64, 64) to hydrostatic equilibrium equation (40), one gets the correction terms to the density and temperature

$$n(r)=\frac{n_0}{r^{7/2}}\left[1+\frac{n_1}{r^{1/2}}+\mathrm{}\right],T(r)=\frac{T_0}{r}\left[1+\frac{T_1}{r^{1/2}}+\mathrm{}\right],$$
(67)
$$5T_0(n_1+T_1)=GM(n_1+8\pi Gn_0)$$
(68)
with constant $`n_1`$, $`T_1`$.

The case with $`u0`$

This case is capable to describe the evaporation phenomenon  as we indicated already in section 3.2. Having this in mind, we will look for a solution, which has the Maxwellian structure for very long distances, since there it is supposed to describe free escaped stars. In other words, for those scales only first two moments $`u=u_0>0`$ and $`T=T_0>0`$ are nonzero, and should be viewed as boundary conditions. Possible solutions of this kind should be then matched with the overdamped solutions having finite $`u`$. Notice that the equation (41) for energy can be written in a simplified form

$`\nu _3(r)^{}+{\displaystyle \frac{J(r)}{u^2(r)}}[u^2(r)T(r)]^{}=0,`$ (69)
where $`M_3=\nu _3(r)r^2`$, and where $`J=r^2nu`$ is the constant current of particles.

Searching again for a self-consistent solution, one gets

$`r^2n(r)n_0+{\displaystyle \frac{n_{10}+n_{11}\mathrm{ln}r}{r}}+{\displaystyle \frac{n_{20}+n_{21}\mathrm{ln}r+n_{22}(\mathrm{ln}r)^2}{r^2}},`$ (70)
$`T(r)T_0+{\displaystyle \frac{T_{10}+T_{11}\mathrm{ln}r}{r}}+{\displaystyle \frac{T_{20}+T_{21}\mathrm{ln}r+T_{22}(\mathrm{ln}r)^2}{r^2}},`$ (71)
$`r\varphi ^{}\varphi _0+{\displaystyle \frac{\varphi _{10}+\varphi _{11}\mathrm{ln}r}{r}}+{\displaystyle \frac{\varphi _{20}+\varphi _{21}\mathrm{ln}r+\varphi _{22}(\mathrm{ln}r)^2}{r^2}},`$ (72)
$`u(r)u_0+{\displaystyle \frac{u_{10}+u_{11}\mathrm{ln}r}{r}}+{\displaystyle \frac{u_{20}+u_{21}\mathrm{ln}r+u_{22}(\mathrm{ln}r)^2}{r^2}},`$ (73)
$`\nu _3(r)\nu _{30}+{\displaystyle \frac{\nu _{310}+\nu _{311}\mathrm{ln}r}{r}}+{\displaystyle \frac{\nu _{320}+\nu _{321}\mathrm{ln}\nu _{322}(\mathrm{ln}r)^2}{r^2}},`$ (74)
$`\nu (r)\gamma {\displaystyle \frac{n_0T_0^{3/2}}{r^2}}.`$ (75)
Higher-order terms can be presented analogously. Notice that density does not have to be integrable at infinity when considering stationary regimes of evaporation. Indeed, the only way to have a truly stationary evaporation regime is to have an infinite amount of stars, the most of them being located at infinity. The appearance of logarithms in the asymptotic expansions is mathematically connected with non-integrable character of the density.

It is seen as well that the leading-order behavior of $`\nu (r)`$ is the same for both underdamped regimes (with and without evaporation).

At the present stage of our asymptotic analysis we able to obtain only certain relations between the coefficients of (70-75).

$`(r^2\varphi ^{})^{}=4\pi Gr^2n(r)\varphi _0=4\pi Gn_0,n_{11}=u_{11}=0,\varphi _{11}=4\pi Gn_{10};`$ (79)
$`unr^2=J,J=u_0n_0,n_0u_{10}+n_{10}u_0=0,n_0u_{11}+n_{11}u_0=0;`$
$`(nT)^{}+\varphi ^{}+{\displaystyle \frac{n}{2}}(u^2)^{}=0,\varphi _0=2T_0,`$
$`n_0\varphi _{10}2n_0u_0u_{10}=3T_{10}n_03T_{11}n_0+2n_{10}T_0,3T_{11}=\varphi _{11};`$
$`\nu _3^{}+{\displaystyle \frac{J}{u^2}}[u^2T]^{}=0,`$ (80)
$`J(T_{11}T_{10})2T_0u_{10}+\nu _{311}+\nu _{310}=0,JT_{11}=\nu _{311}.`$
The equations support an additional simplification $`T_{10}=\varphi _{10}=0`$, and consequently $`4\pi Gn_0=T_02u_0^2`$. This equation connects the energetic characteristics of the considered solution.

## 4 Conclusion

This contribution was devoted to the kinetic theory of many-body self-gravitating systems in the framework of the model kinetic equation developed by Bhatnager, Gross and Krook . The main advantage of using this kinetic model (as well as other, more refined kinetic models) is that it produces rather tractable schemes in contrast to the full kinetic equations, which are typically quite difficult to investigate. On the other hand, model kinetic equations offer rather sensible results, which are usually at least in a qualitative agreement with those obtained from the full Boltzmann kinetic equation. The physical reason of this success is that a good amount of statistical phenomena are quite insensitive to the details of the collision integral. In order to explain them it is enough to use the simplest relaxation mechanism which only takes into account necessary conservation laws of energy, momentum and probability (number of particles) and the correct characteristic relaxation time.

The assumptions on the spherically symmetric character of the considered system and the assumed radial character of non-equilibrium led us to the kinetic equation (1), which was then studied by the method of moments. The straightforward analytic investigation allowed us to identify several self-similar regimes of a steady but non-equilibrium behavior. These solutions corresponds to the presence of the central black hole in the system (section 3.1), the influence of binary-formation and tightening processes to the system (section 3.2) and evaporation phenomena (section 3.3). In this way we reproduced certain known results  and also obtained a number of new regimes of behavior.

We believe that the model kinetic approach will be potentially quite powerful in application to many-body, self-gravitating systems. We plan in future to turn to the dynamical aspects of this approach, as well as construct more realistic (from the viewpoint of applications in astrophysics) kinetic models, that e.g account for the strong anisotropy of the velocity distribution in the outer part of the cluster.