# Orthogonality graphs of matrices
over commutative rings

O. G. Styrt

The paper is devoted to studying the orthogonality graph of the matrix ring
over a commutative ring. It is proved that the orthogonality graph of the ring
of matrices with size greater than $1$ over a commutative ring with zero-
divisors is connected and has diameter $3$ or $4$; a criterion for each value
is obtained. It is also shown that each of its vertices has distance at most
$2$ from some scalar matrix.

Key words: associative ring with identity, commutative ring, zero-divisor,
matrix ring, zero-divisor graph, orthogonality graph.

##  1 Introduction

Researching properties of associative rings in terms of graphs of some
naturally occurring algebraic binary relations takes an important place in
modern mathematics. Thus, a zero-divisor graph was first defined in 1986 by
Beck [1] for a commutative ring. Its vertices were all zero-divisors, and
edges connected exactly all pairs of distinct elements giving zero in product.
But since 1999 one uses its more convenient interpretation introduced by
Anderson and Livingston in [2] via excluding the zero element of the ring from
its vertex set. It is also proved in [2] that the zero-divisor graph of a
commutative ring is connected and has diameter at most three; in the former
treatment of the graph these statements would be trivial. A number of further
papers also studies various characteristics of the zero-divisor graph: center
and radius [8], concepts of planarity [4] and uniqueness of determining the
ring by the graph up to an isomorphism [3, 5]. For non-commutative rings,
there are several types of graphs defined by zero-divisors:

Table 1: *

|  |  |  |  |
---|---|---|---|---|---

|
Name
|
Edge orientation
|
Vertices
|
Edge
from $x$ to $y$
|
See

1)
|
Directed zero-divisor graph
|
Yes
|
One- and two-sided
zero-divisors
|
$xy=0$
|
[6, 7]

2)
|
(Undirected) zero-divisor graph
|
No
|
Nonzero
one- and two-sided
zero-divisors
|
$\left[\begin{array}[]{@{}l@{\quad}l@{}}xy=0\\\ yx=0\end{array}\right.$
|
[7]

3)
|
Orthogonality graph
|
No
|
Nonzero two-sided
zero-divisors
|
$\begin{cases}xy=0\\\ yx=0\end{cases}$
|
[9, 10]

The main results for orthogonality graphs of non-commutative rings found by
now concern primarily matrix rings. Thus, in the case of the basic ring being
a skew field, the following properties of the orthogonality graph of the
$(n\times n)$-matrix ring are obtained: once $n=2$, it is disconnected and all
its connected components have diameters at most $2$, and, once $n\geqslant 3$,
it is connected and has diameter $4$. These statements are proved in 2014 for
a field [9] and later, in 2017 — for an arbitrary skew field [10]; they can
also be easily generalized to integral domains (by reducing to the field of
fractions).

In this paper, there will be the orthogonality graph of the matrix ring over a
commutative ring with zero-divisors studied and the following main result
proved.

######  Theorem 1.1.

Let $R$ be a commutative ring with zero-divisor set $Z_{R}\neq\\{0\\}$. Then,
for any $n>1$, the orthogonality graph of the ring of $(n\times n)$-matrices
over $R$ is connected and has diameter $3$ or $4$, the value $3$ being
equivalent to the relation

$\,\forall\,a_{0}\in Z_{R}\quad\,\exists\,a_{1},a_{2}\in
R\setminus\\{0\\}\quad\quad\,\forall\,i,j\in\\{0,1,2\\},\ i\neq j\colon\quad
a_{i}a_{j}=0,$ (1.1)

and each of its vertices has distance at most $2$ from some scalar matrix.

######  Theorem 1.2.

Let $r$ be the radius of the graph under conditions of Theorem 1.1. Then

  1. 1)

$2\leqslant r\leqslant 4$;

  2. 2)

if (1.1) holds, then $r\in\\{2;3\\}$;

  3. 3)

$r=2$ if and only if there exists an element $c\in R\setminus\\{0\\}$ such
that

$\,\forall\,a\in Z_{R}\quad\operatorname{Ann}(c)\cap\operatorname{Ann}(a)\neq
0.$ (1.2)

##  2 Auxiliary agreements

In the paper, the following notations and agreements will be used.

  1. 1)

Set-theoretical:

     * •

While listing elements of a disordered set, figured brackets are used. As for
elements of an ordered tuple, they are listed in round brackets and can be
repeated.

     * •

$D^{n}:=\underbrace{D\times\nobreak\ldots\times\nobreak D}_{n}$ is the $n$-ary
Cartesian power of a set $D$.

  2. 2)

General algebraic:

     * •

All rings considered are supposed to be associative and with identity.

     * •

$R$ is an arbitrary ring.

     * •

For any subset $D\subset R$, define $D^{*}:=D\setminus\\{0\\}$. In particular,
by $R^{*}$ denote the subset of all nonzero (not necessarily invertible as in
standard interpretation) elements of $R$.

     * •

An ideal in $R$ is proper if it does not equal $R$.

     * •

$M_{m\times n}(R)$ is the $R$-module of $(m\times n)$-matrices over $R$;
$M_{n}(R)$ is the ring $M_{n\times n}(R)$. If in the brackets the ring is
replaced with some of its subsets $D$, then the subset of all matrices with
entries from $D$ is meant.

     * •

$0^{m}_{n}$ is the zero $(m\times n)$-matrix; $0_{n}:=0^{n}_{n}$; $E_{n}$ is
the identity $(n\times n)$-matrix; $J_{r}$ is the Jordan cell of size $r$ with
eigenvalue $0$. If the matrix sizes are clear from the context, then the
indices can be omitted.

     * •

$E_{kl}$ is the matrix unit $(a_{ij})$, $a_{ij}:=\delta_{ki}\delta_{lj}$.

     * •

For a square matrix $A$ over a commutative ring: $\widetilde{A}$ is its
cofactor matrix; $\widehat{A}:=\bigl{(}\widetilde{A}\bigr{)}^{T}$.

     * •

If $A=(a_{k_{1},k_{2}})\in M_{n_{1}\times n_{2}}(R)$,
$P_{i}\in\\{1,\ldots,n_{i}\\}^{m_{i}}$ ($i=1,2$), then $A^{P_{1}}_{P_{2}}$ is
the matrix $(b_{l_{1},l_{2}})\in M_{m_{1}\times m_{2}}(R)$,
$b_{l_{1},l_{2}}:=a_{k_{1}(l_{1}),k_{2}(l_{2})}$, where $k_{i}(l_{i})$ is the
$l_{i}$-th element of $P_{i}$. If numbers are repeated neither in $P_{1}$, nor
in $P_{2}$, then $A^{P_{1}}_{P_{2}}$ is the submatrix of $A$ with row and
column numbers from $P_{1}$ and $P_{2}$ respectively.

  3. 3)

On zero-divisor types:

     * •

An element $a\in R$ is called

       * –

a left (resp. right) zero-divisor if there exists an element $b\in R^{*}$ such
that $ab=0$ (resp. $ba=0$);

       * –

a zero-divisor if it is either left or right zero-divisor;

       * –

a two-sided zero-divisor if it is both left and right zero-divisor.

At that,

       * –

in a commutative ring, the concepts of all zero-divisor types are equivalent;

       * –

zero is a two-sided zero-divisor; if there are no other zero-divisors, then
$R$ is called a ring without zero-divisors.

     * •

An integral domain is a commutative ring without zero-divisors.

  4. 4)

From general graph theory:

     * •

All graphs considered are assumed to be undirected.

     * •

$\Gamma=(V,E)$ is an arbitrary graph; $V$ and $E$ are its vertex and edge sets
respectively. In doing so, one can (usually with more convenience) define $E$
via a symmetric binary relation on $V$.

     * •

Two vertices are adjacent if they are connected with an edge.

     * •

A subgraph is a graph with vertex set $V^{\prime}\subset V$ and, unless
otherwise stated, with the same binary relation restricted on $V^{\prime}$.

     * •

A path is a sequence of vertices where any two neighbor ones are adjacent.

     * •

The length of a path is the number of its edges.

     * •

The distance between vertices $v$ and $w$ (not. $d(v,w)$) is the minimum of
lengths of paths between them; if they do not exist, then set
$d(v,w):=+\infty$; the sign is obvious in this context and therefore will be
omitted. Clearly, $\bigl{(}d(v,w)=0\bigr{)}\Leftrightarrow(v=w)$.

     * •

The distance from a vertex $v$ to a subset $W\subset V$ (not. $d(v,W)$) is the
number***Possibly $\infty$.

$\min\bigl{\\{}d(v,w)\colon w\in W\bigr{\\}}.$

     * •

$d(v):=\sup\bigl{\\{}d(v,w)\colon w\in W\bigr{\\}}$ ($v\in V$).

     * •

The diameter of $\Gamma$ is the number**footnotemark: *

$\operatorname{diam}(\Gamma):=\sup\bigl{\\{}d(v,w)\colon v,w\in
V\bigr{\\}}=\max\bigl{\\{}d(v)\colon v\in W\bigr{\\}}.$

     * •
     * •

The radius of $\Gamma$ is the number**footnotemark: *
$\operatorname{rad}(\Gamma):=\min\bigl{\\{}d(v)\colon v\in W\bigr{\\}}$.
Clearly,

$\operatorname{rad}(\Gamma)\leqslant\operatorname{diam}(\Gamma)\leqslant
2\cdot\operatorname{rad}(\Gamma).$ (2.1)

     * •

A graph is connected if there exists a path between any two of its vertices.

Remark. It is easy to see that a graph with finite diameter is connected. The
converse fails; an example is the set of positive integers with the
neighborhood relation.

  5. 5)

On special graphs in algebraic structures:

     * •

$O(R)$ is the orthogonality graph of the ring $R$ (for a commutative ring it
is the same as the zero-divisor graph).

     * •

Vertices of $O(R)$ are all nonzero two-sided zero-divisors of $R$; the
orthogonality relation ($xy=yx=0$) is written as ($x\perp y$); $O_{R}(x)$ is
the set of all vertices orthogonal to $x$.

##  3 Proofs of the results

Consider an arbitrary commutative ring $R$. Denote by $\operatorname{Ann}(a)$
($a\in R$) the ideal $\\{x\in\nobreak R\colon ax=0\\}$ and by $Z_{R}$ the set
$\bigl{\\{}a\in R\colon\operatorname{Ann}(a)\neq 0\bigr{\\}}$ of all zero-
divisors. Further, let $S$ be the ring $M_{n}(R)$ ($n>1$). Via the natural
ring embedding $R\hookrightarrow S,\,a\to aE$, identify $R$ with the subring
$RE\subset S$ (and, thus, $O(R)$ — with a subgraph of the graph $O(S)$). For
$A\in S$, set $I_{A}:=\operatorname{Ann}(\det A)\lhd R$.

The graph $O(R)$ is connected and has diameter at most $3$ (see Theorem 2.3 in
[2, § 2]). Besides, if $R$ is a skew body, then

  1. 1)

once $n=2$, the graph $O(S)$ is disconnected and all its connected components
have diameters $\leqslant 2$;

  2. 2)

once $n\geqslant 3$, the graph $O(S)$ is connected and has diameter $4$.

These results are obtained in [9, § 4] for fields (Lemma 4.1 and Theorem 4.5
respectively), and in [10, § 2] are generalized to arbitrary skew-fields
(Lemma 2.2 and Theorem 2.1 respectively). They are also shifted to integral
domains (by reducing to the field of fractions).

######  Theorem 3.1.

For any matrix $A\in S$ and proper ideal $I\lhd R$ containing $\det A$, there
exists a matrix $B\in S\setminus\bigl{(}M_{n}(I)\bigr{)}$ such that $AB,BA\in
M_{n}(I)$.

$\square\quad$ For $m\in\mathbb{N}$, set $Q_{m}:=\\{1,\ldots,m\\}$ and
$P_{m}:=(1,\ldots,m)\in\mathbb{N}^{m}$.

Consider all triples $(k,P^{\prime},P^{\prime\prime})$ ($k\geqslant 0$,
$P^{\prime},P^{\prime\prime}\in(Q_{n})^{k}$) satisfying the relation
$\det(A^{P^{\prime}}_{P^{\prime\prime}})\notin\nobreak I$. For each of them,
numbers are repeated neither in $P_{1}$, nor in $P_{2}$, and, by condition,
$k<n$. Besides, at least one of such triples exists: for $k:=0$ and empty
tuples $P^{\prime},P^{\prime\prime}$, the corresponding $(0\times 0)$-matrix
has determinant $1\notin I$. Hence, we can fix one of these triples with the
largest possible $k$, and then $0\leqslant k<n$, $m:=k+1\in Q_{n}$.

######  Case 1).

$P^{\prime}=P^{\prime\prime}=P_{k}$.

By construction, $\det(A^{P_{k}}_{P_{k}})\notin I$. Further, set
$C:=A^{P_{m}}_{P_{m}}\in M_{m}(R)$,

$B:=\begin{pmatrix}\widehat{C}&0^{m}_{n-m}\\\
0^{n-m}_{m}&0_{n-m}\end{pmatrix}\in S.$

Then $b_{m,m}=\det(A^{P_{k}}_{P_{k}})\notin I$ implying $B\notin M_{n}(I)$.
Show that $AB,A^{T}B^{T}\in M_{n}(I)$, i. e. that, for any $p,q\in Q_{n}$, the
matrix entries $(AB)_{p,q}$ and $(A^{T}B^{T})_{p,q}$ belong to $I$. Assume
that $p\in Q_{n}$ and $q\in Q_{m}$ (otherwise
$(AB)_{p,q}=(A^{T}B^{T})_{p,q}=0$). Let $P\in(Q_{n})^{m}$ be the tuple
obtained form $P_{m}$ by changing the $q$-th element with $p$. Due to
maximality of $k$ and the inequality $m>k$, we have
$\det(A^{P}_{P_{m}}),\det(A^{P_{m}}_{P})\in I$,

$\displaystyle(AB)_{p,q}=\sum\limits_{{i\in
Q_{n}}}(a_{p,i}b_{i,q})=\sum\limits_{{i\in
Q_{m}}}\bigl{(}a_{p,i}(\widehat{C})_{i,q}\bigr{)}=$
$\displaystyle\sum\limits_{{i\in
Q_{m}}}\bigl{(}a_{p,i}(\widetilde{C})_{q,i}\bigr{)}=\det(A^{P}_{P_{m}})\in I;$
$\displaystyle(A^{T}B^{T})_{p,q}=\sum\limits_{{i\in
Q_{n}}}\bigl{(}(A^{T})_{p,i}(B^{T})_{i,q}\bigr{)}=$
$\displaystyle\sum\limits_{{i\in
Q_{m}}}\bigl{(}a_{i,p}(\widetilde{C})_{i,q}\bigr{)}=\det(A^{P_{m}}_{P})\in I.$

Thereby, it is proved that $AB,(BA)^{T}=A^{T}B^{T}\in M_{n}(I)$ implying
$BA\in M_{n}(I)$.

######  Case 2).

$P^{\prime},P^{\prime\prime}\in(Q_{n})^{k}$ are arbitrary tuples.

In each of the tuples $P^{\prime}$ and $P^{\prime\prime}$ all numbers are
distinct. Hence, via suitable permutations of rows and columns, one can obtain
from $A$ a matrix $A_{0}$ satisfying Case 1) with the same $k$. By proved
above, there exists a matrix $B_{0}\in S\setminus\bigl{(}M_{n}(I)\bigr{)}$
such that $A_{0}B_{0},B_{0}A_{0}\in\nobreak M_{n}(I)$. At that, there exist
monomial (therefore, invertible) matrices $C_{1},C_{2}\in S$ such that
$A=\nobreak C_{1}A_{0}C_{2}^{-1}$. Left (resp. right) multiplying a matrix by
a monomial one permutes its rows (resp. columns), and, consequently,
$B:=C_{2}B_{0}C_{1}^{-1}\in S\setminus\bigl{(}M_{n}(I)\bigr{)}$, $AB=\nobreak
C_{1}(A_{0}B_{0})C_{1}^{-1}\in\nobreak M_{n}(I)$,
$BA=C_{2}(B_{0}A_{0})C_{2}^{-1}\in M_{n}(I)$. ∎

######  Corollary 3.1.

If $A\in S$ and $c\in I_{A}^{*}$, then, in the subset $(cS)^{*}\subset S$,
there exists an element orthogonal to $A$.

$\square\quad$ By condition, $I:=\operatorname{Ann}(c)\lhd R$ is a proper
ideal containing $\det A$. According to Theorem 3.1, there exists a matrix
$B\in S\setminus\bigl{(}M_{n}(I)\bigr{)}$ such that $AB,BA\in M_{n}(I)$. Thus,
$C:=cB\neq 0$ and $c(AB)=c(BA)=0$, i. e. $C\in(cS)^{*}$ and $AC=CA=0$. ∎

######  Lemma 3.1.

For any $A\in S$, the following conditions are equivalent:

  1. 1)

$\det A\in Z_{R}$;

  2. 2)

$I_{A}\neq 0$;

  3. 3)

in $S^{*}$, there exists an element orthogonal to $A$;

  4. 4)

$A$ is a two-sided zero-divisor;

  5. 5)

$A$ is a zero-divisor.

$\square\quad$ The implications
$\text{\ref{dedi}}\Leftrightarrow\text{\ref{ine}}$ and
$\text{\ref{orto}}\Rightarrow\text{\ref{zdi2}}\Rightarrow\text{\ref{zdi1}}$
obviously follow from definitions, and the implication
$\text{\ref{ine}}\Rightarrow\text{\ref{orto}}$ — from Corollary 3.1.

Prove the implication $\text{\ref{zdi1}}\Rightarrow\text{\ref{dedi}}$. Suppose
that, without loss of generality, $A$ is a left zero-divisor, i. e. that
$AB=0$ for some $B\in S^{*}$. Then $\widehat{A}A=(\det A)E$ implying $(\det
A)B=\nobreak\widehat{A}AB=0$. It remains to use non-triviality of $B$. ∎

######  Corollary 3.2.

All zero-divisors in $S$ are two-sided.

Let $Z_{S}\subset S$ be the subset of all elements $A\in S$ satisfying each of
the equivalent conditions 1)—5) of Lemma 3.1, i. e. the set of all zero-
divisors of the ring $S$. Then the vertex set of the graph $O(S)$ is
$Z_{S}^{*}$.

Further, we will assume that $Z_{R}^{*}\neq\varnothing$.

######  Statement 3.1.

If $I\lhd R$ and $I\neq 0$, then $Z_{R}\cap I\neq\\{0\\}$.

$\square\quad$ Suppose that $Z_{R}\cap I=\\{0\\}$. There exist elements $b\in
I^{*}$ and $c\in Z_{R}^{*}$; then $bc\in\nobreak Z_{R}\cap I=\\{0\\}$. So,
$bc=0\neq c$ that implies $b\in Z_{R}\cap I^{*}=\varnothing$, a contradiction.
∎

######  Lemma 3.2.

If, for a subset $D\subset S$, the ideal $I:=\bigcap\limits_{{A\in
D}}I_{A}\lhd R$ is nonzero, then there exist elements $b\in Z_{R}^{*}$ and
$C_{A}\in S^{*}$, $A\in D$, such that $bE\perp C_{A}\perp A$ ($A\in D$).

$\square\quad$ According to Statement 3.1, the ideal $I$ contains an element
$c\in Z_{R}^{*}$. Then $bc=0$ where $b\in Z_{R}^{*}$. Further, for any $A\in
D$, we have $c\in I_{A}^{*}$ and, by Corollary 3.1, there exist an element
$C_{A}\in(cS)^{*}$ orthogonal to $A$; at that, $bC_{A}\in bcS=0$, $bE\perp
C_{A}$. ∎

######  Corollary 3.3.

  1. 1)

For any $A\in Z_{S}^{*}$, we have $d\bigl{(}A,O(R)\bigr{)}\leqslant 2$.

  2. 2)

If $A_{1},A_{2}\in Z_{S}^{*}$ and $I_{A_{1}}\cap I_{A_{2}}\neq 0$, then
$d(A_{1},A_{2})\leqslant 4$.

$\square\quad$ It suffices to apply Lemma 3.2 to the subsets
$\\{A\\},\\{A_{1},A_{2}\\}\subset S$. ∎

######  Lemma 3.3.

If $A_{i}\in Z_{S}^{*}$, $c_{i}\in I_{A_{i}}^{*}$ ($i=1,2$) and
$c_{1}c_{2}=0$, then $d(A_{1},A_{2})\leqslant 3$.

$\square\quad$ By Corollary 3.1, for each $i=1,2$, there exists an element
$C_{i}\in(c_{i}S)^{*}$ such that $C_{i}\perp A_{i}$. In this case,
$C_{1}C_{2},C_{2}C_{1}\in c_{1}c_{2}S=0$, $C_{1}\perp C_{2}$. ∎

Definition. We will say that an ideal $I\lhd R$ does not have zero-divisors if
$I^{*}I^{*}\not\ni 0$, i. e. if the ring111In general, without identity. $I$
does not have zero-divisors.

######  Lemma 3.4.

If $A_{1},A_{2}\in Z_{S}^{*}$ and $d(A_{1},A_{2})>3$, then $I_{A_{i}}$
($i=1,2$) is the same ideal without zero-divisors.

$\square\quad$ According to Lemma 3.3, $I_{A_{1}}^{*}I_{A_{2}}^{*}\not\ni 0$.
It remains to prove that $I_{A_{1}}=I_{A_{2}}$.

Suppose that $I_{A_{1}}\neq I_{A_{2}}$. Without loss of generality, assume
that there exists an element $c\in I_{A_{1}}\setminus I_{A_{2}}$. Setting
$a:=\det A_{2}$, we have $I_{A_{2}}=\operatorname{Ann}(a)$, $b:=ca\in
I_{A_{1}}^{*}$ and $bI_{A_{2}}=caI_{A_{2}}=0$ implying
$bI_{A_{2}}^{*}\subset\\{0\\}\cap(I_{A_{1}}^{*}I_{A_{2}}^{*})=\varnothing$,
$I_{A_{2}}^{*}=\varnothing$, $I_{A_{2}}=0$, a contradiction. ∎

######  Theorem 3.2.

The graph $O(S)$ is connected and has diameter at most $4$.

$\square\quad$ Suppose that there exist elements $A_{1},A_{2}\in Z_{S}^{*}$
satisfying the inequality $d(A_{1},A_{2})>\nobreak 4$. By Lemma 3.4, $0\neq
I_{A_{1}}=I_{A_{2}}=I_{A_{1}}\cap I_{A_{2}}$ that contradicts with Corollary
3.3. ∎

######  Theorem 3.3.

We have $\operatorname{diam}\bigl{(}O(S)\bigr{)}\geqslant 3$, the strict
inequality being equivalent to the existence of an ideal
$\operatorname{Ann}(a)\lhd R$ ($a\in Z_{R}$) without zero-divisors.

$\square\quad$ Similarly with examples from [9, 10] giving lower estimates of
the diameter, for an arbitrary $a\in Z_{R}$, set $I:=\operatorname{Ann}(a)\lhd
R$ and $A:=J_{n}+aE_{n1}\in S$. Note that

  * •

$A,A^{T}\in Z_{S}^{*}$, $O_{S}(A)=I^{*}E_{1n}$, $O_{S}(A^{T})=I^{*}E_{n1}$;

  * •

$a_{12}=1\neq a_{21}$, $(AA^{T})_{11}=1$ and $O_{S}(A)\cap
O_{S}(A^{T})=\varnothing$, that implies $d(A,A^{T})\geqslant 3$;

  * •

if $I^{*}I^{*}\not\ni 0$, then
$\bigl{(}O_{S}(A)\bigr{)}\bigl{(}O_{S}(A^{T})\bigr{)}=(I^{*}I^{*})E_{11}\not\ni
0$ and, hence, $d(A,A^{T})\geqslant 4$.

Due to mentioned above, $\operatorname{diam}\bigl{(}O(S)\bigr{)}\geqslant 3$,
the strict inequality following from the existence of an ideal
$\operatorname{Ann}(a)\lhd R$ ($a\in Z_{R}$) without zero-divisors.
Conversely, in the case of the strict inequality, by Lemma 3.4, for some
elements $A\in Z_{S}$ and $a:=\det A\in Z_{R}$, the ideal
$I_{A}=\operatorname{Ann}(a)\lhd R$ does not have zero-divisors. ∎

Now the main Theorem 1.1 follows from Theorems 3.2 and 3.3, and Corollary 3.3.
It implies (see (2.1)) the statements 1) and 2) of Theorem 1.2. Let us prove
3).

Suppose that $\operatorname{rad}\bigl{(}O(S)\bigr{)}=2$. There exist elements
$C\in Z_{S}^{*}$, $c\in R^{*}$ and $k,l\in Q_{n}$ such that $d(C,A)\leqslant
2$ ($A\in Z_{S}^{*}$) and $c_{kl}=c$. Further, there exists a permutation
$\sigma\in S_{n}$ such that $m:=\sigma(k)\neq l$.

Let $a\in Z_{R}$ be an arbitrary element.

Set $I:=\operatorname{Ann}(a)\lhd R$ and $A:=\Bigl{(}\sum\limits_{{i\neq
k}}E_{i,\sigma(i)}\Bigr{)}+aE_{km}\in S$. Note that

  * •

$A\in Z_{S}^{*}$, $O_{S}(A)=I^{*}E_{mk}$;

  * •

$A\neq C$ (otherwise $a_{kl}=c\neq 0$, $m=l$);

  * •

$(m,k)\neq(k,l)$ (otherwise $m=k=l$), that implies $C\notin O_{S}(A)$.

Thus, $d(C,A)=2$, so, there exists an element $B\in Z_{S}^{*}$ orthogonal to
$C$ and $A$. We have $B=bE_{mk}$ where $b\in I^{*}$. Meanwhile, $BC=0$,
$0=(BC)_{ml}=bc$, $b\in\operatorname{Ann}(c)\cap I^{*}$.

Due to arbitrariness of $a\in Z_{R}$, the element $c\in R^{*}$ satisfies
(1.2).

Conversely, assume that (1.2) holds for some $c\in R^{*}$. Show that the
element $C:=cE\in\nobreak S^{*}$ satisfies, for each $A\in Z_{S}^{*}$, the
inequality $d(C,A)\leqslant 2$.

Let $A\in Z_{S}^{*}$ be an arbitrary element. Then $\det A\in Z_{R}$, and, by
(1.2), there exists an element $b\in I_{A}^{*}$ such that $cb=0$. Further,
according to Corollary 3.1, there exists an element $B\in(bS)^{*}$ orthogonal
to $A$; in this case, $cB\in cbS=0$, $C\in Z_{S}^{*}$, $C\perp B\perp A$,
$d(C,A)\leqslant 2$.

So, Theorem 1.2 is completely proved.

## Acknowledgements

The author is grateful to Prof. E. B. Vinberg for exciting interest to
algebra.

The author dedicates the article to E. N. Troshina.

## References

  * [1] Beck I. Coloring of commutative rings // J. Algebra. 1988\. Vol. 116. Pp. 208–226.
  * [2] Anderson D. F., Livingston P. S. The zero-divisor graph of a commutative ring // J. Algebra. 1999\. Vol. 217. Pp. 434–447.
  * [3] Anderson D. F., Frazier A., Lauve A., Livingston P. S. The zero-divisor graph of a commutative ring, II // Lect. Notes Pure Appl. Math. 2001\. Vol. 220, Marcel Dekker, New York. Pp. 61–72.
  * [4] Akbari S., Maimani H. R., Yassemi S. When zero-divisor graph is planar or a complete $r$-partite graph // J. Algebra. 2003\. Vol. 270. Pp. 169–180.
  * [5] Akbari S., Mohammadian A. On the zero-divisor graph of a commutative ring // J. Algebra. 2004\. Vol. 274. Pp. 847–855.
  * [6] Akbari S., Mohammadian A. Zero-divisor graphs of non-commutative rings // J. Algebra. 2006\. Vol. 296. Pp. 462–479.
  * [7] Akbari S., Mohammadian A. On zero-divisor graphs of finite rings // J. Algebra. 2007\. Vol. 314. Pp. 168–184.
  * [8] Redmond S. P. Central sets and radii of the zero-divisor graphs of commutative rings // Comm. in Algebra. 2006\. Vol. 34 (is. 7). Pp. 2389–2401.
  * [9] Bakhadly B. R., Guterman A. E., Markova O. V. Graphs defined by orthogonality // Computational methods and algorithms. Part XXVII, Zap. Nauchn. Sem. POMI, St. Petersburg. 2014\. Vol. 428. Pp. 49-80; J. Math. Sci. 2015\. Vol. 207 (is. 5). Pp. 698–717.
  * [10] Guterman A. E., Markova O. V. Orthogonality graphs of matrices over skew fields // Computational methods and algorithms. Part XXX, Zap. Nauchn. Sem. POMI, St. Petersburg. 2017\. Vol. 463. Pp. 81–93; J. Math. Sci. 2018\. Vol. 232 (is. 6). Pp. 797–804.


# Time-dependent condensate fraction in an analytical model

A. Simon  G. Wolschin<EMAIL_ADDRESS>Institute for
Theoretical Physics, Heidelberg University, Philosophenweg 12–16, D-69120
Heidelberg, Germany, European Union

###### Abstract

We apply analytical solutions of a nonlinear boson diffusion equation (NBDE)
that include boundary conditions at the singularity to calculate the time
evolution of the entropy during evaporative cooling of ultracold atoms, and
the time-dependent condensate fraction. For suitable initial conditions it is
found to agree with available data on 23Na.

##  I Introduction

Shortly after the discovery of Bose-Einstein condensation in ultracold atoms
of rubidium [1], sodium [2] and lithium [3, 4], models and theories to
calculate the time-dependent condensate fraction were developed, see Refs. [5,
6, 7] and related works. These approaches usually consider a quantum treatment
such as a nonlinear Schrödinger equation for the condensate and an
equilibrium- or nonequilibrium-statistical description for the thermal
component.

Discrepancies between early theoretical results [5] and data for the time-
dependent condensate fraction [8] were partly due to the fact that the
noncondensed particles were represented by an equilibrium particle reservoir
with a fixed chemical potential – an assumption that does not correspond to
the nonequilibrium initial state that is generated through evaporative
cooling.

The relaxational dynamics of the nonthermal component that occurs in the
presence of the mean field of the condensate was later taken into account
through kinetic equations. In particular, a semiclassical Boltzmann equation
[9, 10, 11, 12, 13, 7] was solved numerically in the ergodic approximation,
which has been used by many authors [14, 15, 16, 17, 18, 19]. In this
approximation, the population of a given state depends only on its energy.
Such an approach provides a detailed description of the growth of the
condensate based on numerical solutions of the equations for the dynamics of
the thermal cloud and its mean-field interaction with the condensate [7],
yielding satisfactory – though not perfect – agreement with the available data
on the time-dependent condensate fraction. These numerical approaches to
thermalization and condensate formation usually rely on an artificial seed
condensate [7].

In view of such involved numerical approaches it is interesting to have a
simple, exactly solvable model for the physics of condensate formation in cold
atoms. For bosons, such a model has been proposed in Ref. [20] and adapted to
cold quantum gases in Refs. [21, 22]. The model is based on a nonlinear boson
diffusion equation (NBDE) which has been derived from the quantum Boltzmann
collision term. For energy-independent transport coefficients, it is
structurally simple, but still complicated to solve exactly due to the
nonlinearity in the drift term, which corresponds to Bose stimulation and
causes the system to reach the Bose-Einstein equilibrium distribution for
sufficiently large times, with a singularity in the infrared. The solutions
represent, in particular, evaporative cooling from an initial temperature
$T_{\text{i}}$ to a final temperature $T_{\text{f}}$ that may be below the
critical value $T_{\text{c}}$. We explore the solutions in this work to
calculate the time-dependent condensate fraction based on the analytic
solutions and on particle-number conservation in the cloud and the condensate,
together with the time evolution of the entropy in successive and single
cooling steps.

In the next section, the NBDE and its exact solutions [21, 22] with the
consideration of boundary conditions at the singularity are briefly reviewed.
In Section 3, the analytic time-dependent solutions of the nonlinear equation
with constant transport coefficients are evaluated and compared to numerical
results. In Section 4, the time evolution of the entropy in evaporative
cooling is considered for several sequential cooling steps, as well as for a
single step that leads to condensate formation. In Section 5, the time-
dependent chemical potential is evaluated from the condition of particle-
number conservation and the condensate fraction is discussed based on the
analytic NBDE solutions. Due to the quantum-statistical properties of the
bosonic system that are encoded in the NBDE, a condensate forms and grows if
the final temperature is below the critical value. With parameters adapted to
the evaporative cooling of 23Na, a comparison of the time-dependent condensate
fraction with MIT data is presented. The conclusions are drawn in Section 6.

##  II Nonlinear boson diffusion equation and analytic solutions

The system of ultracold atoms is viewed as a time-dependent mean-field
supplemented by a collision term. The atomic density is usually so low that
interactions involving more than two particles at a time can be ignored. Only
$s$-wave scattering contributes for bosons at ultracold temperatures, such
that the interaction is a contact interaction with a coupling strength that is
proportional to the $s$-wave scattering length.

The $N$-body density operator $\hat{\rho}_{N}(t)$ can be written based on $N$
single-particle wave functions of the atoms which are solutions of the time-
dependent Hartree-Fock equations supplemented by a time-irreversible collision
term $K_{N}(t)$ that accounts for the thermalization through random two-body
collisions ($\hbar=c=1$)

$i\,\frac{\partial\hat{\rho}_{N}(t)}{\partial
t}=\big{[}\hat{H}_{\text{HF}}(t),\hat{\rho}_{N}(t)\big{]}+i\hat{K}_{N}(t)\;$
(1)

where $\hat{H}_{\text{HF}}(t)$ is the self-consistent Hartree-Fock mean field
of the atoms, with an external potential for ultracold atoms in a trap.

In this work we do not consider the full many-body problem and its reduction
to the one-body level, but rather an approximate version that starts from the
reduced ensemble-averaged single-particle density operator
$\bar{\rho}_{1}(t)$. Its diagonal elements can be interpreted as the
probability for a particle to be in a state $|\alpha\rangle$ with energy
$\epsilon_{\alpha}$

$\big{(}\bar{\rho}_{1}(t)\big{)}_{\alpha,\alpha}=\langle
n(\epsilon_{\alpha},t)\rangle\equiv n_{\alpha}(\epsilon,t)\;.$ (2)

The total number of particles is $N=\sum_{\alpha}n_{\alpha}$, and we neglect
here the off-diagonal terms of the density matrix. The occupation-number
distribution $n_{\alpha}(t)$ in a finite Bose system obeys a Boltzmann-like
collision term [20]

$\displaystyle\frac{\partial n_{\alpha}}{\partial
t}=\sum_{\beta,\gamma,\delta}\overline{V^{2}_{\alpha\beta\gamma\delta}}\,\mathcal{G}_{\alpha\beta\gamma\delta}\big{[}(1+n_{\alpha})(1+n_{\beta})n_{\gamma}n_{\delta}$
(3) $\displaystyle-(1+n_{\gamma})(1+n_{\delta})n_{\alpha}n_{\beta}\big{]}$

where $\overline{V^{2}_{\alpha\beta\gamma\delta}}$ is the second moment of the
pairwise atom-atom interaction. $\mathcal{G}_{\alpha\beta\gamma\delta}$ is an
energy-conserving function which has a finite width, because it is the total
energy of mean field plus collision term that must be conserved. It becomes a
$\delta$-function as in the Boltzmann-Nordheim equation [23] only if the mean-
field energy content stays constant in time – which is, however, usually not
the case in cold-atom systems. Whereas in an infinite homogeneous system that
is described using a quantum Boltzmann equation with an energy-conserving
$\delta$-function an infinite amount of time is required to create a
condensate, this is not expected to be the case for an energy-conserving
function that has a finite width. Hence, the theoretical framework is
consistent with a finite condensate formation time as required from
experiment.

This equation has been transformed into a nonlinear partial differential
equation in Refs. [20, 21]. One starts by defining transition probabilities
from state $\gamma$ to $\alpha$ as

$W_{\gamma\rightarrow\alpha}=\sum_{\beta,\delta}\overline{V^{2}_{\alpha\beta\gamma\delta}}\,\mathcal{G}_{\alpha\beta\gamma\delta}(1+n_{\beta})\,n_{\delta}\;.$
(4)

The collision term can be written as a master equation where the transitions
into and out of the state with energy $\epsilon_{\alpha}$ are more explicit

$\frac{\partial n_{\alpha}}{\partial
t}=(1+n_{\alpha})\sum_{\gamma}W_{\gamma\rightarrow\alpha}n_{\gamma}-n_{\alpha}\sum_{\gamma}W_{\alpha\rightarrow\gamma}(1+n_{\gamma})\;.$
(5)

Introducing the density of states $g_{\alpha}=g(\epsilon_{\alpha})$, such that
$W_{\gamma\rightarrow\alpha}=W_{\gamma,\alpha}g_{\alpha}$ and using the fact
that quantum particles are indistinguishable
$W_{\alpha,\gamma}=W_{\gamma,\alpha}=W(\epsilon_{\alpha},\epsilon_{\gamma})$
one arrives at

$\frac{\partial n_{\alpha}}{\partial
t}=\int_{0}^{\infty}W_{\alpha,\gamma}\left[g_{\alpha}(1+n_{\alpha})n_{\gamma}-g_{\gamma}(1+n_{\gamma})n_{\alpha}\right]\text{d}\epsilon_{\gamma}\;.$
(6)

Defining the transport coefficients $v\equiv v\,(\epsilon_{\alpha},t)$ and
$D\equiv D\,(\epsilon_{\alpha},t)$ as first and second moments of the
transition probability [20], we have obtained the nonlinear boson diffusion
equation (NBDE) for the single-particle expectation-value occupation-number
distribution of the energy eigenstates $\epsilon_{\alpha}$, $n\equiv
n_{\alpha}\equiv\langle n(\epsilon_{\alpha},t)\rangle$, in Refs. [20, 21] as
111The derivative-term of the diffusion coefficient has been modified as
compared to Refs. [20, 21, 22] in order to secure the correct stationary
solution

$\frac{\partial n}{\partial
t}=-\frac{\partial}{\partial\epsilon}\Bigl{[}v\,n\,(1+n)+n\frac{\partial
D}{\partial\epsilon}\Bigr{]}+\frac{\partial^{2}}{\partial\epsilon^{2}}\bigl{[}D\,n\bigr{]}\,.$
(7)

The drift term $v(\epsilon,t)$ accounts for dissipative effects, the term
$D(\epsilon,t)$ for diffusion of particles in the energy space. The many-body
physics is contained in these transport coefficients, which depend on energy,
time, and the second moment of the interaction.

A prerequisite for the reduction to $1+1$ dimensions in the above formulation
is spatial and momentum isotropy. This corresponds to the assumption of
sufficient ergodicity, which has been widely discussed in the literature [14,
15, 16, 17, 18, 19]. For the thermal cloud of cold atoms around a Bose-
Einstein condensate (BEC), it is expected to be a reasonable assumption, even
though the condensate in a trap is spatially anisotropic. Regarding the role
of different spatial dimensions in view of BEC formation, this enters our
present formulation only through the density of states, which differs
according to the number of spatial dimensions, and the type of confinement.
The model calculations in this work are for a 3d system, and we shall
investigate results for the density of states of a free Bose gas, and bosonic
atoms confined in a harmonic trap. One-dimensional systems where no BEC should
be formed have not yet been studied.

To derive the stationary solution $n_{\infty}(\epsilon)$ for the equation with
variable transport coefficients, we rewrite Eq. (1) and set the time
derivative to zero

$0=\frac{\partial}{\partial\epsilon}{\left[v\,n_{\infty}\,(1+n_{\infty})-D\,\frac{\partial
n_{\infty}}{\partial\epsilon}\right]}\,,$ (8)

such that

$v\,n_{\infty}(1+n_{\infty})-D\,\frac{\partial
n_{\infty}}{\partial\epsilon}=c_{1}\;.$ (9)

Dividing by $n_{\infty}(1+n_{\infty})D$ and integrating over $\epsilon$ yields

$\int\frac{\text{d}n_{\infty}/\text{d}\epsilon}{n_{\infty}(1+n_{\infty})}\text{d}\epsilon=\int\left(\frac{v}{D}-\frac{c_{1}}{n_{\infty}(1+n_{\infty})D}\right)\text{d}\epsilon\;.$
(10)

We integrate the l.h.s, resulting in
$[\ln(n_{\infty})-\ln(1+n_{\infty})]=\ln(1-\frac{1}{1+n_{\infty}})$ plus an
integration constant $c_{2}$. Solving for $n_{\infty}$ one obtains

$n_{\infty}=\left[\exp\left(\int\left(-\frac{v}{D}+\frac{c_{1}}{n_{\infty}(1+n_{\infty})D}\right)\text{d}\epsilon+c_{2}\right)-1\right]^{-1}.\;$
(11)

In order to reduce to a Bose-Einstein distribution, $c_{1}=0$ is required, and
the ratio $v/D$ must have no energy dependence for ${t\rightarrow\infty}$ so
that it can be pulled out of the integral. It follows that
$\lim\limits_{t\to\infty}[-v(\epsilon,t)/D(\epsilon,t)]\equiv 1/T$ and
$c_{2}\equiv-\mu/T$ such that the stationary distribution equals the Bose-
Einstein equilibrium distribution ($k_{\text{B}}=1$)

$n_{\infty}(\epsilon)=n_{\text{eq}}(\epsilon)=\frac{1}{e^{(\epsilon-\mu)/T}-1}$
(12)

with the chemical potential $\mu<0$ in a finite Bose system.

In the limit of energy-independent transport coefficients the nonlinear boson
diffusion equation for the occupation-number expectation-value distribution
$n(\epsilon,t)$ becomes

$\frac{\partial n}{\partial
t}=-v\,\frac{\partial}{\partial\epsilon}\Bigl{[}n\,(1+n)\Bigr{]}+D\,\frac{\partial^{2}n}{\partial\epsilon^{2}}\,.$
(13)

As in case of Eq. (7), the thermal equilibrium distribution $n_{\text{eq}}$ is
a stationary solution with $\mu<0$ and $T=-D/v$.

In spite of its simple structure, the NBDE with constant transport
coefficients thus preserves the essential features of Bose-Einstein statistics
which are contained in the bosonic Boltzmann equation. For a given initial
condition $n_{\text{i}}(\epsilon)$, it can be solved exactly using the
nonlinear transformation outlined in Ref. [21]. The resulting solution is

$\displaystyle
n(\epsilon,t)=-\frac{D}{v}\frac{\partial}{\partial\epsilon}\ln{\mathcal{Z}(\epsilon,t)}-\frac{1}{2}=-\frac{D}{v}\frac{1}{\mathcal{Z}}\frac{\partial\mathcal{Z}}{\partial\epsilon}-\frac{1}{2}$
(14)

where the time-dependent partition function ${\mathcal{Z}(\epsilon,t)}$ obeys
a linear diffusion equation

$\displaystyle\frac{\partial}{\partial
t}{\mathcal{Z}}(\epsilon,t)=D\frac{\partial^{2}}{\partial\epsilon^{2}}{\mathcal{Z}}(\epsilon,t)\,.$
(15)

The partition function $\mathcal{Z}(\epsilon,t)$

$\displaystyle\mathcal{Z}(\epsilon,t)=a(t)\int_{-\infty}^{+\infty}G(\epsilon,x,t)\,F(x)\,\text{d}x$
(16)

is an integral over Green’s function $G(\epsilon,x,t)/\sqrt{4\pi\,Dt}$ of Eq.
(15), and an exponential function $F(x)$ that depends on the initial
occupation-number distribution $n_{\text{i}}$

$\displaystyle
F(x)=\exp\Bigl{[}-\frac{1}{2D}\bigl{(}vx+2v\int_{0}^{x}n_{\mathrm{i}}(y)\,\text{d}y\bigr{)}\Bigr{]}\,.$
(17)

The time-dependent prefactor $a(t)=1/\sqrt{4\pi\,Dt}$ in Eq. (16) cancels out
when taking the logarithmic derivative in Eq. (14), and for the same reason
the definite integral over the initial conditions taken at the lower limit in
Eq. (17) drops out in the calculation of $n(\epsilon,t)$. Hence, the definite
integral in Eq. (17) can be replaced [25] by the indefinite integral
$A_{\mathrm{i}}(x)$ over the initial distribution with
$\partial_{x}A_{\mathrm{i}}(x)=n_{\mathrm{i}}(y)$.

With these modifications that do not affect the accuracy, it becomes possible
to compute the partition function and the overall solution for the occupation-
number distribution function Eq. (14) analytically.

When solving the problem with the free Green’s function of Eq. (15), the
physically correct solution with the Bose-Einstein equilibrium limit is
attained in the UV region, but not in the IR [20]. To solve this problem, one
has to consider the boundary conditions at the singularity $\epsilon=\mu<0$
[22]. For fixed chemical potential $\mu$, an initial temperature
$T_{\text{i}}$, and a final temperature $T_{\text{f}}=-D/v$, the combined
initial- and boundary value problem has been solved exactly in Refs. [25, 22]
using the above nonlinear transformation from Eq. (14) [20, 22], and an
infinite series expansion. The analytic solutions are briefly reconsidered
before we use them for further calculations and in particular, to obtain the
time-dependent condensate fraction.

With the boundary condition $\lim_{\epsilon\downarrow\mu}n(\epsilon,t)=\infty$
$\forall$ $t$ at the singularity $\epsilon=\mu$ one obtains
$\mathcal{Z}(\mu,t)=0$, and the energy range is restricted to
$\epsilon\geq\mu$. This requires a new Green’s function [22] that equals zero
at $\epsilon=\mu$ $\forall\,t$. It can be written as

$\displaystyle{G}(\epsilon,x,t)=G_{\text{free}}(\epsilon-\mu,x,t)-G_{\text{free}}(\epsilon-\mu,-x,t)\,,$
(18)

and the partition function with this boundary condition becomes

$\displaystyle{\mathcal{Z}}(\epsilon,t)=\int_{0}^{\infty}{G}(\epsilon,x,t)\,F(x+\mu)\,\text{d}x\,.$
(19)

The function $F$ remains unaltered with respect to Eq. (17), but its argument
is shifted by the chemical potential.

Figure 1: The argument $\ln{F(x)}$ of the exponential function $F(x)$ in Eq.
(17) (solid curve) with a singularity at $x=\mu<0$ in the definite integral of
an initial distribution $n_{\text{i}}$ given by a Bose-Einstein distribution
with $T=240$ nK and $\mu=-8$ nK that is truncated at $30$ nK (dotted vertical
line).

As an initial condition that is appropriate for a schematic description of
evaporative cooling [1, 26, 17] to demonstrate our method, a truncated thermal
equilibrium distribution that is cut off at a maximum energy
$\epsilon_{\text{i}}$ beyond which high-velocity atoms are removed has been
chosen in our previous works

$n_{\text{i}}(\epsilon)=\frac{1}{e^{(\epsilon-\mu)/T}-1}\,\theta(1-\epsilon/\epsilon_{\text{i}})\,.$
(20)

The distribution is truncated at a maximum energy $\epsilon_{\text{i}}$ beyond
which high-velocity atoms are removed. In a trapped ultracold gas, these atoms
leave the trap: The system of atoms in the trap is thus not isolated and its
entropy decreases in the course of cooling. In the subsequent re-
thermalization (and for $T_{\text{f}}<T_{\text{c}}$, condensate formation),
the entropy of the atoms in the trap rises again, but remains below the
initial thermal value. This will be discussed in more detail in Section 5.

For such a truncated initial thermal distribution, the logarithm of the
function $F(x)$ with a singularity at $x=\mu$
(before shifting it by $\mu$) is shown in Fig. 1 with the following set of
parameters: $T_{\text{i}}=240$ nK, $\epsilon_{\text{i}}=30$ nK, $D=100\,$(nK)2
ms-1, $v=-1\,$nK ms-1, $T_{\text{f}}=-D/v=100$ nK. Due to the singularity,
$F(x)$ vanishes at $x=\mu=-8$ nK. The logarithm of $F(x)$ is continuous, but
not differentiable at $x=\epsilon_{\text{i}}$, which holds the key to the
equilibration in the UV region.

Using Eq. (20) as initial distribution, the occupation-number distribution can
be evaluated exactly for $T_{\text{i}}\neq T_{\text{f}}$ with the time-
dependent partition function [25]

$\displaystyle{\mathcal{Z}}(\epsilon,t)=\sqrt{4Dt}\,\exp\Bigl{(}-\frac{\mu}{2T_{\mathrm{f}}}\Bigr{)}\sum_{k=0}^{\infty}\binom{\frac{T_{\mathrm{i}}}{T_{\mathrm{f}}}}{k}\left(-1\right)^{k}\times$
$\displaystyle\Bigg{(}\text{e}^{\alpha_{k}^{2}Dt}\left[\text{e}^{\alpha_{k}(\epsilon-\mu)}\Lambda_{1}^{k}(\epsilon,t)-\text{e}^{\alpha_{k}(\mu-\epsilon)}\Lambda_{2}^{k}(\epsilon,t)\right]$
$\displaystyle+\exp\Bigl{(}\frac{(\mu-\epsilon_{i})k}{T_{\mathrm{i}}}\Bigr{)}\exp\Bigl{(}{\frac{Dt}{4T_{\mathrm{f}}^{2}}}\Bigr{)}\times$
$\displaystyle\Big{[}\exp\Bigl{(}\frac{\epsilon-\mu}{2T_{\mathrm{f}}}\Bigr{)}\Lambda_{3}(\epsilon,t)-\exp\Bigl{(}\frac{\mu-\epsilon}{2T_{\mathrm{f}}}\Bigr{)}\Lambda_{4}(\epsilon,t)\Big{]}\Bigg{)}\,$
(21)

where $\alpha_{k}=1/T_{\text{f}}-k/T_{\text{i}}$ and the auxiliary functions
are

$\displaystyle\Lambda_{1}^{k}(\epsilon,t)=$
$\displaystyle\,\text{erf}\,\Bigl{(}\frac{\epsilon-\mu+2Dt\alpha_{k}}{\sqrt{4Dt}}\Bigr{)}$
$\displaystyle\qquad-\,\text{erf}\,\Bigl{(}\frac{\epsilon-\epsilon_{i}+2Dt\alpha_{k}}{\sqrt{4Dt}}\Bigr{)}\,,$
(22) $\displaystyle\Lambda_{2}^{k}(\epsilon,t)=$
$\displaystyle\,\text{erf}\,\Bigl{(}\frac{\mu-\epsilon+2Dt\alpha_{k}}{\sqrt{4Dt}}\Bigr{)}$
$\displaystyle\qquad-\,\text{erf}\,\Bigl{(}\frac{2\mu-\epsilon-\epsilon_{i}+2Dt\alpha_{k}}{\sqrt{4Dt}}\Bigr{)}\,,$
(23) $\displaystyle\Lambda_{3}(\epsilon,t)=$
$\displaystyle\,\text{erfc}\,\Bigl{(}\frac{\epsilon_{i}-\epsilon+tv}{\sqrt{4Dt}}\Bigr{)}\,,$
(24) $\displaystyle\Lambda_{4}(\epsilon,t)=$
$\displaystyle\,\text{erfc}\,\Bigl{(}\frac{\epsilon-2\mu+\epsilon_{i}+tv}{\sqrt{4Dt}}\Bigr{)}\,.$
(25)

Figure 2:  Schematic representation of evaporative cooling in a a finite Bose
system based on the nonlinear evolution according to Eq. (13) starting from a
truncated Bose-Einstein distribution Eq. (20), upper solid curve with cutoff
at $\epsilon_{\text{i}}$ = 30 nK and $\mu=-8$ nK. The transport coefficients
are $D=100$ (nK)2 ms${}^{-1},~{}\,v=-1$ nK ms-1. The initial temperature is
$T=240$ nK, the final temperature $T_{\text{f}}=T=-D/v=100$ nK. The time
evolution of the single-particle occupation-number distributions is shown at
$t=0.001,0.01,0.1,0.4$, and $0.8$ ms (top to bottom in the IR). For
comparison, the numerical results using Matlab are shown as open circles.

To obtain the time-dependent occupation-number distribution function from Eq.
(14), the derivative $\partial{\mathcal{Z}}/\partial\epsilon$ is also
required. It can be calculated analytically [25], such that an exact
expression for the time-dependent occupation-number distribution function for
evaporative cooling results that can be directly compared with numerical
solutions of the NBDE, see the following section.

##  III Occupation-number distributions

The time-dependent analytic distribution functions that solve the NBDE exactly
with boundary conditions at the singularity are displayed in Fig. 2 for the
same initial conditions as in Fig. 1. The analytic results agree precisely
with numerical solutions (open circles) of the basic equation, and predict the
time-dependent cooling from a thermal distribution with temperature
$T_{\text{i}}$ that is truncated at $\epsilon_{\text{i}}$ to a Bose-Einstein
distribution with $T_{\text{f}}<T_{\text{i}}$, which is the thermal
distribution for $t\rightarrow\infty$. The approach is similar to the kinetic
theory of evaporative cooling in works such as Ref. [17], but now an
analytically solvable model is formulated.

The corresponding numerical solutions (circles in Fig. 2) have been obtained
using Matlab’s routine pdepe for the solution of partial differential
equations with given initial and boundary values [27]. Basic finite-difference
algorithms like the Crank-Nicolson method did not provide sufficient accuracy
for this nonlinear problem. The implementation of pdepe is based on the
algorithm described in Ref. [27] and is suited for (nonlinear) parabolic
partial differential equations. Due to the singularity at $\epsilon=\mu$, the
integration was started at $\mu+\delta<0$. Regarding the associate boundary
condition, we use the known stationary solution $n_{\infty}(\epsilon)$ to
determine the boundary values in the IR as $n_{\infty}(\mu+\delta)$, and
accordingly, in the UV. This minimizes the numerical errors at the boundaries.
The value of $\delta$ is chosen as small as possible, but large enough to
prevent numerical inaccuracies. In this particular case $\delta=0.2$ and an
upper boundary of $\epsilon_{\text{u}}=3\times T_{\text{i}}$ yields
satisfactory agreement with the exact solutions.

The values of the transport coefficients $v,D$ in this specific model
calculation have been derived from their relations to the equilibrium
temperature $T=-D/v$ and the equilibration time $\tau_{\text{eq}}=4D/(9v^{2})$
[20], with $T=100$ nK and $\tau_{\text{eq}}=44$ ms. They do not yet correspond
to a specific experimental situation, and the cut-off temperature
$\epsilon_{\text{i}}=30$ nK is articifially low to better demonstrate the
behaviour of the analytic solutions. For direct comparisons with data, these
values shall be adapted to the corresponding experiment, see Section 5.

In the infrared, thermalization occurs faster [22] than in case of a linear
relaxation ansatz: With the parameter set of Fig. 2, the thermal distribution
is reached within $t\simeq 1$ ms in the IR. The buildup of the thermal slope
in the UV is, however, slower in the nonlinear model as compared to the linear
relaxation ansatz, because the latter enforces an exponential approach to the
Boltzmann-like tail even for $\epsilon\rightarrow\infty$ [22].

Figure 3:  Convergence properties of the analytic NBDE solution with constant
transport coefficients according to Eqs. (14) and (21). The initial condition
is as in Fig. 2, upper solid curve with cutoff at $\epsilon_{\text{i}}$ = 30
nK. The time-dependent nonequilibrium solution is shown for $t=0.4$ ms, lower
solid curve. The convergence is displayed for maximum expansion parameters
$k_{\text{max}}=10,12,14,20$, and $200$ (from bottom to top) in the series
expansion of the partition function.

The convergence of our analytical method is displayed in Fig. 3 using the
above NBDE solution from Fig. 2 for $t=0.4$ ms as an example. Results
according to Eq. (14) are shown for $k_{\text{max}}=10,12,14,20$, and $200$
(from bottom to top) in the series expansion of the partition function Eq.
(21). Here, $k_{\text{max}}$ is the maximum expansion coefficient that is used
in the calculation. A solution for $k_{\text{max}}=40$ is already
indistinguishable from the plotted one for $k_{\text{max}}=200$. It is thus
sufficiently close to the exact solution, which would require
$k_{\text{max}}=\infty$. We shall later check the convergence properties also
for derived quantities such as the time-dependent chemical potential $\mu(t)$
that is needed to maintain particle-number conservation in the course of the
time evolution, and for the condensate fraction, in Section 5.

##  IV Evaporative cooling and time-dependent entropy

To achieve the phase transition to the condensate, successive evaporative
cooling is used, thereby removing high-velocity atoms. In the course of the
subsequent equilibration, the number of condensed particles rises due to a
transfer from the nonlinear kinetic region into the coherent region [28, 15,
29] and an isotropic tail [1, 30] develops in the ultraviolet, smearing out
the sharp cut that corresponds to evaporative cooling.

Solutions of the nonlinear boson diffusion equation fulfil the physically
reasonable condition [22] that the entropy drops in the course of an
evaporative cooling step and then increases with time towards the equilibrium
value that is determined by the final Bose-Einstein distribution, such that
there is an interplay of cooling and re-thermalization.

To compute the entropy, we first consider a bosonic system in statistical
equilibrium that will be the limit of the time-dependent case for
$t\rightarrow\infty$,

$\displaystyle S_{\text{eq}}=\int
g(\epsilon)\Bigl{[}\bigl{(}1+n_{\text{eq}}(\epsilon)\bigr{)}\ln\bigl{(}1+n_{\text{eq}}(\epsilon)\bigr{)}$
(26) $\displaystyle-n_{\text{eq}}(\epsilon)\ln
n_{\text{eq}}(\epsilon)\Bigr{]}\,\text{d}\epsilon\,.$

Here the density of states $g(\epsilon)$ for a three-dimensional isotropic
Bose gas without external potential is given by

$\displaystyle g(\epsilon)=g_{0}\,\sqrt{\epsilon}$ (27)

with

$\displaystyle g_{0}=(2m)^{3/2}\,V/(4\pi^{2})\,,$ (28)

as obtained from the substitution of a summation over the quantum numbers of
the associated states with an energy integration [30]. For a harmonic
oscillator potential the dependence is
$g(\epsilon)=g^{\text{HO}}_{0}\epsilon^{2}$.

In order to check the numerical integration in Eq. (26) against an exact
result, we derive $S_{\text{eq}}^{\mu\rightarrow 0}$ for a free three-
dimensional Bose gas as

$\displaystyle S_{\text{eq}}^{\mu\rightarrow
0}/g_{0}=T^{3/2}\Biggl{(}\zeta\Bigl{(}\frac{5}{2}\Bigr{)}\Gamma\Bigl{(}\frac{5}{2}\Bigr{)}\qquad\qquad\qquad$
(29)
$\displaystyle+\frac{\sqrt{\pi}}{2}\sum_{k=1}^{\infty}\frac{1}{k}\Bigl{[}\zeta(3/2,k)-\zeta(3/2,k+1)\Bigr{]}\Biggr{)}$

with the Hurwitz zeta function $\zeta(s,q)$, yielding
$S_{\text{eq}}^{\mu\rightarrow 0}/(g_{0}\,T^{3/2})=2.9721553$ as compared to
the result $2.97216$ of a numerical integration from zero to infinity.

If the final temperature $T_{\text{f}}$ is below the initial one as is the
case for cooling, the entropy drops below the initial equilibrium value and
then rises again towards the final value in the course of the equilibration
process. For $T_{\text{f}}$ below critical value for condensate formation
$T_{\text{c}}$, particles also occupy the condensed state, but the total
entropy still equals the entropy of the atoms in the thermal cloud. As has
been emphasized in Ref. [31], this is the case even though the entropy of the
particles in the ground state is nonzero, because the latter is cancelled by
the so-called correlation entropy due to the fixed number of particles
distributed among the quantum states.

It is therefore sufficient for a calculation of the total time-dependent
entropy $S(t)$ to consider only the thermal cloud. The analytical solutions
based on Eqs. (14) and (21) for constant chemical potential $\mu$ and boundary
conditions at the singularity $\epsilon=\mu$ are used. The entropy in a
bosonic system for an average number of particles $n(\epsilon,t)$ per single-
particle state can be reformulated [32, 22] as

$\displaystyle
S(t)=\int_{0}^{\infty}g(\epsilon)\Bigl{[}\ln\bigl{(}1+n(\epsilon,t)\bigr{)}\qquad\qquad$
$\displaystyle+n(\epsilon,t)\ln\bigl{(}1+1/n(\epsilon,t)\bigr{)}\Bigr{]}\text{d}\epsilon\,.$
(30)

Again, the spatial dimensionality and the external confinement enter the
present formulation only through the density of states. The properties of the
trapping potential and its effect on the density of states in evaporative
cooling have been previously discussed in Ref. [26].

Figure 4:  Time evolution of the entropy $S(t)/g_{0}$ in an equilibrating Bose
gas for $k=0-6$ evaporative cooling steps starting with $T_{\text{i}}=1000$ nK
and $T_{k}=1000\times(0.8)^{k}$ nK as calculated from the analytical solution
of the NBDE Eq. (13) with constant transport coefficients and constant
chemical potential $\mu=-8$ nK, solid curves. In each step, $S(t)$ drops
instantaneously due to cooling and then rises in the course of re-
thermalization.  Figure 5:  Time evolution of the entropy $S(t)/g_{0}$ in an
equilibrating Bose gas in the course of single-step evaporative cooling from
$T_{\text{i}}=240$ nK with a cut at $\epsilon_{\text{i}}=80$ nK to
$T_{\text{f}}=100$ nK for $\mu=-8$ nK as calculated from the analytical
solution of the NBDE Eq. (13), solid curve. $S(t=0)$ is the entropy following
evaporative cooling. The dashed curve is the wave entropy, the dotted curve
the particle entropy, and the dot-dashed curve the result from the linear
relaxation ansatz [22]. The dotted horizontal line indicates the equilibrium
value at $T_{\text{f}}=100$ nK calculated from Eq. (26).

Results for successive cooling as displayed in Fig. 4 show that the entropy
decreases in each cooling step and subsequently rises due to re-
thermalization. Correspondingly, for $T_{\text{f}}<T_{\text{c}}$ the number of
condensed particles rises due to a transfer from the nonlinear kinetic region
into the condensate [28, 29], as will be discussed in the next section. In the
example shown here with $T_{\text{i}}\simeq 1000$ nK and
$T_{k}=1000\times(0.8)^{k}$ nK,
$\epsilon_{\text{cut},k}=1000\times(0.7)^{k+1}$ nK, $k=0\ldots 6$. The entropy
of the initial thermal distribution – without the cut – is about four times
larger than the one of the final distribution.

The first term in the entropy Eq. (30) is referred to as wave entropy, it
yields the largest contribution when the single-particle state is occupied by
many particles, as in the IR. The second term is the particle entropy, which
is more relevant in case of low occupation $n(\epsilon,t)<1$, as in the UV
[32, 22].

Results for the time-dependent contributions to the entropy are shown in Fig.
5, calculated with the analytical solutions of the NBDE. Here, the equilibrium
value of the entropy at the initial temperature $T_{\text{i}}\simeq 240$ nK
and for $\mu=-8$ nK is $S_{\text{eq}}/g_{0}\,(T_{\text{i}})\simeq 10666$
(nK)3/2. With a sharp cutoff at $\epsilon_{\text{i}}=80$ nK to account for
evaporative cooling, the entropy of the initial (cooled) nonequilibrium
distribution is reduced to $S_{\text{i}}/g_{0}\,(T_{\text{i}})\simeq 1218$
(nK)3/2.

Thermalization during the time evolution then occurs through the analytical
solutions of the NBDE, and the result for the rising entropy $S(t)/g_{0}$ is
shown in the solid curve in Fig. 5. The new equilibrium value of the entropy
at the final temperature $T_{\text{f}}\simeq 100$ nK following evaporative
cooling and thermalization is $S_{\text{eq}}/g_{0}\,(T_{\text{f}})\simeq 2740$
(nK)3/2, dotted horizontal line.

The dashed curve is the wave entropy. It is most relevant for large occupation
numbers, which are present at all times in the IR region
$\epsilon<\epsilon_{\text{i}}$, and therefore, this contribution shows a
rather weak time dependence. The dotted curve is the particle entropy, which
is initially smaller than the wave entropy, because the occupation in the UV
beyond the cut is negligible at small times. In the course of thermalization,
however, it exceeds the wave entropy at $t\simeq 10$ ms for the parameter set
used in this calculation, and rises subsequently. Hence, the entropy at large
times is mostly determined by the slowly rising contribution of the thermal
tail.

In contrast, results for the corresponding linear relaxation-time
approximation [22] (dot-dashed curve in Fig. 5) reach the thermal limit that
is given by Eq. (26) much faster because an exponential buildup of the thermal
tail towards $n_{\text{eq}}$ is enforced. The entropy approaches the same
equilibrium value, but at shorter times than the nonlinear solution.

It is noted that the final state does not change if both transport
coefficients are scaled by the same amount because $T_{\text{f}}=-D/v$, but
the time scale varies since $\tau_{\text{eq}}=4D/(9v^{2})$. By measuring the
time it takes for a given system to thermalize we can thus determine the
absolute magnitude of the transport coefficients. To obtain the time scale
from first principles, a microscopic calculation of the transport coefficients
from a many-body theory is required, which is beyond the scope of our present
approach.

Up to now, the results for the time-dependent occupation-number distributions
and the entropy have been discussed for constant chemical potential. When
computing the condensate fraction, however, particle-number conservation must
be considered and $\mu$ becomes time-dependent, see the following section.

Figure 6:  Chemical potential $\mu(t)$ in the course of evaporative cooling
from $T_{\text{i}}=240$ nK to $T_{\text{f}}=100$ nK as calculated based on
particle-number conservation from the analytical solution of the NBDE Eq. (13)
with $k_{\text{max}}\equiv K=5,10,20,200$ in the series expansion of the exact
solution with $\epsilon_{\text{i}}=100$ nK and $\mu_{\text{i}}=-8$ nK.

##  V Time-dependent condensate fraction

We proceed to calculate the rate of increase of atoms in the condensate
$N_{\text{c}}(t)$ for conserved total particle number
$N=N_{\text{c}}(t)+N_{\text{th}}(t)$ and $T_{\text{f}}<T_{\text{c}}$ based on
the exact analytic NBDE-solutions Eq. (14) for any given $\mu$. Particle-
number conservation – which is a necessary condition for condensate formation
to occur – requires a time-dependent chemical potential, with
$\mu(t)\rightarrow 0$ corresponding to the condensed state.

To compute $\mu(t)$ in accordance with overall particle-number conservation,
we note that the particle number of the thermal cloud is given by

$N_{\text{th}}(t)=\int_{0}^{\infty}g(\epsilon)\,\tilde{n}(\epsilon,t)\,\text{d}\epsilon\;,$
(31)

where $\tilde{n}(\epsilon,t)$ is the distribution function with the time-
dependent value of $\mu\equiv\mu(t)$. For each time step, $\mu(t)$ is
calculated such that the particle number is conserved, $N=N_{\text{th}}$,
until $\mu$ reaches the value of zero and the largest possible value of
$N_{\text{th}}$ is attained. This instant marks the moment of condensation and
any future difference between the initial particle number and the thermal
cloud particle number is accounted for by the condensate

$N-N_{\text{th}}(t)\equiv N_{c}(t)\,,\quad\text{for }\mu=0\;.$ (32)

A typical time evolution of $\mu(t)$ is displayed in Fig. 6, for different
values of $k_{\text{max}}\equiv K$ to show the convergence. Here the initial
chemical potential is $\mu_{\text{i}}=-8$ nK, the final one
$\mu_{\text{f}}=0$. With increasing $k_{\text{max}}$ the accuracy increases
rapidly, precise results are obtained already for $k_{\text{max}}=20$. For the
parameters used in Fig. 6, condensation sets in at around $400$ ms. As is
obvious from Eq. (31), $\mu(t)$ and hence, the onset of condensation, also
depends on the density of states.

Due to the statistical properties of the bosonic system that are encoded in
the NBDE, no condensate forms if the final temperature remains above the
critical value. We have confirmed this in calculations for different final
temperatures $T_{\text{f}}$, with the condensate fraction
$N_{\text{c}}(t)/N=0\,\forall\,t\land T_{\text{f}}>T_{\text{c}}$. It is also
self-evident because $T_{\text{c}}$ is derived from the Bose-Einstein
equilibrium solution, which is the stationary limit of the NBDE for
$t\rightarrow\infty$.

Figure 7:  Condensate fraction $N_{\text{c}}(t)/N$ in an equilibrating Bose
gas of 23Na subsequent to fast evaporative cooling in a single step from
$T_{\text{i}}=876$ nK to $T_{\text{f}}=750$ nK as calculated from the
analytical solution of the NBDE Eq. (13) with $k_{\text{max}}\equiv
K=5,10,20,40$ in the series expansion of the exact solution, cutoff energy
$\epsilon_{\text{i}}=2190$ nK, $\mu_{\text{i}}=-8$ nK, and the density of
states for a free Bose gas. The transport coefficients are $D=3750$ (nK)2
ms${}^{-1},~{}\,v=-5$ nK ms-1. The MIT data for the condensate fraction
(crosses, no error bars) are from Ref. [8]. Figure 8:  Condensate fraction
$N_{\text{c}}(t)/N$ as in Fig. 7, following evaporative cooling (RF-sweep) in
a single step from $T_{\text{i}}=876$ nK to $T_{\text{f}}=660$ nK as
calculated from the analytical solution of the NBDE Eq. (13) with
$k_{\text{max}}\equiv K=5,10,20,40$, $\epsilon_{\text{i}}=1810$ nK,
$\mu_{\text{i}}=-18$ nK, and the density of states in a harmonic trap. Due to
the constant transport coefficients, the diffusion into the condensate starts
instantaneously. The transport coefficients are $D=3300$ (nK)2
ms${}^{-1},~{}\,v=-5$ nK ms-1. The MIT data for 23Na are from Ref. [8], with
an enlarged time scale.

In order to compare our results for $T_{\text{f}}<T_{\text{c}}$ to data we use
the MIT measurements from Ref. [8] with parameters from Ref. [7], where a
nonlinear Schrödinger equation coupled to a Boltzmann-like quantum collision
term has been solved numerically.

In particular, the final number of condensate atoms is taken from Ref. [8].
The initial number of atoms in the trap is unknown but in Ref. [7], $N_{0}$ =
$60\text{\times}{10}^{6}$ has been used in order to achieve agreement of a
numerical simulation and the data. We adopt their values of $T_{\text{i}}=876$
nK and $\epsilon_{\text{i}}/T_{\text{i}}$ = $2.5$. In contrast to Ref. [7] our
treatment also requires the specification of the final temperature which we
infer from energy conservation in a numerical calculation, yielding
$T_{\text{f}}=750$ nK for the free gas, and $T_{\text{f}}=660$ nK with the
density of states in a harmonic oscillator potential as in the MIT trap.

After solving the NBDE for $n(\epsilon,t)$ at fixed initial chemical potential
we calculate $\mu(t)$ at all times and finally compute
$N_{c}(t)/N_{\text{tot}}$ which is shown in Fig. 7 together with the
convergence of the analytic solution for different values of $k_{\text{max}}$
in the series expansion of the exact analytic solution. In this result we have
used the density of states for a free gas, to be compared with the one for a
harmonic potential in Fig. 8 (here the timescale is enlarged to emphasize the
small-time region).

In both cases, the time axis has been shifted to give the best possible fit.
The transport coefficients were adapted such that the observed time scale was
reproduced. With these adjustments, data and model result are found to agree.
The solutions are steeper for the free gas at small times compared to the
harmonic-oscillator density of states.

One feature both solutions fail to account for, however, is the measured
initial slow increase of the condensate fraction in the first 30 ms, which
cannot be reproduced in our model with constant transport coefficients. Should
the small-time gradual rise of the condensate persist in future measurements,
we can adapt the model using time-dependent transport coefficients, such that
the diffusion into the condensate builds up slowly with time: For constant
drift and diffusion coefficients, the diffusion into the condensate sets in
instantaneously with full strength and without a seed condensate.

It is recognized that the experiments deal with an inhomogeneous trapped Bose
gas, whereas the theory in the present form is homogeneous. In particular, we
consider the difference between a finite system with and without a trapping
potential only through the difference in the densities of states. This is
certainly a limitation of the model, it is not meant to present a full
description of the dynamics of condensate formation. But it yields physically
reasonable results on a purely nonequilibrium-statistical basis.

##  VI Conclusion

Exact solutions of the nonlinear boson diffusion equation which take account
of the singularity in the initial conditions at $\epsilon=\mu<0$ and the
necessary boundary conditions at the singularity have been explored in this
work. The analytical solutions for constant transport coefficients are
compared in detail with numerical solutions of the NBDE obtained with Matlab
and agreement is found. The exact solutions are used to calculate the time-
dependent entropy in a cooling Bose gas. Whereas cooling drastically reduces
the entropy, the subsequent re-thermalization causes a gradual increase
towards the equilibrium value, which coincides with the Bose-Einstein result
and is significantly below the initial value before cooling. Wave and particle
entropy are discussed. Together with particle-number conservation and a time-
dependent chemical potential, the condensate fraction in a trap has been
calculated as a function of time.

The solutions of the NBDE are thus shown to properly describe equilibration
processes that occur in quantum gases in the course of evaporative cooling,
and the associated condensate formation for $T_{\text{f}}<T_{\text{c}}$. They
are specifically applied to evaporative cooling and condensate formation in
23Na. Agreement with MIT data for the time-dependent condensate fraction is
found. The convergence properties of the infinite series expansion in the
analytic solution of the NBDE have been tested.

In this work, we have not aimed at a detailed treatment of the quantum-
mechanical properties of the condensate, which is coupled to the time-
dependent evolution of the thermal cloud. Instead, the characteristic
nonequilibrium-statistical features of the bosonic system and in particular,
the boson stimulation that is encoded in the NBDE, have been used together
with particle-number conservation to infer the time-dependent evolution of the
condensate fraction.

Further refinements of the model such as variable transport coefficients to
describe a gradual buildup of the condensate, and the coupling to the quantum
fluctuations or excitations of the condensate are conceivable, but are
unlikely to allow for analytic solutions. Extensions of the NBDE itself to
higher dimensions in order to better account for anisotropic systems should
also be tested.

More extensive comparisons of the results to new data from cold-atom
experiments with an improved signal-to-noise ratio and systematic as well as
statistical error bars would be most welcome. In particular, time-dependent
measurements for other bosonic alkali atoms such as rubidium or lithium with
substantially different timescales for condensate formation and equilibration
are very desirable. This could, in particular, offer the opportunity to learn
more about the initial BEC formation phase in the first 50 ms. Here, the
quantum Boltzmann-type approaches often use an artificial seed condensate.
Although this is needed to trigger condensate growth and achieve agreement
with the data, it seems difficult to develop a sound physical justification
for the initial condensate occupation. In contrast, in the NBDE approach no
seed is needed.

###### Acknowledgements.

We thank both referees for questions and comments, and Joseph Indekeu of
Leuven University for the editorial work.

## References

  * Anderson  _et al._ [1995] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Observation of Bose-Einstein condensation in a dilute atomic vapor, Science 269, 198 (1995).
  * Davis  _et al._ [1995a] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Bose-Einstein condensation in a gas of sodium atoms, Phys. Rev. Lett. 75, 3969 (1995a).
  * Bradley  _et al._ [1995] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Evidence of Bose-Einstein condensation in an atomic gas with attractive interactions, Phys. Rev. Lett. 75, 1687 (1995).
  * Bradley  _et al._ [1997] C. C. Bradley, C. A. Sackett, and R. G. Hulet, Bose-Einstein condensation of lithium: Observation of limited condensate number, Phys. Rev. Lett. 78, 985 (1997).
  * Gardiner  _et al._ [1997] C. W. Gardiner, P. Zoller, R. J. Ballagh, and M. J. Davis, Kinetics of Bose-Einstein condensation in a trap, Phys. Rev. Lett 79, 1793 (1997).
  * Gardiner and Zoller [1997] C. W. Gardiner and P. Zoller, Quantum kinetic theory: A quantum kinetic master equation for condensation of a weakly interacting Bose gas without a trapping potential, Phys. Rev. A 55, 2902 (1997).
  * Bijlsma  _et al._ [2000] M. J. Bijlsma, E. Zaremba, and H. T. C. Stoof, Condensate growth in trapped Bose gases, Phys. Rev. A 62, 063609 (2000).
  * Miesner  _et al._ [1998] H.-J. Miesner, D. M. Stamper-Kurn, M. R. Andrews, D. S. Durfee, S. Inouye, and W. Ketterle, Bosonic stimulation in the formation of a Bose-Einstein condensate, Science 279, 1005 (1998).
  * Eckern [1984] U. Eckern, Relaxation processes in a condensed Bose gas, J. Low Temp. Phys. 54, 333 (1984).
  * Kirkpatrick and Dorfman [1985] T. R. Kirkpatrick and J. R. Dorfman, Transport coefficients in a dilute but condensed Bose gas, J. Low Temp. Phys. 58, 301 (1985).
  * Stoof [1997] H. T. C. Stoof, Initial stages of Bose-Einstein condensation, Phys. Rev. Lett. 78, 768 (1997).
  * Stoof [1999] H. T. C. Stoof, Coherent versus incoherent dynamics during Bose-Einstein condensation in atomic gases, J. Low Temp. Phys. 114, 11 (1999).
  * Zaremba  _et al._ [1999] E. Zaremba, T. Nikuni, A. Griffin, and J. Low, Dynamics of trapped Bose gases at finite temperatures, J. Low Temp. Phys. 116, 277 (1999).
  * Snoke and Wolfe [1989] D. W. Snoke and J. P. Wolfe, Population dynamics of a Bose gas near saturation, Phys. Rev. B 39, 4030 (1989).
  * Kagan  _et al._ [1992] Y. M. Kagan, B. V. Svistunov, and G. V. Shlyapnikov, Kinetics of Bose condensation in an interacting Bose gas, Sov. Phys. JETP 74, 279 (1992).
  * Semikoz and Tkachev [1995] D. V. Semikoz and I. I. Tkachev, Kinetics of Bose condensation, Phys. Rev. Lett. 74, 3093 (1995).
  * Luiten  _et al._ [1996] O. J. Luiten, M. W. Reynolds, and J. T. M. Walraven, Kinetic theory of the evaporative cooling of a trapped gas, Phys. Rev. A 53, 381 (1996).
  * Holland  _et al._ [1997] M. Holland, J. Williams, and J. Cooper, Bose-Einstein condensation: Kinetic evolution obtained from simulated trajectories, Phys. Rev. A 55, 3670 (1997).
  * Jaksch  _et al._ [1997] D. Jaksch, C. W. Gardiner, and P. Zoller, Quantum kinetic theory. II. Simulation of the quantum Boltzmann master equation, Phys. Rev. A 56, 575 (1997).
  * Wolschin [2018a] G. Wolschin, Equilibration in finite Bose systems, Physica A 499, 1 (2018a).
  * Wolschin [2018b] G. Wolschin, An exactly solvable model for equilibration in bosonic systems, EPL 123, 20009 (2018b).
  * Wolschin [2020] G. Wolschin, Time-dependent entropy of a cooling Bose gas, EPL 129, 40006 (2020).
  * Nordheim [1928] L. W. Nordheim, On the kinetic method in the new statistics and application in the electron theory of conductivity, Proc. R. Soc. Lond. A 119, 689 (1928).
  * Note [1] The derivative-term of the diffusion coefficient has been modified as compared to Refs.+.1667em[20, 21, 22] in order to secure the correct stationary solution.
  * Rasch and Wolschin [2020] N. Rasch and G. Wolschin, Solving a nonlinear analytical model for bosonic equilibration, Physics Open 2, 100013 (2020).
  * Davis  _et al._ [1995b] K. B. Davis, M.-O. Mewes, and W. Ketterle, An analytical model for evaporative cooling of atoms, Appl. Phys. B 60, 155 (1995b).
  * Skeel and Berzins [1990] R. D. Skeel and M. Berzins, A method for the spatial discretization of parabolic equations in one space variable, SIAM J. Sci. Statist. Comput. 11, 1 (1990).
  * Svistunov [1991] B. V. Svistunov, Highly nonequilibrium Bose condensation in a weakly interacting gas, J. Mosc. Phys. Soc. 1, 373 (1991).
  * Kagan and Svistunov [1997] Y. Kagan and B. V. Svistunov, Evolution of correlation properties and appearance of broken symmetry in the process of Bose-Einstein condensation, Phys. Rev. Lett. 79, 3331 (1997).
  * Pitaevskii and Stringari [2003] L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation_ , International series of monographs on physics (Clarendon Press, Oxford, 2003).
  * Kim  _et al._ [2018] M. B. Kim, A. Svidzinsky, G. S. Agrawal, and M. O. Scully, Entropy of the Bose-Einstein-condensate ground state: Correlation versus ground-state entropy, Phys. Rev. A 97, 013605 (2018).
  * Yamamoto and Haus [1986] Y. Yamamoto and H. A. Haus, Preparation, measurement and information capacity of optical quantum states, Rev. Mod. Phys. 58, 1001 (1986).