text:
# Maximal work extraction from quantum systems

## Abstract

Thermodynamics teaches that if a system initially off-equilibrium is coupled to work sources, the maximum work that it may yield is governed by its energy and entropy. For finite systems this bound is usually not reachable. The maximum extractable work compatible with quantum mechanics (“ergotropy”) is derived and expressed in terms of the density matrix and the Hamiltonian. It is related to the property of majorization: more major states can provide more work. Scenarios of work extraction that contrast the thermodynamic intuition are discussed, e.g. a state with larger entropy than another may produce more work, while correlations may increase or reduce the ergotropy.

The generality of the laws of thermodynamics for macroscopic bodies led to a long-lasting effort to derive them from microphysics landau ; balian ; thirring . This program is by now completed. In contrast, for finite systems, the application and even the formulation of the laws of thermodynamics are still the subject of studies AN ; QL2L . The origin of qualitative differences between large and finite systems was recognized long ago old1 ; old2 . A quantum system submitted to time-dependent external potentials that describe work sources undergoes a unitary transformation. During such an evolution, the density matrix has constant eigenvalues, thus it cannot become gibbsian when starting from an arbitrary initial state. In contrast, most macroscopic systems have a thermodynamic behavior: they evolve close to a Gibbs state under the effect of slowly varying external potentials, keeping nearly constant entropy old2 ; woron . This holds in general, although their evolution is Hamiltonian at the microscopic level, owing to the large and smooth density of states of the spectrum old2 . Accordingly, the responses of finite and of infinite systems to external perturbations differ qualitatively old1 . As an example, processes were studied which would be reversible in thermodynamics, but entail a specific irreversibility due to the finite size of the involved system sato ; satotasaki .

The problem we treat here is well known in thermodynamics. It initiated its birth in 1824: What is the maximal amount of work that can be extracted from a system S by means of an external source of work acting cyclically in a thermally isolated process? And what are the criteria for comparing different states with respect to their work-providing ability? We find answers to both questions, and show that for finite systems they qualitatively disagree with the standard thermodynamical ones.

A specific case in which a finite system S behaves thermodynamically in regard to its work production is well known. If S is initially prepared in a gibbsian state, it cannot produce work when coupled to a source of work through an external variable which returns to its initial value woron ; ANThomson ; thirring . Such a state is thus called passive. We shall therefore be interested in off-equilibrium initial states. A scenario for preparing such states is to take S consisting of two non-interacting (or weakly interacting) subsystems, to set separately the two parts in thermal contact with heat baths at different temperatures, and then to decouple them from the baths. Subsequent coupling with a source of work may yield an amount of work that we wish to evaluate.

The maximal work-extraction problem is thus posed in the following way landau ; callen . Consider a system S which can exchange work with external macroscopic sources. The evolution of its density operator $`\rho (t)`$ is then generated by a Hamiltonian $`H(t)=H+V(t)`$, where the time-dependence of $`V(t)`$ accounts for work transfer. Following Refs. landau ; callen , we call “cyclic” a process in which S, originally isolated, is coupled at the time $`t=0`$ to external sources of work, and decouples from them at the time $`\tau `$. Thus, the driving variables of the sources are cyclic, and the potential $`V(t)`$ vanishes before $`t=0`$ and after $`\tau `$: $`V(0)=V(\tau )=0`$. However, S need not return to its initial state at the time $`\tau `$. The system S is thermally isolated but may involve energy exchanges between its parts. The initial state $`\rho (0)=\rho _0`$ and the Hamiltonian $`H`$ being given, we look for the maximum work $`𝒲`$ that may be extracted from S for arbitrary $`V(t)`$. According to the dynamics

$`\mathrm{i}\mathrm{}\dot{\rho }=[H(t),\rho (t)],`$ (1)
the work $`\mathrm{d}W=\mathrm{tr}[\rho (t)\dot{V}(t)]\mathrm{d}t`$ done on S is landau ; balian ; callen

$`W=\mathrm{tr}[\rho (\tau )H]E(\rho _0),E(\rho _0)\mathrm{tr}[\rho _0H]E_\mathrm{i}.`$ (2)
Among all final states $`\rho (\tau )`$ reached from $`\rho _0`$ under the action of any potential $`V(t)`$, we are therefore looking for the one with lowest final energy $`E_\mathrm{f}=\mathrm{tr}[\rho (\tau )H]`$.

The standard answer landau ; callen  relies on the idea that the final state $`\rho (\tau )`$ is Gibbsian and that its von Neumann entropy $`S=\mathrm{tr}\rho \mathrm{ln}\rho `$ cannot decrease between the times $`0`$ and $`\tau `$. The maximum value of the work (2) is reached when the final state has the equilibrium form

$`\rho (\tau )=\rho _{\mathrm{eq}}={\displaystyle \frac{e^{\beta H}}{Z}},Z=\mathrm{tr}e^{\beta H},`$ (3)
with $`\beta >0`$ determined by the equality of the initial and final entropies: $`\mathrm{ln}Z\beta \mathrm{ln}Z/\beta =\mathrm{tr}\rho _0\mathrm{ln}\rho _0`$. The largest amount of work $`W`$ extractable from S is thus

$`𝒲_{\mathrm{th}}=E(\rho _0)TS(\rho _0)+T\mathrm{ln}Z,`$ (4)
the familiar difference of free energy between initial and final state, both evaluated with the final temperature $`T`$.

Finite systems. The above derivation involves two arguments which call for some discussion. Following thermodynamical intuition, we have first stated that the entropy cannot decrease. In fact, the von Neumann entropy $`S(\rho )`$ remains constant during the evolution (1). If S is macroscopic, it is the coarse-grained entropy which can increase. Anyhow, this point is harmless since the bound (4) corresponds to constant-entropy processes of the total system. We also implicitly assumed that S may be brought into an equilibrium state $`\rho _{\mathrm{eq}}`$ by means of some evolution (1). For macroscopic systems this is usually allowed since dissipative processes within S may occur while it evolves under the influence of the coupling $`V(t)`$.

For finite systems, the sole action of $`V(t)`$ is in general not sufficient to allow reaching at the time $`\tau `$ a gibbsian state of the form (3). Indeed, not only is the entropy $`S(\rho )`$ conserved during the evolution generated by (1), but all the eigenvalues of $`\rho `$old1 ; old2 . In contrast to thermodynamic systems, finite systems keep memory of their initial state and do not involve any relaxation mechanism. One may therefore expect that the maximal amount of work $`𝒲`$ extracted from S is generally smaller than $`𝒲_{\mathrm{th}}`$.

More precisely, the evolution (1) of $`\rho `$ is unitary, so that $`\rho (\tau )=U\rho _0U^{}`$. We look for the minimum of the final energy $`E_\mathrm{f}=\mathrm{tr}U\rho _0U^{}H`$ over all unitary operators $`U`$. We can parameterize the variations $`\delta U`$ of $`U`$ as $`\delta U=XU`$, where $`X`$ is an arbitrary infinitesimal antihermitian operator. Hence, we find $`\delta E_\mathrm{f}=\mathrm{tr}(XU\rho _0U^{}HU\rho _0U^{}XH)=\mathrm{tr}X[\rho (\tau ),H]`$. The stationarity of $`E_\mathrm{f}`$ thus implies that $`\rho (\tau )`$ should commute with $`H`$ and have the same eigenvalues as $`\rho _0`$, a condition which replaces (3). In the spectral resolutions of $`\rho _0`$ and $`H`$,

$`\rho _0={\displaystyle \underset{j1}{}}r_j|r_jr_j|,H={\displaystyle \underset{k1}{}}\epsilon _k|\epsilon _k\epsilon _k|,`$ (5)
we order the eigenvalues as

$`r_1r_2\mathrm{},\epsilon _1\epsilon _2\mathrm{}`$ (6)
The minimum of $`E_\mathrm{f}`$ is then $`_jr_j\epsilon _j`$, and it is reached for

$`\rho (\tau )={\displaystyle \underset{j}{}}r_j|\epsilon _j\epsilon _j|,`$ (7)
which is stationary since it commutes with $`H`$. This result is consistent with the extension for finite systems of the second law in Thomson’s formulation woron ; ANThomson ; thirring : the two conditions which characterize the state $`\rho (\tau )`$ (commutation with $`H`$ and ordering (6)) are the ones which ensure that this state is passive: no further work can be extracted from S after time $`\tau `$ by means of cyclic processes.

Note that, if the spectrum of $`H`$ involves degeneracies which have no counterpart in $`\rho _0`$, the final state $`\rho (\tau )`$ is not unique (in contrast to $`\rho _{\mathrm{eq}}`$ associated with $`𝒲_{\mathrm{th}}`$).

Altogether, the maximum of the amount of work (2) that can be extracted from S is

$`𝒲={\displaystyle \underset{j,k}{}}r_j\epsilon _k(|r_j|\epsilon _k|^2\delta _{jk}).`$ (8)
For $`𝒲`$, which depends only on the initial state and Hamiltonian, we coin the name ergotropy ($`{}_{}{}^{}\epsilon \rho \gamma o\nu `$: work; $`\tau \rho o\pi \eta ^{}`$: transformation, turn). By construction, we have $`𝒲_{\mathrm{th}}𝒲0`$. The ergotropy $`𝒲`$ vanishes if $`\rho _0`$ is passive. It equals the thermodynamical upper bound $`𝒲_{\mathrm{th}}`$ only if there exist two numbers $`\beta `$ and $`Z`$ such that the eigenvalues (5, 6) of $`\rho _0`$ and $`H`$ satisfy $`\mathrm{ln}r_j=\beta \epsilon _j\mathrm{ln}Z`$, so as to allow $`\rho (\tau )=U\rho _0U^{}`$ to reach a gibbsian form (3) in spite of the lack of thermalization mechanism. Noticeable examples include: (i) a pure initial state of S; (ii) two-level systems; (iii) harmonic oscillators in case $`\rho _0`$ is a gaussian state, since both the sequences $`\epsilon _k`$ and $`\mathrm{ln}r_j`$ are then equidistant.

For macroscopic systems the difference $`𝒲_{\mathrm{th}}𝒲`$ is typically relatively small, since the final state (7) may lie close to an equilibrium state if $`\tau `$ is large; the spectra of $`\mathrm{ln}\rho _0`$ and $`H`$ are dense, and a linear relation between them is approximately satisfied in the relevant region. However, for finite systems $`𝒲_{\mathrm{th}}𝒲`$ can be significant.

It remains to show that the bound $`𝒲`$ can actually be reached by coupling S with some source of work that realizes a cyclic process. We thus want to find a time $`\tau `$ and an interaction $`V(t)`$ which vanishes at $`t=0`$ and at $`t=\tau `$, so that, when added to the Hamiltonian $`H`$, $`V(t)`$ leads from the initial state $`\rho _0`$ to $`\rho (\tau )`$ defined in (7). An evolution operator which realizes this goal is

$`U={\displaystyle \underset{j}{}}|\epsilon _jr_j|.`$ (9)
A Hamiltonian $`H+V(t)`$ generates in the interaction representation an evolution operator $`U_\mathrm{I}`$ which satisfies

$`i\mathrm{}{\displaystyle \frac{\mathrm{d}U_\mathrm{I}(t)}{\mathrm{d}t}}=e^{iHt/\mathrm{}}V(t)e^{iHt/\mathrm{}}U_\mathrm{I}(t),`$ (10)
with $`U_\mathrm{I}(0)=1`$. We define $`U_\mathrm{I}(\tau )e^{iH\tau /\mathrm{}}U=e^{i\mathrm{\Lambda }\tau /\mathrm{}}`$, where $`\mathrm{\Lambda }`$ is obtained by diagonalization. We choose for $`U_\mathrm{I}(t)`$ the simple form $`U_\mathrm{I}(t)=e^{i\mathrm{\Lambda }\phi (t)/\mathrm{}}`$, where $`\phi (0)=\dot{\phi }(0)=\dot{\phi }(\tau )=0`$ and $`\phi (\tau )=\tau `$. Then the potential $`V(t)=\dot{\phi }(t)e^{iHt/\mathrm{}}\mathrm{\Lambda }e^{iHt/\mathrm{}}`$ describes according to (10) a source of work that extracts during the time $`\tau `$ the work (8) from S. This potential is far from unique.

In case not only $`\rho (\tau )`$ but also the initial state is stationary, $`[\rho _0,H]=0`$, we can choose the same eigenbases for $`\rho _0`$ and $`H`$, but (6) implies that $`|r_1,|r_2,\mathrm{}`$ in (5) are deduced from $`|\epsilon _1,|\epsilon _2,\mathrm{}`$ by some permutation. Then the matrix $`|r_j|\epsilon _k|^2`$ in (8) is a permutation matrix. For instance, if the lowest two levels $`\epsilon _1<\epsilon _2`$ have initially the inverted populations $`r_2<r_1`$, respectively, so that $`|r_1=|\epsilon _2`$ and $`|r_2=|\epsilon _1`$, we may easily implement the transformation $`U`$, a permutation that interchanges $`r_1`$ and $`r_2`$, either in a rapid or in a quasistatic regime.

Comparison of activities. We wish to compare two states $`\rho _0`$ and $`\sigma _0`$ of a system S as regards the maximum work that they may provide. To make such a comparison meaningful, we assume the initial energies to be the same, $`E(\rho _0)=E(\sigma _0)`$. If S is macroscopic and can reach equilibrium at the end of the process, $`𝒲_{\mathrm{th}}`$ depends only on the entropy $`S`$ of the initial state, and it decreases when $`S`$ increases, since $`\mathrm{d}𝒲_{\mathrm{th}}/\mathrm{d}S`$ is the temperature of the final state reached. However, the situation is different for finite systems. Consider, for instance, a three-level system with eigenenergies $`\epsilon _{1,3}=1`$, $`\epsilon _2=0`$. The eigenstates of $`\rho _0`$ and $`\sigma _0`$ are taken as $`|r_{1,3}=|s_{1,3}=(|\epsilon _1|\epsilon _3)/\sqrt{2}`$, $`|r_2=|s_2=|\epsilon _2`$, so that $`E(\rho _0)=E(\sigma _0)=0`$. If their eigenvalues are $`\{r_j\}=\{0.90,0.08,0.02\}`$ and $`\{s_j\}=\{0.91,0.05,0.04\}`$, the entropy $`S(\rho _0)0.375`$ exceeds $`S(\sigma _0)0.364`$. Accordingly, the thermodynamic bound $`𝒲_{\mathrm{th}}(\rho _0)0.882`$ for the work is smaller than $`𝒲_{\mathrm{th}}(\sigma _0)0.887`$. Nevertheless, the ergotropy $`𝒲(\rho _0)=0.88`$ of $`\rho _0`$ is larger than the ergotropy $`𝒲(\sigma _0)=0.87`$ of $`\sigma _0`$. The actually reachable bounds are, as expected, lower than the corresponding $`𝒲_{\mathrm{th}}`$’s, but they are reversed in order: the entropically more disordered state $`\rho _0`$ may provide more work.

Thus, the entropy criterion fails for comparing the ergotropies. The theory of majorization major , that we briefly recall now, provides another criterion which may be helpful. In quantum statistical mechanics thirring ; fabio  a density operator $`\rho `$ is said to majorize $`\sigma `$ if their eigenvalues $`r_j`$ and $`s_j`$, set in the decreasing order (6), satisfy

$`{\displaystyle \underset{j=1}{\overset{k}{}}}r_j{\displaystyle \underset{j=1}{\overset{k}{}}}s_j,\mathrm{for}\mathrm{any}k1.`$ (11)
This property, denoted as $`\rho \sigma `$, is transitive ($`\rho \sigma `$ and $`\sigma \tau `$ imply $`\rho \tau `$). It characterizes order, but in a stronger way than entropy since $`\rho \sigma `$ implies not only $`S(\rho )S(\sigma )`$, but also $`\mathrm{tr}f(\rho )\mathrm{tr}f(\sigma )`$ for any concave function $`f(x)`$. Pure states majorize all states, while in a Hilbert space of dimension $`n`$ all states majorize $`\rho =1/n`$. If we have both $`\rho \sigma `$ and $`\sigma \rho `$, then $`\rho `$ and $`\sigma `$ are unitarily equivalent. However, the order defined by majorization is incomplete since, for $`n3`$, pairs of states $`\rho `$ and $`\sigma `$ exist of which neither majorizes the other.

Returning to the comparison of the activities of $`\rho _0`$ and $`\sigma _0`$ with $`E(\rho _0)=E(\sigma _0)`$, we find after summation by parts the ergotropy difference $`\delta 𝒲𝒲(\rho _0)𝒲(\sigma _0)`$ as

$`\delta 𝒲={\displaystyle \underset{j1}{}}(s_jr_j)\epsilon _j={\displaystyle \underset{k1}{}}(\epsilon _{k+1}\epsilon _k){\displaystyle \underset{j=1}{\overset{k}{}}}(r_js_j).`$ (12)
Hence a sufficient condition for $`\rho _0`$ to be more active than $`\sigma _0`$ is $`\rho _0\sigma _0`$ ($`S(\rho _0)S(\sigma _0)`$ alone is neither necessary nor sufficient). There exists a wide class of non-unitary evolutions fabio  which lead from a state $`\rho _0`$ to $`\sigma _0`$ such that $`\rho _0\sigma _0`$. If $`\rho _0`$ and $`\sigma _0`$ have the same energy, we have $`𝒲(\rho _0)𝒲(\sigma _0)`$. For instance, the diagonal part $`\overline{\rho }_0=_k|\epsilon _k\epsilon _k|\rho _0|\epsilon _k\epsilon _k|`$ of $`\rho _0`$ has the same energy as $`\rho _0`$ itself, but $`𝒲(\rho _0)𝒲(\overline{\rho }_0)`$ since $`\rho _0\overline{\rho }_0`$ major . In order to find opposite behaviors in the comparison of ergotropies and of free energies (or entropies), we have to search for cases when $`\rho _0\sigma _0`$ and $`\rho _0\sigma _0`$. All pairs of inequalities $`𝒲(\rho _0)𝒲(\sigma _0)`$ and $`S(\rho _0)S(\sigma _0)`$ may then occur, as illustrated by examples given above and below.

Auxiliary system. If S is supplemented with an auxiliary system $`\mathrm{\Omega }`$ with Hamiltonian $`H_\mathrm{\Omega }`$ and initial state $`\omega _0`$, the overall Hamiltonian $`H+H_\mathrm{\Omega }+V(t)`$, where $`V(t)`$ couples $`\mathrm{S}+\mathrm{\Omega }`$ with external sources, generates a unitary transformation in the product Hilbert space, which is more general than for work sources coupled separately to S and $`\mathrm{\Omega }`$. The initial state $`\rho _0\omega _0`$ is uncorrelated and the evolution conserves its (factorized) eigenvalues. The ergotropy satisfies $`𝒲(\rho _0\omega _0)𝒲(\rho _0)+𝒲(\omega _0)`$, the same inequality as for $`𝒲_{\mathrm{th}}`$. Consider, however, again two states $`\rho _0`$ and $`\sigma _0`$ with the same energy. If they are macroscopic and satisfy $`𝒲_{\mathrm{th}}(\rho _0)>𝒲_{\mathrm{th}}(\sigma _0)`$, this ordering of activities is not changed by the introduction of the auxiliary system, since the additivity of entropy for the initial states implies $`𝒲_{\mathrm{th}}(\rho _0\omega _0)>𝒲_{\mathrm{th}}(\sigma _0\omega _0)`$. But for finite systems the order of ergotropies can be reversed. As a first example consider for S the same three-level system as above, with eigenvalues of $`\rho _0`$ and $`\sigma _0`$ now equal to $`\{r_j\}=\{0.8,0.1,0.1\}`$ and $`\{s_j\}=\{0.5,0.5,0\}`$. For $`\mathrm{\Omega }`$ we take a two-level system with eigenenergies $`0`$ and $`\mathrm{\Delta }>0`$, initially in a pure state. Although $`\rho _0\sigma _0`$, we have both $`S(\rho _0)0.639<S(\sigma _0)0.693`$ and $`𝒲(\rho _0)=0.7>𝒲(\sigma _0)=0.5`$. However, coupling with $`\mathrm{\Omega }`$, which does not change the entropies, reverses the inequality for the ergotropies if $`\mathrm{\Delta }<1/4`$, since (for $`\mathrm{\Delta }<1`$) $`𝒲(\rho _0\omega _0)𝒲(\sigma _0\omega _0)=0.4\mathrm{\Delta }0.1`$. For $`\mathrm{\Delta }=0`$, the auxiliary system $`\mathrm{\Omega }`$ does not contribute to the energy balance but nevertheless raises $`𝒲`$ owing to its order.

The opposite situation is also possible, provided S has at least four levels. Consider for instance a system S with eigenenergies $`\epsilon _{1,4}=1`$, $`\epsilon _{2,3}=(1x)`$, with $`0<x<1`$. As eigenstates of $`\rho _0`$ and $`\sigma _0`$ we take $`|r_{1,4}=|s_{1,4}=(|\epsilon _1\pm |\epsilon _4)/\sqrt{2}`$, $`|r_{2,3}=|s_{2,3}=(|\epsilon _2\pm |\epsilon _3)/\sqrt{2}`$, which ensures $`E(\rho _0)=E(\sigma _0)=0`$, and as eigenvalues:

$`\{r_j\}={\displaystyle \frac{1}{(1+w)^2}}\{w(w+3),{\displaystyle \frac{1w}{2}},{\displaystyle \frac{1w}{2}},\mathrm{\hspace{0.17em}0}\},`$ (13)
$`\{s_j\}={\displaystyle \frac{1}{(1+w)^2}}\{2w,\mathrm{\hspace{0.17em}2}w,w(1w)^2,(1w)^3\},`$ (14)
which are ordered according to (6) provided $`1>w>1/2`$. The fact that $`\rho _0\sigma _0`$ and $`\sigma _0\rho _0`$ allows to violate the thermodynamical ordering, since we have simultaneously $`S(\rho _0)<S(\sigma _0)`$ and $`𝒲(\rho _0)<𝒲(\sigma _0)`$ for sufficiently small $`x`$, as seen from $`(1+w)^2[𝒲(\sigma _0)𝒲(\rho _0)]=(1w)(2w1)xw(1+2ww^2)`$. Take $`\mathrm{\Omega }`$ as a two-level system; its only relevant feature will be the eigenvalues $`\{w,1w\}`$ of its initial state $`\omega _0`$. Provided $`1/4<w^2<1/2`$, the equations (11) for $`k=1,2,\mathrm{},8`$ are satisfied and hence $`\rho _0\omega _0\sigma _0\omega _0`$. This implies $`𝒲(\rho _0\omega _0)>𝒲(\sigma _0\omega _0)`$: $`\mathrm{\Omega }`$ restores for the ergotropies the order inferred from the thermodynamic relation $`S(\rho _0)<S(\sigma _0)`$.

All these contradictions between the predictions of thermodynamics and the behavior of finite systems do not occur for the subset of states that may ordered in the sense of majorization: then $`\rho _0\sigma _0`$ implies both $`S(\rho _0)S(\sigma _0)`$ and $`𝒲(\rho _0)𝒲(\sigma _0)`$ and $`𝒲(\rho _0\omega _0)𝒲(\sigma _0\omega _0)`$ for an arbitrary $`\omega _0`$.

Correlations. In the above examples the initial state of S$`+\mathrm{\Omega }`$ was uncorrelated, but not the final state, because evolution permutes its eigenvectors although its eigenvalues remain factorized. Thus, contrary to thermodynamic intuition, the maximum work can be achieved owing to creation of correlations.

Conversely, if the initial state $`Q_0`$ of S$`+\mathrm{\Omega }`$ is correlated, thermodynamics predicts that the available work is increased, due to the subadditivity of entropy, $`S(Q_0)S(\rho _0)+S(\sigma _0)`$, where $`\rho _0=\mathrm{tr}_\mathrm{\Omega }Q_0`$ and $`\sigma _0=\mathrm{tr}_\mathrm{S}Q_0`$ are the marginal states of S and $`\mathrm{\Omega }`$, respectively. However, for finite quantum systems, we have to compare the ergotropies $`𝒲(Q_0)`$ and $`𝒲(\rho _0\sigma _0)`$. Take for S a system with three energy levels $`\epsilon _i`$, $`i=1,2,3`$ and for $`\mathrm{\Omega }`$ a two-level system $`k=1,2`$ with energy levels $`0,\epsilon `$ such that $`0<\epsilon <\epsilon _2\epsilon _1`$, $`\epsilon <\epsilon _3\epsilon _2`$, and for $`\mathrm{S}+\mathrm{\Omega }`$ a stationary initial state $`Q_0`$. Denoting the common eigenstates of $`H+H_\mathrm{\Omega }`$ and $`Q_0`$ as $`|i,k`$, we assume that the eigenvalues $`q_{ik}`$ of $`Q_0`$ are ordered as $`q_{11}>q_{12}>q_{21}>q_{22}>q_{31}>q_{32}`$. Then $`Q_0`$ is both correlated and passive. Suppressing correlations leads to a factorized state $`\rho _0\omega _0`$ with the same energy, and eigenvalues $`r_i=_{k=1}^2q_{ik}`$ for $`\rho _0`$ and $`p_k=_{i=1}^3q_{ik}`$ for $`\omega _0`$. The ordering of the set $`r_ip_k`$ may now differ from that of $`q_{ik}`$; for instance, if $`q_{11}`$ is close to one and all other $`q_{ik}`$’s are small with the same order of magnitude, we have $`r_2p_2<r_3p_1`$, and $`𝒲(\rho _0\omega _0)=(\epsilon _3\epsilon _2\epsilon )(r_3p_1r_2p_2)>0`$. Suppressing the correlations has thus let the system $`\mathrm{S}+\mathrm{\Omega }`$ become active. Altogether, the order carried by correlations, although it manifests itself directly in the entropy, has no relation with increase or decrease of ergotropy.

Fluctuations of work. The work $`W`$ extracted from S by means of external sources acting in a cycle was defined by (2) as the expectation value of the difference between the initial and final energies of S. This expectation value can also be observed as an average satotasaki  by coupling the macroscopic sources of work to $`N1`$ identical systems S operating in parallel, in such a way that their total Hamiltonian is the sum of the Hamiltonians $`H+V(t)`$ of each one. However, for a single system S, the work is expected to fluctuate. In order to evaluate such fluctuations, we need to define work as a random variable for the dynamics described by the optimal unitary transformation (9). To do this, we need to determine, for each single realization in the statistical ensemble described by $`\rho _0`$, the initial and final energy of S. This can be done unambiguously, without perturbing the overall ensemble $`\rho _0`$, only if both $`\rho _0`$ and $`\rho (\tau )`$ commute with $`H`$. Then the measurement of energy at the time $`t=0`$ allows to filter systems which have the same initial energy $`\epsilon _k^{(\mathrm{i})}`$ after the measurement. The evolution operator (9) then determines for each initial $`\epsilon _k^{(\mathrm{i})}`$ the final energy $`\epsilon _l^{(\mathrm{f})}`$, hence the work $`\epsilon _k^{(\mathrm{i})}\epsilon _l^{(\mathrm{f})}`$ is given by S to the sources. We can thus regard $`\epsilon _k^{(\mathrm{i})}\epsilon _l^{(\mathrm{f})}`$ as a classical random variable describing work, and governed by the diagonal elements $`r_j`$ of the density matrix $`\rho _0`$. In the same way as $`𝒲=_{jk}r_jP_{jk}(\epsilon _k\epsilon _j)`$, where $`P_{jk}=|\epsilon _j|U|\epsilon _k|^2`$ is a permutation matrix, we get the ergotropy dispersion:

$`\mathrm{\Delta }𝒲^2={\displaystyle \underset{jk}{}}r_jP_{jk}(\epsilon _k\epsilon _j)^2𝒲^2.`$ (15)
It vanishes for passive states, for which $`P_{jk}=\delta _{jk}`$ and $`𝒲\mathrm{max}W=0`$: one cannot hope to put no work in average, and get some work out by making use of fluctuations. For active states, $`\mathrm{\Delta }𝒲`$ may become sizeable. In the example of a two level-system with the initial state: $`\rho _0=r_1|\epsilon _2\epsilon _2|+r_2|\epsilon _1\epsilon _1|`$, the ergotropy is $`𝒲=(r_1r_2)(\epsilon _2\epsilon _1)`$ and its fluctuation from (15) is $`\mathrm{\Delta }𝒲`$$`=2\sqrt{r_1r_2}(\epsilon _2\epsilon _1)`$. The relative fluctuation $`\mathrm{\Delta }𝒲/𝒲`$ vanishes for $`r_1=1`$, $`r_2=0`$, when $`𝒲`$ is largest; it becomes large when $`𝒲`$ is small, i.e., $`r_1r_2\frac{1}{2}`$. Many systems S should then operate in parallel to produce on average the amount $`𝒲`$ in a predictable fashion.

Conclusion. Maximal work extraction is one of the basic problems of thermodynamics and has applications in various processes of energy conversion landau ; callen ; QL2L . In macroscopic physics, the answer (4) is governed by the non-decrease of entropy. We have shown that finite devices are less efficient in this respect: any evolution of a thermally (but not mechanically) isolated quantum system must leave unchanged not only this entropy, but all the eigenvalues of the density operator, which prevents in most situations the thermodynamic bound from being attainable. We have given a general explicit expression (8) for the “ergotropy”, the quantity which, for finite systems, replaces the free energy: it is the upper bound of the work that a finite system S in a non-equilibrium initial state $`\rho _0`$ may yield if it is coupled to external sources of work undergoing a cyclic transformation. Many interaction Hamiltonians $`V(t)`$ allow to reach this bound.

The proper measure of order for comparing the abilities of work production of finite systems is thus ergotropy and not free energy. Several consequences of this result contradict thermodynamic intuition. Consider, for instance, a state $`\sigma _0`$ of S having the same energy as $`\rho _0`$ and lower entropy. Thermodynamics suggests that more work might be extracted from $`\sigma _0`$ than from $`\rho _0`$, and moreover that the presence of an auxiliary system $`\mathrm{\Omega }`$, in a state $`\omega _0`$ initially uncorrelated with S, preserves this property. Such statements can be violated in finite quantum systems. However, even for finite systems, there is a domain of states defined by the majorization relation, where predictions of thermodynamics are qualitatively correct.

This work is supported by the Dutch Science Foundation FOM/NWO. R. B. acknowledges visits at the UvA.


format: text

text:
# Calculation of electron density of periodic systems using non-orthogonal localised orbitals

## 1 Introduction

Understanding of the electronic structure of extended systems with a local perturbation, e.g. point defects in the crystal bulk  or adsorption of molecules at crystal surfaces  is of fundamental importance in solid state physics and chemistry. One way of calculating the electronic structure of the mentioned systems is based on the usage of so-called cluster methods in which a finite fragment of an extended system (a quantum cluster) is considered in detail while the rest of the system is treated at a lower level of theory . The main problem of any existing cluster based scheme is in choosing an appropriate termination of the cluster. Usually, the quantum cluster is surrounded by point charges , pseudoatoms (see, e.g. ), link atoms  or pseudopotentials . In more sophisticated methods the environment region is described by an electronic wavefunction which could be either frozen  or recalculated self-consistently with that of the quantum cluster region  (a general theory of cluster embedding which comprises most of the existing cluster schemes is considered in ).

A rather general cluster method based on overlapping (not orthogonal) localised orbitals is presently being developed in our laboratory. Our method which is similar in spirit to some one-electron methods  is based on a construction of strongly localised orbitals which are designed to represent the true electronic density of the entire system via a combination of elementary densities associated in simple cases with atoms, ions and/or bonds; these are called *regions* . Our intention is to create a rather general technique which can be valid for systems of different chemical character, ranging from purely ionic to strongly covalent (excluding metals). Therefore, the proper choice of the localisation technique as well as a general method of calculating electron density out of strongly localised non-orthogonal orbitals localised within corresponding regions are crucial for our method to work for a wide range of systems.

The issue of calculating orbitals localised in appropriate regions for extreme cases of strongly ionic and covalent crystals has been considered separately . It is the main objective of this paper to discuss methods of calculating the electron density of periodic systems described via localised non-orthogonal orbitals.

It should be mentioned that the literature on this topic is quite scarce which is probably explained by the lack of interest (until recently) to non-orthogonal (non-canonical) molecular orbitals: in most techniques used in the solid state community orthogonal Bloch functions are employed in practical calculations. There are only a few exceptions (see e.g. ). If a set of non-orthogonal orbitals is used, the expression for the electron density is much more complicated since it contains an inverse of an infinite overlap matrix constructed out of the non-orthogonal orbitals of the whole system under consideration .

As far as we are aware, there have only been two methods developed which enable calculation of the electron density of a periodic system from non-orthogonal orbitals. Both methods are based on a series expansion of the density: while the first method  relies on the so-called cluster expansion of the density, the second one  is based on the power expansion of the inverse overlap matrix. In this paper we analyse only the second of the methods in detail since the first one is very complicated and much more difficult to implement. In section 2 we reexamine the second method from the point of view of the correct density normalisation. Then, in section 3 we suggest an alternative technique which does not require any series expansion. Both methods are compared in section 4 using a very simple model system. The paper is finished with a short discussion and conclusions in section 5.

## 2 Electron density of a periodic system

Let Capital letters $`A`$, $`B`$, etc. be used to indicate regions, while the corresponding small letters $`a`$, $`b`$, etc. - localised orbitals associated with them, i.e. $`aA`$, $`bB`$, etc., see Fig. 1. Each region may have several localised orbitals. We assume that the orbitals are real. They are expanded over atomic orbitals centred only on atoms which are inside the region border. Two localised orbitals belonging to different regions are not orthogonal either because they have common atomic orbitals or, if they do not, then due to their exponential tails.

Each region $`A`$ is prescribed with an even number of $`N_A`$ electrons. Thus, there is a finite number $`n_A=N_A/2`$ of double occupied orbitals associated with the given region $`A`$. The localised orbitals $`\phi _{Aa}(𝐫)`$ belonging to the same region are assumed to be orthonormal; orbitals belonging to different regions are not orthogonal, i.e the corresponding overlap integral $`S_{Aa,Bb}=\phi _{Aa}|\phi _{Bb}`$ is assumed to be not zero in general. Note that there might be several regions within every primitive unit cell in the crystal. Localised orbitals belonging to physically equivalent regions in different unit cells are obtained by appropriate lattice translations, i.e. $`\phi _{Ba}(𝐫)=\phi _{Aa}(𝐫𝐋)`$, where physically equivalent regions $`B`$ and $`A`$ are separated by the translation vector $`𝐋`$.

In general, the spinless electron density takes on the following general form

$$\stackrel{~}{\rho }(𝐫,𝐫^{})=2\underset{Aa}{}\underset{Bb}{}\phi _{Aa}(𝐫)(𝐒^1)_{Aa,Bb}\phi _{Bb}(𝐫^{})$$
(1)
which contains the inverse of the overlap matrix, $`𝐒=S_{Aa,Bb}`$. The overlap matrix can also be written as a set of finite matrix blocks $`𝐒_{AB}=S_{Aa,Bb}`$ associated with every pair of regions. Note that for an infinite crystal the matrix $`𝐒`$ has an infinite size. As usual, the factor of two is due to the fact that each orbital is occupied by two electrons with opposite spins.

In both summations above localised orbitals from all unit cells are taken into account. To stress the periodic symmetry of the crystal, it is useful to rewrite the density in a slightly different form:

$$\stackrel{~}{\rho }(𝐫,𝐫^{})=\underset{𝐋}{}\rho (𝐫𝐋,𝐫^{}𝐋)$$
(2)
where we introduced a *periodic image* of the density (“density image” for short):

$$\rho (𝐫,𝐫^{})=2\underset{Aa}{\overset{}{}}\underset{Bb}{}\phi _{Aa}(𝐫)(𝐒^1)_{Aa,Bb}\phi _{Bb}(𝐫^{})$$
(3)
where in the first sum (indicated by a prime) the summation is run only over localised orbitals within the single primitive cell associated with the zero lattice translation; the other summation runs over all orbitals in the whole infinite system. Note that the density image is normalised on the number of electrons in the unit cell only:

$$\rho (𝐫,𝐫)\text{d}𝐫=\underset{A}{\overset{}{}}N_A$$
(4)
### 2.1 Method based on the expansion of the inverse of the overlap matrix

Following the original prescription by Löwdin , we present the overlap matrix as $`𝐒=\mathrm{𝟏}+𝚫`$, where the matrix $`𝚫=𝚫_{Aa,Bb}`$ is the same as the original overlap matrix except for its elements when $`A=B`$ which are all equal to zero, $`𝚫_{Aa,Aa^{}}=0`$. Then, one can formally write a matrix expansion:

$$𝐒^1=(\mathrm{𝟏}+𝚫)^1=\mathrm{𝟏}𝚫+𝚫^2𝚫^3+\mathrm{}$$
(5)
One can show (using diagonalisation of the matrix $`𝐒`$ or its expansion over the eigenstates) that the expansion (5) can only be used if absolute values of *all* eigenvalues of the matrix $`𝚫`$ are less than unity.

Using expansion of Eq. (5), one obtains the following expansion for the image density (3):

$$\rho (𝐫,𝐫^{})=\underset{n=0}{\overset{\mathrm{}}{}}\rho ^{(n)}(𝐫,𝐫^{})=\underset{n=0}{\overset{\mathrm{}}{}}(1)^n\left[2\underset{Aa}{\overset{}{}}\underset{Bb}{}\phi _{Aa}(𝐫)(𝚫^n)_{Aa,Bb}\phi _{Bb}(𝐫^{})\right]$$
(6)
Note that a general $`n`$-th order term (for $`n2`$) contains additional $`n1`$ summations over all regions due to matrix multiplications in $`𝚫^n`$.

In principle, formulae (2) and (6) allow an approximate calculation of the electron density by terminating the infinite expansion. Care should be taken, however, in doing so in order to preserve the correct normalisation of the density.

The zero order contribution,

$$\rho ^{(0)}(𝐫,𝐫^{})=2\underset{Aa}{\overset{}{}}\phi _{Aa}(𝐫)\phi _{Aa}(𝐫^{})$$
(7)
does not contain overlap integrals at all and is normalised to the total number of electrons in the unit cell. Therefore, if any higher order terms are kept in the terminated expansion (6), they should be integrated (normalised) to zero. Consider this point in more detail.

The first order contribution to the image density,

$$\rho ^{(1)}(𝐫,𝐫^{})=2\underset{Aa}{\overset{}{}}\underset{Bb}{}\phi _{Aa}(𝐫)\mathrm{\Delta }_{Aa,Bb}\phi _{Bb}(𝐫^{})$$
(8)
contains different regions $`AB`$ and thus its contribution to the charge (or normalisation) becomes:

$$\mathrm{\Delta }N^{(1)}=\rho ^{(1)}(𝐫,𝐫)\text{d}𝐫=2\underset{Aa}{\overset{}{}}\underset{Bb}{}\mathrm{\Delta }_{Bb,Aa}\mathrm{\Delta }_{Aa,Bb}=2\underset{A}{\overset{}{}}\text{Tr}_A\left(𝚫^2\right)$$
(9)
where the trace $`\text{Tr}_A(\mathrm{})`$ here is calculated with respect to the localised orbitals belonging to region $`A`$ only. We see that the first order term has a finite nonzero charge (in fact, it is negative).

Any higher order contributions in Eq. (6) for $`n2`$ contains additional summations over regions so that equal regions $`A=B`$ in the double summation there are also possible. Therefore, every such contribution, $`\rho ^{(n)}(𝐫,𝐫^{})`$, will be split into two terms: a diagonal term,

$$\rho _d^{(n)}(𝐫,𝐫^{})=2(1)^n\underset{Aa,a^{}}{\overset{}{}}\phi _{Aa}(𝐫)(𝚫^n)_{Aa,Aa^{}}\phi _{Aa^{}}(𝐫^{})$$
(10)
in which $`A=B`$, and a non-diagonal term,

$$\rho _{nd}^{(n)}(𝐫,𝐫^{})=2(1)^n\underset{Aa}{\overset{}{}}\underset{B(A),b}{}\phi _{Aa}(𝐫)(𝚫^n)_{Aa,Bb}\phi _{Bb}(𝐫^{})$$
(11)
associated with $`AB`$ in Eq. (6). Correspondingly, we obtain the following contributions to the charge:

$$\mathrm{\Delta }N_d^{(n)}=2(1)^n\underset{Aaa^{}}{\overset{}{}}(𝚫^n)_{Aa,Aa^{}}S_{Aa,Aa^{}}=2(1)^n\underset{A}{\overset{}{}}\text{Tr}_A\left(𝚫^n\right)$$
(12)
$$\mathrm{\Delta }N_{nd}^{(n)}=2(1)^n\underset{Aa}{\overset{}{}}\underset{B(A),b}{}(𝚫^n)_{Aa,Bb}\mathrm{\Delta }_{Bb,Aa}=2(1)^n\underset{A}{\overset{}{}}\text{Tr}_A\left(𝚫^{n+1}\right)$$
(13)
Thus, we see that in any order $`n2`$ we have $`\mathrm{\Delta }N_{nd}^{(n)}\mathrm{\Delta }N_d^{(n+1)}`$. This means that the non-diagonal contribution to the density (11) is compensated exactly by the diagonal one (10) of the next order. For instance, the non-zero charge (9) is to be exactly eliminated by a charge due to the diagonal second order density; in turn, a nonzero charge due to non-diagonal second order density is compensated exactly by the diagonal third order density contribution, and so on.

This result is very useful since it allows one to balance properly a terminated expansion for the image density so that it would correspond (in any order!) to the correct total charge. To do this, the final expression for the density of any $`n`$-th order should also include the diagonal ($`A=B`$) term from the contribution of the next order. We stress that this fact was ignored in the previous applications of this method . We obtain, that the correct $`n`$th order expansion for the image density in the notations of Eqs. (10), (11) should have the form:

$$\rho (𝐫,𝐫^{})\rho ^{[n]}(𝐫,𝐫^{})\underset{i=0}{\overset{n}{}}\left[\rho _d^{(i)}(𝐫,𝐫^{})+\rho _{nd}^{(i)}(𝐫,𝐫^{})\right]+\rho _d^{(n+1)}(𝐫,𝐫^{})$$
(14)
By employing this termination of the series, the normalisation condition (4) is satisfied exactly.

Thus, in order to calculate the density up to the $`n`$th order, one has to calculate the matrix elements $`\left(𝚫^k\right)_{Aa,Bb}`$ of the matrix $`\mathrm{\Delta }`$ for all powers $`k=1,\mathrm{},n`$; in addition, one also need diagonal $`A=B`$ elements of $`𝚫^{n+1}`$. Then, the contributions from all density images corresponding to all lattice translations, Eq. (2), are to be added together to get the final electron density.

The method described here relies on the convergence of the density expansion (6). The better localisation of the orbitals $`\phi _{Aa}(𝐫)`$, the faster convergence and thus smaller number of terms is needed. We shall demonstrate in section 3 that in some cases of not very well localised orbitals one has to consider the density expansion up to a very high order which makes the calculation extremely time-consuming. Moreover, if the orbitals localisation becomes worse than a certain criteria (to be also discussed in section 3), then this method fails altogether as the expansion diverges. A general and an extremely efficient technique which is not based on a perturbative expansion of any kind and can be used for localised orbitals of practically any degree of localisation is suggested in the next subsection.

### 2.2 Method based on the Fourier transform of localised orbitals

In Eq. (1) for the electron density, regions $`A`$ and $`B`$ are to be chosen from all unit cells of the infinite periodic system. It is convenient in this section to identify explicitly the lattice vector for every localised orbital in its index. Therefore, in the following we shall use letters $`A`$, $`B`$, etc. only for regions within the zeroth unit cell; in particular, the orbital $`\phi _{Aa}(𝐫)`$ is assumed to be from the zeroth cell. Localised orbitals from other cells are characterised by the combined index $`(𝐋Aa)`$, i.e. $`\phi _{𝐋Aa}(𝐫)=\phi _{Aa}(𝐫𝐋)`$ is the $`a`$-th localised orbital from region $`A`$ in the unit cell separated from the zeroth cell by the lattice translation $`𝐋`$.

Correspondingly, Eq. (1) is rewritten in the following way:

$$\stackrel{~}{\rho }(𝐫,𝐫^{})=2\underset{𝐋Aa}{}\underset{𝐌Bb}{}\phi _{Aa}(𝐫𝐋)(𝐒^1)_{𝐋Aa,𝐌Bb}\phi _{Bb}(𝐫^{}𝐌)$$
(15)
where $`𝐋`$ and $`𝐌`$ are two lattice translations. A further transformation is possible here since the overlap integral $`S_{𝐋Aa,𝐌Bb}`$ depends in fact only on the difference $`𝐌𝐋`$ of the lattice translations. This allows expansion of the overlap integral into the Fourier integral

$$S_{𝐋Aa,𝐌Bb}=\frac{1}{N}\underset{𝐤}{}S_{Aa,Bb}(𝐤)e^{\text{i}𝐤(𝐋𝐌)}$$
(16)
where the summation is performed over $`N`$ points $`𝐤`$ in the first Brillouin zone (BZ) and

$$S_{Aa,Bb}(𝐤)=\underset{𝐋}{}S_{\mathrm{𝟎}Aa,𝐋Bb}e^{\text{i}\mathrm{𝐤𝐋}}$$
(17)
is the corresponding Fourier image. The direct lattice summation in the last formula is easily terminated due to (usually) exponential decay of the overlap integrals between localised orbitals.

Using the Fourier representation of the overlap matrix, one can exactly calculate its inverse as follows:

$$\left(𝐒^1\right)_{𝐋Aa,𝐌Bb}=\frac{1}{N}\underset{𝐤}{}\left[𝐒^1(𝐤)\right]_{Aa,Bb}e^{\text{i}𝐤(𝐋𝐌)}$$
(18)
Note that the matrix $`𝐒(𝐤)=S_{Aa,Bb}(𝐤)`$ has a finite size of the number of localised orbitals per unit cell. Therefore, in order to calculate the inverse of the overlap matrix in direct space, one has to perform the calculation of $`𝐒^1(𝐤)`$ for finite size matrices for every $`𝐤`$ point necessary to sample the BZ. Substituting Eq. (18) into Eq. (15), we arrive at the following final expression for the electron density:

$$\stackrel{~}{\rho }(𝐫,𝐫^{})=\frac{2}{N}\underset{𝐤}{}\left\{\underset{Aa}{}\underset{Bb}{}\phi _{Aa}(𝐫,𝐤)\left[𝐒^1(𝐤)\right]_{Aa,Bb}\phi _{Bb}^{}(𝐫^{},𝐤)\right\}$$
(19)
where

$$\phi _{Aa}(𝐫,𝐤)=\underset{𝐋}{}\phi _{Aa}(𝐫𝐋)e^{\text{i}\mathrm{𝐤𝐋}}$$
(20)
is the Fourier expansion of the localised orbital. Due to exponential decay of the localised orbitals, the summation over lattice vectors $`𝐋`$ in the last expression is in fact finite.

The obtained formula for the density is exact. In particular, it contains the periodicity of the lattice built in. It is also extremely convenient for numerical implementation. Indeed, what is needed is the calculation of the Fourier images, according to Eq. (20), of every localised orbital in the primitive unit cell for every $`𝐤`$ point. The summations in the curly brackets in Eq. (19) are finite (limited to the orbitals within the zeroth cell only) and are thus easily performed. The extend to which the orbitals $`\phi _{Aa}(𝐫)`$ are localised is reflected by the number of cells to be taken into account while performing the lattice summations in Eqs. (17) and (20). Even for orbitals which are not very well localised, the amount of work needed to perform these lattice summations is incomparable with the cost of the first method (section 2.1) which requires including more terms in the perturbation expansion if the localisation is not good enough.

## 3 Results

Atomic units are used throughout this section. The application of the two methods considered in the previous sections is illustrated here on a simple cubic lattice model containing a single region in every unit cell. The lattice constant $`a`$ will be assumed to be equal to 1 a.u. for simplicity. Each region is represented by a single localised orbital in a form of a normalised $`s`$ type Gaussian

$$\phi _{𝐋Aa}(𝐫)\phi _𝐋(𝐫)\phi (𝐫𝐋),\phi (r)=\left(\frac{2\alpha }{\pi }\right)^{3/4}e^{\alpha 𝐫^2}.$$
(21)
By choosing various values for the exponent $`\alpha `$, one can vary the degree of localisation of the orbitals. Indeed, the size of the *spatial extent* of the orbital can be measured in terms of $`r_{eff}=\sqrt{\frac{\mathrm{ln}10}{\alpha }}1.52\alpha ^{1/2}`$, which corresponds to $`e^{\alpha r_{eff}^2}=`$0.1. We found this approach more convenient in our particular case than the application of the existing localisation criteria (see, e.g. ).

For this model system it is possible to do some preliminary analytical estimates of the convergence of the series (5). We know from section 2.1 that the series will converge if *all* eigenvalues $`\mathrm{\Delta }_\lambda `$ of the matrix $`𝚫=𝐒\mathrm{𝟏}`$ are between -1 and 1. It is easy to notice that the eigenvalues are in fact given by the Fourier transforms $`\mathrm{\Delta }_𝐤`$ of the matrix $`𝚫`$ which is introduced much in the same way as $`𝐒(𝐤)`$ in Eq. (17). Indeed, because $`\mathrm{\Delta }_{𝐋,𝐌}=\mathrm{\Delta }_{\mathrm{𝟎},𝐌𝐋}`$, one can write:

$$\underset{𝐌}{}\mathrm{\Delta }_{𝐋,𝐌}e^{\text{i}\mathrm{𝐤𝐌}}=\left(\underset{𝐌}{}\mathrm{\Delta }_{\mathrm{𝟎},𝐌𝐋}e^{\text{i}𝐤(𝐌𝐋)}\right)e^{\text{i}\mathrm{𝐤𝐋}}=\mathrm{\Delta }_𝐤e^{\text{i}\mathrm{𝐤𝐋}}$$
(22)
This is nothing but the eigenproblem for the matrix $`𝚫`$ with $`\mathrm{\Delta }_𝐤`$ being its eigenvalues (numbered by vectors $`𝐤`$ from the BZ) and $`e^{\text{i}\mathrm{𝐤𝐋}}`$ \- eigenvectors. Therefore, the convergence criteria for the series (5) reduces to the inequalities $`\left|\mathrm{\Delta }_𝐤\right|<1`$ which should be valid for *any* $`𝐤`$. Taking into account the overlap only between nearest neighbours, we obtain:

$$\left|\mathrm{\Delta }_𝐤\right|=\left|2\delta (\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)+\mathrm{cos}(k_za))\right|6\delta <1$$

with the overlap between neighbouring orbitals being $`\delta =e^{\alpha a^2/2}`$. This results in the following criterion for the convergence of the Löwdin expansion (for $`a=`$1 a.u.):

$$\alpha \alpha _1^{}=2\mathrm{l}\mathrm{n}(6)3.6$$
(23)
Similar analysis which takes into account the next nearest neighbours gives a very similar estimate of $`\alpha _1^{}4.05`$. These estimates correspond to the maximum spatial extent of the orbitals (21) of the order of $`r_{eff}`$0.76 a.u., i.e. there is very small overlap between neighbouring orbitals which, we recall, are separated by 1 a.u. in the lattice.

The other method based on the Fourier transform of the orbitals has also its limits which are hidden in the formulae (17) and (20): if a certain cut-off $`\left|𝐋\right|r_c`$ for the direct lattice summation $`𝐋`$ is assumed in the calculation of $`S_{Aa,Bb}(𝐤)`$ and $`\phi _{Aa}(𝐫,𝐤)`$, then there will be some limitations on the allowed degree of localisation of the orbitals. The required criterion can be worked out e.g. by analysing the Fourier transform (20) of the orbital at its maximum in the centre of the BZ (i.e. of $`\phi (𝐫=\mathrm{𝟎},𝐤=\mathrm{𝟎})`$) as follows:

$$\underset{|𝐋|>r_c}{}\phi (𝐋)\underset{|𝐋|<r_c}{}\phi (𝐋)$$

Replacing the sums by the corresponding volume integrals, we obtain the following criterion:

$$xe^{x^2}+\frac{\sqrt{\pi }}{2}\text{erfc}(x)\frac{\sqrt{\pi }}{4}$$

where $`x=r_c\sqrt{\alpha }`$. The inequality above is satisfied if $`x2`$, i.e. $`\alpha 4/r_c^2`$. Assuming that $`r_c`$ is equal to 4$`÷`$5 lattice constants, we obtain the necessary condition for the exponent of the localised orbitals,

$$\alpha \alpha _2^{}0.2$$
(24)
for which our Fourier transform method should work. The obtained critical value of $`\alpha _2^{}`$ results in the maximum spatial extent of the orbitals of the order of $`r_{eff}`$3.4 a.u. which corresponds to very diffuse orbitals spreading over more than six unit cells.

Similar criteria is obtained for the ovelrap integrals as well. Thus, the method we suggest should have a much wider range of applicability than the Löwdin method as far as the degree of localisation of the non-orthogonal orbitals is concerned since $`\alpha _1^{}\alpha _2^{}`$. This conclusion is also supported by our numerical calculations which we now describe.

Numerical calculations of the necessary powers of the $`𝚫`$ matrix needed for the Löwdin method were done in the following way. Since the density is calculated in the same point $`𝐫=𝐫^{}`$ in Eqs. (10), (11) and (14), the regions $`A`$ and $`B`$ in these equations are either the same or not far away from each other. Therefore, to calculate $`\left(𝚫^n\right)_{Aa,Bb}`$ one can simply choose a sufficiently big finite cluster of atoms (in fact, the cluster radius should be at least of the order of $`\frac{n}{2}r_c^{}`$, where $`r_c^{}`$ is the decay length of the overlap integral) with regions $`A`$ and $`B`$ somewhere in its centre and then calculate the complete ovelrap matrix for it, $`\stackrel{~}{𝚫}`$. Then, by performing the necessary $`n1`$ matrix multiplications, one can calculate $`\left(𝚫^n\right)_{Aa,Bb}`$ as $`\left(\stackrel{~}{𝚫}^n\right)_{Aa,Bb}`$.

When using the Fourier transform method, we employed the Monkhorst-Pack (MP) method  for the $`𝐤`$ point sampling and the same cut-off distance for the direct lattice summations in Eqs. (17) and (20) as in the previous method. In all our calculations we used the 4x4x4 MP set which was found to be sufficient in all cases.

Results of our calculations for a large value of the exponent ($`\alpha \alpha _1^{}\alpha _2^{}`$) are shown in Fig.2. This case corresponds to strongly localised orbitals as is the case in ionic systems such as MgO and NaCl. Overlap between orbitals is negligible and even zero order approximation Löwdin method, Eq. (7), was found sufficient to give the correct density. Density curves for both methods are indistinguishable from each other.

The calculated densities in the intermediate case ($`\alpha \alpha _1^{}`$) are shown in Fig. 3. This value of $`\alpha `$ may correspond to ion-covalent and covalent systems. One can see that high order approximations (up to $`n=`$8) of the Löwdin method, Eq. (14), are needed here to converge the density and thus the calculation is quite time consuming.

Finally, we show in Fig. 4 the densities calculated using both methods for orbitals which are least localised when $`\alpha _1^{}>\alpha >\alpha _2^{}`$. The density, obtained using the Fourier transform method, Eq. (19), is spread almost uniformly in the crystal volume and thus may correspond to a metallic band. At the same time, the Löwdin expansion method, Eq. (14), does not converge at all and the density is clearly diverges.

One can expect that the latter situation can happen only for metallic systems. Interestingly, our calculations (not reported here) for such a realistic covalent system as a crystalline Si show that the Löwdin approach also fails in some cases when the orbitals are not sufficiently localised. Note that various degree of localisation of the orbitals can be obtained using different localisation techniques and different choice of regions, see  for more details.

## 4 Conclusions

In summary, we have considered two numerical methods which allow calculation of the electron density of a 3D periodic system constructed via a set of non-orthogonal localised orbitals. The first, so-called Löwdin, method based on the power expansion of the inverse of the overlap matrix has been found to be efficient only for strongly localised orbitals. For an intermediate degree of orbitals localisation this method has been found to be quite computationally demanding since many terms in the series are to be retained. However, if orbitals are not sufficiently localised (the exact criterion has also been suggested), the method fails altogether and the power expansion has been shown to be divergent.

Then, we have suggested another method based on the Fourier transform of the localised orbitals which involves calculations of inverse of only finite matrices and a $`𝐤`$ point summation over the Brillouin zone. This method is computationally much less demanding and does not have any convergency problems. Using a simple model for the crystal electron density represented via a set of Gaussian $`s`$ type orbitals in a simple cubic lattice (one orbital per unit cell), we have shown that our method works equally well within a rather wide range of orbitals having different localisation, whereas the first method fails for a relatively weakly localised orbitals. The application of the Fourier transform method to realistic systems such as MgO and Si perfect crystals is published elsewhere .

## Acknowledgements

We are extremely grateful to I. V. Abarenkov and I. Tupitsin for useful and stimulating discussions. O.D. would also like to acknowledge the financial support from the Leverhulme Trust (grant F/07134/S) which made this work possible.


text:
# Collisional oscillations of trapped boson-fermion mixtures approaching collapse

## I Introduction

Boson-fermion mixtures have been produced with atomic gases in several experiments DeMarco and Jin (1999); Granade et al. (2002); Truscott et al. (2001); Schreck et al. (2001); Roati et al. (2002); Hadzibabic et al. (2002) and the deep quantum-degeneracy regime has been reached for both the bosonic and the fermionic component Strecker et al. (2003); Hadzibabic et al. (2003). The aim of fermion cooling in these mixtures has been the realization of a superfluid state Stoof et al. (1996); Ohashi and Griffin (2002). Even though the mixtures are very dilute, the interactions between bosons and fermions play important roles and in particular the $`s`$-wave collisions between the two species are exploited in the evaporative cooling process to obtain rethermalization. Mean-field interactions affect the expansion rate of the cloud and its thermodynamic properties Roati et al. (2002); Hu and Liu (2003).

For a 87Rb-40K mixture the boson-fermion interactions are large and attractive Ferrari et al. (2002), while in other systems they could be tuned to become attractive by exploiting Feshbach resonances Feshbach (1962); Stwalley (1976); Tiesinga et al. (1993). Boson-fermion attractions enhance the overlap between the two species, favoring in principle a boson-mediated fermion pairing Heiselberg et al. (2000); Bijlsma et al. (2000), and ultimately lead to the collapse of the cloud as the boson-fermion coupling strength overcomes the Fermi pressure. Collapse has been observed experimentally in the 87Rb-40K mixture as a sudden disappearance of the fermion cloud when the number of bosons, and hence the coupling between the two species, is increased Modugno et al. (2002). The equilibrium properties and the phase diagram of a mixture with attractive interactions have been investigated by several authors Miyakawa et al. (2001); Roth and Feldmeier (2002); Roth (2002); Liu et al. (2003), finding good agreement between the experimental data and a mean-field theoretical description Modugno et al. (2003); Capuzzi et al. (2003a).

The study of the collective modes of a boson-fermion mixture with repulsive interactions has proved to be a good indicator of the quantum transition to a spatially demixed state Capuzzi et al. (2003b). For the case of attractive interactions some analyses using semi-analytical methods Miyakawa et al. (2000); Liu and Hu (2003) have suggested that the lowest-lying monopolar mode should show a softening in the approach to collapse. A similar frequency softening is found for a pure condensate with increasingly large boson attractions Dodd et al. (1996); Pérez-García et al. (1996). In contrast, the numerical solution of the equations of motion for boson-fermion mixtures in the collisionless Random-Phase-Approximation (RPA) does not give indications of frequency softening for the choice of parameters adopted in the calculations Sogo et al. (2002). Moreover, it is found that when the number of bosons is much larger than that of fermions some families of modes are strongly blue-shifted as a consequence of the increase of the particle densities in the cloud approaching collapse Capuzzi et al. (2003c). However, the RPA calculations could not be brought to the very proximity of collapse because of very severe numerical difficulties.

In this paper we study the collective modes of a trapped boson-fermion mixture with attractive interactions using the equations of generalized hydrodynamics. These allow us to span the entire range of boson-fermion interaction strength up to the collapse point. Although the hydrodynamic equations are strictly valid only in the collisional regime, their predictions also serve as a guideline to interpret the RPA spectra. By evaluating the frequencies of the first ten low-lying modes we find that some modes undergo a blue shift which can become quite large at large numbers of bosons, while some other modes display a softening which becomes pronounced in a very narrow window in the proximity to collapse.

## II Equations of generalized hydrodynamics

We consider a model for a dilute mixture of two species of alkali atoms inside harmonic traps, one species being spin-polarized fermions of mass $`m_F`$ and the other bosons of mass $`m_B`$. The boson-boson and boson-fermion interactions are described by contact potentials involving the coupling constants $`g_{BB}=4\pi \mathrm{}^2a_{BB}/m_B`$ and $`g_{BF}=2\pi \mathrm{}^2a_{BF}/m_r`$, where $`a_{BB}`$ and $`a_{BF}`$ are the corresponding $`s`$-wave scattering lengths and $`m_r=m_Bm_F/(m_F+m_B)`$. The fermion-fermion interactions are negligible, since the Pauli principle forbids collisions in the $`s`$-wave channel between fermions of the same spin. In the following we treat the case of boson-fermion attractions ($`a_{BF}<0`$) and boson-boson repulsions ($`a_{BB}>0`$).

The dynamics of the mixture in the collisional regime is described by the equations of generalized hydrodynamics as in Ref. Capuzzi et al. (2003b). These are (i) the equations for particle conservation,

$$_t\rho _\sigma +𝐣_\sigma =0,$$
(1)
relating the particle density $`\rho _\sigma (𝐫,t)`$ to the current density $`𝐣_\sigma (𝐫,t)`$ for each component ($`\sigma =B,F`$); and (ii) the quantum Navier-Stokes equations for momentum conservation, which are

$$m_\sigma _t𝐣_\sigma =\rho _\sigma (𝐅_\sigma V_\sigma g_{BF}\rho _{\overline{\sigma }})$$
(2)
to linear order in the velocity fields. Here $`\overline{\sigma }`$ denotes the component different from $`\sigma `$, $`V_\sigma (𝐫)=m_\sigma \omega _\sigma ^2r^2/2`$ are isotropic trapping potentials, and $`𝐅_\sigma (𝐫,t)`$ are self-consistent internal forces. We take for the latter the expressions

$$𝐅_B=\left[g_{BB}\rho _B\frac{\mathrm{}^2}{2m_B}\frac{^2\sqrt{\rho _B}}{\sqrt{\rho _B}}\right]$$
(3)
and

$$𝐅_F=\left[A\rho _F^{2/3}\frac{\mathrm{}^2}{6m_F}\frac{^2\sqrt{\rho _F}}{\sqrt{\rho _F}}\right],$$
(4)
with $`A=\mathrm{}^2(6\pi ^2)^{2/3}/2m_F`$. The equations of generalized hydrodynamics are thereby closed by means of a local-density approximation on the kinetic stress tensors, which also includes surface kinetic contributions through the last term on the rhs of Eqs. (3) and (4) Capuzzi et al. (2003b). These surface correction terms transcend the standard Thomas-Fermi approximation and for bosons the present approach is in fact equivalent to the full time-dependent Gross-Pitaevskii equation.

The equilibrium state of the mixture is given by the stationary solutions of Eq. (2). These also correspond to a minimum of the energy functional

$`E[\rho _B,\rho _F]`$ $`=`$ $`{\displaystyle d^3r\left(V_B\rho _B+\frac{g_{BB}}{2}\rho _B^2+\xi _B\right)}`$ (5)
$`+{\displaystyle d^3r\left(V_F\rho _F+\frac{3}{5}A\rho _F^{5/3}+\xi _F\right)}`$
$`+g_{BF}{\displaystyle d^3r\rho _F\rho _B},`$
where the surface kinetic-energy terms are $`\xi _B=\mathrm{}^2|\sqrt{\rho _B}|^2/2m_B`$ and $`\xi _F=\mathrm{}^2|\sqrt{\rho _F}|^2/6m_F`$.

At increasingly large boson-fermion attractions the densities of the two species increase in their overlap region and collapse occurs when these attractions overcome the Fermi kinetic pressure and the boson-boson repulsions. Collapse is identified in the numerical search for the equilibrium state as the point where it is no longer possible to find a stable minimum for the energy functional in Eq. (5). We have verified that the location of collapse found in this way is well approximated by the estimate Mølmer (1998); Miyakawa et al. (2001)

$$|a_{BF}^c|=\left(\frac{a_{BB}}{\alpha k_F}\right)^{1/2},$$
(6)
where $`k_F=(6\pi ^2\rho _F(0))^{1/3}`$ is the Fermi wave number estimated from the fermion density at the center of the trap and $`\alpha =[3^{1/3}/(2\pi )^{2/3}](m_F+m_B)^2/(4m_Fm_B)`$.

In the following we characterize the approach of the mixture to the collapse instability by following the behavior of its low-lying collective oscillation frequencies. The normal modes of Eqs. (1) and (2) are found by expanding $`\rho _\sigma (𝐫,t)`$ around the equilibrium state $`\rho _\sigma ^0(r)`$ as $`\rho _\sigma (𝐫,t)=\rho _\sigma ^0(r)+\delta \rho _\sigma (𝐫)e^{i\omega t}`$, linearizing Eqs. (1) and (2), and Fourier transforming with respect to the time variable. This yields the coupled eigenvalue equations

$$m_\sigma \omega ^2\delta \rho _\sigma =(\rho _\sigma ^0\delta 𝐅_\sigma )g_{BF}(\rho _\sigma ^0\delta \rho _{\overline{\sigma }}),$$
(7)
where $`\delta 𝐅_\sigma `$ are the forces obtained as $`\delta 𝐅_\sigma =𝐅_\sigma [\rho _\sigma ^0+\delta \rho _\sigma ]𝐅_\sigma [\rho _\sigma ^0]`$ to linear order in the density fluctuations (for their expressions see Ref. Capuzzi et al. (2003b)). We also determine the solutions of the eigenvalue equations in the case where the dynamical coupling between the two components is neglected. Hereafter we shall call these solutions the uncoupled boson and fermion modes and calculate them by solving Eq. (7) for $`\sigma =B`$ and $`F`$ with the second term in the rhs set to zero. This calculation will help us in attributing bosonic or fermionic character to the coupled modes of the mixture by comparison with the frequencies of the uncoupled modes. Even though the bosonic and fermionic oscillations that are obtained in this way are dynamically uncoupled, the effect of the mutual boson-fermion interaction is included in their calculation through the use of the density profiles of the mixture at equilibrium.

The numerical procedure that we have used can be summarized as follows: (i) we find the equilibrium density profiles by minimizing the energy functional (5) using the steepest-descent method; (ii) we project Eq. (7) into subspaces of different angular momentum $`l`$ by factorizing the density fluctuations as $`\delta \rho _\sigma (𝐫)=\delta \rho _\sigma ^l(r)Y_{lm}(\widehat{r})`$ where $`Y_{lm}(\widehat{r})`$ are the spherical harmonics; and (iii) we solve the coupled equations in a given $`l`$-subspace by means of standard linear-algebra routines Anderson et al. (1999).

## III Collective modes

We have solved Eq. (7) for two sets of experimentally relevant system parameters, corresponding to a 7Li-6Li mixture Schreck et al. (2001) and to a 87Rb-40K mixture Modugno et al. (2002). In both cases we have set the frequency of the isotropic trap at the geometric average of the frequencies in the experimental setup and have adopted the experimental value of the boson-boson scattering length. In our calculations we follow two different routes to reach collapse. We drive the 7Li-6Li mixture to collapse by varying the value of $`a_{BF}`$ from positive to strongly negative, as could be achieved by means of a Feshbach resonance. In the 87Rb-40K mixture instead we keep $`a_{BF}`$ fixed and increase the number $`N_B`$ of bosons, as is done in the experiments of Modugno et al. Modugno et al. (2002).

### III.1 The 7Li-6Li mixture

Figure 1 shows the frequencies of the low-lying monopole ($`l=0`$) modes in the 7Li-6Li mixture as functions of $`a_{BF}`$, with the choice of parameters $`\omega _F=\omega _B=2\pi \times 1000\mathrm{s}^1`$, $`a_{BB}=0.27`$nm, and particle numbers $`N_F=10^4`$ and $`N_B=10^6`$. Pronounced blue shifts are observed at intermediate values of the boson-fermion scattering length both in a set of eigenfrequencies of Eq. (7) (circles) and in the frequencies of the uncoupled fermionic modes (dashed lines). This effect is a consequence of the increase of fermion density in the central part of the trap as the strength of the coupling is increased. For our choice of parameters the bosonic density profile is almost unaffected in this range of boson-fermion coupling and essentially acts as an effective attractive well for the fermions. The proximity of these modes to the uncoupled fermionic modes indicates that the effect of the dynamical coupling is negligible, thus confirming the interpretation in terms of static mean-field effects. At larger values of $`|a_{BF}|`$ the other eigenfrequencies of Eq. (7) (circles) show an increasing departure from the uncoupled bosonic modes (solid lines) and tend to rapidly soften as collapse is approached. We conclude that, in contrast to the blue-shift of the fermionic modes, the softening of “bosonic” modes is a truly dynamical signature of the impending collapse.

Similar trends of the mode frequencies are also found for the dipole ($`l=1`$) and quadrupole ($`l=2`$) oscillations, except that the lowest ($`n=0`$) mode is in both cases essentially unaffected by the interactions. For the dipole mode this behavior is an exact consequence of the generalized Kohn theorem Kohn (1961); Dobson (1994), while for modes of higher $`l`$ it is strictly valid only in the Thomas-Fermi limit. This can be explicitly verified from Eq. (7) by making the Ansatz $`\delta \rho _\sigma (𝐫)r^lY_{lm}(\widehat{r})`$ for the density fluctuations and by exploiting the property $`^2\delta \rho _\sigma =0`$ for these modes.

### III.2 The 87Rb-40K mixture

We evaluate the collective modes of a 87Rb-40K mixture with scattering lengths $`a_{BB}=5.5`$nm and $`a_{BF}=21.7`$nm. Although both components are inside the same magnetic trap, their trapping frequencies differ considerably as a consequence of the large difference in atomic masses. In the numerical calculations we take $`\omega _B=2\pi \times 90.9\mathrm{s}^1`$ and $`\omega _F=2\pi \times 134\mathrm{s}^1`$. Collapse is approached by varying $`N_B`$ from $`4\times 10^3`$ to approximately $`8\times 10^4`$, while keeping $`N_F=2\times 10^4`$ fixed.

In Fig. 2 we show the frequencies of the low-lying monopolar modes for this mixture. At increasing values of the boson number we find that the frequencies of the bosonic modes, identified as those that are closer to the uncoupled bosonic fluctuations, show a monotonic decrease which becomes very pronounced very close to collapse. On the other hand, the fermionic modes only show a modest blue shift. This is not as large as in our results for the 7Li-6Li mixture, since the numbers of particles of the two species are here more similar.

## IV Spectral functions

We have also calculated the whole spectral functions under the effect of the perturbing fields $`\delta U_\sigma (𝐫)e^{i\omega t}`$, in order to make a comparison of the collisional spectra with those of the mixtures in the RPA collisionless regime. Peaks in the spectral functions indicate the location of collective modes, while their height measures the probability of a transition from the equilibrium state to an excited state with energy $`\mathrm{}\omega `$.

In the collisional regime the spectral peaks are located in correspondence of the eigenfrequencies $`\omega _i`$ of Eq. (7), while the strength of a transition at a frequency $`\omega `$ is estimated from the corresponding density fluctuations $`\delta \rho _\sigma ^i(𝐫)`$ as

$$S_{\sigma \sigma ^{}}^{\text{coll}}(\omega )=\frac{1}{\pi }\mathrm{Im}\underset{i}{}\frac{𝒵_\sigma ^{}^i}{\omega _i^2\omega ^2}d^3r\delta U_\sigma ^{}(𝐫)\delta \rho _\sigma ^{}^i(𝐫)$$
(8)
(see Appendix). Here the sum is performed over the whole spectrum of Eq. (7) and $`𝒵_\sigma ^i`$ is given by

$$𝒵_F^i=\frac{2A}{3m_F}d^3r(\rho _F^0)^{1/3}\delta \rho _F^i\left(\rho _F^0\delta U_F\right)$$
(9)
and

$$𝒵_B^i=\frac{g_{BB}}{m_B}d^3r\delta \rho _B^i\left(\rho _B^0\delta U_B\right).$$
(10)
Since we are using in Eq. (8) the equilibrium profiles and the mode frequencies of the fully coupled dynamical equations, we expect it to yield the spectral functions quite accurately.

More generally, the spectral functions are related to the density-density responses $`\chi _{\sigma \sigma ^{}}(𝐫,𝐫^{},\omega )`$ by the fluctuation-dissipation theorem,

$$S_{\sigma \sigma ^{}}(\omega )=\frac{1}{\pi }\mathrm{Im}d^3rd^3r^{}\delta U_\sigma ^{}(𝐫)\chi _{\sigma \sigma ^{}}(𝐫,𝐫^{},\omega )\delta U_\sigma ^{}(𝐫^{}).$$
(11)
We evaluate the response functions in the collisionless regime within the RPA, following the approach described in Refs. Capuzzi and Hernández (2001); Capuzzi et al. (2003c). For a monopolar drive ($`\delta U_\sigma (𝐫)r^2`$) we calculate both the fermionic response $`S_{FF}^{l=0}(\omega )`$ and the bosonic response $`S_{BB}^{l=0}(\omega )`$, corresponding to excitations induced by applying the drive separately on the fermions or on the bosons.

In Figs. 3 and 4 we report the spectral functions of the 7Li-6Li mixture as obtained from the RPA formalism (top panels) and from Eq. (8) for the collisional regime (bottom panels) at various values of $`a_{BF}`$ and with the same choice of parameters as in Fig. 1. In the absence of boson-fermion coupling a fermionic perturbing field excites the pure $`n=0`$ mode of the fermions, which gives the only contribution to the FF spectra. In contrast, the BB spectra contain in the same limit several weak peaks at higher frequencies above the pure $`n=0`$ mode of the bosons, owing to the boson-boson interactions. For $`a_{BF}0`$ the RPA spectra are very complex, from the lifting of degeneracy in the discrete-level structure which is due to the interactions with the bosons. The fermionic peaks can be recognized as broad fragmented contributions of largest height in the fermionic spectral functions, while the bosonic peaks are usually isolated and are most marked in the bosonic spectra Capuzzi et al. (2003c).

Quite interestingly, the spectra in the two regimes display similar gross features. At small values of $`|a_{BF}|`$ the peaks are located near the frequencies of the higher-order monopolar modes, but at larger values they tend to spread over the whole range of frequency. This phenomenon occurs for the same choice of parameters for which the blue shift of the collisional fermionic modes becomes pronounced and several mode crossings occur. For $`a_{BF}=3`$ nm the first fermionic peak is strongly blue-shifted while the low-frequency part of the spectrum is dominated by bosonic modes having considerable oscillator strength in both bosonic and fermionic spectral functions. This indicates a strong dynamical coupling of the density fluctuations of the two species in a regime of parameters where the uncoupled frequencies are close to the coupled solutions for the collisional fermionic modes (see again Fig. 1). A careful analysis also shows that for $`a_{BF}=3`$ nm the frequencies of the bosonic modes in the RPA spectra are slightly red-shifted, showing the first indications of a trend towards softening on the approach to collapse.

Finally, in Fig. 5 we compare the collisional results for the monopolar spectral function (bottom panel) with those of the RPA (top panel) for a 87Rb-40K mixture with the same choice of parameters as in Fig. 2. The RPA spectra show a fragmented fermionic contribution around the frequency $`2\omega _F`$ of the bare monopolar oscillation as well as two main bosonic peaks at its sides. The frequency of these two peaks tends to decrease with increasing $`N_B`$, similarly to what is found for the same modes in the collisional regime. This suggests that such frequency softening is a signature of the incipient collapse.

## V Summary and concluding remarks

In summary, we have used the equations of generalized hydrodynamics to study the collective modes of trapped boson-fermion mixtures with mutual attractive interactions driving the mixture towards the collapse instability. We have focused on two specific systems of experimental interest, for which we have chosen two different routes to collapse.

In both cases we find that the collective spectra show a frequency softening of a family of modes of bosonic nature as a signature of the incipient collapse. This softening becomes most pronounced in a very narrow region of parameters near collapse. A second effect of increasing the attractive coupling between the two species is a blue shift of the modes of fermionic origin, which becomes more evident for large numbers of bosons and reflects the compression of the fermionic cloud from the interactions with the bosons.

A comparison of the hydrodynamic spectra with the spectra calculated in the collisionless regime within the random-phase approximation suggests that both the blue shift of the “fermionic” modes and the softening of the “bosonic” modes are general features of the dynamics of the mixture. By comparing the height of the peaks in the bosonic and fermionic spectral functions we also conclude that the two families of modes can be most efficiently excited by applying the driving fields separately on the species.

Our analysis has been restricted to the linear regime, neglecting nonlinear and beyond-mean-field effects which may start to play a role close to collapse. We hope to address these issues in future work.

###### Acknowledgements.

This work was partially supported by INFM through the PRA-Photonmatter Program.

*

## Appendix A Spectral functions in the collisional regime

The spectral functions of a mixture in the collisional regime can be obtained by adding the perturbing fields $`\delta U_\sigma (𝐫)e^{i\omega t}`$ to the trapping potential $`V_\sigma (𝐫)`$ in Eqs. (1) and (2) and then linearizing the equations of motion in terms of the induced density fluctuation $`\delta \rho _\sigma (𝐫,\omega )`$. This leads in Fourier transform with respect to time to the equation

$$m_\sigma \omega ^2\delta \rho _\sigma =(\rho _\sigma ^0\delta 𝐅_\sigma )g_{BF}(\rho _\sigma ^0\delta \rho _{\overline{\sigma }})(\rho _\sigma ^0\delta U_\sigma )$$
(12)
for each component of the mixture. Here, at variance from Eq. (7) $`\omega `$ is fixed by the frequency of the drive. Equation (12) can be solved by expanding $`\delta \rho _\sigma (𝐫,\omega )`$ in the basis of the eigenvectors $`\delta \rho _\sigma ^i(𝐫)`$ of Eq. (7) as

$$\delta \rho _\sigma (𝐫,\omega )=\underset{i}{}C^i(\omega )\delta \rho _\sigma ^i(𝐫).$$
(13)
In turn, the expansion coefficients $`C^i`$ are obtained by projecting Eq. (12) into the subspace generated by a given fluctuation $`\delta \rho _\sigma ^j(𝐫)`$. This procedure requires defining a scalar product between two eigenmodes of Eq. (7), which we shall denote by the braket $`\delta \rho ^i|\delta \rho ^j`$. The projection yields

$$C^i(\omega )=\frac{\delta \rho ^i|\left(\rho _\sigma ^0\delta U_\sigma \right)}{\omega _i^2\omega ^2}$$
(14)
where $`\omega _i`$ is the frequency corresponding to the eigenvector $`\delta \rho _\sigma ^i(𝐫)`$.

A generic scalar product can be written as

$$\delta \rho ^i|\delta \rho ^j=\underset{\sigma ,\sigma ^{}}{}d^3r\delta \rho _\sigma ^i(𝐫)w_{\sigma \sigma ^{}}(𝐫)\delta \rho _\sigma ^{}^j(𝐫)$$
(15)
where $`w_{\sigma \sigma ^{}}(𝐫)`$ are suitable weights to be determined. Here we take for the scalar product (15) that of a bosonic and of a fermionic cloud in the Thomas-Fermi limit with vanishing mutual interactions. We thus set $`w_{FF}=2A(\rho _F^0)^{1/3}/3m_F`$ and $`w_{BB}=g_{BB}/m_B`$ and neglect the off-diagonal weight functions. This choice yields

$$C_\sigma ^i(\omega )=\frac{𝒵_\sigma ^i}{\omega _i^2\omega ^2}$$
(16)
for fermions and bosons, where $`𝒵_\sigma ^i`$ is given in Eqs. (9) and (10). We have taken $`\delta \rho _\sigma ^i`$ as normalized to unity ( $`\delta \rho ^i|\delta \rho ^i=1`$).

The spectral functions are then calculated from the definition

$$S_{\sigma \sigma ^{}}(\omega )=\frac{1}{\pi }\mathrm{Im}d^3r\delta U_\sigma ^{}(𝐫)\delta \rho _\sigma ^{}(𝐫,\omega )$$
(17)
where $`\delta \rho _\sigma ^{}`$ is evaluated by setting $`\delta U_{\overline{\sigma }^{}}=0`$. The use of Eqs. (13) and (16) in Eq. (17) directly leads to Eq. (8) in the main text.


text:
# Firing Time Statistics for Driven Neuron Models: Analytic Expressions versus Numerics

## Abstract

Analytical expressions are put forward to investigate the forced spiking activity of abstract neuron models such as the driven leaky integrate-and-fire model. The method is valid in a wide parameter regime beyond the restraining limits of weak driving (linear response) and/or weak noise. The novel approximation is based on a discrete state Markovian modeling of the full long-time dynamics with time-dependent rates. The scheme yields excellent agreement with numerical Langevin and Fokker-Planck simulations of the full non-stationary dynamics, not only for the first-passage time statistics, but also for the important interspike interval (residence time) distribution.

The detailed modeling of neural behavior presents a prominent challenge on the intriguing vista towards the understanding of neural coding principles. The leaky integrate-and-fire (LIF) model, whose deterministic formulation has been introduced long ago Lap , likely is one of the most studied abstract neuron models Tuckwell . It is characterized by marked simplicity and lack of memory: Whenever the neuron has been excited to fire a pulse it is reset to a predefined state. The beneficial rôle of an appreciable dose of noise has proved to bestow a key part to the interspike statistics of neurons Tuckwell . There exist by now numerous studies and important generalizations of realistic synaptic models, mostly of a numerical nature, that demonstrate the rich behavior of the renewal firing probability, e.g. see in Ref. list ; BulElsDoe96 ; Plesser .

In presence of a time-dependent input-stimulation the stochastic firing process becomes non-stationary, which in turn significantly complicates the stochastic firing statistics. Nevertheless, the signal transmission and its detection can exhibit a remarkable improvement via the phenomenon of Stochastic Resonance GamHaeJun98 . The dynamics of the neuronal firing probability emerges due to a large bombardment of synaptic spike events; consequently, it is customary to employ a diffusion approximation for the stochastic dynamics of the membrane potential $`x(t)`$. The complexity of the driven abstract LIF model thus assumes the archetype, non-stationary Langevin dynamics

$$\dot{x}(t)=\lambda x(t)+\mu +f(t)+\sqrt{2D}\xi (t)$$
(1)
where the process starts at a time $`s`$ at $`x(s)=x_0`$ and fires when it reaches the threshold voltage $`x=a`$. Here, $`f(t)`$ presents a general, time-dependent stimulus which, for example, can be chosen to be oscillatory, and $`\xi (t)`$ is white Gaussian noise. The dynamics of the process $`x(t)`$ is equivalently described by a Fokker-Planck (FP) equation for the conditional probability density $`\rho (x,t|x_0,s)`$ in a time-dependent quadratic potential, $`U(x,t)=\lambda \left[xx_{\text{min}}(t)\right]^2/2`$ with $`x_{\text{min}}(t)=\left(\mu +f(t)\right)/\lambda `$, reading

$$_t\rho =L(t)\rho =_x\left(U^{}(x,t)\rho \right)+D_x^2\rho ,$$
(2)
with the absorbing boundary and initial conditions

$`\rho (a,t|x_0,s)`$ $`=0\text{for all }t\text{}s\text{, and }x_0`$ (3)
$`\rho (x,s|x_0,s)`$ $`=\delta (xx_0).`$ (4)
After firing the process immediately restarts at the instantaneous minimum of the potential.

The set of eqs. (14) defines our starting point for obtaining the firing statistics of this driven neuron model. This is a rather intricate problem because the presence of non-stationarity and multiple time-scales for driving and relaxation, in combination with the absorbing boundary condition prohibits an analytical exact solution RMP . Our main objective is, nevertheless, to develop a most accurate analytical approximation that supersedes all prior attempts known to us. Those attempts, in fact, all involve the use of either of the following limiting approximation schemes such as the limit of linear response theory (i.e. a weak stimulus $`f(t)`$) LR , the limit of asymptotically weak noise GamHaeJun98 ; WN  or the use of the method of images which appears to present an uncontrollable approximation for the case with $`\lambda 0`$ BulElsDoe96 ; Plesser . A most appealing numerical approach is based on an exact integral equation for the first-passage time density of time-dependent Gauss-Markov processes with an absorbing boundary r . Our scheme detailed below yields novel analytic and tractable expressions beyond the linear response and weak noise limit; these are limited solely by the use of a discrete, Markovian stochastic dynamics for the population of the attracting domain and slowly varying (in comparison to intra-well relaxation time-scale) stimuli $`f(t)`$. As demonstrated below, this novel scheme indeed provides analytical formulae that compare very favorably with precise numerical results of the full dynamics in eqs. (1, 24). Different from other approaches, we obtain the distribution not only of the first-passage time but also of the residence time, which is the more interesting variable, concerning neurons.

To start, we approximate the solution to eqs. (24) in the regime where the statistics of times at which the threshold is reached can be characterized by a time-dependent firing rate $`\kappa (t)`$ RMP ; TALU .

This rate then follows from a time-scale separation in the full FP dynamics (2) with boundary condition (3). After a few times $`\lambda ^1`$ of the fast relaxation the probability density $`\rho (x,t)`$ assumes a slowly varying pattern that decays with the rate $`\kappa (t)`$. As in the time-independent case, this slowly varying part of $`\rho (x,t)`$ can be expressed by a product of the (normalized) instantaneous stationary solution $`\rho _0(x,t)\mathrm{exp}\left\{U(x,t)/D\right\}`$ to the FP equation, satisfying $`L(t)\rho _0(x,t)=0`$, and a form function $`\zeta (x,t)`$. Our ansatz thus reads

$$\rho (x,t|x_0,s)\stackrel{\left(ts\right)>\lambda ^1}{\text{}}\zeta (x,t)\rho _0(x,t)\mathrm{exp}\left(\underset{s}{\overset{t}{}}\kappa (s^{})𝑑s^{}\right).$$
(5)
The initial time $`s`$ enters only through the exponential factor. The dependence of the conditional probability on the initial value $`x_0`$ decays exponentially on the timescale $`\lambda ^1`$ and therefore can be neglected for long times.

Deep inside the attracting well, i.e. for $`xa`$ there is no sensible difference between $`\rho (x,t)`$ and $`\rho _0(x,t)`$, and consequently $`\zeta (x,t)`$ approaches one. At the absorbing boundary, however, $`\rho (a,t)`$ and consequently $`\zeta (a,t)`$ both must vanish. The quantitative form of $`\zeta (x,t)`$ in the crossover region follows from

$$L^+(t)\zeta (x,t)=0,$$
(6)
where the potential in the backward operator $`L^+(t)=U^{}(x,t)_x+D_x^2`$ can be linearized about the threshold $`x=a`$. Eq. (6) then yields for the form function the result

$$\zeta (x,t)=1\mathrm{exp}\left\{(xa)\frac{U^{}(a,t)}{D}\right\}.$$
(7)
The rate $`\kappa (t)`$ is determined by multiplying the FP equation (2) in the long-time limit (5) by the form function $`\zeta (x,t)`$ and integrating over $`x`$ from $`\mathrm{}`$ to the threshold voltage $`a`$. In doing so, we account for prominent finite barrier corrections, yielding

$$\kappa (t)=\frac{_{\mathrm{}}^a𝑑x\zeta (x,t)L(t)\zeta (x,t)\rho _0(x,t)}{_{\mathrm{}}^a𝑑x\zeta ^2(x,t)\rho _0(x,t)}.$$
(8)
Upon insertion of eq. (7) for the form function one can exactly perform the integrations and obtains for the rate

$$\kappa (t)=\lambda \frac{\mathrm{\Delta }U(t)}{D}\frac{1\mathrm{erf}\left(\sqrt{\mathrm{\Delta }U(t)/D}\right)}{1\mathrm{exp}(\mathrm{\Delta }U(t)/D)},$$
(9)
where $`\mathrm{\Delta }U(t)`$ denotes the instantaneous potential height at the threshold as seen from the minimum, and $`\mathrm{erf}(z)`$ is the error function. For very small $`D`$ an expansion of the error function leads to the well-known weak noise result for the time-dependent rate GamHaeJun98 ; RMP , i.e.

$$\kappa ^{\text{wn}}(t)=\lambda \sqrt{\mathrm{\Delta }U(t)/(\pi D)}\mathrm{exp}(\mathrm{\Delta }U(t)/D).$$
(10)
Firing time distributions.—With the expression for the exit rate $`\kappa (t)`$ in (9) we can calculate the properties of interest, namely the densities for the first-passage time and the residence time TD  of the attracting ”integrating” state that covers the domain $`\mathrm{}<x(t)<a`$.

The first-passage time distribution is given by the negative rate of change of probability finding the process at time $`t`$ in the ”integrating” state, i.e.

$`g(t|s)`$ $`=_t{\displaystyle _{\mathrm{}}^a}\rho (x,t|x_0,s)𝑑x`$ (11)
$`=\kappa (t)\mathrm{exp}\left(\underset{s}{\overset{t}{}}\kappa (s^{})𝑑s^{}\right),`$ (12)
Here, the integral over the spatial part of $`\rho (x,t)`$ yields unity and the exponential factor in (5) is obtained. It gives the probability for the process to stay in the ”integrating” state from time $`s`$ until $`t`$ without interruption.

The distribution of the residence times $`h(\tau )`$, also termed the interspike interval density, follows as the average of the first-passage time density over the density of resetting times $`s`$, which coincides with the firing rate $`\kappa (s)`$. It thus reads:

$$h(\tau )=\underset{T\mathrm{}}{lim}\frac{_T^Tg(\tau +s|s)\kappa (s)𝑑s}{_T^T\kappa (s)𝑑s}.$$
(13)
Equations (12) and (13), together with the expression for the rate (9) constitute the main results of this work. Their quantitative validity for an extended parameter regime will be checked next.

Numerical comparison.—We have employed three methods for the numerical analysis. A first one is based on the Langevin equation (1) where the position $`x(t)`$ is updated sequentially. For the second we have solved the FP equation (2) numerically, using a Chebychev collocation method to reduce the problem to a coupled system of ordinary differential equations, see also Lehmann . The third method solves an integral equation for the first-passage time probability density and is described in r . All three methods have provided practically identical results.

In the prior stimulating work BulElsDoe96  it has been left open in what relevant parameter regime the employed approximation possesses validity Plesser . Here, using a periodic modulation $`f(t)=A\mathrm{cos}(\omega t)`$ we like to determine a minimal set of relevant parameters that can be taken for comparison with physiological measurements. The most general Langevin equation (1) considered has seven constant parameters, including the threshold $`a`$ and the reset-position $`x_0`$. Three of them, $`(\lambda ,\mu ,a)`$, can be chosen to be $`(1,0,1)`$ by transforming to dimensionless coordinates, using the time-unit $`\lambda ^1`$, the space-unit $`(a\mu /\lambda )`$, and the coordinate-origin at $`\mu /\lambda `$. The resulting parameters thus read (the bars indicate dimensionless coordinates)

$$\overline{A}=A/(\lambda a\mu ),\overline{\omega }=\omega /\lambda ,\overline{D}=D/\left(\lambda (a\mu \lambda )^2\right).$$
(14)
For our purposes it is advantageous to use instead the equivalent set

$$\overline{U}_+/\overline{D},\overline{U}_{}/\overline{D},\text{ and }\overline{\omega },$$
(15)
where $`\overline{U}_+`$ and $`\overline{U}_{}`$ are the maximum and the minimum of $`\overline{U}(1,\overline{t})`$ during a whole period of modulation. These so chosen parameters have the benefit that they provide an on-hand estimate for the validity of our approximations and can be evaluated directly from (1).

The equation we have used in our simulations thus reads (omitting the overbars)

$$\dot{x}(t)=x(t)+A\mathrm{cos}(\omega t)+\sqrt{2D}\xi (t)$$
(16)
with the threshold located at $`x=1`$. For obtaining the residence time, $`x`$ has been reset into the minimum of $`U(x,t)`$ immediately after firing. The Figs. 2 and 2 depict the probability densities of the first-passage and the residence times, respectively. The residence time distribution exhibits a less pronounced modulation in comparison with the first-passage time distribution. Both analytical expressions in (12) and (13) compare very favorably with the numerical results obtained by iteration of the Langevin equation (16). The remaining deviations in the two figures are of purely statistical nature (see the inset in Fig. 2) and can be diminished further by increasing the number of events in the simulations.

In order to further test the range of validity of our novel approximation scheme we – on purpose – have chosen extreme values for the lower barrier height $`U_{}/D`$ and angular driving frequency $`\omega `$, respectively, see Figs. 4 and 4. Here, deviations from the numerical results are not the result of statistics but are systematic. For the low potential barrier in Fig. 4, the time-scales in the process are not separated sufficiently. The fast intra-well fluctuations begin to influence the behavior of the modulated firing dynamics. Moreover, the difference between the moderate noise result for the time-dependent rate $`\kappa (t)`$ in (9) and its weak noise approximation in (10) becomes visibly increased, as expected. Figure 4 depicts the other extreme situation with a modulation time-scale that is not slow enough. Because the system cannot follow the driving instantaneously we find a shift in the maxima of the first-passage time density. This shift is not reproduced by our approximation in (5) and (12); nevertheless, our scheme yields amazingly good results even within this extreme parameter regime. The results based on the numerical evaluation of the FP equation (2) and the integral equation r  virtually collapse into one curve and perfectly coincide with the Langevin simulations. The same result was obtained for the other parameter values.

Conclusions.—The precise theoretical modeling of the neuronal spiking activity under external time-dependent driving presents a challenge of considerable importance in neurophysiology and physics. Due to the presence of non-stationarity, absorbing boundary conditions of the underlying first-passage problem and differing time-scales the task of obtaining reliable analytical estimates for the firing statistics is anything but trivial. By reference to a discrete Markovian dynamics for the corresponding full space-continuous stochastic dynamics we succeeded in obtaining analytical approximations for the time-dependent first-passage time and the residence time statistics that are valid beyond the restraining limits of linear response and asymptotically weak noise. We have tested our findings for the case of a periodically driven LIF model. The obtained agreement with precise numerical simulations of either the Langevin type in (16) or, equivalently, of the FP type in (2) turns out to be very good. Our method is not restricted to an oscillatory forcing but applies as well to arbitrary drive functions $`f(t)`$ such as an exponentially decaying drive (e.g. simulating a decaying threshold). Our scheme even yields good results in extreme parameter regimes where agreement cannot be expected a priori.

Our method, primarily aimed at describing first-passage time and residence time probabilities of driven dynamical systems, is also readily extended to more realistic neuron models such as e.g. the two-dimensional driven FitzHugh-Nagumo model FN  for neuronal spiking activity, whose multiple attractors may be considered as discrete states. Likewise, the scheme can also be employed to study yet other time-dependent switching dynamics and synchronization phenomena such as the paradigm of Stochastic Resonance GamHaeJun98  and discrete or continuous Brownian motor transport BM .

This work has been supported by the Deutsche Forschungsgemeinschaft via Projects No. HA1517/13-4 and SFB-486, Projects No. A10 and No. B13.