# One-body reduced density-matrix functional theory for the canonical ensemble

S.M. Sutter<EMAIL_ADDRESS>K.J.H. Giesbertz<EMAIL_ADDRESS>Theoretical Chemistry, Faculty of Exact Sciences, VU University, De Boelelaan
1083, 1081 HV Amsterdam, The Netherlands

(August 27, 2024)

###### Abstract

We establish one-body reduced density-matrix functional theory for the
canonical ensemble in a finite basis set at an elevated temperature. Including
temperature guarantees differentiability of the universal functional by
occupying all states and additionally not fully occupying the states in a
fermionic system. We use convexity of the universal functional and
invertibility of the potential-to-1RDM map to show that the subgradient
contains only one element which is equivalent to differentiability. This
allows us to show that all 1RDMs with a purely fractional occupation number
spectrum ($0<n_{i}<1\;\forall_{i}$) are uniquely $v$-representable up to a
constant.

##  I Introduction

Quantum chemistry and physics deal with the description of many interacting
particles. Often we limit ourselves to a single particle species. In quantum
chemistry these are usually electrons, but in physics also bosonic particles
are of interest. Though the many-body Schrödinger equation involves only
linear operators, the daunting dimensionality of the many-body wave function
renders a direct solution intractable, but for a few particles. This is one of
the prime reasons to aim directly for reduced quantities.

In 1964 Hohenberg and Kohn presented their revolutionary work about density
functional theory (DFT) [1]. They showed that any observable can be regarded
as a functional of the density. Especially the Kohn–Sham (KS) formulation [2]
has been important to the success of DFT. Their idea was to approximate the
true kinetic energy by the kinetic energy of the KS system: a non-interacting
system with the same density as the interacting system. The KS kinetic energy
turns out to be a decent approximation to the true kinetic energy of the
interacting system. The difference in the kinetic energy is then lumped
together with the interaction beyond Hartree (classical Coulomb) in the
exchange-correlation energy functional. Though formally exact, in practice KS-
DFT has some weaknesses, since the exchange-correlation functional needs to be
approximated. A famous example is the stretching of the H2 bond [3, 4, 5].

One way to bypass some of these problems in constructing an approximate
exchange-correlation energy functional is (one-body) reduced density matrix
(1RDM) functional theory. One advantage over DFT is that we have also an
explicit expression for the kinetic energy while still having the total energy
as a functional of the 1RDM [6]. However, in the zero temperature setting,
mapping back from 1RDMs to (non-local) potentials is problematic, as already
noted by Gilbert [6] and others [7, 8, 9, 10]. This is most clear in the case
of non-interacting particles, since typically ground state 1RDMs are
idempotent. It therefore seems that non-idempotent 1RDMs cannot be
$v$-representable in the absence of interactions. There is the possibility for
orbital energies to be degenerate, however, which allows fractionally occupied
orbitals and hence non-idempotent 1RDMs [11, 12]. But one quickly realizes
that the scaled identity operator is the one-body Hamiltonian which has all
1RDMs as ground state 1RDM. It is clear that the interaction should play a
crucial role in the back mapping, but there has been no progress in this
direction.

An alternative to regularize the theory is to introduce entropy, i.e., work at
finite temperature, as proposed more than a decade ago [8, 9, 13, 14] Though
this is a theoretical motivation to introduce temperature, also physically
this is a well justified choice, since most experiments are conducted at
$T>0$. Important examples where temperature plays an important role are metal-
insulator transitions in transition metal oxides [15], high $T_{c}$ super
conductors [16], hot plasmas [17], etc.

In Ref. [10] 1RDM functional theory (1RDMFT) was presented for the grand
canonical ensemble within a finite basis set. However, the use of a grand
canonical ensemble is inappropriate if the number of particles is relatively
low as in ultra cold atom experiments [18], but also in the low temperature
limit the grand canonical ensemble can lead to unphysical results [19, 20]. A
canonical formulation of 1RDMFT is therefore desirable and will be the goal of
this article.

In classical thermodynamics the grand potential can be reached by a Legendre
transformation of the Helmholtz free energy with respect to the number of
particles. In the quantum mechanical setting we can not do that. The reason is
that the grand potential and the Helmholtz functional act on different spaces:
the Fock space and the $N$-particle Hilbert spaces, respectively. Thus, we can
not simply transform it back to obtain the canonical case. Another major
difference w.r.t. the grand canonical ensemble is that in the non-interacting
case the occupation numbers are not explicitly given by either the Fermi/Bose
function for fermions/bosons. Instead, they need to be calculated recursively,
using auxiliary partition functions [21]. In Ref. [14] it therefore remained
an open question whether every thermal 1RDM (only fractional occupation
numbers) would be non-interacting $v$-representable and not much progress
could be made. However, we do not rely on such an explicit relation and we are
able to prove a one-to-one correspondence between thermal 1RDMs (all
occupation numbers fractional) and (non-local) potentials for any interaction.
The non-interacting system is just a particular case. This result justifies
the existence of an algorithm which finds for any thermal 1RDM the
corresponding non-interacting Hamiltonian as published recently by Kooi [22].

In this work we present 1RDMFT in a rigorous way for a fixed number of
particles, finite basis set and elevated temperature. We show that the
universal functional is differentiable and it holds

$\displaystyle\frac{\partial F_{N,\pm}}{\partial\gamma}=-v,$ (1)

where $\gamma$ is the ground state 1RDM for the potential $v$. Here and in the
following $+$ and $-$ stand for the bosonic and the fermionic case
respectively. If we have a handy expression for $F_{N,\pm}[\gamma]$ then we
can circumvent the handling of the density-matrix operator to compute the free
energy and instead we only need to deal with the reduced quantity $\gamma$.
The minimizer for the Helmholtz functional can then be determined through the
above relation (1).

This work is built up in the following way. In Section II we introduce all the
relevant spaces, then, in Section III, we present the Helmholtz functional,
its minimizing density-matrix operator and the general approach for 1RDMFT.
For this task we make use of the universal functional $F_{N,\pm}$. To show
differentiability of $F_{N,\pm}$ we utilize results from convex analysis. In
Section IV we show that all the relevant functionals are convex. Additionally,
we show that two potentials differing by more than a constant can not generate
the same density-matrix operator. The proof of differentiability of
$F_{N,\pm}$ is finalized in Section V.

##  II Setting

We build our $N$-particle space from a finite number of single particle states
$\ket{i}$, for $i\in\\{1,\dotsc,N_{b}\\}$ and $N_{b}<\infty$. We require the
states to be orthonormal. The one-particle Hilbert space $\mathcal{H}$ is now
the ${\mathbb{C}}$-vector space generated by the states $\ket{i}$, i.e.,
$\mathcal{H}\coloneqq\mathrm{span}\\{\ket{1},\dotsc,\ket{N_{b}}\\}\cong{\mathbb{C}}^{N_{b}}$.
To build the $N$-particle space we need to distinguish between bosons and
fermions.

Bosons

A system with $N$ bosons is described by a symmetric wave function. Therefore,
the bosonic $N$-particle Hilbert state, denoted by $\mathcal{H}^{N}_{+}$,
consists of all symmetric tensors of order $N$, i.e.,
$\mathcal{H}_{+}^{N}\coloneqq\mathrm{Sym}^{N}(\mathcal{H})$. The dimension of
$\mathcal{H}_{+}^{N}$ is $\binom{N_{b}+N-1}{N}$.

Fermions

Fermionic systems are described by anti-symmetric wave functions. Thus, we
consider the space of anti-symmetric tensors of order $N$, i.e.
$\mathcal{H}_{-}^{N}\coloneqq\wedge^{N}\mathcal{H}$. The dimension is given by
$\binom{N_{b}}{N}$. Note that we need to have $N_{b}\geq N$. The case
$N_{b}=N$ is trivial since we have only one possible state. Thus, we will only
consider $N_{b}>N$.

The set of density-matrix operators on the $N$-particle space
$\mathcal{H}_{\pm}^{N}$ is defined as

$\overline{\mathcal{P}}_{N,\pm}\coloneqq\Set{\\!\hat{\rho}:\mathcal{H}_{\pm}^{N}\\!\to\mathcal{H}_{\pm}^{N}}{\hat{\rho}=\hat{\rho}^{\dagger},\hat{\rho}\geq
0,\operatorname{Tr}\\{\hat{\rho}\\}=1\\!\\!},$ (2)

which we endow with the norm

$\lVert\hat{\rho}\rVert_{2}=\big{(}\operatorname{Tr}\\{\lvert\hat{\rho}\rvert^{2}\\}\big{)}^{1/2}.$
(3)

A density-matrix operator $\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}$ has a
spectral decomposition

$\hat{\rho}=\sum\limits_{l}\lambda_{l}\ket{\psi_{l}}\bra{\psi_{l}},$ (4)

and its kernel is given by

$\rho(x_{1},\ldots,x_{N};y_{1},\ldots,y_{N})\\\
=\sum\limits_{l}\lambda_{l}\psi_{l}(x_{1},\ldots,x_{N})\psi_{l}^{*}(y_{1},\ldots,y_{N}).$
(5)

We can define the 1RDM $\gamma$ by 111Creation and annihilation operator are
only defined pairwise for a $N$ particle space.

$\displaystyle\gamma_{ij}[\hat{\rho}]=\operatorname{Tr}\\{\hat{\rho}\,\hat{a}^{\dagger}_{j}\hat{a}^{\vphantom{\dagger}}_{i}\\}.$

It turns out (see Appendix) that the relevant spaces for the 1RDMs are subsets
of the space of all Hermitian $N_{b}\times N_{b}$ matrices denoted by
$\mathbb{H}(N_{b}),$

$\displaystyle\overline{\mathscr{N}}_{N,+}\\!$
$\displaystyle\coloneqq\\!\Set{\\!\gamma\in\mathbb{H}(N_{b})}{\\!\gamma\geq
0,\operatorname{tr}\\{\gamma\\}=N},$ (6a)
$\displaystyle\overline{\mathscr{N}}_{N,-}\\!$
$\displaystyle\coloneqq\\!\Set{\\!\gamma\in\mathbb{H}(N_{b})}{\\!\gamma\geq
0,\gamma^{2}\\!\leq\gamma,\operatorname{tr}\\{\gamma\\}=N\\!\\!}.$ (6b)

We have used $\operatorname{tr}\\{\cdot\\}$ to emphasize that the trace is
over the one-particle Hilbert space $\mathcal{H}$ as opposed to the
$\operatorname{Tr}\\{\cdot\\}$ which is over a $\mathcal{H}^{N}_{\pm}$ Hilbert
space. By convention, the eigenvalues and eigenstates of the 1RDM $\gamma$ are
called natural occupation numbers and natural orbitals (NO) respectively.
Coleman has shown that all elements of $\overline{\mathscr{N}}_{N,\pm}$ can be
obtained from a density-matrix operator in $\overline{\mathcal{P}}_{N,\pm}$,
so it is a true 1RDM [24].

######  Theorem II.1 (Coleman).

For any $\gamma\in\overline{\mathscr{N}}_{N,\pm}$ there is a density matrix
$\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}$ which generates $\gamma$.

The proof can be found in the Appendix.

##  III General Approach

The Helmholtz functional for the canonical ensemble is defined as

$\Omega_{v}[\hat{\rho}]\coloneqq E_{v}[\hat{\rho}]-\beta^{-1}S[\hat{\rho}],$
(7)

where

$E_{v}[\hat{\rho}]\coloneqq\operatorname{Tr}[\hat{\rho}\hat{H}_{v}]$ (8)

is the energy of a system with Hamiltonian
$\hat{H}_{v}\coloneqq\hat{H}_{0}+\hat{V}_{v}$ ($\hat{H}_{0}$ contains the
kinetic and interaction part and $\hat{V}_{v}$ is the potential with kernel
$v(x,x^{\prime})$). The second term contains the entropy

$S[\hat{\rho}]\coloneqq-\operatorname{Tr}\\{\hat{\rho}\log(\hat{\rho})\\}$ (9)

and the inverse temperature $\beta=1/T$. With $\log$ we mean the natural
logarithm. The minimizer $\hat{\rho}_{v}$ of the Helmholtz functional can be
found by variations in the density-matrix operator which yields the equation

$\operatorname{Tr}\\{\delta\hat{\rho}(\hat{H}_{v}+\beta^{-1}\log(\hat{\rho}_{v}))\\}+\beta^{-1}\operatorname{Tr}\big{\\{}\delta\hat{\rho}\big{\\}}=0.$
(10)

From the unit trace condition and (10) it follows that

$\displaystyle\hat{\rho}_{v}$ $\displaystyle=e^{-\beta\hat{H}_{v}}/Z[v],$
where $\displaystyle Z[v]$
$\displaystyle\coloneqq\operatorname{Tr}\bigl{\\{}e^{-\beta\hat{H}_{v}}\bigr{\\}}.$
(11)

The minimizer $\hat{\rho}_{v}$ is called Gibbs state. Note that we only have a
proper solution for $0<Z[v]<\infty$. This is always the case since we work in
a finite basis setting with a fixed number of particles, so the trace only
runs over a finite number of elements.

One aim is to show that the map from the potential $v$ to the density-matrix
operator $\hat{\rho}_{v}$ is invertible. However, this is only doable up to a
constant since adding a constant to the potential does not change the density-
matrix operator. To achieve a one-to-one correspondence we allow only
potentials from the following set,

$\mathcal{V}\coloneqq\Set{v\in\mathbb{H}(N_{b})}{\operatorname{tr}\\{v\\}=0}.$
(12)

We can also think of $v\in\mathcal{V}$ being a representative of the
equivalence class containing potentials differing by a constant.

In Theorem II.1 we have seen that all
$\gamma\in\overline{\mathscr{N}}_{N,\pm}$ are $N$-representable. However,
physically relevant are only the 1RDMs that are associated with a Gibbs state
$\hat{\rho}_{v}$. Thus, we denote the set of all $v$-representable 1RDM by

$\mathscr{V}_{N,\pm}\coloneqq\Set{\gamma\in\overline{\mathscr{N}}_{N,\pm}}{\exists\;v\in\mathcal{V}\mapsto\gamma}.$
(13)

The approach is to partition the minimization in the Helmholtz functional as

$\Omega_{N,\pm}[v]\coloneqq\inf_{\mathclap{\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}}}\Omega_{v}[\hat{\rho}]=\inf_{\mathclap{\gamma\in\overline{\mathscr{N}}_{N,\pm}}}\big{(}F_{N,\pm}[\gamma]+\operatorname{tr}\\{v\gamma\\}\big{)}$
(14)

where

$\displaystyle F_{N,\pm}[\gamma]\coloneqq{}$
$\displaystyle\inf_{\begin{subarray}{c}\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}\\\
\hat{\rho}\to\gamma\end{subarray}}\Omega_{0}[\hat{\rho}]$ $\displaystyle{}={}$
$\displaystyle\inf_{\begin{subarray}{c}\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}\\\
\hat{\rho}\to\gamma\end{subarray}}\operatorname{Tr}\big{\\{}\hat{\rho}\big{(}\hat{H}_{0}+\beta^{-1}\log(\hat{\rho})\big{)}\big{\\}}$
(15)

is called the universal functional which takes the value $\infty$ in case no
$\hat{\rho}\to\gamma$ exists. Here and in the following,
$\Omega_{0}[\hat{\rho}]=\Omega_{v=0}[\hat{\rho}]$. The aim is to show that
$F_{N,\pm}$ is differentiable. Then the minimizer can be found through the
relation

$\frac{\partial F_{N,\pm}}{\partial\gamma}=-v$

and we know that $\gamma$ is a canonical eq-1RDM (equilibrium 1-RDM) which was
an open question in [14].

##  IV General Properties of the Helmholtz Functional and Implications on the
Universal Functional

######  Theorem IV.1.

The mapping $\hat{H}_{v}\mapsto\hat{\rho}_{v}$ with $v\in\mathcal{V}$ is
invertible up to a constant in the Hamiltonian.

###### Proof.

Assume that two Hamiltonians $\hat{H}_{v}$ and $\hat{H}_{v}^{\prime}$
differing in their potential yield the same density-matrix operator
$\hat{\rho}_{v}$. From (10) it follows that $\hat{\rho}_{v}$ fulfills

$\displaystyle\frac{1}{\beta}\log(\hat{\rho}_{v})+\hat{H}_{v}$
$\displaystyle=C,$
$\displaystyle\frac{1}{\beta}\log(\hat{\rho}_{v})+\hat{H}_{v}^{\prime}$
$\displaystyle=C^{\prime}.$

Subtracting these equations gives
$\hat{H}_{v}-\hat{H}_{v}^{\prime}=C-C^{\prime}$. ∎

######  Remark 1.

Since we have a fixed number of particles in the Hilbert space, the constant
in Theorem IV.1 can be of the form $f(\hat{N})$ where
$f:{\mathbb{R}}\to{\mathbb{R}}$, so this includes the arbitrary constant shift
in the potential.

######  Corollary 1.

The map $v\mapsto\hat{\rho}_{v}$ with $v\in\mathcal{V}$ is invertible.

Note that we only have a one-to-one correspondence because we require
$\operatorname{tr}\\{v\\}=0$. Otherwise a constant shift in the potential
would lead to the same density-matrix operator.

At this point we want to mention that the density-matrix operator
$\hat{\rho}_{v}$ is positive definite, $\hat{\rho}_{v}>0$, and lies in the
following subspace of $\overline{\mathcal{P}}_{N,\pm}$,

$\mathcal{P}_{N,\pm}\coloneqq\Set{\hat{\rho}:\mathcal{H}^{N}_{\pm}\to\mathcal{H}_{\pm}^{N}}{\hat{\rho}=\hat{\rho}^{\dagger},\hat{\rho}>0,\operatorname{Tr}\\{\hat{\rho}\\}=1}.$
(16)

It follows that the natural occupation numbers $n_{i}$ are positive and in the
fermionic case additionally $n_{i}<1$. To see this let
$\phi_{1},\dotsc,\phi_{N_{b}}$ be the NO basis and
$\hat{\rho}_{v}=\sum_{j}\lambda_{j}\ket{\psi_{j}}\bra{\psi_{j}}$ be the
spectral decomposition of the density-matrix operator. Then, as the
$\psi_{j}$’s build a basis of $\mathcal{H}_{\pm}^{N}$, each NO $\phi_{i}$
contributes to at least one of the eigenstates. So,

$n_{i}=\sum\limits_{j}\lambda_{j}\int\mathrm{d}x\mathrm{d}y\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\phi_{i}^{*}(x)\phi_{i}(y)\\\
\psi_{j}(x,x_{2},\dotsc,x_{N})\psi_{j}^{*}(y,x_{2},\dotsc,x_{N})>0,$ (17)

where we used the fact that all weights $\lambda_{j}=e^{-\beta E_{j}}/Z$ are
positive. In case of fermions we have already showed that $n_{i}\leq 1$. The
$i$th NO can not be present in all $\psi_{j}$’s (in case $N\neq N_{b})$, so

$n_{i}=\sum\limits_{j}\lambda_{j}\int\mathrm{d}x\mathrm{d}y\mathrm{x_{2}}\dotsi\mathrm{d}x_{N}\phi_{i}^{*}(x)\phi_{i}(y)\\\
\psi_{j}(x,x_{2},\dotsc,x_{N})\psi_{j}^{*}(y,x_{2},\dotsc,x_{N})<1,$ (18)

because the integral is at least for one $j$ not equal to 1. Summarized, the
1RDMs produced by a potential are contained in

$\displaystyle\mathscr{N}_{N,+}$
$\displaystyle\coloneqq\Set{\gamma\in\mathbb{H}(N_{b})}{\gamma>0},$ (19)
$\displaystyle\mathscr{N}_{N,-}$
$\displaystyle\coloneqq\Set{\gamma\in\mathbb{H}(N_{b})}{\gamma>0,\gamma^{2}<\gamma}.$
(20)

We want to show that the most important functionals are either convex or
concave to be able to use results from convex analysis. We start with the
functional $\Omega_{N,\pm}[v]$ which is achieved through a minimization and
thus turns out to be concave [25].

######  Theorem IV.2.

The functional
$\Omega_{N,\pm}[v]=\min\limits_{\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}}\Omega_{v}[\hat{\rho}]=\min\limits_{\hat{\rho}\in\mathcal{P}_{N,\pm}}\Omega_{v}[\hat{\rho}]=-\beta^{-1}\log\big{(}Z[v]\big{)}$
is strictly concave in $v$.

###### Proof.

Let $v_{1}\neq v_{2}$ be two potentials in $\mathcal{V}$ and let $0<t<1$. Then
we have

$\displaystyle\Omega_{N,\pm}[tv_{1}+(1-t)v_{2}]$
$\displaystyle=\min_{\hat{\rho}\in\mathcal{P}_{N,\pm}}\operatorname{Tr}\Big{\\{}\hat{\rho}\big{(}t\hat{H}_{v_{1}}+t\frac{1}{\beta}\log(\hat{\rho})\big{)}\Big{\\}}+\operatorname{Tr}\Big{\\{}\hat{\rho}\big{(}(1-t)\hat{H}_{v_{2}}+(1-t)\frac{1}{\beta}\log(\hat{\rho})\big{)}\Big{\\}}$
$\displaystyle>t\min_{\hat{\rho}_{1}\in\mathcal{P}_{N,\pm}}\operatorname{Tr}\Big{\\{}\hat{\rho}_{1}\big{(}\hat{H}_{v_{1}}+\frac{1}{\beta}\log(\hat{\rho}_{1})\big{)}\Big{\\}}+(1-t)\min_{\hat{\rho}_{2}\in\mathcal{P}_{N,\pm}}\operatorname{Tr}\Big{\\{}\hat{\rho}_{2}\big{(}\hat{H}_{v_{2}}+\frac{1}{\beta}\log(\hat{\rho}_{2})\big{)}\Big{\\}}$
$\displaystyle=t\Omega_{N,\pm}[v_{1}]+(1-t)\Omega_{N,\pm}[v_{2}]$

where the strict inequality follows from Corollary 1. ∎

With Corollary 1 it is possible to show a generalization of the Hohenberg–Kohn
theorem for 1RDMs and non-local potentials [26].

######  Theorem IV.3.

The map $v\mapsto\gamma_{v}$ for $v\in\mathcal{V}$ is invertible.

###### Proof.

Assume there are two potentials $v_{1}\neq v_{2}\in\mathcal{V}$ yielding
different density matrix operators $\hat{\rho}_{1}\neq\hat{\rho}_{2}$ but the
same 1RDM $\gamma$. Then we get

$\displaystyle\Omega_{N,\pm}[v_{1}]$
$\displaystyle=\Omega_{v_{1}}[\hat{\rho}_{1}]=\Omega_{v_{2}}[\hat{\rho}_{1}]+\operatorname{tr}\big{\\{}\gamma(v_{1}-v_{2})\big{\\}}$
$\displaystyle>\Omega_{v_{2}}[\hat{\rho}_{2}]+\operatorname{tr}\big{\\{}\gamma(v_{1}-v_{2})\big{\\}}$
$\displaystyle=\Omega_{N,\pm}[v_{2}]+\operatorname{tr}\big{\\{}\gamma(v_{1}-v_{2})\big{\\}}.$

Changing the role of $v_{1}$ and $v_{2}$ and adding the two inequalities gives

$\displaystyle\Omega_{N,\pm}[v_{1}]+\Omega_{N,\pm}[v_{2}]<\Omega_{N,\pm}[v_{2}]+\Omega_{N,\pm}[v_{1}],$

which is a contradiction. ∎

One aim is to show that the universal functional is convex. For this purpose
we first show that the entropy is strictly concave ([27], [28], [29]).

######  Theorem IV.4.

The entropy is strictly concave, i.e., for any
$\hat{\rho}_{0},\hat{\rho}_{1}\in\overline{\mathcal{P}}_{N,\pm}$ and
$\lambda\in(0,1)$ we have
$S[\lambda\hat{\rho}_{0}+(1-\lambda)\hat{\rho}_{1}]>\lambda
S[\hat{\rho}_{0}]+(1-\lambda)S[\hat{\rho}_{1}]$.

###### Proof.

Let
$\hat{\rho}_{\lambda}=\lambda\hat{\rho}_{0}+(1-\lambda)\hat{\rho}_{1}=\sum_{k}w_{k}\ket{\psi_{k}}\bra{\psi_{k}}$.
We use strict concavity of the function $s(x)=-x\log(x)$ and we get

$\displaystyle S[\hat{\rho}_{\lambda}]$
$\displaystyle=-\sum_{k}w_{k}\log(w_{k})=\sum_{k}s\big{(}\braket{\psi_{k}}{\hat{\rho}_{\lambda}}{\psi_{k}}\big{)}$
$\displaystyle=\sum_{k}s\big{(}\lambda\braket{\psi_{k}}{\hat{\rho}_{0}}{\psi_{k}}+(1-\lambda)\braket{\psi_{k}}{\hat{\rho}_{1}}{\psi_{k}}\big{)}$
$\displaystyle>\lambda\sum_{k}s\big{(}\braket{\psi_{k}}{\hat{\rho}_{0}}{\psi_{k}}\big{)}+(1-\lambda)\sum_{k}s\big{(}\braket{\psi_{k}}{\hat{\rho}_{1}}{\psi_{k}}\big{)}$
$\displaystyle\geq\lambda\sum_{k}\braket{\psi_{k}}{s(\hat{\rho}_{0})}{\psi_{k}}+(1-\lambda)\sum_{k}\braket{\psi_{k}}{s(\hat{\rho}_{1})}{\psi_{k}}$
$\displaystyle=\lambda S[\hat{\rho}_{0}]+(1-\lambda)S[\hat{\rho}_{1}]$

where we used Jensen’s inequality for the last inequality. ∎

######  Corollary 2.

The Helmholtz functional $\Omega_{v}[\hat{\rho}]$ is strictly convex in the
density operator $\hat{\rho}$.

###### Proof.

It follows directly from the fact that the Helmholtz functional is the sum of
a linear and a strictly convex functional. ∎

######  Theorem IV.5.

The universal functional $F_{N,\pm}[\gamma]$ is convex on
$\overline{\mathscr{N}}_{N,\pm}$.

###### Proof.

Let $\gamma_{0},\gamma_{1}\in\overline{\mathscr{N}}_{N,\pm}$,
$\lambda\in[0,1]$ and
$\gamma_{\lambda}=\lambda\gamma_{0}+(1-\lambda)\gamma_{1}$ and taking
$\hat{\rho},\hat{\rho}_{0},\hat{\rho}_{1}\in\overline{\mathcal{P}}_{N,\pm}$ we
get

$\lambda F_{N,\pm}[\gamma_{0}]+(1-\lambda)F_{N,\pm}[\gamma_{1}]\\\
\begin{aligned}
&=\lambda\inf_{\hat{\rho}_{0}\to\gamma_{0}}\Omega_{0}[\hat{\rho}_{0}]+(1-\lambda)\inf_{\hat{\rho}_{1}\to\gamma_{1}}\Omega_{0}[\hat{\rho}_{1}]\\\
&=\inf_{\hat{\rho}_{0}\to\gamma_{0}}\inf_{\hat{\rho}_{1}\to\gamma_{1}}\lambda\Omega_{0}[\hat{\rho}_{0}]+(1-\lambda)\Omega_{0}[\hat{\rho}_{1}]\\\
&\geq\inf_{\hat{\rho}_{0}\to\gamma_{1}}\inf_{\hat{\rho}_{1}\to\gamma_{1}}\Omega_{0}[\lambda\hat{\rho}_{0}+(1-\lambda)\hat{\rho}_{1}]\\\
&=\inf_{\hat{\rho}\to\gamma_{\lambda}}\Omega_{0}[\hat{\rho}]=F_{N,\pm}[\gamma_{\lambda}].&\qed\end{aligned}$

##  V Final Result

Now we want to show that the universal functional $F_{N,\pm}$ is
differentiable. Differentiability is only defined on an open set. However, the
set $\overline{\mathscr{N}}_{N,\pm}$ has empty interior in
$\mathbb{H}(N_{b})$. Thus, we need to embed $\overline{\mathscr{N}}_{N,\pm}$
in a topological space where $\mathscr{N}_{N,\pm}$ is the interior of
$\overline{\mathscr{N}}_{N,\pm}$. The idea is to use the following result
about subgradients and subdifferentials.

######  Theorem V.1.

Let $X$ be a finite dimensional vector space and let
$f:X\to{\mathbb{R}}\cup\\{\infty\\}$ be a convex function with domain $M$.
Assume $M$ is contained in $a+\mathfrak{L}$ such that $\mathfrak{L}$ is a
subspace with the lowest dimension such that there exists $a\in X$ with
$M\subset a+\mathfrak{L}$ 222This means $M$ is contained in an affine space.
Let $\partial_{\mathfrak{L}}f(x)\coloneqq\partial f(x)\cap\mathfrak{L}$ where
$\partial f(x)$ is the subdifferential of $f$ at a point $x$ in the set $M$.
Then the following properties hold for $\partial_{\mathfrak{L}}f(x)$

  1. (i)

the set $\partial_{\mathfrak{L}}f(x)$ is nonempty,

  2. (ii)

$f$ is differentiable at $x$ if and only if $\partial_{\mathfrak{L}}f(x)$
contains only one element. In that case this element equals the usual
gradient. (With differentiable we mean that there is a linear map
$J:\mathfrak{L}\to{\mathbb{R}}$ such that for all $h\in\mathfrak{L}$ we have
$\lim_{h\to 0}\frac{1}{\lVert h\rVert_{\mathfrak{L}}}\lvert
f(x+h)-f(x)-J(h)\rvert=0$.

$n_{1}$$n_{2}$$n_{3}$$a$$2$$2$$2$ Figure 1: Representation of
$\overline{\mathscr{N}}_{N,+}$ in terms of the occupation numbers for
$N_{b}=3$ and $N=2$ (See proof of Theorem II.1). The vector $a$ shifts a $2$
dimensional subspace to an affine space which contains
$\overline{\mathscr{N}}_{N,+}$.

The universal functional has domain $\overline{\mathscr{N}}_{N,\pm}$ which is
contained in $\mathfrak{L}+a$ with
$\mathfrak{L}=\Set{x\in\mathbb{H}(N_{b})}{\operatorname{tr}\\{x\\}=0}$ and
$a=N/N_{b}\cdot\mathds{1}$ 333We could have taken a different point in
$\mathscr{N}_{N,\pm}$, but the current choice is the most natural in this
setting. and its (relative) interior is $\mathscr{N}_{N,\pm}$. Fig. 1 shows
$\overline{\mathscr{N}}_{N,+}$ for $2$ particles and $3$ basis functions. This
also justifies the choice of the potential gauge. The aim is to get the
relation $\partial F_{N,\pm}/\partial\gamma=-v$. But as mentioned in Theorem
(V.1), the differential is a map $J:\mathfrak{L}\to{\mathbb{R}}$, i.e. it is
contained in $\mathcal{V}$. We can now apply the above theorem for all 1RDMs
contained in $\mathscr{N}_{N,\pm}$.

######  Theorem V.2.

If the infimum in (III) is attained, then

  1. (i)

$\mathscr{N}_{N,\pm}=\mathscr{V}_{N,\pm}$

  2. (ii)

the universal functional $F_{N,\pm}[\gamma]$ is differentiable on
$\mathscr{N}_{N,\pm}$.

###### Proof.

Convexity of $F_{N,\pm}$ implies that for any $\gamma\in\mathscr{N}_{N,\pm}$
there exists at least one subgradient $h\in\mathcal{V}$. So for all
$\tilde{\gamma}\in\overline{\mathscr{N}}_{N,\pm}$ it holds
$F_{N,\pm}[\tilde{\gamma}]+\braket{-h}{\tilde{\gamma}}\geq
F_{N,\pm}[\gamma]+\braket{-h}{\gamma}$, which implies

$\displaystyle
F_{N,\pm}[\gamma]+\braket{-h}{\gamma}\leq\min_{\mathclap{\tilde{\gamma}\in\overline{\mathscr{N}}_{N,\pm}}}F_{N,\pm}[\tilde{\gamma}]+\braket{-h}{\tilde{\gamma}}=\Omega_{N,\pm}[-h].$

Thus, the negative of the subgradient, $-h$, yields a potential generating
$\gamma$ and hence $\mathscr{N}_{N,\pm}=\mathscr{V}_{N,\pm}$. By Theorem IV.3,
we get that there is only one such potential. Hence the subgradient is unique
and $F_{N,\pm}[\gamma]$ is differentiable for all
$\gamma\in\mathscr{N}_{N,\pm}$ by Theorem V.1. ∎

We proved $v$-representability under the assumption that the minimum in (III)
is attained. To finish the proof we still need to show that this is indeed the
case. The idea is to show that the relevant functions are continuous and then
use the fact that continuous functions attain their minima (and maxima) over
compact sets.

######  Proposition V.3.

The energy $E_{v}[\hat{\rho}]$ is Lipschitz continuous on
$\overline{\mathcal{P}}_{N,\pm}$.

###### Proof.

The Hamiltonian acts on a finite dimensional space and thus it has a maximum
eigenvalue, $\lVert\hat{H}_{v}\rVert_{\infty}<\infty$. For two density-matrix
operators $\hat{\rho}_{0},\hat{\rho}_{1}$ we have

$\displaystyle\lvert E_{v}[\hat{\rho}_{0}]-E_{v}[\hat{\rho}_{1}]\rvert$
$\displaystyle\leq\lVert\hat{H}_{v}\rVert_{\infty}\lVert\hat{\rho}_{0}-\hat{\rho}_{1}\rVert.\qed$

######  Proposition V.4.

The entropy is continuous on $\overline{\mathcal{P}}_{N,\pm}$.

###### Proof.

The function $\hat{\rho}\mapsto(a_{1},\dotsc,a_{m})$ where
$(a_{1},\dotsc,a_{m})$ contains the eigenvalues of $\hat{\rho}$ in an ordered
fashion and with multiplicity ($a_{i}\leq a_{i+1}$ for all $i=1,\dotsc,m-1$)
is continuous. The statement follows from continuity of the function
$(a_{1},\dotsc,a_{m})\mapsto\sum_{j=1}^{m}a_{j}\log(a_{j})$. ∎

To finish the proof we make use of the following theorem.

######  Theorem V.5.

Let $X$ be a compact metric space and let $f:X\to{\mathbb{R}}$ be a continuous
function. Then $f$ is bounded and it attains its maximum and minimum.

Compact sets in finite dimensional affine spaces (with the usual metric) are
fully characterized by closedness and boundedness. Thus,
$\overline{\mathcal{P}}_{N,\pm}$ and
$\set{\hat{\rho}\in\overline{\mathcal{P}}_{N,\pm}}{\hat{\rho}\to\gamma}$ are
compact and we get the following corollary.

######  Corollary 3.

The infima in the the Helmholtz functional $\Omega_{N,\pm}[v]$ and the
universal functional $F_{N,\pm}[\gamma]$ are achieved in the fermionic and
bosonic case, so the infima in (7) and (14) can be replaced by minima.

##  VI Conclusion

For a fixed number of particles, finite numbers of basis functions and
elevated temperature the universal functional in 1RDM functional theory is
differentiable with $\partial F_{N,\pm}/\partial\gamma=-v$ for all 1RDMs
$\gamma$ in $\mathscr{N}_{N,\pm}$. This relation holds for potential with
$\operatorname{tr}\\{v\\}=0$. However, all potentials differing from $v$ by
only a constant lead to the same 1RDM. From this relation it follows directly
that the map $v\mapsto\gamma$ is bijective up to a constant in the potential
and it gives a characterization of the set of $v$-representable 1RDMs.
Additionally, for every potential the Helmholtz functional and the universal
functional attain a minimum. One might think of reaching the $T=0$ case by
taking the limit $T\to 0$. The Gibbs state will just be an equi-ensemble of
the ground states, but a difficulty is that it changes discontinuously when
the potential is varied such that other states become ground states.
Additionally, we can no longer guarantee that the 1RDM will be in the
(relative) interior of $\overline{\mathscr{N}}_{N,\pm}$ which prevents
invertibility of $v\mapsto\gamma_{v}$ and also differentiability will probably
no longer be in the cards.

###### Acknowledgements.

The authors acknowledge support by the Netherlands Organisation for Scientific
Research (NWO) under Vici grant 724.017.001.

##  Appendix A Properties of the 1RDM

We are still left with showing that the set of 1RDMs is
$\overline{\mathscr{N}}_{N,\pm}$. To show the properties (6a) and (6b) for the
1RDM we use the following equivalent definition.

The kernel of $\gamma$ is given by tracing out $N-1$ particles in the density
matrix operator,

$\gamma(x,y)=N\int\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\rho(x,x_{2},\dotsc,x_{N};y,x_{2},\dotsc,x_{N}).$
(21)

The 1RDM can be worked out in a $1$-particle orthonormal basis
$\phi_{1},\dotsc,\phi_{N_{b}}$ for a matrix representation with elements

$\gamma_{ij}=N\int\mathrm{d}x\mathrm{d}y\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\phi_{i}^{*}(x)\phi_{j}(y)\\\
\rho(x,x_{2},\ldots,x_{N};y,x_{2},\ldots,x_{N}).$ (22)

Note that we have the following inequality,

$\gamma_{ii}=N\sum_{l}\lambda_{l}\int\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\\\
\biggl{(}\int\mathrm{d}x\phi_{i}^{*}(x)\psi_{l}(x,x_{2},\dotsc,x_{N})\biggr{)}^{2}\geq
0.$ (23)

Next we want to show that for the fermionic 1RDM the diagonal entries are
bounded from above by 1. For this, note that the entries of a statistical
ensemble are bounded by the maximum value of the eigenstates of the
corresponding density-matrix operator, i.e.,

$\displaystyle\gamma_{ij}$
$\displaystyle=N\sum_{l}\lambda_{l}\int\mathrm{d}x\mathrm{d}y\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\phi_{i}^{*}(x)\phi_{j}(y)$
$\displaystyle\qquad\psi_{l}(x,x_{2},\dotsc,x_{N})\psi_{l}^{*}(y,x_{2},\dotsc,x_{N})$
$\displaystyle\leq\max_{l}N\int\mathrm{d}x\mathrm{d}y\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\phi_{i}^{*}(x)\phi_{j}(y)$
$\displaystyle\qquad\psi_{l}(x,x_{2},\dotsc,x_{N})\psi_{l}^{*}(y,x_{2},\dotsc,x_{N})$
$\displaystyle=\max\limits_{l}(\gamma_{\psi_{l}})_{ij},$ (24)

where $\gamma_{\psi_{l}}$ means the 1RDM generated from the wave function
$\psi_{l}$. Thus, we need to show the desired upper bound only for pure
states. With a similar argument, it suffices to show the bound for Slater
determinants built from $1$-particle orthonormal states $f_{1},\dotsc,f_{N}$,

$\psi(x_{1},\dotsc,x_{N})\\\ =\frac{1}{\sqrt{N!}}\sum_{\sigma\in
S_{N}}(-1)^{\operatorname{sgn}(\sigma)}f_{\sigma(1)}(x_{1})\dotsb
f_{\sigma(N)}(x_{N}).$ (25)

The kernel of the 1RDM can be worked out as

$\gamma_{\psi}(x,y)=\sum\limits_{j=1}^{N}f_{j}(x)f_{j}^{*}(y)$ (26)

and its diagonal elements are

$(\gamma_{\psi})_{ii}=\sum\limits_{j=1}^{N}|\langle
f_{j}|\phi_{i}\rangle|^{2}\leq\langle\phi_{i}|\phi_{i}\rangle=1.$ (27)

It is easy to see that the 1RDM is Hermitian. Thus, it has a spectral
decomposition

$\gamma=\sum_{l=1}^{N_{b}}\lambda_{l}\ket{\varphi_{l}}\bra{\varphi_{l}}.$ (28)

Note that since the diagonal elements are non-negative for any basis, it
follows that the eigenvalues $\lambda_{l}$ are non-negative. Therefore,
$\gamma\geq 0$ and for fermions we have additionally $\gamma\leq 1$. The trace
of $\gamma$ can be calculated through its integral kernel (21),

$\displaystyle\operatorname{tr}\\{\gamma\\}$
$\displaystyle=\int\mathrm{d}x\gamma(x,x)$
$\displaystyle=N\sum_{l}\lambda_{l}\int\mathrm{d}x\mathrm{d}x_{2}\dotsi\mathrm{d}x_{N}\lvert\psi_{l}(x,x_{2},\dotsc
x_{N})\rvert^{2}$ $\displaystyle=N.$ (29)

All these properties together show that the set of 1RDMs is contained in
$\overline{\mathscr{N}}_{N,\pm}$.

Next, we want to prove Theorem II.1.

###### Proof.

For $N=1$ we can simply take $\hat{\rho}=\gamma$. So let us consider the case
$N\geq 2$. We represent $\gamma$ in the NO basis
$\gamma=\sum_{j=1}^{N_{b}}\lambda_{j}\ket{\varphi_{j}}\bra{\varphi_{j}}$. We
need to distinguish between bosons and fermions.

Bosonic case

We define the $N$-particle wave function

$\psi(x_{1},\ldots,x_{N})\coloneqq\frac{1}{\sqrt{N}}\sum_{j=1}^{N_{b}}\lambda_{j}^{1/2}\prod_{i=1}^{N}\varphi_{j}(x_{i}).$

It is now easy to see that $\psi$ is symmetric, normalized and that it
generates $\gamma$.

Fermionic case

We work with a polytope. The 1RDM $\gamma$ can be expressed as a vector of
length $N_{b}$ containing its occupation numbers
$\textbf{n}=(\lambda_{1},\dotsc,\lambda_{N_{b}})$. The extreme points of the
polytope are all possible permutations of $N$ occupation numbers set to one
and all other set to zero

$\overline{\gamma}_{I}\coloneqq\overline{\gamma}_{i_{1}\ldots
i_{N}}\coloneqq\textbf{e}_{i_{1}}+\dotsb+\textbf{e}_{i_{N}},$

for $1\leq i_{1}<\ldots<i_{N}\leq N_{b}$ and where the $\textbf{e}_{i}$’s are
unit vectors. The index $I$ is a renumeration of $i_{1}\ldots i_{N}$ and can
take $K=\binom{N_{b}}{N}$ values. The vector n is an element of the polytope

$\Gamma\coloneqq\Set{\sum_{I=1}^{K}\mu_{I}\overline{\gamma}_{I}}{\mu_{I}\geq
0,\;\sum\limits_{I=1}^{K}\mu_{I}=1}.$

The extreme points $\overline{\gamma}_{i_{1}\dotsc i_{N}}$ can now be
identified with
$\ket{\varphi_{i_{1}}\dotsc\varphi_{i_{N}}}\bra{\varphi_{i_{1}}\dotsc\varphi_{i_{N}}}$.
Since the mapping $\hat{\rho}\to\gamma$ is linear we find that $\gamma$ is
generated from a linear combination of the Slater determinants
$\ket{\varphi_{i_{1}}\dotsc\varphi_{i_{N}}}\bra{\varphi_{i_{1}}\dotsc\varphi_{i_{N}}}$.
∎

## References

  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Physical review 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Physical review 140, A1133 (1965).
  * Cohen  _et al._ [2012] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Challenges for density functional theory, Chemical reviews 112, 289 (2012).
  * Becke [2014] A. Becke, Perspective: Fifty years of density-functional theory in chemical physics, J. Chem. Phys. 140, 18A301 (2014).
  * Vuckovic  _et al._ [2015] S. Vuckovic, L. O. Wagner, A. Mirtschink, and P. Gori-Giorgi, Hydrogen molecule dissociation curve with functionals based on the strictly correlated regime, J. Chem. Theory Comput. 11, 3153 (2015).
  * Gilbert [1975] T. L. Gilbert, Hohenberg-kohn theorem for nonlocal external potentials, Physical Review B 12, 2111 (1975).
  * Pernal [2005] K. Pernal, Effective potential for natural spin orbitals, Phys. Rev. Lett. 94, 233002 (2005).
  * van Leeuwen [2007] R. van Leeuwen, Functionals of the one-particle density matrix: Uniqueness theorems and $v$-representability in equilibrium and time-dependent systems (2007).
  * Baldsiefen [2012] T. Baldsiefen,  _Reduced Density Matrix Functional Theory at Finite Temperature_ , Ph.D. thesis, Institut für Theoretische Physik Freie Universität Berlin (2012).
  * Giesbertz and Ruggenthaler [2019] K. J. Giesbertz and M. Ruggenthaler, One-body reduced density-matrix functional theory in finite basis sets at elevated temperatures, Physics Reports 806, 1 (2019).
  * Requist and Pankratov [2008] R. Requist and O. Pankratov, Generalized Kohn-Sham system in one-matrix functional theory, Phys. Rev. B 77, 235121 (2008).
  * Giesbertz and Baerends [2010] K. J. H. Giesbertz and E. J. Baerends, Aufbau derived from a unified treatment of occupation numbers in Hartree–Fock, Kohn–Sham, and natural orbital theories with the Karush–Kuhn–Tucker conditions for the inequality constraints $n_{i}\leq 1$ and $n_{i}\geq 0$, J. Chem. Phys. 132, 194108 (2010).
  * Baldsiefen and Gross [2013] T. Baldsiefen and E. Gross, Minimization procedure in reduced density matrix functional theory by means of an effective noninteracting system, Computational and Theoretical Chemistry 1003, 114 (2013).
  * Baldsiefen  _et al._ [2015] T. Baldsiefen, A. Cangi, and E. K. U. Gross, Reduced-density-matrix-functional theory at finite temperature: Theoretical foundations, Physical Review A 92, 052514 (2015).
  * Yoo  _et al._ [2005] C. S. Yoo, B. Maddox, J. H. P. Klepeis, V. Iota, W. Evans, A. McMahan, M. Y. Hu, P. Chow, M. Somayazulu, D. Häusermann, R. T. Scalettar, and W. E. Pickett, First-order isostructural Mott transition in highly compressed MnO, Phys. Rev. Lett. 94, 115502 (2005).
  * Nagamatsu  _et al._ [2001] J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani, and J. Akimitsu, Superconductivity at 39 K in magnesium diboride, Nature 410, 63 (2001).
  * Dharma-wardana and Perrot [1982] M. W. C. Dharma-wardana and F. Perrot, Density-functional theory of hydrogen plasmas, Phys. Rev. A 26, 2096 (1982).
  * Ewald  _et al._ [2019] N. V. Ewald, T. Feldker, H. Hirzler, H. A. Fürst, and R. Gerritsma, Observation of interactions between trapped ions and ultracold rydberg atoms, Phys. Rev. Lett. 122, 253401 (2019).
  * Bedingham [2003] D. J. Bedingham, Bose-einstein condensation in the canonical ensemble, Phys. Rev. D 68, 105007 (2003).
  * Mullin and Fernández [2003] W. J. Mullin and J. P. Fernández, Bose–Einstein condensation, fluctuations, and recurrence relations in statistical mechanics, Am. J. Phys. 71, 661 (2003).
  * Barghathi  _et al._ [2020] H. Barghathi, J. Yu, and A. Del Maestro, Theory of noninteracting fermions and bosons in the canonical ensemble, Physical Review Research 2, 043206 (2020).
  * Kooi [2022] D. P. Kooi, Efficient bosonic and fermionic sinkhorn algorithms for non-interacting ensembles in one-body reduced density matrix functional theory in the canonical ensemble, arXiv:2205.15058 [physics.chem-ph] (2022).
  * Note [1] Creation and annihilation operator are only defined pairwise for a $N$ particle space.
  * Coleman [1963] A. J. Coleman, Structure of fermion density matrices, Reviews of modern Physics 35, 668 (1963).
  * Eschrig [2010] H. Eschrig, T$>0$ ensemble-state density functional theory via legendre transform, Physical Review B 82, 205120 (2010).
  * Mermin [1965] N. D. Mermin, Thermal properties of the inhomogeneous electron gas, Physical Review 137, A1441 (1965).
  * Ruelle [1969] D. Ruelle,  _Statistical Mechanics_ , Mathematical Physics Monograph Series (W.A. Benjamin, Inc., 1969).
  * Lieb [1975] E. H. Lieb, Some convexity and subadditivity properties of entropy, Bull. Amer. Math. Soc 81, 1 (1975).
  * Wehrl [1978] A. Wehrl, General properties of entropy, Reviews of Modern Physics 50, 221 (1978).
  * Note [2] This means $M$ is contained in an affine space.
  * Note [3] We could have taken a different point in $\mathscr{N}_{N,\pm}$, but the current choice is the most natural in this setting.


# Fueling limits in a cylindrical viscosity-limited reactor

T. Rubin<EMAIL_ADDRESS>Department of Astrophysical Sciences, Princeton
University, Princeton, New Jersey 08544, USA  E. J. Kolmes  Department of
Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544, USA
I. E. Ochs  Department of Astrophysical Sciences, Princeton University,
Princeton, New Jersey 08544, USA  M. E. Mlodik  Department of Astrophysical
Sciences, Princeton University, Princeton, New Jersey 08544, USA  N. J. Fisch
Department of Astrophysical Sciences, Princeton University, Princeton, New
Jersey 08544, USA

###### Abstract

Recently, a method to achieve a “natural hot-ion mode” was suggested, by
utilizing ion viscous heating in a rotating plasma with a fixed boundary. We
explore the steady-state solution to the Braginskii equations and find the
parameter regime in which a significant temperature difference between ions
and electrons can be sustained in a driven steady state. The threshold for
this effect occurs at $\rho_{i}\gtrsim 0.1R$. An analytic, leading order low
flow solution is obtained, and a numerical, moderate Mach number $M\lesssim 2$
is investigated. The limitation is found to be at moderate Mach numbers.

††preprint:

##  I Introduction

Magnetic plasma confinement assisted by rotation has been explored in several
configurations, such as mirrors [1, 2, 3] and toroidal devices[4, 5]. Rotating
mirrors in particular are receiving renewed interest, leading to new
experimental devices in the near future[6, 7, 8]. Plasma mass filters[9, 10,
11, 12, 13, 14, 15, 16, 17] are another rotating plasma application in which
density gradients are of particular importance, and are similar to rotating
mirrors in many respects.

Sufficiently long mirrors may be analyzed using classical transport theory.
Radial cross field ion currents [18, 19] in such devices appear to be an
attractive fueling method, as they can induce rotation in the plasma due to
their interaction with the magnetic field. The hydrodynamic variables -
densities, momenta and pressures - in such configurations can be
asymptotically solved for[20] or numerically integrated using a variety of
tools, such as the MITNS: Multiple-Ion Transport Numerical Solver code[21].

Nuclear fusion in magnetic devices is realized by confinement of hot ions for
sufficient time[22]. In these devices, the plasma confinement is often limited
by the total plasma pressure, which sums the electron and ion pressures. As
such, a hot-ion mode is preferable, as it can produce more fusion power for
the same magnetic field strength[23], in addition to a decrease in energy
radiation losses through electrons.

Plasma heating can be accomplished using a variety of methods. One proposed
method is to use viscous heating due to sheared rotation. Kolmes  _et al._
[24] showed how heat dissipation channels in an axisymmetric cylindrical
plasma could preferentially heat the ion population. Such a configuration may
be realized by a radial flow of fuel ions into the hot center of the cylinder,
and the removal of ash ions by a fast process other than classical transport -
such as $\alpha$-channeling[25, 26]. The heating is the result of the viscous
dissipation of the ion fluid, i.e., the rate of work done by the viscous
stress times strain-rate. The electron fluid viscous stress, and the resultant
heating rate, is smaller by a factor of $(m_{e}/m_{i})^{3/2}$.

Kolmes  _et al._ [24] predicted that the ion temperature could be higher than
the electron temperature, and suggested that the temperature difference could
be large in cases with sufficiently large radial influxes of fuel ions and
high radial electric fields. However, that paper left open the question of
what kinds of radial influxes and fields could be self-consistently supported,
and of precisely what are the density, velocity, and pressure profiles in this
case. The present paper addresses that question. One of the key results that
follows from this calculation is that there are nontrivial limitations on
these hot-ion-mode solutions.

In this paper, we explore a particular solution to the proposed concept. We
consider a one-dimensional Braginskii fluid model[27] to explore the nonlinear
effects due to the density and temperature dependence of the viscosity and
heat conductivities for a uniform volumetric charge extraction, i.e. a radial
ion current. We compare analytic solutions using constant coefficients to the
full numerical non-linear solution.

Because the Braginskii fluid model is the most usual model, it serves well our
purposes here. We do note that several authors published corrections to the
transport coefficients[28, 29, 30, 31, 32]. In addition, a magnetic mirror
machine would have a distribution function that contains voids where particles
are not confined. In the magnetic mirror case, it is not clear that
Braginskii, or other closures that expand the distribution function in a
polynomial basis, produce the precise transport coefficients. However, the
effects discussed in this work do not rely on the specific transport
coefficients.

The limitation on the viscous heating arises because the radial ion current
into the center of the cylinder leads to a large rotation of the plasma, and
the centrifugal force caused by this rotation empties out the density at the
core. The quadratic dependence of the viscosity on the density means that
large angular velocity gradients are required to produce the viscous shear
needed to balance the torque produced by the radial current and the magnetic
field - leading to progressively larger rotation. This effect limits the
amount of charge extraction possible in this configuration. The viscous
heating itself, which increases the plasma temperature, further reduces the
viscosity coefficient. Of course, diverging angular velocities are not
physical, and are the result of an attempt to balance a finite torque when the
viscosity coefficient approaches 0. In a physical system, the torque would be
limited.

As a result, the proposed “natural hot-ion mode” is limited. Beyond a certain
radial ion flux, the nonlinearities in the viscosity coefficient would cut off
the shear stress in the plasma. This limitation in shear limits the viscous
heating; the limit depends on the magnetization of the plasma. We calculate a
magnetization threshold above which the viscous heating is small compared to
the rate of temperature equilibration between species.

This paper is organized as follows: In Sec. II we present the
nondimensionalized two-fluid equations. In Sec. III we present the low-flux
approximate solution, assuming constant coefficients. In Sec. IV we discuss
the deviation of the full nonlinear (variable coefficients) solution from the
linear approximation.

##  II Model

Consider an axisymmetric, infinitely long plasma cylinder, in equilibrium,
such that $\frac{\partial}{\partial t}=\frac{\partial}{\partial
z}=\frac{\partial}{\partial\theta}=0$, with a constant axial magnetic field.
The plasma species to be considered, for simplicity, are fuel ions and
electrons. This setting might be realizable in a steady-state rotating mirror
machine, which is fueled radially rather than axially, and in which fusion ash
is removed quickly radially using $\alpha$-channeling, before it can interact
with the fuel ions or the electrons. It is assumed the $\alpha$-channeling
does not affect electrons.

In steady state, fuel ions fuse, the ash is removed, and more fuel is supplied
continuously from the outer edge, while the electrons have no average radial
velocity. The ion sink term produces in steady state an inward-flowing ion
current. The $\mathbf{r}\times(\mathbf{j}\times\mathbf{B})$ torque due to this
current induces rotation in the plasma. The plasma rotation, in addition to
the other radial forces acting on it, determines the steady-state density
profile.

The radial expulsion of the fusion ash produces the opposite
$\mathbf{r}\times(\mathbf{j}\times\mathbf{B})$ torque on the ash ions.
However, the fusion ash is kept at such a low density, using
$\alpha$-channeling, that its collisional interactions with other plasma
constituents can be ordered out of the momentum and energy equations. A
density ratio of $n_{a}/n_{i}\sim\rho_{*}^{2}$ is sufficient, with the small
parameter $\rho_{*}$ being the normalized ion Larmor radius, also defined
later.

In this work, we consider only classical transport effects in order to
determine the temperatures. In many real-life plasmas, other effects
contribute to the energy balance - examples include radiative cooling, and RF
heating. These, and energy exchange with fusion ash may even be the dominant
mechanisms, over and above classical transport effects. However, one purpose
of this paper is to determine the merit of the proposed natural ion mode,
which cannot be separated from the other effects considered here

For each fluid, the continuity, radial momentum, angular momentum, and
pressure equations are,

$\displaystyle\frac{\partial n_{s}}{\partial
t}+\frac{1}{r}\frac{\partial}{\partial r}rn_{s}v_{rs}=s_{s},~{}$ (1)
$\frac{\partial}{\partial
t}(m_{s}n_{s}v_{rs})+\frac{1}{r}\frac{\partial}{\partial
r}rm_{s}n_{s}v_{rs}^{2}+\frac{\partial p_{s}}{\partial
r}+(\nabla\cdot\pi_{s})_{r}=\\\
s_{s}m_{s}v_{sr}^{src}+m_{s}n_{s}r\omega_{s}^{2}+Z_{s}en_{s}\left(E_{r}+r\omega_{s}B_{z}\right)+\sum_{s^{\prime}}R_{ss^{\prime}r},~{}$
(2) $\frac{\partial}{\partial
t}(m_{s}n_{s}r^{2}\omega_{s})+\frac{1}{r}\frac{\partial}{\partial
r}(r^{3}m_{s}n_{s}\omega_{s}v_{rs})=rs_{s}m_{s}v_{s\theta}^{src}\\\
-Z_{s}ern_{s}v_{rs}B_{z}-r(\nabla\cdot\pi_{s})_{\theta}+r\sum_{s^{\prime}}\left(R_{ss^{\prime}\theta}+f_{ss^{\prime}\theta}\right),~{}$
(3) $\frac{3}{2}\frac{\partial p_{s}}{\partial
t}+\frac{1}{r}\frac{\partial}{\partial
r}r\left(q_{rs}+\frac{5}{2}p_{s}v_{rs}\right)+\pi_{s}:\nabla\mathbf{v}_{s}\\\
=v_{rs}\frac{\partial p_{s}}{\partial
r}+\sum_{s^{\prime}}3\frac{m_{s}n_{s}\nu_{ss^{\prime}}}{m_{s}+m_{s^{\prime}}}(T_{s^{\prime}}-T_{s})+\frac{3}{2}s_{s}T^{src}_{s}\\\
+\frac{1}{2}m_{s}s_{s}\left[\left(v^{src}_{rs}-v_{rs}\right)^{2}+\left(v^{src}_{\theta
s}-v_{\theta s}\right)^{2}\right]\\\
+\sum_{s^{\prime}}\Bigg{[}\frac{m_{s^{\prime}}}{m_{s}+m_{s^{\prime}}}(\mathbf{v}_{s^{\prime}}-\mathbf{v}_{s})\cdot\left(\mathbf{R}_{ss^{\prime}}+\mathbf{f}_{ss^{\prime}}\right)\Bigg{]}.~{}$
(4)

Quantities with a subscript $s$ represent a species-dependent quantity, with
the index $s$ representing ions, $s=i$, or electrons, $s=e$.

Number density, velocity, pressure, and temperature are denoted by $n$, $v$,
$p$ and $T$ respectively, with $p_{s}=n_{s}T_{s}$. The quantity $s_{s}$ is a
particle source or sink for species $s$. Time and radius (spatial coordinate)
are denoted by $t$ and $r$. Particle mass and charge number are denoted by $m$
and $Z$, while the elementary charge is denoted by $e$. The (radial) electric
field and (axial) magnetic field are denoted by $E$ and $B$. The constant
magnetic field assumption can be construed to stem from a low plasma
$\beta\doteq\frac{2\mu_{0}p}{B^{2}}$, with $\mu_{0}$ being the permeability of
vacuum. The value used for the numerical simulations of the nonlinear
equations is $\beta=0.002$.

If the particle source is negative, as a sink term, the source temperature
$T_{s}^{src}=T_{s}$ and velocity $\mathbf{v}_{s}^{src}=\mathbf{v}_{s}$. If the
source term is positive, $T_{s}^{src}$ and $\mathbf{v}_{s}^{src}$ need to be
specified.

The friction body force $\mathbf{R}_{ss^{\prime}}$ and the thermal friction
(“Nernst”) body force $\mathbf{f}_{ss^{\prime}}$ between species $s$ and
$s^{\prime}$ are expressed as

$\displaystyle\mathbf{R}_{ss^{\prime}}$
$\displaystyle=m_{s}n_{s}\nu_{ss^{\prime}}(\mathbf{v}_{s^{\prime}}-\mathbf{v}_{s}),$
(5) $\displaystyle\mathbf{f}_{ss^{\prime}}$
$\displaystyle=\frac{3}{2}\frac{m_{s}n_{s}\nu_{ss^{\prime}}}{Z_{s}Z_{s^{\prime}}eB}\hat{b}\times\frac{Z_{s^{\prime}}m_{s^{\prime}}T_{s}\nabla
T_{s}-Z_{s}m_{s}T_{s^{\prime}}\nabla
T_{s^{\prime}}}{m_{s}T_{s^{\prime}}+m_{s^{\prime}}T_{s}},$ (6)

where $\hat{b}$ is a unit vector in the direction of the magnetic field. In
the case of an electron-ion plasma, the thermal friction can be approximated
as

$\displaystyle\mathbf{f}_{ie}=\frac{3}{2}\frac{Z_{i}n_{i}\nu_{ie}}{\Omega_{i}}\frac{\partial
T_{e}}{\partial r}\hat{\theta}.$ (7)

The resultant viscous stress tensor divergence, using the Braginskii [27]
closure, is [24]

$\displaystyle\nabla\cdot\pi_{s}=$ $\displaystyle\frac{\partial}{\partial
r}\left[\frac{\eta_{0s}}{3r}\frac{\partial}{\partial
r}\left(rv_{rs}\right)\right]\hat{r}$ $\displaystyle-\frac{1}{r^{2}}$
$\displaystyle\frac{\partial}{\partial
r}r^{3}\left[\eta_{1s}\frac{\partial}{\partial
r}\left(\frac{v_{rs}}{r}\right)+\eta_{3s}\frac{\partial\omega_{s}}{\partial
r}\right]\hat{r}$ $\displaystyle-\frac{1}{r^{2}}$
$\displaystyle\frac{\partial}{\partial
r}r^{3}\left[\eta_{1s}\frac{\partial\omega_{s}}{\partial
r}-\eta_{3s}\frac{\partial}{\partial
r}\left(\frac{v_{rs}}{r}\right)\right]\hat{\theta}.$ (8)

The viscosity coefficients for an electron-ion plasma, with $Z=1$ are

$\displaystyle\eta_{0i}=0.96\sqrt{2}\frac{p_{i}}{\nu_{ii}},$ (9)
$\displaystyle\eta_{1i}=\frac{3}{10\sqrt{2}}\frac{p_{i}\nu_{ii}}{\Omega_{i}^{2}},$
(10) $\displaystyle\eta_{3i}=\frac{p_{i}}{2\Omega_{i}},$ (11)
$\displaystyle\eta_{0e}=0.73\frac{p_{e}}{\nu_{ei}},$ (12)
$\displaystyle\eta_{1e}=0.51\frac{p_{e}\nu_{ei}}{\Omega_{e}^{2}},$ (13)
$\displaystyle\eta_{3e}=\frac{p_{e}}{2\Omega_{e}},$ (14)

where $\Omega_{s}$ denotes the signed Larmor frequency for species $s$. The
viscous heating term is expressed as

$\pi_{s}:\nabla\mathbf{v}_{s}=-\frac{\eta_{0s}}{3}\left(\frac{1}{r}\frac{\partial}{\partial
r}\left(rv_{rs}\right)\right)^{2}\\\
-\eta_{1s}\left(\left(r\frac{\partial}{\partial
r}\left(\frac{v_{rs}}{r}\right)\right)^{2}+\left(r\frac{\partial\omega_{s}}{\partial
r}\right)^{2}\right).$ (15)

The heat conduction $q_{r}$ is

$\displaystyle q_{rs}=-\kappa_{s}\frac{\partial T_{s}}{\partial
r}-\kappa_{ss^{\prime}}^{\omega}r\left(\omega_{s^{\prime}}-\omega_{s}\right),$
(16)

with heat diffusivities,

$\displaystyle\kappa_{i}=\sqrt{2}\frac{p_{i}\nu_{ii}}{m_{i}\Omega_{i}^{2}},$
(17) $\displaystyle\kappa_{e}=4.66\frac{p_{e}\nu_{ei}}{m_{e}\Omega_{e}^{2}},$
(18) $\displaystyle\kappa_{ie}^{\omega}\approx 0,$ (19)
$\displaystyle\kappa_{ei}^{\omega}=\frac{3}{2}\frac{p_{e}\nu_{ei}}{\Omega_{e}}.$
(20)

The Coulomb collision frequency between particles of species $s$ and
$s^{\prime}$ is expressed as[33]

$\nu_{ss^{\prime}}=\frac{\sqrt{2}e^{4}\log\Lambda_{ss^{\prime}}}{12\pi^{3/2}\varepsilon_{0}^{2}}\cdot\\\
Z_{s}^{2}Z_{s^{\prime}}^{2}\sqrt{\frac{m_{s^{\prime}}}{m_{s}}\frac{1}{m_{s}+m_{s^{\prime}}}}\left(\frac{m_{s}+m_{s^{\prime}}}{m_{s}T_{s^{\prime}}+m_{s^{\prime}}T_{s}}\right)^{3/2}n_{s^{\prime}}.$
(21)

with $\log\Lambda_{ss^{\prime}}$ being the Coulomb logarithm and
$\varepsilon_{0}$ the permittivity of free space. Note that this definition of
collision frequencies for like-species collisions differs from its expression
in Braginskii by a factor of $1/\sqrt{2}$. The corresponding transport
coefficients incorporate a $\sqrt{2}$ factor, so the transport itself is
unchanged.

In steady state, the particle flux

$\displaystyle\Gamma_{s}\doteq
rn_{s}v_{rs}=\Gamma_{s}(R)+\int_{R}^{r}rs_{s}\mathrm{d}r,$ (22)

is the quantity driving the system. Here we use $R$ as the outer radius of the
cylinder.

###  II.1 Source terms and boundary conditions

The source terms and fluxes for the ions and electrons,

$\displaystyle\Gamma_{e}\equiv 0,$ (23) $\displaystyle
s_{i}=\mathrm{const.}<0,$ (24) $\displaystyle\Gamma_{i}=\frac{s_{i}}{2}r^{2},$
(25)

correspond to a constant ion sink $s_{i}<0$, understood to be a fusion
process, and no radial motion of electrons.

We choose to use a uniform particle sink because we would like to compare a
constant coefficient solution to the same case with variable coefficients. In
the comparison, we want to be clear that the nonlinear behavior originates in
the nonlinear transport coefficients rather than from the nonlinear source
term.

The boundary conditions for the plasma are:

$\displaystyle\omega_{i}(R)=0,$ (26) $\displaystyle\omega_{e}(R)=0,~{}$ (27)
$\displaystyle Z_{i}n_{i}(R)=n_{e}(R)=Z_{i}n_{0},$ (28) $\displaystyle
T_{i}(R)=T_{e}(R)=T_{0},$ (29)
$\displaystyle\omega_{i}^{\prime}(0)=\omega_{e}^{\prime}(0)=0,$ (30)
$\displaystyle n_{i}^{\prime}(0)=n_{e}^{\prime}(0)=0,$ (31) $\displaystyle
T_{i}^{\prime}(0)=T_{e}^{\prime}(0)=0.$ (32)

The integrated viscous heating of the plasma, to leading order[24], is

$\displaystyle\int_{0}^{R}\pi:\nabla\mathbf{v}r\mathrm{d}r=\int_{0}^{R}v_{\theta}(\nabla\cdot\pi)_{\theta}r\mathrm{d}r+\left[rv_{\theta}\pi_{r\theta}\right]_{r=R}.$
(33)

The choice of the boundary conditions (26) and (27) drops the boundary term in
(33), which is the work done by the boundary on the plasma. Boundary
conditions (28) and (29) are a choice of normalization, where the equality of
the ion and electron temperatures might require further justification (see
appendix). The boundary conditions (30), (31), and (32) are the result of the
cylindrical geometry.

###  II.2 Nondimensionalization

Nondimensionalizing the equations of motion allows us to factor out small
parameters for use in an asymptotic expansion. Denoting $X=X_{0}\tilde{X}$,
with $X_{0}$ being a reference quantity:

$\displaystyle m_{0}\doteq m_{p},$ (34) $\displaystyle
v_{0}\doteq\sqrt{T_{0}/m_{p}},$ (35) $\displaystyle r_{0}\doteq R,$ (36)
$\displaystyle\Gamma_{0}\doteq Rn_{0}v_{0},$ (37)
$\displaystyle\nu_{0}\doteq\frac{\sqrt{2}e^{4}\log\Lambda}{12\pi^{3/2}\epsilon_{0}}\frac{n_{0}}{T_{0}^{3/2}m_{0}^{1/2}}$
(38) $\displaystyle\eta_{00}\doteq\frac{n_{0}T_{0}}{\nu_{0}}$ (39)
$\displaystyle\eta_{10}\doteq\frac{n_{0}T_{0}\nu_{0}}{\Omega_{p0}^{2}}=\epsilon^{2}\eta_{00},$
(40)
$\displaystyle\eta_{30}\doteq\frac{n_{0}T_{0}}{\Omega_{p0}}=\epsilon\eta_{00},$
(41) $\displaystyle R_{ss^{\prime}0}\doteq m_{0}n_{0}\nu_{0}v_{0},$ (42)
$\displaystyle
f_{ss^{\prime}0}\doteq\frac{n_{0}\nu_{0}T_{0}}{\Omega_{p0}R}=m_{0}n_{0}\nu_{0}v_{0}\rho_{*},$
(43)
$\displaystyle\kappa_{s0}\doteq\frac{n_{0}T_{0}\nu_{0}}{m_{0}\Omega_{p0}^{2}}=\frac{\epsilon^{2}}{m_{0}}\eta_{00},$
(44)
$\displaystyle\kappa_{ss^{\prime}0}^{\omega}\doteq\frac{n_{0}T_{0}\nu_{0}}{\Omega_{p0}}.$
(45)

In this paper, reference quantities are chosen as their value at the outer
radius.

This fluid closure features two small parameters, the normalized ion Larmor
radius, $\rho_{*}\doteq v_{0}/\Omega_{p0}R$, and the ratio of collision
frequency to Larmor frequency $\epsilon\doteq\nu_{0}/\Omega_{p0}$, which is
the inverse Hall parameter $\epsilon=1/C_{H}$. We use the reference quantities
to define the values of these constants, i.e., we use them as constants rather
than as functions of radius. When dealing with electron-ion plasma, a third
small parameter is present, the square root of the electron-to-proton mass
ratio, $\sqrt{\tilde{m}_{e}}$. An asymptotic expansion for a parameters
$\tilde{X}$, in powers of the small parameters $\rho_{*}$ and $\epsilon$ would
be denoted by
$\tilde{X}=\sum_{\alpha,\beta}\rho_{*}^{\alpha}\epsilon^{\beta}\tilde{X}^{(\alpha,\beta)}$.

The steady-state dimensionless angular momentum equation for a single fluid
species is,

$\tilde{m}_{s}\frac{d}{d\tilde{r}}\tilde{\Gamma}_{s}\tilde{r}^{2}\tilde{\omega}_{s}+\frac{1}{\rho_{*}}Z_{s}\tilde{\Gamma}_{s}\tilde{B}_{z}\tilde{r}=\tilde{r}^{2}\tilde{s}_{s}\tilde{m}_{s}\tilde{\omega}_{s}\\\
\rho_{*}\frac{d}{d\tilde{r}}\tilde{r}^{3}\left[\epsilon\tilde{\eta}_{1s}\frac{d\tilde{\omega}_{s}}{d\tilde{r}}-\frac{\tilde{\eta}_{s3}\tilde{\Gamma}_{s}}{\tilde{r}^{2}\tilde{n}_{s}}\left(\frac{d\ln(\tilde{\Gamma}_{s})}{d\tilde{r}}-\frac{2}{\tilde{r}}-\frac{d\ln(\tilde{n}_{s})}{d\tilde{r}}\right)\right]\\\
+\epsilon\tilde{r}^{2}\sum_{s^{\prime}}\left(\frac{1}{\rho_{*}}\tilde{R}_{ss^{\prime}\theta}+\tilde{f}_{ss^{\prime}\theta}\right).~{}$
(46)

The steady-state dimensionless radial force-balance equation for a single
fluid species is,

$\frac{d\tilde{p}_{s}}{d\tilde{r}}=\frac{1}{\rho_{*}}Z_{s}\tilde{n}_{s}\left(\tilde{E}_{r}+\tilde{r}\tilde{\omega}_{s}\tilde{B}\right)+\frac{\epsilon}{\rho_{*}}\sum_{s^{\prime}}\tilde{R}_{ss^{\prime}r}\\\
+\tilde{m}_{s}\tilde{s}_{s}\frac{\tilde{\Gamma}_{s}}{\tilde{r}\tilde{n}_{s}}+\tilde{m}_{s}\tilde{n}_{s}\tilde{r}\tilde{\omega}_{s}^{2}+\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\left(\frac{\tilde{m}_{s}\tilde{\Gamma}_{s}^{2}}{\tilde{r}\tilde{n}_{s}}\right)\\\
+\frac{\rho_{*}}{\epsilon}\frac{d}{d\tilde{r}}\left(\frac{\tilde{\eta}_{s0}}{3\tilde{r}}\frac{d}{d\tilde{r}}\left(\frac{\tilde{\Gamma}_{s}}{\tilde{n}_{s}}\right)\right)\\\
+\frac{\rho_{*}}{\tilde{r}^{2}}\frac{d}{d\tilde{r}}\tilde{r}^{3}\left(\epsilon\tilde{\eta}_{s1}\frac{d}{d\tilde{r}}\left(\frac{\tilde{\Gamma}_{s}}{\tilde{r}^{2}\tilde{n}_{s}}\right)+\tilde{\eta}_{s3}\frac{d\tilde{\omega}_{s}}{d\tilde{r}}\right).$
(47)

The steady-state dimensionless temperature equation for a single fluid species
is,

$-\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\tilde{r}\left[\tilde{\kappa}_{s}\frac{d\tilde{T}_{s}}{d\tilde{r}}+\frac{1}{\rho_{*}}\tilde{\kappa}_{ss^{\prime}}^{\omega}\tilde{r}\left(\tilde{\omega}_{s^{\prime}}-\tilde{\omega}_{s}\right)\right]+\frac{5}{2}\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\tilde{T}_{s}\frac{\tilde{\Gamma}_{s}}{\rho_{*}\epsilon}\\\
=\frac{d\tilde{p}_{s}}{d\tilde{r}}\frac{\tilde{\Gamma}_{s}}{\rho_{*}\epsilon\tilde{r}\tilde{n}_{s}}+\frac{1}{\rho_{*}^{2}}\sum_{s^{\prime}}3\frac{\tilde{m}_{s}\tilde{n}_{s}\tilde{\nu}_{ss^{\prime}}}{\tilde{m}_{s}+\tilde{m}_{s^{\prime}}}(\tilde{T}_{s^{\prime}}-\tilde{T}_{s})+\frac{3}{2}\frac{\tilde{s}_{s}\tilde{T}_{s}}{\rho_{*}\epsilon}\\\
+\sum_{s^{\prime}}\Bigg{[}\frac{\tilde{m}_{s^{\prime}}}{\tilde{m}_{s}+\tilde{m}_{s^{\prime}}}\frac{1}{\rho_{*}}(\tilde{\mathbf{v}}_{s^{\prime}}-\tilde{\mathbf{v}}_{s})\cdot\left(\frac{1}{\rho_{*}}\tilde{\mathbf{R}}_{ss^{\prime}}+\tilde{\mathbf{f}}_{ss^{\prime}}\right)\Bigg{]}\\\
+\frac{1}{\epsilon^{2}}\frac{\tilde{\eta}_{0}}{3}\left[\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\frac{\tilde{\Gamma}_{s}}{\tilde{n}_{s}}\right]^{2}+\tilde{\eta}_{1}\tilde{r}^{2}\left[\left[\frac{d}{d\tilde{r}}\left(\frac{\tilde{\Gamma}_{s}}{\tilde{r}^{2}\tilde{n}_{s}}\right)\right]^{2}+\left[\frac{d\tilde{\omega}_{s}}{d\tilde{r}}\right]^{2}\right]~{}$
(48)

##  III Constant coefficients solution

In this section, the leading order solution for the angular velocities,
density and temperatures are derived, assuming $\eta_{1}$, $\kappa$ and
$\nu_{ie}$ are constants, and do not depend on the variations in density or
temperature within the domain. This is an approximation that holds well for
slow rotation ($\tilde{v}_{\theta}\ll 1$), when the density and temperature
are indeed nearly uniform.

###  III.1 Electron angular velocity

For $\Gamma_{e}\equiv 0$ the electron angular velocity equation (46) reads,

$\displaystyle\rho_{*}\frac{d}{d\tilde{r}}\left(\tilde{r}^{3}\tilde{\eta}_{1e}\frac{d\tilde{\omega}_{e}}{d\tilde{r}}\right)=\frac{1}{\rho_{*}}\tilde{m}_{i}\tilde{n}_{i}\tilde{\nu}_{ie}\tilde{r}^{3}(\tilde{\omega}_{e}-\tilde{\omega}_{i})+\tilde{r}^{2}\tilde{f}_{ie\theta}.$
(49)

The electron angular velocity can be solved asymptotically, to leading and
first order in $\rho_{*}$,

$\displaystyle\tilde{\omega}_{e}=\tilde{\omega}_{i}-\frac{3}{2}\frac{\rho_{*}}{\tilde{r}\tilde{B}}\frac{d\tilde{T}_{e}}{d\tilde{r}}\left(1-\frac{I_{1}\left(\frac{\tilde{r}}{\rho_{*}}\sqrt{\frac{\tilde{m}_{e}\tilde{n}_{e}\tilde{\nu}_{ei}}{\tilde{\eta}_{1e}}}\right)}{I_{1}\left(\frac{1}{\rho_{*}}\sqrt{\frac{\tilde{m}_{e}\tilde{n}_{e}\tilde{\nu}_{ei}}{\tilde{\eta}_{1e}}}\Big{|}_{\tilde{r}=1}\right)}\right)~{},$
(50)

where $I_{1}$ is the modified Bessel function of the first kind. This solution
sets the azimuthal component of the total friction force
$R_{ie\theta}+f_{ie\theta}=0$, except for a boundary layer at
$\tilde{r}\approx 1$.

The electron viscosity coefficient is smaller by a factor of
$\tilde{m}_{e}^{3/2}$ relative to the ion viscosity coefficient, and the
contribution of the viscosity is of $\mathcal{O}(\rho_{*}^{2}\tilde{m}_{e})$
relative to the ion angular velocity. The ion angular velocity equation can be
exactly solved if the sink term
$\tilde{r}^{2}\tilde{s}_{s}\tilde{m}_{s}\tilde{\omega}_{s}$ is dropped. This
term would turn out to be of $\mathcal{O}(\rho_{*}^{2}\epsilon)$ later. The
exact solution to $\tilde{\omega}_{i}$ (as a function of $\tilde{\omega}_{e}$)
without the sink term, with boundary conditions $\tilde{\omega}_{i}(1)=0$ is

$\tilde{\omega}_{i}=e^{\int_{1}^{\tilde{r}}\frac{\tilde{m}_{i}\tilde{\Gamma}_{i}}{\rho_{*}\epsilon\tilde{\eta}_{1i}}\frac{d\tilde{r}^{\prime}}{\tilde{r}^{\prime}}}\int_{1}^{r}d\tilde{r}^{\prime}e^{-\int_{1}^{\tilde{r}^{\prime}}\frac{\tilde{m}_{i}\tilde{\Gamma}_{i}}{\rho_{*}\epsilon\tilde{\eta}_{1i}}\frac{d\tilde{r}^{\prime\prime}}{\tilde{r}^{\prime\prime}}}\cdot\\\
\left[\frac{Z_{i}}{\rho_{*}^{2}\epsilon}\frac{\int_{0}^{\tilde{r}^{\prime}}\tilde{\Gamma}_{i}\tilde{B}\tilde{r}^{\prime\prime}d\tilde{r}^{\prime\prime}}{\tilde{r}^{\prime
3}\tilde{\eta}_{1i}}+\frac{\tilde{\eta}_{1e}\tilde{\omega}_{e}^{\prime}}{\tilde{\eta}_{1i}}\right.\\\
\left.-\frac{\tilde{\Gamma}_{i}\tilde{\eta}_{3i}}{\epsilon\tilde{\eta}_{1i}\tilde{r}^{\prime
2}\tilde{n}_{i}}\left[\frac{2}{\tilde{r}^{\prime}}+\frac{d\ln(\tilde{n}_{i})}{d\tilde{r}}\right]\right].~{}$
(51)

Notice that each term in (51) is proportional to

$\displaystyle\frac{1}{\tilde{\eta}_{1i}}\propto\frac{\sqrt{\tilde{T}_{i}}}{\tilde{n}_{i}^{2}}.$
(52)

If $\tilde{n}_{i}$ were to continuously decrease from $1$ to near $0$ when
moving from $\tilde{r}=1$ inwards or if the ion temperature should diverge,
$\tilde{\omega}_{i}$ would diverge. This is due to the decrease in viscosity
while the magnetic field torque remains the same.

The leading order solution for $\tilde{\omega}_{i}$, taking
$\tilde{\Gamma}_{i}=\frac{1}{2}\tilde{s}_{i}\tilde{r}^{2}$, which is the
solution to the continuity equation with a uniform steady source term
$\tilde{s}_{i}=\mathrm{const.}$, and $\tilde{B}=\mathrm{const.}$, in addition
to $\tilde{\eta}_{1i}=\mathrm{const.}$ is,

$\tilde{\omega}_{i}^{(0,0)}=\frac{Z_{i}\tilde{B}\tilde{s}_{i}}{8\rho_{*}^{2}\epsilon\tilde{\eta}_{1i}}e^{\frac{\tilde{s}_{i}\tilde{m}_{i}}{4\rho_{*}\epsilon\tilde{\eta}_{1i}}(\tilde{r}^{2}-1)}\cdot\int_{1}^{r}e^{-\frac{\tilde{s}_{i}\tilde{m}_{i}}{4\rho_{*}\epsilon\tilde{\eta}_{1i}}(\tilde{r}^{\prime
2}-1)}\tilde{r}^{\prime}\mathrm{d}\tilde{r}^{\prime}~{}\\\
=\frac{\tilde{\Omega}_{i}}{4\rho_{*}}\left(e^{F_{i}(\tilde{r}^{2}-1)}-1\right),$
(53)

It is useful to define the variables
$F_{i}=\frac{\tilde{s}_{i}\tilde{m}_{i}}{4\rho_{*}\epsilon\tilde{\eta}_{1i}}$,
which is the mass flux over viscosity, and
$\tilde{\Omega}_{i}\doteq\frac{Z_{i}\tilde{B}}{\tilde{m}_{i}}$ which is the
dimensionless gyro-frequency. The solution is exponentially dependent on the
strength of the source term $\tilde{s}_{i}$.

The azimuthal velocity in this case is

$\displaystyle\tilde{v}_{\theta
i}^{(0,0)}=\frac{\tilde{\Omega}_{i}}{4\rho_{*}}\tilde{r}\left(e^{F_{i}(\tilde{r}^{2}-1)}-1\right).$
(54)

The angular velocity becomes $\mathcal{O}(1)$ if
$F_{i}\sim\mathcal{O}(\rho_{*})$. An azimuthal Mach 1 ($\tilde{v}_{\theta
i}=1$) is obtained when $F_{i}\approx-6\sqrt{3}\rho_{*}/\tilde{\Omega}_{i}$,
corresponding to a source term
$\tilde{s}_{i}\sim\mathcal{O}(\rho_{*}^{2}\epsilon)$. The next corrections are
of $\mathcal{O}(\tilde{m}_{e}^{3/2})$ due to the electron viscosity and
$\mathcal{O}(\rho_{*}^{2})$ due to the $\tilde{\eta}_{3i}$ rotation term.

###  III.2 Temperatures

The ion temperature equation, when substituting
$\tilde{\Gamma}_{i}=\rho_{*}^{2}\epsilon\tilde{\Gamma}_{i}^{(2,1)}$, and
remembering $\tilde{\omega}_{i}\sim\mathcal{O}(1)$

$-\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\left(\tilde{r}\tilde{\kappa}_{i}\frac{d\tilde{T}_{i}}{d\tilde{r}}\right)+\rho_{*}\frac{5}{2}\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\tilde{T}_{i}\tilde{\Gamma}_{i}^{(2,1)}\\\
=3\frac{1}{\rho_{*}^{2}}\tilde{n}_{i}\tilde{\nu}_{ie}(\tilde{T}_{e}-\tilde{T}_{i})+\tilde{\eta}_{1i}\left(\tilde{r}\frac{d\tilde{\omega}_{i}}{d\tilde{r}}\right)^{2}\\\
+\rho_{*}\frac{d\tilde{p}_{i}}{d\tilde{r}}\frac{\tilde{\Gamma}_{i}^{(2,1)}}{\tilde{r}\tilde{n}_{i}}+\frac{3}{2}\rho_{*}\tilde{s}_{i}^{(2,1)}\tilde{T}_{i}\\\
+\rho_{*}^{4}\frac{\tilde{\eta}_{0i}}{3}\left(\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\frac{\tilde{\Gamma}_{i}^{(2,1)}}{\tilde{n}_{i}}\right)^{2}+\tilde{\eta}_{1i}\left(\rho_{*}^{2}\epsilon\tilde{r}\frac{d}{d\tilde{r}}\left(\frac{\tilde{\Gamma}_{i}^{(2,1)}}{\tilde{r}^{2}\tilde{n}_{i}}\right)\right)^{2},~{}$
(55)

and the electron equation,

$-\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\left(\tilde{r}\tilde{\kappa}_{e_{eff}}\frac{d\tilde{T}_{e}}{d\tilde{r}}\right)=-3\frac{1}{\rho_{*}^{2}}\tilde{n}_{i}\tilde{\nu}_{ie}(\tilde{T}_{e}-\tilde{T}_{i})\\\
+\rho_{*}^{2}\epsilon^{2}\frac{\tilde{m}_{i}\tilde{\nu}_{ie}}{\tilde{r}^{2}\tilde{n}_{i}}\tilde{\Gamma}_{i}^{(2,1)^{2}}+\tilde{\eta}_{1e}\left(\tilde{r}\frac{d\tilde{\omega}_{e}}{d\tilde{r}}\right)^{2}$
(56)

with

$\displaystyle\tilde{\kappa}_{e_{eff}}=\left[4.66+\frac{9}{4}\left[1-\frac{I_{1}\left(\frac{\tilde{r}}{\rho_{*}}\sqrt{\frac{\tilde{m}_{e}\tilde{n}_{e}\tilde{\nu}_{ei}}{\tilde{\eta}_{1e}}}\right)}{I_{1}\left(\frac{1}{\rho_{*}}\sqrt{\frac{\tilde{m}_{e}\tilde{n}_{e}\tilde{\nu}_{ei}}{\tilde{\eta}_{1e}}}\Big{|}_{\tilde{r}=1}\right)}\right]\right]\frac{\tilde{p}_{e}\tilde{\nu}_{ei}}{\tilde{m}_{e}\tilde{\Omega}_{e}^{2}}~{},$
(57)

being the effective heat transfer coefficient for the electrons, due to the
contribution of the Ettingshausen, $\tilde{\kappa}^{\omega}$ term. The plasma
heat transfer coefficient, for $Z_{i}=1$, $\tilde{m}_{i}=1$, is
$\tilde{\kappa}_{tot}=\tilde{\kappa}_{i}+\tilde{\kappa}_{e_{eff}}\approx
1.16\tilde{\kappa}_{i}$, outside of the boundary layer.

This is a linear ODE for $\tilde{T}_{i}$, with a single term - the collisional
equilibration between ions and electrons $\sim\rho_{*}^{-2}$. Thus, in the
low-flow case, $F_{i}\sim\mathcal{O}(\rho_{*})$,
$\Gamma_{i}\sim\mathcal{O}(\rho_{*}^{2}\epsilon)$, the temperature difference
between the ion and electron fluids is of
$\mathcal{O}(\rho_{*}^{2}\tilde{\omega}_{i}^{2}/\sqrt{\tilde{m}_{e}})\sim\mathcal{O}(\rho_{*}^{2}/\sqrt{\tilde{m}_{e}})$.

To leading order, with constant coefficients, the sum of the two equations is,

$\displaystyle-\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\left(\tilde{r}\tilde{\kappa}_{tot}\frac{d\tilde{T}_{i}^{(0,0)}}{d\tilde{r}}\right)=\tilde{\eta}_{1i}\left(\tilde{r}\frac{d\tilde{\omega}_{i}}{d\tilde{r}}\right)^{2},$
(58)

yielding

$\tilde{T}_{i}^{(0,0)}=-\frac{\tilde{\eta}_{1i}}{\tilde{\kappa}_{tot}}\frac{\tilde{\Omega}_{i}^{2}e^{-2F_{i}}}{32\rho_{*}^{2}F_{i}}\Bigg{(}\frac{1}{2}(Ei(2F_{i}\tilde{r}^{2})-Ei(2F_{i}))-\ln(\tilde{r})\\\
+\frac{1}{4}(e^{2F_{i}\tilde{r}^{2}}(2F_{i}\tilde{r}^{2}-3)-e^{2F_{i}}(2F_{i}-3))\Bigg{)}+\tilde{T}_{i}(1)$
(59)

with $\tilde{\eta}_{1i}/\tilde{\kappa}_{i}=0.15\tilde{m}_{i}$ and
$\tilde{\eta}_{1i}/\tilde{\kappa}_{tot}\approx 0.13\tilde{m}_{i}$. These
ratios are related to the Prandtl number as
$\tilde{\eta}/\tilde{\kappa}=\frac{2}{3}\tilde{m}\mathrm{Pr}$. The function
$Ei(x)=\int_{-\infty}^{x}\frac{e^{t}}{t}dt$ is the exponential integral. For
small values of $F_{i}$, such as $F_{i}\approx-\rho_{*}$, the boundary term
$\tilde{T}_{i}(1)=1$ is the dominant term in the leading-order temperature
solution.

The temperature difference between electrons and ions $\Delta
T_{ei}=\tilde{T}_{e}-\tilde{T}_{i}$, is determined, to leading order, by

$-\frac{1}{\tilde{r}}\frac{d}{d\tilde{r}}\left(\tilde{r}\frac{d\Delta
T_{ei}}{d\tilde{r}}\right)=-\frac{3\tilde{n}_{i}\tilde{\nu}_{ie}}{\rho_{*}^{2}}\left(\frac{1}{\tilde{\kappa}_{e_{eff}}}+\frac{1}{\tilde{\kappa}_{i}}\right)\Delta
T_{ei}\\\
+\frac{2}{3}\left(\tilde{m}_{e}\mathrm{Pr}_{e}-\tilde{m}_{i}\mathrm{Pr}_{i}\right)\left(\tilde{r}\frac{d\tilde{\omega}_{i}}{d\tilde{r}}\right)^{2}.~{}$
(60)

The leading order solution to this equation is

$\Delta
T_{ei}=\left.-\frac{2\rho_{*}^{2}}{9\tilde{n}_{i}\tilde{\nu}_{ie}}\frac{\tilde{m}_{i}\mathrm{Pr}_{i}-\tilde{m}_{e}\mathrm{Pr}_{e}}{\frac{1}{\tilde{\kappa}_{i}}+\frac{1}{\tilde{\kappa}_{e_{eff}}}}\right(\tilde{r}^{2}\left[\frac{d\tilde{\omega}_{i}}{d\tilde{r}}\right]^{2}\\\
\left.-\left[\frac{d\tilde{\omega}_{i}}{d\tilde{r}}\right]^{2}\Bigg{|}_{\tilde{r}=1}\frac{I_{0}\left(\frac{\tilde{r}}{\rho_{*}}\sqrt{3\tilde{n}_{i}\tilde{\nu}_{ie}\left(\tilde{\kappa}_{e_{eff}}^{-1}+\tilde{\kappa}_{i}^{-1}\right)}\right)}{I_{0}\left(\frac{1}{\rho_{*}}\sqrt{3\tilde{n}_{i}\tilde{\nu}_{ie}\left(\tilde{\kappa}_{e_{eff}}^{-1}+\tilde{\kappa}_{i}^{-1}\right)\Big{|}_{\tilde{r}=1}}\right)}\right)~{}$
(61)

with $I_{0}$ being the modified Bessel function of the first kind.

Since $\tilde{m}_{e}\mathrm{Pr}_{e}<\tilde{m}_{i}\mathrm{Pr}_{i}$, the
solution is negative, indeed yielding a hot-ion mode, but it also contains a
pre-factor of $\rho_{*}^{2}/\sqrt{\tilde{m}_{e}}$. In order to achieve a
temperature difference of $\mathcal{O}(1)$, we identify a course of action,
that does not suffer from the fueling limit discussed later. In any case, one
has to push $\tilde{\omega}$ to be as large as possible. However, the
$\rho_{*}^{2}/\sqrt{\tilde{m}_{e}}$ pre-factor has to be dealt with as well:
Use low magnetization, such that $\rho_{*}\gtrsim\sqrt[4]{\tilde{m}_{e}}\sim
0.1$. This would bring the pre-factor to be $\mathcal{O}(1)$. In this case,
$\tilde{s}_{i}^{(2,1)}$ can remain small such that the angular velocity is
still linear in $\tilde{s}_{i}$. The source term magnitude must then be larger
than the square root of thecombined value of the reminder of (61), which is of
$\mathcal{O}(0.001)$, when substituting $\tilde{\omega}_{i}$ as a function of
$\tilde{s}_{i}$. In case of Braginskii coefficients for
$Z_{i}=\tilde{m}_{i}=1$, this source term is
$\tilde{s}_{i}\approx-20\rho_{*}^{2}\epsilon$. It is perhaps easier to
increase $\rho_{*}$ rather than $\tilde{s}_{i}$, as the maximal temperature
difference increases as $\rho_{*}^{3}\left(\tilde{s}_{i}^{(2,1)}\right)^{2}$
at large $\rho_{*}$, see Figure 3. In the figure, the temperature difference
is $\sim 0.01$ because the source term magnitude kept relatively small for
numerical reasons. This Larmor radius is becoming somewhat large for a small
parameter, required for the fluid approximation. Some magnetic traps operate
with this Larmor radius / machine size scale.

###  III.3 Density

The density profile is determined from the radial force balance. Summing
electrons and ions, taking $\tilde{n}_{e}=Z_{i}\tilde{n}_{i}$, and neglecting
electron viscosity, the radial force balance of the plasma as a whole is

$\frac{d\tilde{p}_{i}}{d\tilde{r}}+\frac{d\tilde{p}_{e}}{d\tilde{r}}=\frac{3}{2}Z_{i}\tilde{n}_{i}\frac{d\tilde{T}_{e}}{d\tilde{r}}+\tilde{n}_{i}\tilde{r}\left(\tilde{m}_{i}\tilde{\omega}_{i}^{2}+Z_{i}\tilde{m}_{e}\tilde{\omega}_{e}^{2}\right)\\\
+\frac{\tilde{m}_{i}\tilde{\Gamma}_{i}^{2}}{\tilde{r}^{2}\tilde{n}_{i}}\left(\frac{1}{\tilde{r}}+\frac{d\ln(\tilde{n}_{i})}{d\tilde{r}}-\frac{d\ln(\tilde{\Gamma}_{i}^{2})}{d\tilde{r}}\right)\\\
+\frac{\rho_{*}}{\epsilon}\frac{d}{d\tilde{r}}\left(\frac{\tilde{\eta}_{0i}\tilde{\Gamma}_{i}}{3\tilde{r}\tilde{n}_{i}}\left[\frac{d\ln(\tilde{\Gamma}_{i})}{d\tilde{r}}-\frac{d\ln(\tilde{n}_{i})}{d\tilde{r}}\right]\right)\\\
+\frac{\rho_{*}}{\tilde{r}^{2}}\frac{d}{d\tilde{r}}\tilde{r}^{3}\left(\tilde{\eta}_{3i}\frac{d\tilde{\omega}_{i}}{d\tilde{r}}+\tilde{m}_{e}\tilde{\eta}_{3e}\frac{d\tilde{\omega}_{e}}{d\tilde{r}}\right.\\\
\left.+\epsilon\frac{\tilde{\eta}_{1i}\tilde{\Gamma}_{i}}{\tilde{r}^{2}\tilde{n}_{i}}\left[\frac{d\ln(\tilde{\Gamma}_{i})}{d\tilde{r}}-\frac{d\ln(\tilde{n}_{i})}{d\tilde{r}}-\frac{2}{\tilde{r}}\right]\right).~{}$
(62)

This is a second-order nonlinear ODE in $\tilde{n}_{i}$, and solving the
leading-order terms in it would both remove the non-linearity, and reduce it
to a first-order differential equation. To leading order,

$\displaystyle\frac{d\ln(\tilde{n}_{i})}{d\tilde{r}}+\frac{1-\frac{1}{2}Z_{i}}{1+Z_{i}}\frac{d\ln(\tilde{T}_{i})}{d\tilde{r}}=\frac{\tilde{r}\left(\tilde{m}_{i}\tilde{\omega}_{i}^{2}+Z_{i}\tilde{m}_{e}\tilde{\omega}_{e}^{2}\right)}{\tilde{T}_{i}(1+Z_{i})},$
(63)

or for $Z_{i}=1$,

$\displaystyle\tilde{n}_{i}^{(0,0)}=\frac{1}{\sqrt[4]{\tilde{T}_{i}}}e^{\frac{\tilde{m}_{i}\tilde{\Omega}_{i}^{2}}{32\rho_{*}^{2}}\int_{1}^{\tilde{r}}\frac{\tilde{r}^{\prime}\left(e^{F_{i}(\tilde{r}^{\prime
2}-1)}-1\right)^{2}}{\tilde{T}_{i}}d\tilde{r}^{\prime}}.$ (64)

This density profile is hollow - matter is pushed to outer radii from the
center.

Figure 1: Comparison between solutions to the constant-coefficient leading
order Braginskii equations in steady state (full line), and MITNS nonlinear
results (dashed line with markers). The source term used for the ODE set and
MITNS were the same. The nonlinear solution starts diverging from the leading
order linear solution at Mach number of $M\approx 0.5$. At a Mach number of
$M\approx 1$ in the linear solution, the difference between the linear and
nonlinear solutions is quite significant.

###  III.4 Validity of the constant coefficients solution

The viscosity profile, to leading order, is

$\displaystyle\tilde{\eta}_{1i}\propto\frac{\tilde{n}_{i}^{2}}{\sqrt{\tilde{T}_{i}}}=\frac{1}{\tilde{T}_{i}}e^{\frac{\tilde{m}_{i}\tilde{\Omega}_{i}^{2}}{16\rho_{*}^{2}}\int_{1}^{\tilde{r}}\frac{\tilde{r}^{\prime}\left(e^{F_{i}(\tilde{r}^{\prime
2}-1)}-1\right)^{2}}{\tilde{T}_{i}}\mathrm{d}\tilde{r}^{\prime}}.$ (65)

Some solutions to the constant coefficient equations are presented in Figure
1. Notice how $\tilde{\eta}_{1i}\propto\tilde{\kappa}_{i}$ drops rapidly as
the source term is increased. The disagreements between the linear and
nonlinear solutions are visible at Mach number of $M\approx 0.5$, and become
increasingly severe.

##  IV Variable coefficients solution

Braginskii’s fluid model is inherently nonlinear, with collision frequencies
that depend on the densities and temperatures. The diffusion coefficients
(viscosities and heat transfer coefficients), present additional nonlinearity.
From the linear, constant coefficients solution, it is evident that the
viscosity and heat transfer coefficients drop as the magnitude of the particle
flux increases. The shorthand $F_{i}$ would become larger at smaller radii,
and the angular velocity increase beyond its constant coefficients solution
values.

There are two ways in which a particle flux would fail to produce a physical
solution:

  1. 1.

The solution breaks the ordering $\rho_{i}\ll R$.

  2. 2.

The solution produces a negative pressure, or is unable to satisfy both
boundary conditions.

For the first case, even a constant coefficient solution with
$F_{i}\sim\mathcal{O}({\rho_{*}})$ would produce rotations that are
$\tilde{\omega}\sim\mathcal{O}(\rho_{*}^{-1})$, and temperature that is
$\mathcal{O}(\rho_{*}^{-2})$. The Larmor radius $\rho_{i}\propto\sqrt{T_{i}}$,
and increasing the temperature above its reference value by a factor of
$\rho_{*}^{-2}$ puts us firmly in the kinetic regime, where the fluid model is
inapplicable. Using a variable-coefficients solution would bring that
threshold to smaller values of particle flux.

The second case is a feature of non-linear ODEs, where there is no guarantee
for the existence of a boundary-value-problem solution in all cases. We can
attempt solving initial value problems, and using a shooting method to
pinpoint the correct values and derivatives at one boundary to hit the correct
values at the other, or transform the ODEs into PDEs, and attempt to relax the
solution to a steady state at finite times.

Even at Mach numbers in which a solution for the nonlinear case exists, the
solution might develop an angular velocity boundary layer at $\tilde{r}=0$
which might be nonphysical. This boundary layer enforces the axisymmetry
condition $\tilde{\omega}_{i}^{\prime}(0)=0$, which appears in (30). The
boundary layer width, shrinks rapidly around $M=2$, while the temperature
rises, as can be seen in Figure 2. When the boundary layer width becomes
smaller than an ion Larmor radius, the solution becomes nonphysical. See
appendix. Steady state solutions for $M\gtrsim 2$ don’t exist in the full
nonlinear case, as the viscosity drops quite significantly, and it could not
balance the torque produced by the magnetic field.

Figure 2 presents the values of the density, temperature, and angular velocity
at the center of the cylinder, and the Mach number as a function of the
particle sink magnitude.

Figure 2: Values of angular velocity at the center of the cylinder and Mach
number (maximal azimuthal velocity) (top), and values of density and
temperature at the center of the cylinder, as a function of $|\tilde{s}_{i}|$
(bottom). These nonlinear calculations were performed using MITNS, with
$\rho_{*}=0.01$ and $\epsilon=0.1$.

Figure 3 shows the temperature difference between ions and electrons at the
center of the cylinder, and the maximal value of the temperature difference
between the species, as a function of $\rho_{*}$. The sink term was chosen
such that $F_{i}/\rho_{*}$ would remain constant at the boundary, and the
azimuthal Mach number would approximately be constant at $M\approx 0.93$, with
small variation due to the nonlinear nature of the system. This logarithmic
plot shows a near-perfect power law, with the two curves crossing at
$\rho_{*}=0.11$, when the center of the cylinder becomes the point of maximal
temperature difference between species. The maximal temperature difference
depends very nearly on $\rho_{*}^{2}$, even in the nonlinear case.

Figure 3: Temperature difference between ions and electrons at the center of
the cylinder, and the maximal temperature difference as a function of
$\rho_{*}$ (top), and the power-law dependence of the same curves (bottom).
Results from the nonlinear solution in MITNS. The ion sink term was
$\tilde{s}_{i}=-5\rho_{*}^{2}\epsilon$, such that $F_{i}\sim\rho_{*}$, with
$\epsilon=0.1$.

##  V Conclusion

A solution to the flow and temperature profiles in a highly magnetized
rotating cylindrical two-fluid (ion-electron) plasma, driven by constant ion
charge extraction following Braginskii’s fluid was investigated. First, a
leading-order solution in the low-flow limit was presented. Second, the
physical validity of the asymptotic solution was considered, and the existence
of a hot-ion mode was evaluated. It was shown that the collisional temperature
coupling between the fluids is stronger than the difference in viscous heating
between the fluids by a factor of $\sqrt{\tilde{m}_{e}}/\rho_{*}^{2}$, and
this limits the temperature difference to be of
$\mathcal{O}(\rho_{*}^{2}/\sqrt{\tilde{m}_{e}})$.

The ratio of ion to electron heating is the same as the ratio of viscosity
coefficients $\tilde{m}_{e}^{3/2}$. Picking a
$\rho_{*}\sim\sqrt[4]{\tilde{m}_{e}}\sim 0.1$ would bring the temperature
difference to $\mathcal{O}(0.1-1)$, depending on the exact sink term
magnitude. This would be a large $\rho_{*}$, but not impossibly so - some
magnetic traps or FRCs[34] operate in this regime. Finally, the departure of
the non-linear solution from the linear approximation at moderate Mach numbers
($M\gtrsim 0.5$) was demonstrated and is explained by the hollowing-out of the
viscosity and heat conductivity profiles.

We have shown the inherent difficulties in achieving a significant hot-ion
mode, even in the absence of electron heating, due to the limitation on ion
heating. We have shown that low magnetization devices present an easier avenue
for a hot-ion mode. Paradoxically, this is a point in favor of devices that
contain fewer ion Larmor radii - which can be accomplished using a smaller
magnetic field or a smaller device size - in view of the large advantage a
hot-ion mode might present for a fusion plasma.

There are other possible ion sink profiles, $\tilde{s}_{i}(\tilde{r})$, and
the details of such solutions may differ from the solution presented here. The
main effect discussed here, however, the hollowing out of the density and the
viscosity profiles, is weakly dependant on the specific
$\tilde{s}_{i}(\tilde{r})$ profile, and requires only the functional
dependence $\tilde{\eta}_{1}\propto\tilde{n}^{2}$.

### Acknowledgments

This work was supported by Cornell NNSA 83228-10966 [Prime No. DOE (NNSA) DE-
NA0003764] and by NSF-PHY-1805316. IEO also acknowledges support from the DOE
FES Fellowship [DOE Contract DE-SC0014664].

## Author Declarations

### Conflict of interest

The authors have no conflicts to disclose.

### Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Boundary layers

The above solution includes several boundary layers, in the electron angular
velocity, effective heat transfer coefficient, and the temperature difference
between electrons and ions. The arguments of the boundary layers are expressed
in terms of the electron viscosity and heat transfer coefficients. Using
Braginskii’s classical transport coefficients, the argument of the electron
angular velocity boundary layer term is

$\displaystyle\frac{\tilde{r}}{\rho_{*}}\sqrt{\frac{\tilde{m}_{e}\tilde{n}_{e}\tilde{\nu}_{ei}}{\tilde{\eta}_{1e}}}=\frac{\tilde{r}}{\rho_{*}}\frac{\tilde{B}}{\sqrt{0.51\tilde{T}_{e}\tilde{m}_{e}}}=1.4\frac{r\Omega_{e}}{v_{the}}=1.4\frac{r}{\rho_{e}}.$
(66)

Here, $\rho_{e}$ is the electron Larmor radius. This boundary layer is thinner
than a Larmor orbit. As such, the classical transport model breaks down at
these length scales, hence, we must abandon this boundary layer solution and
the boundary condition that produced it, and resort to a full slip condition
for the electron fluid at the outer edge of the cylinder. This same boundary
layer appears in the effective electron heat conduction, and must be discarded
there as well.

As such, in the limit of Braginskii’s classical transport, equations (27),
(50), and (57) should simply read

$\displaystyle(\tilde{\nabla}\cdot\tilde{\pi}_{e})_{\theta}(1)=0,$ (67)
$\displaystyle\tilde{\omega}_{e}=\tilde{\omega}_{i}-\frac{3}{2}\frac{\rho_{*}}{\tilde{r}\tilde{B}}\frac{d\tilde{T}_{e}}{d\tilde{r}},$
(68)
$\displaystyle\tilde{\kappa}_{e_{eff}}=6.91\frac{\tilde{p}_{e}\tilde{\nu}_{ei}}{\tilde{m}_{e}\tilde{\Omega}_{e}^{2}}.$
(69)

However, some non-classical transport effects [35], such as turbulent
anomalous transport, or the effects of perturbations to the magnetic field
[36], might enhance the electron viscosity above Braginskii’s values, while
still maintaining the validity of the fluid approximation. In these cases, if
the effective electron viscosity might be large enough such that the boundary
layer thickness might encompass several Larmor orbits, the boundary layer
solution in the electron angular velocity might be a true physical effect.

The argument of the temperature-difference boundary layer, using Braginskii’s
classical transport coefficients, is

$\displaystyle\frac{\tilde{r}}{\rho_{*}}\sqrt{3\tilde{n}_{i}\tilde{\nu}_{ie}\left(\tilde{\kappa}_{e_{eff}}^{-1}+\tilde{\kappa}_{i}^{-1}\right)}=\frac{\tilde{r}}{\rho_{*}}\frac{0.71}{\sqrt{\tilde{m}_{i}}}.$
(70)

This boundary layer does encompass many electron gyro-orbits, and is plausible
because of it. Physically, it must be enforced by some effect on the boundary.
We are not concerned with modeling the plasma-surface interactions that might
lead to this electron temperature boundary layer, and conclude that the fluid
approximation is valid for a boundary layer, with, possibly a different
boundary value for the electrons other than the ions.

Note that the magnitude of $\Delta T_{ei}$, for the constant coefficient
solution, is controlled by the scaling of the flux ($\tilde{s}_{i}$ or
equivalently $F_{i}$), which directly affects $\tilde{\omega}_{i}^{\prime}$,
and is not dependent on the exact form of the viscosity or heat transfer
coefficients. Trying to maximize $\tilde{T}_{i}-\tilde{T}_{e}$ in equation
(61), by changing $\tilde{\kappa}_{e_{eff}}$ only, or in conjunction with
$\tilde{\eta}_{1e}$, using an anomalous electron transport, for example, would
yield limited results.

In a fusion reactor, when $\alpha$-channeling is missing or is insufficiently
effective at removing the fusion ash, the ash would slow down on the electrons
and introduce significant electron heating. Indeed this mechanism typically
causes a hot electron mode. Alternatively, Bremsstrahlung radiation introduces
a heat sink in the electrons, which would help maintain a hot-ion mode, when
dealing with mildly relativistic plasmas such as in $p-B^{11}$ reactors.

## References

## References

  * Lehnert [1971] B. Lehnert, Nucl. Fusion 11, 485 (1971).
  * Bekhtenev  _et al._ [1980] A. Bekhtenev, V. Volosov, V. Pal'chikov, M. Pekker, and Y. Yudin, Nuclear Fusion 20, 579 (1980).
  * Hassam [1999] A. B. Hassam, Physics of Plasmas 6, 3738 (1999).
  * Rax  _et al._ [2017] J.-M. Rax, R. Gueroult, and N. J. Fisch, Phys. Plasmas 24, 032504 (2017).
  * Ochs and Fisch [2017] I. E. Ochs and N. J. Fisch, Phys. Plasmas 24, 092513 (2017).
  * Ellis  _et al._ [2001] R. F. Ellis, A. B. Hassam, S. Messer, and B. R. Osborn, Physics of Plasmas 8, 2057 (2001).
  * Ellis  _et al._ [2005] R. F. Ellis, A. Case, R. Elton, J. Ghosh, H. Griem, A. Hassam, R. Lunsford, S. Messer, and C. Teodorescu, Phys. Plasmas 12, 055704 (2005).
  * Teodorescu  _et al._ [2010] C. Teodorescu, W. C. Young, G. W. S. Swan, R. F. Ellis, A. B. Hassam, and C. A. Romero-Talamas, Phys. Rev. Lett. 105, 085003 (2010).
  * Rax and Gueroult [2016] J.-M. Rax and R. Gueroult, J. Plasma Phys. 82, 595820504 (2016).
  * Zweben  _et al._ [2018] S. J. Zweben, R. Gueroult, and N. J. Fisch, Phys. Plasmas 25, 090901 (2018).
  * Ochs  _et al._ [2017] I. E. Ochs, R. Gueroult, N. J. Fisch, and S. J. Zweben, Phys. Plasmas 24, 043503 (2017).
  * Gueroult and Fisch [2012] R. Gueroult and N. J. Fisch, Phys. Plasmas 19, 122503 (2012).
  * Gueroult  _et al._ [2015] R. Gueroult, D. T. Hobbs, and N. J. Fisch, J. Hazard. Mater. 297, 153 (2015).
  * Gueroult  _et al._ [2019] R. Gueroult, J.-M. Rax, and N. J. Fisch, Phys. Plasmas 26, 122106 (2019).
  * Fetterman and Fisch [2011] A. J. Fetterman and N. J. Fisch, Phys. Plasmas 18, 094503 (2011).
  * Bonnevier [1966] B. Bonnevier, Ark. Fys. 33, 255 (1966).
  * Bonnevier [1971] B. Bonnevier, Plasma Phys. 13, 763 (1971).
  * Kolmes  _et al._ [2019] E. J. Kolmes, I. E. Ochs, M. E. Mlodik, J.-M. Rax, R. Gueroult, and N. J. Fisch, Physics of Plasmas 26, 082309 (2019).
  * Rax  _et al._ [2019] J.-M. Rax, E. J. Kolmes, I. E. Ochs, N. J. Fisch, and R. Gueroult, Phys. Plasmas 26, 012303 (2019).
  * Rubin  _et al._ [2021] T. Rubin, E. J. Kolmes, I. E. Ochs, M. E. Mlodik, and N. J. Fisch, Physics of Plasmas 28, 122303 (2021).
  * Kolmes  _et al._ [2021a] E. J. Kolmes, I. E. Ochs, and N. J. Fisch, Comp. Phys. Communications 258, 107511 (2021a).
  * Lawson [1957] J. D. Lawson, Proceedings of the Physical Society. Section B 70, 6 (1957).
  * Clarke [1980] J. Clarke, Nuclear Fusion 20, 563 (1980).
  * Kolmes  _et al._ [2021b] E. J. Kolmes, I. E. Ochs, M. E. Mlodik, and N. J. Fisch, Phys. Rev. E 104, 015209 (2021b).
  * Fisch and Rax [1992] N. J. Fisch and J.-M. Rax, Phys. Rev. Lett. 69, 612 (1992).
  * Fisch and Herrmann [1994] N. J. Fisch and M. C. Herrmann, Nucl. Fusion 34, 1541 (1994).
  * Braginskii [1965] S. I. Braginskii, Transport processes in a plasma, in  _Reviews of Plasma Physics_ , Vol. 1, edited by M. A. Leontovich (Consultants Bureau, New York, 1965) p. 205.
  * Epperlein and Haines [1986] E. M. Epperlein and M. G. Haines, Phys. Fluids 29, 1029 (1986).
  * Simakov [2022] A. N. Simakov, Physics of Plasmas 29, 022304 (2022).
  * Ji and Held [2013] J.-Y. Ji and E. D. Held, Physics of Plasmas 20, 042114 (2013).
  * Sadler  _et al._ [2021] J. D. Sadler, C. A. Walsh, and H. Li, Phys. Rev. Lett. 126, 075001 (2021).
  * Davies  _et al._ [2021] J. R. Davies, H. Wen, J.-Y. Ji, and E. D. Held, Physics of Plasmas 28, 012305 (2021).
  * Fundamenski and Garcia [2007] W. Fundamenski and O. Garcia, Report No. EFDA-JET-R (07) 01 (2007).
  * Steinhauer [2011] L. C. Steinhauer, Physics of Plasmas 18, 070501 (2011).
  * Rognlien and Ryutov [1999] T. Rognlien and D. Ryutov, Plasma Physics Reports 25, 943 (1999).
  * Finn  _et al._ [1992] J. M. Finn, P. N. Guzdar, and A. A. Chernikov, Physics of Fluids B: Plasma Physics 4, 1152 (1992).