$\displaystyle\sum_{g\geq 0}N^{-2g}\sum_{\alpha\in
S_{n}|\mathrm{df}(\alpha,(1,\ldots,n))\leq
2g}\kappa^{(2g-\mathrm{df}(\alpha,(1,\ldots,n)))}_{\alpha}(A_{N})\operatorname{tr}_{\alpha^{-1}(1,\ldots,n)}(B_{N}),$
(40)

  2. 2.

both tensors $\mathbb{E}\left[(V_{N}A_{N}V_{N}^{*})^{\otimes n}\right]$ and
$\mathbb{E}\left[(U_{N}B_{N}U_{N}^{*})^{\otimes n}\right]$, one gets that the
moment $\mathbb{E}[\tau_{A_{N},U_{N}B_{N}U_{N}^{*}}(P)]$ is equal to:

$\displaystyle\sum_{g\geq 0}N^{-2g}\sum_{\begin{subarray}{c}a,b\geq
0,\alpha,\beta\in S_{n}\\\
a+b+\mathrm{df}(\alpha,\alpha\circ\beta)+\mathrm{df}(\alpha\circ\beta,(1,\ldots,n))=2g\end{subarray}}\kappa_{\alpha}^{(a)}(A_{N})\kappa_{\beta}^{(b)}(B_{N}).$

Using Corollary 5.7 together with Equations (35) and (37) which allows us to
bound the higher order cumulants, we obtain the expansion stated in the
theorem and the formulae for $\mu_{A_{N}}\star_{g}\mu_{B_{N}}(P)$. At last,
note that when $g=0$, the first term in the expansion in Equation (40) is
simply

$\sum_{\alpha\prec(1,\ldots,n)}\kappa^{(0)}_{\alpha}(A_{N})\operatorname{tr}_{\alpha^{-1}(1,\ldots,n)}(B_{N}).$

Using Lemma 5.9, we can conclude that the first term is simply given by the
free convolution $\tau_{A_{N}}\star\tau_{B_{N}}(P)$. ∎

The bounded coefficients
$\tau_{\mathbf{A}_{N}}\star_{g}\tau_{\mathbf{B}_{N}}(P)$ do depend on $N$. In
concrete cases, we might be given sequences of matrices having distribution
with a $\tfrac{1}{N}$ expansion up to order $K\in\mathbb{N}$ in the following
sense:

$\tau_{\mathbf{A}_{N}}(P)=\sum_{k=0}^{K}\frac{1}{N^{k}}\tau_{\mathbf{A}}^{(k)}(P)+o\left(N^{-K}\right).$

From such expansions, it is possible to expand each term
$\tau_{\mathbf{A}_{N}}\star_{g}\tau_{\mathbf{B}_{N}}(P)$ as a series in
$\tfrac{1}{N}$, and Theorem 5.10 yields an expansion

$\mathbb{E}[\tau_{\mathbf{A}_{N},U_{N}\mathbf{B}_{N}U_{N}^{*}}(P)]=\sum_{k=0}^{K}\frac{1}{N^{k}}\tau^{(k)}(P)+o\left(N^{-K}\right)$

with coefficients $\tau^{(k)}(P)$ which do not depend on $N$. We do not
provide here the formulae for computing the coefficients $\tau^{(k)}(P)$ from
the expansions of $\mathbf{A}_{N}$ and $\mathbf{B}_{N}$, as it has already
been done independently in [Gab15c] and [BCGF+21]. More precisely, the
existence of this expansion is ensured by [Gab15c, Theorem 8.5] and the
computation is possible thanks to the notion of _freeness up to order $K$ of
fluctuations_. This expansion is also a consequence of [BCGF+21, Theorem 4.28]
and the computation is possible thanks to the notion of _$(\tfrac{K}{2},1)$
-freeness_. The difference of Theorem 5.10 from [Gab15c, Theorem 8.5] and
[BCGF+21, Theorem 4.28] is that we only assume the boundedness in n.c.
distribution of $(\mathbf{A}_{N})_{N\geq 1}$ and $(\mathbf{B}_{N})_{N\geq 1}$
instead of assuming the expansion of $\tau_{\mathbf{A}_{N}}$ and
$\tau_{\mathbf{B}_{N}}$ in $\tfrac{1}{N}$.

######  Remark 5.11.

Note that from the computations in the previous proof, the mixed moment
$\mathbb{E}[\tau_{A_{N},U_{N}B_{N}U_{N}^{*}}((A_{N}B_{N})^{n})]$ is equal to
$\sum_{\alpha\in
S_{n}}N^{-\mathrm{df}(\alpha,(1,\ldots,n))}\kappa^{N}_{\alpha}(A_{N})\operatorname{tr}_{\alpha^{-1}(1,\ldots,n)}(B_{N})$
or $\sum_{\alpha,\beta\in
S_{n}}N^{-\mathrm{df}(\alpha,\alpha\circ\beta)-\mathrm{df}(\alpha\circ\beta,(1,\ldots,n))}\kappa^{N}_{\alpha}(A_{N})\kappa^{N}_{\beta}(B_{N})$.
This shows how one can use the mixed matricial cumulants in order to compute
mixed moments (see Example 5.12).

####  5.2.3 Computations of the higher order matricial cumulants

We illustrate the notion of matricial cumulants by first showing how one can
compute them explicitly when $n=2$. As a corollary of this, we obtain the
higher order matricial cumulants. Then, we compute a mixed moment of $A$ and
$UBU^{*}$ using these matricial cumulants. At last, we provide the relation
between these higer order matricial cumulants and the Weingarten function.

######  Example 5.12.

Consider $A\in M_{N}(\mathbb{C})$ fixed. Using Proposition 5.5,

$\displaystyle\operatorname{tr}_{1_{2}}(A,A)$
$\displaystyle=\kappa_{1_{2}}^{N}(A,A)+N^{-2}\kappa_{(1\,2)}^{N}(A,A),$
$\displaystyle\operatorname{tr}_{(1\,2)}(A,A)$
$\displaystyle=\kappa_{1_{2}}^{N}(A,A)+\kappa_{(1\,2)}^{N}(A,A).$

Therefore,

$\displaystyle\kappa_{1_{2}}^{N}(A,A)$
$\displaystyle=\frac{1}{1-N^{-2}}\operatorname{tr}_{N}(A)^{2}-\frac{N^{-2}}{1-N^{-2}}tr_{N}(A^{2})$
$\displaystyle\kappa_{(1\,2)}^{N}(A,A)$
$\displaystyle=\frac{1}{1-N^{-2}}(\operatorname{tr}_{N}(A^{2})-\operatorname{tr}_{N}(A)^{2}).$

Using Corollary 5.7, it follows that

$\kappa_{1_{2}}^{(2g)}(A,A)=\left\\{\begin{array}[]{lc}\operatorname{tr}_{N}(A)^{2}&\text{
if }g=0,\\\ \operatorname{tr}_{N}(A)^{2}-\operatorname{tr}_{N}(A^{2})&\text{
if }g\geq 1,\end{array}\right.$

whereas

$\kappa_{(1\,2)}^{(2g)}(A,A)=\operatorname{tr}_{N}(A^{2})-\operatorname{tr}_{N}(A)^{2},\hskip
14.22636pt\forall g\geq 0.$

Consider now $A,B\in M_{N}(\mathbb{C})$ and $\tilde{B}=UBU^{*}$ where $U$ is a
Haar unitary matrix. Let us compute
$\mathbb{E}[\operatorname{tr}(A\tilde{B}A\tilde{B})]$ using the matricial
cumulants (see also Remark 5.11):

$\displaystyle\mathbb{E}[\operatorname{tr}(A\tilde{B}A\tilde{B})]$
$\displaystyle=\frac{1}{N}\operatorname{Tr}_{2,N}\left[\left(\kappa_{1_{2}}^{N}(A,A)1_{2}+N^{-1}\kappa_{(1\,2)}^{N}(A,A)(1\,2)\right)B^{\otimes
2}(1,2)^{t}\right]$
$\displaystyle=\kappa_{1_{2}}^{N}(A,A)\operatorname{tr}_{N}(B^{2})+\kappa_{(1,2)}^{N}(A,A)\operatorname{tr}_{N}(B)^{2}.$

Note that it can also be expressed as:

$\displaystyle\kappa_{1_{2}}(A,A)\operatorname{tr}_{N}(B^{2})+\kappa_{(1,2)}(A,A)\operatorname{tr}_{N}(B)^{2}-\frac{N^{-2}}{1-N^{-2}}\kappa_{(1\,2)}(A,A)\kappa_{(1\,2)}(B,B).$

In Lemma 5.9, we identified the $0$-order matricial cumulants as the free
cumulants. This identification yields a way to compute explicitly them. All
higher order cumulants $\kappa^{(g)}$ can also be expressed thanks to Lemma
5.8. Let us give an example with $g=1.$

######  Lemma 5.13.

For any $\sigma\in S_{n},$ any n-tuple of matrices
$\mathbf{M}=(M_{1},\ldots,M_{n})\in M_{N}(\mathbb{C})$,

$\displaystyle\kappa^{(2)}_{\sigma}(\mathbf{M})$
$\displaystyle=-\sum_{\begin{subarray}{c}\alpha,\beta\in S_{n}\\\
\mathrm{df}(\beta,\alpha)=2,\alpha\prec\sigma\end{subarray}}\mathrm{Moeb}(\alpha^{-1}\sigma)\kappa_{\beta}(\mathbf{M})$
$\displaystyle=-\sum_{\begin{subarray}{c}\alpha,\beta,\gamma\in S_{n}\\\
\gamma\prec\beta,\mathrm{df}(\beta,\alpha)=2,\alpha\prec\sigma\end{subarray}}\mathrm{Moeb}(\alpha^{-1}\sigma)\mathrm{Moeb}(\gamma^{-1}\beta)\operatorname{tr}_{\gamma}(\mathbf{M}).$

###### Proof.

The function $\mathrm{Moeb}$ is related to the Moebius function for the order
$\prec$. More precisely, two families $(a_{\sigma})_{\sigma\in S_{n}}$ and
$(b_{\sigma})_{\sigma\in S_{n}}$ satisfy $a_{\sigma}=\sum_{\alpha\in
S_{n},\alpha\prec\sigma}b_{\alpha}$ for any $\sigma\in S_{n}$, if and only if

$\displaystyle b_{\sigma}=\sum_{\alpha\in
S_{n},\alpha\prec\sigma}\mathrm{Moeb}(\alpha^{-1}\sigma)a_{\alpha}$

holds for any $\sigma\in S_{n}$. From Lemma 5.8, we know that:

$\sum_{\alpha\prec\sigma}\kappa_{\alpha}(\mathbf{M})=\operatorname{tr}_{\sigma}(\mathbf{M}),\quad\quad\sum_{\alpha\prec\sigma}\kappa_{\alpha}^{(2)}(\mathbf{M})=\zeta_{\sigma},$

where $\zeta_{\sigma}=\sum_{\alpha\in
S_{n},\alpha\neq\sigma\atop\mathrm{df}(\alpha,\sigma)=2}\kappa_{\alpha}(\mathbf{M})$.
This yields the two formulas for $\kappa_{\sigma}^{(2)}(\mathbf{M})$. ∎

Actually, matricial cumulants relate to the Weingarten function. Recall that,
in Section 5.1, we saw that the function
$\mathrm{Moeb}_{N}(\cdot):=N^{n+|\cdot|}\mathrm{Wg}_{N}(\cdot)$ can be
expanded in power of $N^{-2}$ (see [CŚ04]): for any $\sigma\in S_{n}$,
$\mathrm{Moeb}_{N}(\sigma)=\sum_{g\geq 0}\mathrm{Moeb}^{(g)}(\sigma)N^{-2g}.$
The matricial cumulants can then be expressed using $\mathrm{Moeb}_{N}$ and
$\mathrm{Moeb}^{(g)}$.

######  Lemma 5.14.

For any $\sigma\in S_{n},$ any n-tuple of matrices
$\mathbf{M}=(M_{1},\ldots,M_{n})\in M_{N}(\mathbb{C})$,

$\kappa_{\sigma}^{N}(\mathbf{M})=\sum_{\alpha\in
S_{n}}N^{-\mathrm{df}(\alpha,\sigma)}\mathrm{Moeb}_{N}(\alpha)\operatorname{tr}_{\alpha^{-1}\sigma}(\mathbf{M}).$

Hence, for any $g\geq 0$,

$\kappa_{\sigma}^{(g)}(\mathbf{M})=\sum_{\alpha\in
S_{n}:\mathrm{df}(\alpha,\sigma)\leq
2g}\mathrm{Moeb}^{(g-\mathrm{df}(\alpha,\sigma)/2)}(\alpha)\operatorname{tr}_{\alpha^{-1}\sigma}(\mathbf{M}).$
(41)

###### Proof.

Equation (23) can be expressed as:

$\displaystyle\mathbb{E}\left[U^{\otimes n}\otimes(U^{*})^{\otimes n}\
\right]=\sum_{\alpha,\beta\in
S_{n}}\mathrm{Wg}_{N}(\alpha\beta)\tau\circ(\alpha\otimes\beta),$

where $\tau(i)=n+i$ mod. $2n$. This implies that:

$\displaystyle\mathbb{E}\left[UM_{1}U^{*}\otimes\ldots\otimes
UM_{n}U^{*}\right]=\sum_{\alpha,\beta\in
S_{n}}\mathrm{Wg}_{N}(\alpha\beta)Tr_{n,N}(M_{1}\otimes\ldots\otimes
M_{n}\circ\beta)\alpha.$

Using the definition of the matricial cumulants,

$\displaystyle\kappa_{\sigma}^{N}(\mathbf{M})$
$\displaystyle=N^{|\sigma|}\sum_{\beta\in
S_{n}}\mathrm{Wg}_{N}(\sigma\beta)Tr_{n,N}(M_{1}\otimes\ldots\otimes
M_{n}\circ\beta),$ $\displaystyle=\sum_{\alpha\in
S_{n}}N^{-\mathrm{df}(\alpha,\sigma)}\mathrm{Moeb}_{N}(\alpha)\operatorname{tr}_{\alpha^{-1}\sigma}(\mathbf{M})$

which allows us to conclude for the first equation. The second equation is
obtained by using the expansion of $\mathrm{Moeb}_{N}$ and identifying the
coefficients in the expansion. ∎

##  6 The orthogonal case

In this section, we indicate which results remain true and which ones do not
when Haar distributed unitary random matrices are replaced by Haar distributed
orthogonal random matrices. Let us denote by $O(N)$ the orthogonal group of
real $N\times N$ matrices $U_{N}$ such that $U_{N}U_{N}^{t}=I_{N}$. We
consider also the subgroup $SO(N)=\\{U_{N}\in O(N):\det(U_{N})=1\\}$.

The first proof of Theorem 3.3 can be adapted to the orthogonal case. It has
been done by Collins and Śniady in [CŚ04] (see also [Lév17, Chapter 3]). The
main difference is that the Weingarten function is converging up to order
$O(N^{-1})$. It yields to the following asymptotic freeness.

######  Theorem 6.1.

Let $(\mathbf{A}_{N})_{N\geq 1}$ and $(\mathbf{B}_{N})_{N\geq 1}$ be two
deterministic sequences of $k$-tuples and $\ell$-tuples of real symmetric
$N\times N$ matrices. Let $(U_{N})_{N\geq 1}$ be a sequence of orthogonal
random matrices of size $N$, uniformly distributed in $SO(N)$, or $O(N)$.

Then, for all $P\in\mathbb{C}\langle
X_{1},\ldots,X_{k},Y_{1},\ldots,Y_{\ell}\rangle$,

$\mathbb{E}[\tau_{\mathbf{A}_{N},U_{N}\mathbf{B}_{N}U_{N}^{*}}(P)]=\tau_{\mathbf{A}_{N}}\star\tau_{\mathbf{B}_{N}}(P)+N^{-1}R(\mathbf{A}_{N},\mathbf{B}_{N}),$

with
$|R(\mathbf{A}_{N},\mathbf{B}_{N})|=O(\sum_{a+b\leq\mathrm{deg}(P)}\|\mathbf{A}_{N}\|_{*,a}\|\mathbf{B}_{N}\|_{*,b})$.
In other words, $\mathbb{C}\langle X_{1},\ldots,X_{k}\rangle$ and
$\mathbb{C}\langle Y_{1},\ldots,Y_{\ell}\rangle$ are free up to order
$O(N^{-1})$ with respect to
$\mathbb{E}[\tau_{\mathbf{A}_{N},U_{N}\mathbf{B}_{N}U_{N}^{*}}]$ whenever
$(\mathbf{A}_{N})_{N\geq 1}$ and $(\mathbf{B}_{N})_{N\geq 1}$ are bounded in
n.c. distribution.

Note that the rate of convergence is not sufficient here to conclude about
infinitesimal freeness. In fact, the invariance in law by orthogonal
conjugation does not imply infinitesimal freeness, as shown by Mingo in
[Min18]. As well, the rate of convergence seems not to be sufficient to
conclude about freeness of type $B$. However, note that if
$\|\mathbf{A}_{N}\|_{*,a}\|\mathbf{B}_{N}\|_{*,b}$ is $O(N^{-1})$, we have
$N^{-1}R(\mathbf{A}_{N},\mathbf{B}_{N})=O(N^{-2})$. This simple observation
allows to get the asymptotic freeness of type $B$ as follows.

######  Theorem 6.2.

In Theorem 3.5, the matrices $(U_{N})_{N\geq 1}$ can be taken uniformly
distributed in $O(N)$ or in $SO(N)$ whenever $(\mathbf{A}_{N})_{N\geq 1}$ and
$(\mathbf{B}_{N})_{N\geq 1}$ are real symmetric matrices.

###### Proof.

The proof follows the one of Theorem 6.2. It is not possible to apply directly
Proposition 2.9 here because the freeness is not of order $O(N^{-2})$.
However, we have the rate of convergence $O(N^{-2})$ on $\mathcal{I}$, which
yields the convergence of
$\mathbb{E}[\tau_{\mathbf{A}_{N},U_{N}\mathbf{B}_{N}U_{N}^{*}}]$ to $\tau$ and
$\tau^{\prime}$, and the freeness of type $B$ as explained in Remark 2.10.

The almost sure convergence is a consequence of the concentration phenomenon.
The analogue of Theorem 3.6 is true for the special orthogonal group $SO(N)$
(see [AGZ10, Corollary 4.4.28]). Unfortunately, this concentration fail to
apply to $O(N)$ because it is not connected. Following [MM13b], we condition
on the two possible values of $\det(U_{N})$. For this purpose, it is useful to
consider the uniform measure on the homogeneous space $SO^{-}(N)=\\{U_{N}\in
O(N):\det(U_{N})=-1\\}$ (it is the invariant measure under the action of
$SO(N)$). Theorem 6.1 and Theorem 6.2 remain true if $(U_{N})_{N\geq 1}$ is
uniformly distributed in $SO^{-}(N)$ because the concentration phenomenon
holds for $SO^{-}(N)$ (see Proposition 2.2 of [MM13b]). The almost sure result
for $O(N)$ is obtained by conditioning on the two possible values of the
determinant, and using the result for $SO(N)$ and $SO^{-}(N)$. ∎

As for the unitary case, we obtain (cyclic) monotone independence and
conditional independence as consequences of Theorem 6.2.

######  Theorem 6.3.

In Theorem 3.8, Theorem 3.10 and Corollary 3.11, the matrices $(U_{N})_{N\geq
1}$ can be taken uniformly distributed in $O(N)$ or in $SO(N)$ whenever
$(\mathbf{A}_{N})_{N\geq 1}$ and $(\mathbf{B}_{N})_{N\geq 1}$ are real
symmetric matrices, and $v_{N}\in\mathbb{R}^{N}$.

## References

  * [AB14] Romain Allez and Jean-Philippe Bouchaud.  Eigenvector dynamics under free addition.  Random Matrices: Theory and Applications, 3(03):1450010, 2014.
  * [AC21] Octavio Arizmendi and Adrian Celestino.  Polynomial with cyclic monotone elements with applications to random matrices with discrete spectrum.  Random Matrices: Theory and Applications, 10(02):2150020, 2021.
  * [AGZ10] Greg W Anderson, Alice Guionnet, and Ofer Zeitouni.  An introduction to random matrices.  Number 118. Cambridge university press, 2010.
  * [AHL22] Octavio Arizmendi, Takahiro Hasebe, and Franz Lehner.  Cyclic independence: Boolean and monotone.  arXiv preprint arXiv:2204.00072, 2022.
  * [Au21] Benson Au.  Finite-rank perturbations of random band matrices via infinitesimal free probability.  Communications on Pure and Applied Mathematics, 74(9):1855–1895, 2021.
  * [BABP16] Joël Bun, Romain Allez, Jean-Philippe Bouchaud, and Marc Potters.  Rotational invariant estimator for general noisy matrices.  IEEE Transactions on Information Theory, 62(12):7475–7490, 2016\.
  * [BBAP05] Jinho Baik, Gérard Ben Arous, and Sandrine Péché.  Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices.  Ann. Probab., 33(5):1643–1697, 2005.
  * [BBCF17] Serban T. Belinschi, Hari Bercovici, Mireille Capitaine, and Maxime Février.  Outliers in the spectrum of large deformed unitarily invariant models.  Ann. Probab., 45(6A):3571–3625, 11 2017.
  * [BBP17] Joël Bun, Jean-Philippe Bouchaud, and Marc Potters.  Cleaning large correlation matrices: tools from random matrix theory.  Physics Reports, 666:1–109, 2017.
  * [BBP18] Joël Bun, Jean-Philippe Bouchaud, and Marc Potters.  Overlaps between eigenvectors of correlated random matrices.  Physical Review E, 98(5):052145, 2018.
  * [BCGF+21] Gaëtan Borot, Séverin Charbonnier, Elba Garcia-Failde, Felix Leid, and Sergey Shadrin.  Analytic theory of higher order free cumulants.  arXiv preprint arXiv:2112.12184, 2021.
  * [Bel08] Serban Belinschi.  C-free convolution of measures with unbounded support.  von Neumann algebras in Sibiu, Theta Ser. Adv. Math., (10):1–7, 2008.
  * [BES16] Zhigang Bao, László Erdős, and Kevin Schnelli.  Local stability of the free additive convolution.  Journal of Functional Analysis, 271(3):672–719, 2016.
  * [BES17a] Zhigang Bao, László Erdős, and Kevin Schnelli.  Convergence rate for spectral distribution of addition of random matrices.  Advances in Mathematics, 319:251–291, 2017.
  * [BES17b] Zhigang Bao, László Erdős, and Kevin Schnelli.  Local law of addition of random matrices on optimal scale.  Communications in Mathematical Physics, 349(3):947–990, 2017.
  * [BES20] Zhigang Bao, László Erdős, and Kevin Schnelli.  Spectral rigidity for addition of random matrices at the regular edge.  Journal of Functional Analysis, 279(7):108639, 2020.
  * [BG22] Florent Benaych-Georges.  A short proof of Ledoit-Péché’s RIE formula for covariance matrices.  arXiv preprint arXiv:2201.05690, 2022.
  * [BGEM19] Florent Benaych-Georges, Nathanaël Enriquez, and Alkéos Michaïl.  Empirical spectral distribution of a matrix under perturbation.  Journal of Theoretical Probability, 32(3):1220–1251, 2019.
  * [BGGM11] Florent Benaych-Georges, Alice Guionnet, and Mylène Maida.  Fluctuations of the Extreme Eigenvalues of Finite Rank Deformations of Random Matrices.  Electronic Journal of Probability, 16(none):1621 – 1662, 2011.
  * [BGN03] Philippe Biane, Frederick Goodman, and Alexandru Nica.  Non-crossing cumulants of type B.  Transactions of the American Mathematical Society, 355(6):2263–2303, 2003.
  * [BGN11] Florent Benaych-Georges and Raj Rao Nadakuditi.  The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices.  Adv. Math., 227(1):494–521, 2011.
  * [BGN12] Florent Benaych-Georges and Raj Rao Nadakuditi.  The singular values and vectors of low rank perturbations of large rectangular random matrices.  J. Multivariate Anal., 111:120–135, 2012.
  * [Bia97] Philippe Biane.  Some properties of crossings and partitions.  Discrete Mathematics, 175(1):41–53, 1997.
  * [Bia98a] Philippe Biane.  Processes with free increments.  Mathematische Zeitschrift, 227(1):143–174, 1998.
  * [Bia98b] Philippe Biane.  Representations of Symmetric Groups and Free Probability.  Advances in Mathematics, 138(1):126–181, 1998.
  * [BJ21] Zdzislaw Burda and Andrzej Jarosz.  Cleaning large-dimensional covariance matrices for correlated samples.  arXiv preprint arXiv:2107.01352, 2021.
  * [BK18] Joël Bun and A Knowles.  An optimal rotational invariant estimator for general covariance matrices: The outliers.  Preprint, 2018.
  * [BLS96] Marek Bożejko, Michael Leinert, and Roland Speicher.  Convolution and limit theorems for conditionally free random variables.  Pacific Journal of Mathematics, 175(2):357–388, 1996.
  * [BMP07] ZD Bai, BQ Miao, and GM Pan.  On asymptotics of eigenvectors of large sample covariance matrix.  The Annals of Probability, 35(4):1532–1572, 2007.
  * [BS12] Serban T. Belinschi and Dimitri Shlyakhtenko.  Free probability of type B: analytic interpretation and applications.  American Journal of Mathematics, 134(1):193–234, 2012.
  * [BV93] Hari Bercovici and Dan Voiculescu.  Free convolution of measures with unbounded support.  Indiana University Mathematics Journal, 42(3):733–773, 1993.
  * [BV98] Hari Bercovici and Dan Voiculescu.  Regularity questions for free convolution.  In Nonselfadjoint operator algebras, operator theory, and related topics, pages 37–47. Springer, 1998.
  * [BW08] Hari Bercovici and Jiun-Chau Wang.  On freely indecomposable measures.  Indiana University Mathematics Journal, pages 2601–2610, 2008.
  * [CC06] M. Capitaine and M. Casalis.  Cumulants for random matrices as convolutions on the symmetric group.  Probability Theory and Related Fields, 136(1):19–36, 2006.
  * [CC08] M. Capitaine and M. Casalis.  Geometric interpretation of the cumulants for random matrices previously defined as convolutions on the symmetric group, pages 93–119.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.
  * [CDM18] M. Capitaine and C. Donati-Martin.  Spectrum of deformed random matrices and free probability, volume 53.  Panoramas et synthèses, 2018.
  * [CHS18] Benoît Collins, Takahiro Hasebe, and Noriyoshi Sakuma.  Free probability for purely discrete eigenvalues of random matrices.  Journal of the Mathematical Society of Japan, 70(3):1111–1150, 2018\.
  * [CLS22] Benoît Collins, Felix Leid, and Noriyoshi Sakuma.  Matrix models for cyclic monotone and monotone independences.  arXiv preprint arXiv:2202.11666, 2022.
  * [CMN21] Benoit Collins, Sho Matsumoto, and Jonathan Novak.  The Weingarten Calculus.  arXiv e-prints, page arXiv:2109.14890, September 2021.
  * [Col03] Benoît Collins.  Moments and cumulants of polynomial random variables on unitary groups, the Itzykson-Zuber integral, and free probability.  Int. Math. Res. Not., (17):953–982, 2003.
  * [CŚ04] Benoît Collins and Piotr Śniady.  Integration with respect to the Haar measure on unitary, orthogonal and symplectic group.  Communications in Mathematical Physics, 264(3):773–795, 2004.
  * [CS11] Stephen Curran and Roland Speicher.  Asymptotic Infinitesimal Freeness with Amalgamation for Haar Quantum Unitary Random Matrices.  Communications in Mathematical Physics, 301(3):627–659, 2011.
  * [DF19] Sandrine Dallaporta and Maxime Fevrier.  Fluctuations Of Linear Spectral Statistics Of Deformed Wigner Matrices.  mar 2019.
  * [DJ20] Xiucai Ding and Hong Chang Ji.  Local laws for multiplication of random matrices and spiked invariant model.  arXiv preprint arXiv:2010.16083, 2020.
  * [dS18] Ricardo Correa da Silva.  Lecture notes on noncommutative $L_{p}$-spaces.  arXiv preprint arXiv:1803.02390, 2018.
  * [FK81] Zoltán Füredi and János Komlós.  The eigenvalues of random symmetric matrices.  Combinatorica, 1(3):233–241, 1981.
  * [FMNS19] Maxime Février, Mitja Mastnak, Alexandru Nica, and Kamil Szpojankowski.  A construction which relates c-freeness to infinitesimal freeness.  Advances in Applied Mathematics, 110:299–341, 2019.
  * [FN10] Maxime Février and Alexandru Nica.  Infinitesimal non-crossing cumulants and free probability of type B.  Journal of Functional Analysis, 258(9):2983–3023, 2010.
  * [Fra06] Uwe Franz.  Multiplicative monotone convolutions, in ”quantum probability”.  Banach Center Publications, 73:153–166, 2006.
  * [Gab15a] Franck Gabriel.  A combinatorial theory of random matrices III: random walks on $\mathfrak{S}(N)$, ramified coverings and the $\mathfrak{S}(\infty)$ Yang-Mills measure, 2015.
  * [Gab15b] Franck Gabriel.  Combinatorial theory of permutation-invariant random matrices I: partitions, geometry and renormalization, 2015.
  * [Gab15c] Franck Gabriel.  Combinatorial theory of permutation-invariant random matrices II: cumulants, freeness and Lévy processes, 2015.
  * [GNR16] Fabrice Gamboa, Jan Nagel, and Alain Rouault.  Sum rules via large deviations.  Journal of Functional Analysis, 270(2):509–559, 2016.
  * [GNR19] Fabrice Gamboa, Jan Nagel, and Alain Rouault.  Sum rules and large deviations for spectral matrix measures.  Bernoulli, 25(1):712–741, 2019.
  * [GNR22] Fabrice Gamboa, Jan Nagel, and Alain Rouault.  Sum rules via large deviations: extension to polynomial potentials and the multi-cut regime.  Journal of Functional Analysis, 282(3):109307, 2022.
  * [GR11] Fabrice Gamboa and Alain Rouault.  Large deviations for random spectral measures and sum rules.  Applied Mathematics Research eXpress, 2011(2):281–307, 2011.
  * [Has10] Takahiro Hasebe.  New associative product of three states generalizing free, monotone, anti-monotone, Boolean, conditionally free and conditionally monotone products.  arXiv preprint arXiv:1009.1505, 2010.
  * [Ho22] Ching-Wei Ho.  A local limit theorem and delocalization of eigenvectors for polynomials in two matrices.  International Mathematics Research Notices, 2022(3):1734–1769, 2022\.
  * [HP00] Fumio Hiai and Dénes Petz.  The semicircle law, free random variables and entropy.  Number 77. American Mathematical Soc., 2000.
  * [Kar15] Vladislav Kargin.  Subordination for the sum of two random matrices.  The Annals of Probability, 43(4):2119–2150, 2015.
  * [KY14] Antti Knowles and Jun Yin.  The outliers of a deformed Wigner matrix.  The Annals of Probability, 42(5):1980–2031, 2014.
  * [Lév17] Thierry Lévy.  The master field on the plane.  Astérisque, 388, dec 2017.
  * [LP11] Olivier Ledoit and Sandrine Péché.  Eigenvectors of some large sample covariance matrix ensembles.  Probability Theory and Related Fields, 151(1):233–264, 2011.
  * [Min18] James A Mingo.  Non-crossing annular pairings and the infinitesimal distribution of the goe.  Journal of the London Mathematical Society, 2018.
  * [MM13a] Elizabeth Meckes and Mark Meckes.  Spectral measures of powers of random matrices.  Electronic communications in probability, 18:1–13, 2013.
  * [MM13b] Elizabeth S Meckes and Mark W Meckes.  Concentration and convergence rates for spectral measures of random matrices.  Probability Theory and Related Fields, 156(1):145–164, 2013.
  * [Mur00] Naofumi Muraki.  Monotonic convolution and monotonic Lévy-Hincin formula.  preprint, 2000.
  * [Mur01] Naofumi Muraki.  Monotonic independence, monotonic central limit theorem and monotonic law of small numbers.  Infinite Dimensional Analysis, Quantum Probability and Related Topics, 4(01):39–58, 2001.
  * [Noi21] Nathan Noiry.  Spectral measures of spiked random matrices.  Journal of Theoretical Probability, 34(2):923–952, 2021.
  * [NR20] Nathan Noiry and Alain Rouault.  Large deviations for spectral measures of some spiked matrices.  arXiv preprint arXiv:2009.12541, 2020.
  * [NS06] Alexandru Nica and Roland Speicher.  Lectures on the Combinatorics of Free Probability, volume 335 of London Mathematical Society Lecture Note Series.  Cambridge University Press, 2006.
  * [PB20] Marc Potters and Jean-Philippe Bouchaud.  A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists.  Cambridge University Press, 2020.
  * [Péc06] S. Péché.  The largest eigenvalue of small rank perturbations of Hermitian random matrices.  Probab. Theory Related Fields, 134(1):127–173, 2006.
  * [PW11] Mihai Popa and Jiun-Chau Wang.  On multiplicative conditionally free convolution.  Transactions of the American Mathematical Society, 363(12):6309–6335, 2011.
  * [Sam80] Stuart Samuel.  $U(N)$ Integrals, $1/N$, and the De Wit–’t Hooft anomalies.  Journal of Mathematical Physics, 21(12):2695–2703, dec 1980.
  * [Shl18] D. Shlyakhtenko.  Free probability of type B and asymptotics of finite-rank perturbations of random matrices.  Indiana University Mathematics Journal, 67(7), 2018.
  * [Spe93] Roland Speicher.  Free convolution and the random sum of matrices.  Publications of the Research Institute for Mathematical Sciences, 29(5):731–744, 1993.
  * [Tho00] S Thorbjørnsen.  Mixed Moments of Voiculescu’s Gaussian Random Matrices.  Journal of Functional Analysis, 176(2):213–246, 2000.
  * [Voi85] Dan Voiculescu.  Symmetries of some reduced free product $C^{*}$-algebras.  In Operator algebras and their connections with topology and ergodic theory, pages 556–588. Springer, 1985.
  * [Voi91a] Dan Voiculescu.  Limit laws for random matrices and free products.  Invent. Math., 104(1):201–220, 1991.
  * [Voi91b] Dan-Virgil Voiculescu.  Limit laws for random matrices and free products.  Invent. Math., 104(1):201–220, 1991.
  * [XB19] Ningning Xia and Zhidong Bai.  Convergence rate of eigenvector empirical spectral distribution of large wigner matrices.  Statistical Papers, 60(3):983–1015, 2019.
  * [XQB13] Ningning Xia, Yingli Qin, and Zhidong Bai.  Convergence rates of eigenvector empirical spectral distribution of large dimensional sample covariance matrix.  The Annals of Statistics, 41(5):2572–2607, 2013.
  * [Xu97] Feng Xu.  A random matrix model from two-dimensional Yang-Mills theory.  Comm. Math. Phys., 190(2):287–307, 1997.
  * [XYY20] Haokai Xi, Fan Yang, and Jun Yin.  Convergence of eigenvector empirical spectral distribution of sample covariance matrices.  The Annals of Statistics, 48(2):953–982, 2020.


11institutetext: Beihang University, China 22institutetext: College of
Computing and Data Science, Nanyang Technological University, Singapore
33institutetext: Lee Kong Chian School of Medicine, Nanyang Technological
University, Singapore 44institutetext: Sun Yat-sen University, China
55institutetext: National University of Singapore, Singapore 66institutetext:
South China University of Technology, China

# ETSCL: An Evidence Theory-Based Supervised Contrastive Learning Framework
for Multi-modal Glaucoma Grading

Zhiyuan Yang Zhiyuan Yang and Bo Zhang contributed equally.11
0009-0007-8277-5772 Bo Zhang0 22 0009-0002-3736-4070 Yufei Shi 33
0000-0003-1822-1621 Ningze Zhong 44 0009-0007-6758-058X Johnathan Loh 55
0000-0001-8176-0498 Huihui Fang 22 0000-0003-3380-7970 Yanwu Xu Yanwu Xu and
Si Yong Yeo are the corresponding authors. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>0000-0002-1779-931X Si Yong Yeo0 33
0000-0001-6403-6019

###### Abstract

Glaucoma is one of the leading causes of vision impairment. Digital imaging
techniques, such as color fundus photography (CFP) and optical coherence
tomography (OCT), provide quantitative and noninvasive methods for glaucoma
diagnosis. Recently, in the field of computer-aided glaucoma diagnosis, multi-
modality methods that integrate the CFP and OCT modalities have achieved
greater diagnostic accuracy compared to single-modality methods. However, it
remains challenging to extract reliable features due to the high similarity of
medical images and the unbalanced multi-modal data distribution. Moreover,
existing methods overlook the uncertainty estimation of different modalities,
leading to unreliable predictions. To address these challenges, we propose a
novel framework, namely ETSCL, which consists of a contrastive feature
extraction stage and a decision-level fusion stage. Specifically, the
supervised contrastive loss is employed to enhance the discriminative power in
the feature extraction process, resulting in more effective features. In
addition, we utilize the Frangi vesselness algorithm as a preprocessing step
to incorporate vessel information to assist in the prediction. In the
decision-level fusion stage, an evidence theory-based multi-modality
classifier is employed to combine multi-source information with uncertainty
estimation. Extensive experiments demonstrate that our method achieves state-
of-the-art performance. The code is available at https://github.com/master-
Shix/ETSCL.

###### Keywords:

Glaucoma grading Multi-modal learning Contrastive learning Evidence theory

##  1 Introduction

Glaucoma, a leading cause of vision impairment, causes progressive damage to
the optic nerve, including the optic nerve head (ONH), the retinal nerve fibre
layer (RNFL) and the ganglion cell-inner plexiform layer (GCIPL) [1]. Digital
imaging techniques, such as color fundus photography (CFP) [2] and optical
coherence tomography (OCT) [3, 4], offer quantitative and noninvasive ways to
evaluate the optic nerve structure for glaucoma diagnosis. In Fig. 1 (top and
middle), we showcase detailed CFP and OCT samples, demonstrating the glaucoma
progression across three stages and highlighting the evolution of key clinical
measures, such as the vertical Cup-to-Disc Ratio (vCDR) and the thickness of
the RNFL and GCIPL.

(a) Normal

(b) Early

(c) Advanced

Figure 1: Multi-modal imaging from the GAMMA dataset demonstrates glaucoma
progression across three stages. Column (a) features images without glaucoma.
Column (b) features early-stage glaucoma. Column (c) features intermediate-
advanced-stage glaucoma. The top row displays CFP images, showing increased
vCDR. The middle row displays OCT samples, showing thinning of the RNFL and
GCIPL. The bottom row displays vessel structures, extracted from CFP images.

In the field of computer-aided glaucoma diagnosis, previous studies have
primarily focused on using either CFP [5, 6, 7] or OCT [8, 9, 10]
individually. However, there is less exploration in multi-modality methods
that utilize both modalities simultaneously. The GAMMA challenge [11], which
introduces the first public multi-modality dataset, aims to advance the
research in multi-modality glaucoma grading using both CFP and OCT. MM-MIL
[12] employs over-sampling to augment CFP data, aiming to balance the multi-
modal data distribution. Corolla [13] introduces the extraction of the retinal
thickness map as an alternative modality, leading to more efficient
calculations and reduced memory usage. MM-RAF [14] develops a feature-level
fusion strategy that facilitates cross-modality correlation, with an emphasis
on leveraging spatial interaction between modalities through co-attention
mechanism. Although existing methods [11, 12, 13, 14] have demonstrated
impressive performance, they still face several challenges: First, current
methods simply concatenate feature embeddings from separate feature encoders
in the final stage, overlooking the unique uncertainty estimation inherent to
each modality. The approach that treats all modalities with equal weight can
result in unreliable predictions. Furthermore, given the inherently high
similarity of medical images and the scarcity of labeled medical datasets,
standard supervised learning methods often produce sub-optimal and
inconsistent feature representations.

To address the challenges discussed above, we propose a novel framework,
namely ETSCL, which consists of a contrastive feature extraction stage and a
decision-level fusion stage. In the contrastive feature extraction stage, we
employ the supervised contrastive loss to enhance the discriminative power,
resulting in more effective feature extraction. In addition, we utilize the
Frangi vesselness algorithm [15] as a preprocessing step to incorporate
previously neglected vessel information as a third complementary branch
alongside the CFP and OCT branches. In the decision-level fusion stage, we
propose a multi-modality classifier based on the evidence theory [16, 17]. The
proposed classifier uses the belief functions to quantitatively estimate the
uncertainty of information from different modalities, and combines the
information based on Dempster’s rule of combination. To the best of our
knowledge, our framework represents the first attempt of applying decision-
level fusion and vessel information in the multi-modality glaucoma grading
problem. Extensive experiments on the public GAMMA dataset [11] demonstrate
that our method achieves state-of-the-art (SOTA) performance.

##  2 Methodology

As shown in Fig. 2, the framework is trained in two stages. In the contrastive
feature extraction stage, alongside the CFP and OCT branches, vessel
information is extracted from the CFP images as a third complementary branch.
Vascular genesis is considered a key triggering factor in the development of
glaucoma [18, 19]. However, vessel information has been overlooked in previous
glaucoma grading methods. As a preprocessing step, we employ the Frangi filter
[15] to extract vessel information from the CFP images. Fig. 1 (bottom)
illustrates the extracted vessel structures. Then the supervised contrastive
loss is employed to train the three feature encoders separately, generating
reliable feature embeddings. In the decision-level fusion stage, we develop an
evidence theory-based classifier to fuse information from multiple sources,
utilizing uncertainty estimation across different modalities.

Figure 2: Overview of the proposed ETSCL framework.

###  2.1 Contrastive Feature Extraction

In clinical practice, extracting glaucoma-indicative features from CFP and OCT
images using traditional convolutional neural networks is challenging due to
the high similarity of medical images and the unbalanced multi-modal data
distribution. Compared to traditional supervised learning, contrastive
learning [20, 21, 22, 23] leverages the semantic similarity between data
points to generate more stable representations for downstream tasks. This
strategy has influenced the development of existing methods such as Corolla
[13] and MM-RAF [14]. As inspired by the work of [22], we implement a
supervised contrastive learning strategy. Given $N$ randomly sampled pairs
$\\{\,x_{h},y_{h}\,\\}_{h=1}^{N}$, the augmented batch for contrastive
learning comprises $2N$ pairs
$\\{\,\tilde{x}_{h},\tilde{y}_{h}\,\\}_{h=1}^{2N}$, where
$y_{h}=\tilde{y}_{2h}=\tilde{y}_{2h-1}$. Let $i\in I=\\{1,\ldots,2N\\}$ be the
index of an augmented sample. The positive is the index of the other augmented
sample derived from the same sample $i$. The embedding vectors $\bm{z}$, with
$\bm{z_{i}}$ representing the embedding vector for sample $i$, are derived via
the encoder and projection head. The supervised contrastive loss [22] is
generalized from the self-supervised contrastive loss [20, 21], and takes the
form of:

$\displaystyle L^{sup}$ $\displaystyle=\sum_{i\in
I}\frac{-1}{\left|P(i)\right|}\sum_{p\in
P(i)}\log\frac{\exp(\bm{z_{i}}\cdot\bm{z_{p}}/\tau)}{\sum_{a\in
A(i)}\exp(\bm{z_{i}}\cdot\bm{z_{a}}/\tau)},$ (1)

where the $\cdot$ symbol denotes the inner dot product, $\tau$ is a scalar
temperature parameter, $A(i)=I\setminus\\{i\\}$, and $P(i)=\\{p\in
A(i):\tilde{y}_{p}=\tilde{y}_{i}\\}$ is the set of indices of all the
positives in the augmented batch. The supervised contrastive loss leverages
label information to generate the positives, making it more effective at
grouping features of the same class closely together and separating those of
different classes.

###  2.2 Decision-Level Fusion

After feature extraction, the information from different modalities is
integrated in the classifier. The evidence theory, also known as the Dempster-
Shafer Theory (DST) [16, 17], introduces the concept of belief to quantify the
uncertainty of evidence. Informed by recent developments in DST-based multi-
view classification [24, 25], we design a multi-modality classifier to
integrate multi-source information more effectively. In the decision-level
fusion stage, the output of the multilayer perceptron (MLP) is regarded as
evidence, denoting the support for each class. For each modality, we employ
the Softplus function on the feature embeddings to obtain the evidence values
of all classes. The Dirichlet distribution is adopted to produce more reliable
class probabilities for each modality, with the loss function designed
accordingly. Subsequently, Dempster’s rule of combination is employed to
integrate multiple Dirichlet distributions from different modalities.

####  2.2.1 Dirichlet distribution

Given $N$ labeled samples $\\{(x_{1},y_{1}),\ldots,(x_{N},y_{N})\\}$,
conventional classifiers employ a Softmax layer to estimate the class
probabilities. Mathematically, the Softmax layer’s outputs can be interpreted
as parameters for a multinomial distribution. For a sample with index $i$, the
probabilities of $K$ classes are $\bm{\mu_{i}}=[\mu_{i1},\ldots,\mu_{iK}]$.
The optimization process of Maximum Likelihood Estimation (MLE), with respect
to the network’s parameters $\theta$, is defined as follows:

$\displaystyle\arg\max_{\theta}\mathrm{Mult}(\mathbf{Y};\theta)=\prod_{i=1}^{N}\prod_{j=1}^{K}\mu_{ij}^{y_{ij}},$
(2)

where $y_{ij}$ is an indicator that is 1 if the sample $i$ belongs to class
$j$ and 0 otherwise. Minimizing the log-likelihood, an approach widely known
as the cross-entropy loss, essentially mirrors the objective of MLE. Given its
limitations as a frequentist technique, MLE often struggles to effectively
handle varying conditions and small data scenarios, which can result in
overconfidence in the estimates. Alternatively, we employ the Dirichlet
distribution as the conjugate prior for the multinomial distribution, allowing
for more reliable and robust statistical inference. For the sample with index
$i$, Dirichlet distribution is defined by $K$ parameters
$\bm{\alpha_{i}}=[\alpha_{i1},\ldots,\alpha_{iK}]$. The probability density
function of the Dirichlet distribution is:

$\displaystyle\mathrm{Dir}(\bm{\mu_{i}}\mid\bm{\alpha_{i}})=\begin{cases}\frac{1}{B(\bm{\alpha_{i}})}\prod_{j=1}^{K}\mu_{ij}^{\alpha_{ij}-1}&\text{for
}\bm{\mu_{i}}\in S_{unit}^{K-1},\\\ 0&\text{otherwise},\end{cases}$ (3)

where $S_{unit}^{K-1}$ denotes the $K-1$ dimensional unit simplex, and
$B(\bm{\alpha_{i}})$ represents the $K$-dimensional multivariate beta
function. Consequently, label prediction for sample $i$ is modeled as a
generative process:

$\displaystyle\bm{\mathbf{y}_{i}}\sim\mathrm{Mult}(\bm{\mathbf{y}_{i}}\mid\bm{\mu_{i}}),\quad\bm{\mu_{i}}\sim\mathrm{Dir}(\bm{\mu_{i}}\mid\bm{\alpha_{i}}).$
(4)

To design the loss function for the Dirichlet distribution of a single sample
with index $i$ and label $\bm{\mathbf{y}_{i}}$, given the Dirichlet
distribution parameters $\bm{\alpha_{i}}$, the Mean Square Error (MSE) loss is
derived as follows:

$\displaystyle\mathrm{MSE}_{i}$ $\displaystyle=\int_{\bm{\mu_{i}}\in
S_{unit}^{K-1}}\|\bm{\mathbf{y}_{i}}-\bm{\mu_{i}}\|_{2}^{2}\cdot\mathrm{Dir}(\bm{\mu_{i}}\mid\bm{\alpha_{i}})d\bm{\mu_{i}}=\mathbb{E}\left[\sum_{j=1}^{K}\left(y_{ij}-\mu_{ij}\right)^{2}\right]$
$\displaystyle=\sum_{j=1}^{K}\left(y_{ij}^{2}-2y_{ij}\mathbb{E}\left[\mu_{ij}\right]+\mathbb{E}\left[\mu_{ij}^{2}\right]\right)=\sum_{j=1}^{K}\left(y_{ij}-\mathbb{E}\left[\mu_{ij}\right]\right)^{2}+\mathrm{Var}(\mu_{ij})$
$\displaystyle=\sum_{j=1}^{K}\left(y_{ij}-\frac{\alpha_{ij}}{S_{i}}\right)^{2}+\frac{\alpha_{ij}\left(S_{i}-\alpha_{ij}\right)}{S_{i}^{2}\left(S_{i}+1\right)},$
(5)

where $S_{i}=\sum_{j=1}^{K}\alpha_{ij}$, known as the Dirichlet strength. To
counteract any drift towards an unrealistic distribution, a penalty term
derived from the Kullback-Leibler (KL) divergence, is incorporated to adjust
the Dirichlet distribution of the sample $i$ accordingly. Consequently, the
loss function for this Dirichlet distribution is:

$\displaystyle
L_{i}=\mathrm{MSE}_{i}+\lambda_{i}\cdot\mathrm{KL}[\mathrm{Dir}(\bm{\mu_{i}}\mid\bm{\alpha_{i}})\,||\,\mathrm{Dir}(\bm{\mu_{i}}\mid\mathbf{1})],$
(6)

where $\mathbf{1}$ represents a $K$-dimensional unit vector, and $\lambda_{i}$
is the annealing coefficient.

####  2.2.2 Belief, uncertainty and combination

Subjective logic (SL) [26] provides a framework that links the Dirichlet
distribution parameters to the belief and uncertainty associated with each
class. For the Dirichlet distribution of a single sample, SL assigns a belief
mass $b_{j}$ to class $j$, reflecting the confidence in the occurrence of each
class. An aggregate uncertainty mass $u$ quantifies the overall uncertainty
across all classes. The summation of belief masses of all classes and the
uncertainty mass equals 1, ensuring a normalized representation of belief and
uncertainty. Given the evidence $e_{j}$ for class $j$, the Dirichlet
distribution parameter, $\alpha_{j}$, is calculated as $\alpha_{j}=e_{j}+1$.
The belief masses $\mathbf{b}$ and the uncertainty mass $u$ are derived from
the Dirichlet distribution parameters as follows:

$\displaystyle\mathbf{b}=\frac{\mathbf{e}}{S}=\frac{\bm{\alpha}-1}{S},\quad
u=\frac{K}{S},$ (7)

where $\bm{\alpha}-1$ indicates the subtraction of 1 from each parameter in
the Dirichlet distribution parameter vector $\bm{\alpha}$, and
$S=\sum_{j=1}^{K}(e_{j}+1)=\sum_{j=1}^{K}\alpha_{j}$, known as the Dirichlet
strength. We then adopt a reduced Dempster’s combination rule [24] to manage
computational complexity by avoiding the large numbers of input and output
masses in the original rule. For the fusion of the CFP and OCT branches, the
mass set of the fusion $M^{F}=\\{b_{1}^{F},\ldots,b_{K}^{F},u^{F}\\}$ is
calculated from the mass sets of the CFP branch
$M^{c}=\\{b_{1}^{c},\ldots,b_{K}^{c},u^{c}\\}$ and the OCT branch
$M^{o}=\\{b_{1}^{o},\ldots,b_{K}^{o},u^{o}\\}$ as follows:

$\displaystyle
b_{j}^{F}=\frac{1}{1-K^{\prime}}(b_{j}^{c}b_{j}^{o}+b_{j}^{c}u^{o}+b_{j}^{o}u^{c}),u^{F}=\frac{1}{1-K^{\prime}}u^{c}u^{o},$
(8)

where $K^{\prime}=\sum_{j\neq i}b_{j}^{c}b_{i}^{o}$ is the conflict measure
between the CFP and OCT branches. Then we apply the same computation to
combine this CFP and OCT fusion $M^{F}$ with the Vessel branch to calculate
the final mass set. Subsequently, based on the final mass set, we can obtain
the corresponding parameters of the joint Dirichlet distribution.

####  2.2.3 Evidence theory-based classifier

After deriving the joint Dirichlet distribution from different modalities, the
class with the highest probability is identified as the final prediction
output. Based on Eq. (6), the loss for each modality is calculated as the
summation of losses for all samples: $\sum_{i=1}^{N}L_{i}$. Consequently, we
have the losses for the CFP, OCT, and Vessel branches, denoted as $L^{c}$,
$L^{o}$, and $L^{v}$, respectively. Additionally, we sum up the loss of the
joint Dirichlet distribution across all samples to calculate the fusion loss,
$L^{F}$. The overall loss function of the classifier is then calculated as the
summation of the losses from all the modalities and the fusion loss, as
illustrated in Fig. 2.

##  3 Experiments

###  3.1 Dataset

The public GAMMA dataset [11] consists of 300 pairs of CFP images and OCT
volumes, with each OCT volume consisting of 256 2D slices. The dataset is
categorized into three labels: non-glaucoma, early-glaucoma, and intermediate-
advanced-glaucoma, presenting a multiclass classification task. We divide the
dataset into 200 samples for training and 100 samples for testing.

###  3.2 Implementation Details

All experiments are conducted on two NVIDIA A40 GPUs. For image augmentation,
we apply normalization, color jittering, random grayscale conversion, center
cropping, and random horizontal flipping. In the feature extraction stage, we
utilize ResNet50, pretrained on ImageNet, as the backbone for all the three
branches. We process the 3D OCT inputs as 256-channel 2D images. For each
feature encoder, we then employ 2 fully-connected layers with ReLU activation
to generate a 128-dimensional tensor, which is subsequently used for the
supervised contrastive loss $L^{sup}$, with the temperature parameter $\tau$
set to 0.05. The network is trained using Adam optimizer, with a learning rate
of 0.001 and a batch size of 14, over 10 epochs. In the decision-level fusion
stage, the embeddings generated by each feature encoder from the previous
stage are independently passed to separate MLP networks. Each MLP, with an
input size of 128 and an output size of 3, is trained using the Adam optimizer
over 200 epochs.

###  3.3 Results and Ablation Study

Following the evaluation approach in the GAMMA challenge, we use Cohen’s kappa
coefficient as the primary metric and accuracy as the secondary metric. For
this ordinal ternary classification task, kappa is quadratically weighted to
reflect the varying degrees of disagreement across ordinal categories [11],
which is particularly effective in assessing models on unbalanced datasets.

We select both single-modality and multi-modality methods as the baselines.
For the single-modality baselines, ResNet50 is selected for CFP, while
ResNet50 and 3D EfficientNet are selected for OCT. For the multi-modality
baselines, we choose a dual-branch ResNet50, and a hybrid of ResNet50 for CFP
and 3D EfficientNet for OCT. Additionally, the GAMMA challenge winner’s method
is selected as another multi-modality baseline. As shown in Table 1, our
method outperforms the best baseline, the GAMMA challenge winner, by over 3%
in kappa and by 3% in accuracy.

Table 1: Performance of the baselines and our ETSCL on the GAMMA dataset. Method | Kappa | Accuracy
---|---|---
Single-modality | CFP - ResNet | 0.6199 | 0.68
OCT - ResNet | 0.5006 | 0.62
OCT - 3D EfficientNet | 0.4930 | 0.60
Multi-modality | CFP+OCT - Dual ResNet | 0.7129 | 0.75
CFP+OCT - ResNet+3D EfficientNet | 0.7200 | 0.75
GAMMA challenge winner [11] | 0.8535 | 0.81
ETSCL | 0.8844 | 0.84
Table 2: Ablation study results on the GAMMA dataset. CFP means the presence of the CFP branch. OCT means the presence of the OCT branch. Vessel denotes the addition of the Vessel branch. SCL indicates the use of the supervised contrastive loss. ET represents the evidence theory-based classifier. CFP | OCT | Vessel | SCL | ET | Kappa | Accuracy
---|---|---|---|---|---|---
✓ |  |  |  |  | 0.6199 | 0.68
✓ | ✓ |  |  |  | 0.7129 | 0.75
✓ | ✓ | ✓ |  |  | 0.7478 | 0.77
✓ |  |  | ✓ |  | 0.8170 | 0.77
✓ | ✓ |  | ✓ |  | 0.8273 | 0.80
✓ | ✓ | ✓ | ✓ |  | 0.8592 | 0.81
✓ | ✓ | ✓ | ✓ | ✓ | 0.8844 | 0.84

As shown in Table 2, the improvements observed when adding the supervised
contrastive loss to conventional CFP, dual branches (CFP and OCT), and triple
branches (CFP, OCT, and Vessel) demonstrate the significant contribution of
the supervised contrastive loss. Furthermore, comparisons between dual
branches (CFP and OCT) and triple branches (CFP, OCT, and Vessel), both with
and without the supervised contrastive loss, validate the benefit of
incorporating the vessel information. Lastly, the advancement from using the
conventional Softmax layer to using the evidence theory-based classifier in
the triple branches (CFP, OCT, and Vessel) highlights the superior performance
of the evidence theory-based classifier.

In Table 3, we present a comparison of our ETSCL with other SOTA methods that
were published after the GAMMA challenge. We include only the kappa values as
reported in their respective publications, since not all studies disclose
their accuracy. Among these published SOTA methods, our ETSCL achieves the
highest kappa. It is worth noting that this comparison might not be strict,
due to variations in the datasets used by each method.

Table 3: Comparison between our ETSCL and other SOTA methods on the GAMMA dataset. $\dagger{}$ Corolla uses 100 samples for training and testing. * MM-RAF trains on a private dataset and tests with 100 GAMMA samples on both its method and MM-MIL. Method | Kappa
---|---
Corolla$\dagger{}$ [13] | 0.8550
MM-MIL* [12] | 0.8562
MM-RAF* [14] | 0.8467
ETSCL | 0.8844

##  4 Conclusion

In this study, we propose an evidence theory-based supervised contrastive
learning framework for multi-modality glaucoma grading. The supervised
contrastive loss is employed to extract effective feature embeddings. To
integrate multi-modality information with uncertainty estimation, we design an
evidence theory-based classifier. In addition, vessel information is utilized
to assist in the prediction. However, the limited size of the GAMMA dataset
constrains our method choices, such as Vision Transformer which often leads to
overfitting. The lack of private datasets restricts our ability to further
test the generalizability and robustness of our method, as the GAMMA dataset
is currently the only high-quality public dataset with multi-stage glaucoma
grading. Furthermore, subsequent generalizability tests can reinforce the
validation of the Vessel branch’s efficacy, considering the novelty of
incorporating vessel information.

####  4.0.1 Acknowledgements

This project is supported by the Ministry of Education, Singapore, under its
Academic Research Fund Tier 1 (RS16/23).

####  4.0.2

The authors have no competing interests to declare that are relevant to the
content of this article.

## References

  * [1] Weinreb, R.N., Leung, C.K., Crowston, J.G., Medeiros, F.A., Friedman, D.S., Wiggs, J.L., Martin, K.R.: Primary open-angle glaucoma. Nature Reviews Disease Primers 2(1), 1–19 (2016)
  * [2] Besenczi, R., Tóth, J., Hajdu, A.: A review on automatic analysis techniques for color fundus photographs. Computational and Structural Biotechnology Journal 14, 371–384 (2016)
  * [3] Huang, D., Swanson, E.A., Lin, C.P., Schuman, J.S., Stinson, W.G., Chang, W., Hee, M.R., Flotte, T., Gregory, K., Puliafito, C.A., et al.: Optical coherence tomography. Science 254(5035), 1178–1181 (1991)
  * [4] De Fauw, J., Ledsam, J.R., Romera-Paredes, B., Nikolov, S., Tomasev, N., Blackwell, S., Askham, H., Glorot, X., O’Donoghue, B., Visentin, D., et al.: Clinically applicable deep learning for diagnosis and referral in retinal disease. Nature Medicine 24(9), 1342–1350 (2018)
  * [5] Al-Bander, B., Al-Nuaimy, W., Al-Taee, M.A., Zheng, Y.: Automated glaucoma diagnosis using deep learning approach. In: Proceedings of the 14th International Multi-Conference on Systems, Signals & Devices (SSD). pp. 207–210 (2017)
  * [6] Fang, H., Li, F., Wu, J., Fu, H., Sun, X., Son, J., Yu, S., Zhang, M., Yuan, C., Bian, C., et al.: Refuge2 challenge: A treasure trove for multi-dimension analysis and evaluation in glaucoma screening. arXiv preprint arXiv:2202.08994 (2022)
  * [7] He, H., Lin, L., Cai, Z., Tang, X.: Joined: Prior guided multi-task learning for joint optic disc/cup segmentation and fovea detection. In: Proceedings of the International Conference on Medical Imaging with Deep Learning (MIDL). pp. 477–492 (2022)
  * [8] Asaoka, R., Hirasawa, K., Iwase, A., Fujino, Y., Murata, H., Shoji, N., Araie, M.: Validating the usefulness of the “random forests” classifier to diagnose early glaucoma with optical coherence tomography. American Journal of Ophthalmology 174, 95–103 (2017)
  * [9] Maetschke, S., Antony, B., Ishikawa, H., Wollstein, G., Schuman, J., Garnavi, R.: A feature agnostic approach for glaucoma detection in oct volumes. PLoS One 14(7), e0219126 (2019)
  * [10] Wang, H., Guo, X., Song, K., Sun, M., Shao, Y., Xue, S., Zhang, H., Zhang, T.: Octformer: An efficient hierarchical transformer network specialized for retinal optical coherence tomography image recognition. IEEE Transactions on Instrumentation and Measurement 72, 1–17 (2023)
  * [11] Wu, J., Fang, H., Li, F., Fu, H., Lin, F., Li, J., Huang, Y., Yu, Q., Song, S., Xu, X., et al.: Gamma challenge: glaucoma grading from multi-modality images. Medical Image Analysis 90, 102938 (2023)
  * [12] Li, X., Zhou, Y., Wang, J., Lin, H., Zhao, J., Ding, D., Yu, W., Chen, Y.: Multi-modal multi-instance learning for retinal disease recognition. In: Proceedings of the 29th ACM International Conference on Multimedia. pp. 2474–2482 (2021)
  * [13] Cai, Z., Lin, L., He, H., Tang, X.: Corolla: An efficient multi-modality fusion framework with supervised contrastive learning for glaucoma grading. In: Proceedings of the 19th International Symposium on Biomedical Imaging (ISBI). pp. 1–4 (2022)
  * [14] Zhou, Y., Yang, G., Zhou, Y., Ding, D., Zhao, J.: Representation, alignment, fusion: A generic transformer-based framework for multi-modal glaucoma recognition. In: Proceedings of the 26th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI). pp. 704–713 (2023)
  * [15] Frangi, A.F., Niessen, W.J., Vincken, K.L., Viergever, M.A.: Multiscale vessel enhancement filtering. In: Proceedings of the First International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI). pp. 130–137 (1998)
  * [16] Shafer, G.: A Mathematical Theory of Evidence, vol. 42. Princeton University Press, Princeton, NJ (1976)
  * [17] Dempster, A.P.: Upper and lower probabilities induced by a multivalued mapping. In: Classic Works of the Dempster-Shafer Theory of Belief Functions, vol. 219, pp. 57–72. Springer Berlin Heidelberg, Berlin, Heidelberg (2008)
  * [18] Cherecheanu, A.P., Garhofer, G., Schmidl, D., Werkmeister, R., Schmetterer, L.: Ocular perfusion pressure and ocular blood flow in glaucoma. Current Opinion in Pharmacology 13(1), 36–42 (2013)
  * [19] Wang, L., Cull, G.A., Fortune, B.: Optic nerve head blood flow response to reduced ocular perfusion pressure by alteration of either the blood pressure or intraocular pressure. Current Eye Research 40(4), 359–367 (2015)
  * [20] Chen, T., Kornblith, S., Norouzi, M., Hinton, G.: A simple framework for contrastive learning of visual representations. In: Proceedings of the 37th International Conference on Machine Learning. pp. 1597–1607 (2020)
  * [21] Henaff, O.: Data-efficient image recognition with contrastive predictive coding. In: Proceedings of the 37th International Conference on Machine Learning. pp. 4182–4192 (2020)
  * [22] Khosla, P., Teterwak, P., Wang, C., Sarna, A., Tian, Y., Isola, P., Maschinot, A., Liu, C., Krishnan, D.: Supervised contrastive learning. Advances in Neural Information Processing Systems 33, 18661–18673 (2020)
  * [23] Chen, X., Xie, S., He, K.: An empirical study of training self-supervised vision transformers. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 9640–9649 (2021)
  * [24] Han, Z., Zhang, C., Fu, H., Zhou, J.T.: Trusted multi-view classification with dynamic evidential fusion. IEEE Transactions on Pattern Analysis and Machine Intelligence 45(2), 2551–2566 (2023)
  * [25] Liu, W., Chen, Y., Yue, X., Zhang, C., Xie, S.: Safe multi-view deep classification. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 37, pp. 8870–8878 (2023)
  * [26] Jøsang, A.: Subjective Logic: A Formalism for Reasoning Under Uncertainty. Springer Cham (2016)


#  Turbulent mixing and a generalized
phase transition in shear-thickening fluids

Helmut Z. Baumert  IAMARIS, Ludwigslust, Germany
<EMAIL_ADDRESS> Bernhard Wessling  LaoWei Chemical Technology
Consulting Inc., Shenzhen, China<EMAIL_ADDRESS>
###### Abstract

This paper presents a new theory of turbulent mixing in stirred reactors. The
degree of homogeneity of a mixed fluid may be characterized by (among other
features) Kolmogorov’s microscale, $\lambda$. The smaller its value, the
better the homogeneity. According to Kolmogorov, $\lambda$ scales inversely
with the fourth root of the energy flux applied in the stirring process,
$\varepsilon$. The higher $\varepsilon$, the smaller $\lambda$, and the better
the homogeneity in the reactor. This is true for Newtonian fluids. In non-
Newtonian fluids the situation is different. For instance, in shear-thickening
fluids it is plausible that high shear rates thicken the fluid and might
strangle the mixing. The internal interactions between different fluid-
mechanical and colloidal variables are subtle, namely due to the (until
recently) very limited understanding of turbulence.
Starting from a qualitatively new turbulence theory for inviscid fluids
[Baumert(2013)], giving e.g. Karman’s constant as $1/\sqrt{2\,\pi}=0.40$ [[,
Princeton’s superpipe gives 0.40$\pm 0.02$,]]baileyetal2014, we generalize
this approach to the case of viscous fluids and derive equations which in the
steady state exhibit two solutions. One solution branch describes a state of
good mixing, the other strangled turbulence. The physical system cannot be in
two different steady states at the same time. But it is physically admissible
to switch between the two steady states by non-stationary transitions, maybe
in a chaotic fashion.

Keywords: Turbulence, kinematic viscosity, eddy viscosity, colloidal systems,
dispersions, shear-thickening, mixing, stirring, reactor, non-Newtonian fluid,
Kolmogorov scale, bifurcation

35pc Corresponding author: Dr. Helmut Z. Baumert,
IAMARIS, Bürgermeister-Jantzen-Str. 3, 19288 Ludwigslust, Germany
<EMAIL_ADDRESS>Submitted for publication to Physica Scripta (IOP): 29
March, 2015

##  1 Introduction

### Plausibility considerations

A mixer is a reactor wherein fluids are mixed with other fluids, with dry or
liquid chemicals, heat or other so-called ‘scalars’. To that aim propellers,
screws, raising bubbles or other technological stirring means and devices are
applied. The near-field turbulence initiated by them is characterized by two
parameters: a characteristic length scale and a characteristic time scale (or,
equivalently, a corresponding wavenumber and a corresponding frequency). These
two parameters may be tuned independently by modifying the geometry and the
rotation rate of the device (the bubble-release frequency in bubble or
convective mixing).

If the operator of a mixer wishes a higher degree of homogeneity then he is
inclined to increase the energy spend for mixing by increasing the rotation
rate, the geometric parameters of the propeller, or both. For Newtonian fluids
this is a reasonable strategy. However, in shear-thickening fluids it is
plausible that with an increase of the forcing one might contribute to
thickening which in the extreme strangles mixing until breakdown of
turbulence. Here much depends on details, i.e. how one increases the stirring
power, by geometric means, by the rotation rate, or by both.

As mixing and dispersion are fundamental technological steps in nearly all
industries and connected with huge financial fluxes if considered for the
globe, this circle of questions deserves attention.

### Historical remarks

Whereas the arts of hydraulics including the phenomenology of eddying flows
and vortices are well known at least for about $2.4\times 10^{3}$ years, back
to e.g. the ‘hydraulic empires’ of Nofretete and Echnaton, the mathematical
science of dynamic fluid motions emerged ’just recently’: in idealized, purely
geometric form by Leonhard Euler (1707 - 1783), supplemented somewhat later
with frictional aspects by Claude-Louis Navier and George Stokes. The theory
of vortices and circulation was more explicitely elaborated about 30 years
later by Herrmann von Helmholtz and William Thomson (Lord Kelvin), and later
came Osborne Reynolds, Ludwig Prandtl, Theodore von Karman, Alexander A.
Friedmann and Lev V. Keller with their specific innovations.

A systematic statistical treatment of turbulence was initiated by Geoffrey I.
Taylor and, in form of spectral considerations, by Andrej N. Kolmogorov. These
attempts were more or less singular pieces, kept loosely together by semi-
empirical scaling rules and similarity arguments. Great contributions in this
respect stem here from Andrej S. Monin, Alexander Obukhov, Rostislav V.
Osmidov, Chuck van Atta and Stephen Thorpe.

With the advent of sufficiently powerful computers the focus moved towards so-
called two-equation turbulence models. Major authors are David C. Wilcox,
George L. Mellor, Tetsuji Yamada, and Wolfgang Rodi. The overwhelming number
of turbulence theories in the above sense were elaborated for Newtonian fluids
(and gases) wherein kinematic viscosity is more or less constant in larger
elements of space and time and not dependent neither on the state of the fluid
motion nor on external forces, nor on time.

### Smart fluids / colloidal systems

Today so-called smart fluids become increasingly relevant in various
industries, from medicine over small-scale bio-technology to construction
engineering, waste-water treatment and defence. The central property of a
smart fluid is its kinematic viscosity, which may explicitely depend on time
or not, which may be tuned by shear or light, by electric or magnetic fields
like in magneto-rheological fluids.

As their kinematic viscosity is not a constant, smart fluids have always non-
Newtonian character. They all are at least 2-phase systems, whereas the – at
least one – dispersed phase is mostly nanoscopic (i.e., less than 100 nm
diameter) in size. Hence these are colloidal systems, either colloidal
dispersions (the dispersed phase being a solid) or emulsions (the dispersions
being a liquid). In contrast to wide-spread assumptions, such colloidal
systems exhibit complex structures. I.e., the dispersed phase is not
statistically evenly distributed, but highly structured [Wessling(1993)]. This
is due to the thermodynamically non-equilibrium character of such colloidal
dispersions / emulsion systems [Wessling(1991), Wessling(1995)].

As long as the flow of a non-Newtonian fluid or gas is laminar, the
consequences of those external controls on the kinematic viscosity remain more
or less predictable and no dispersion would be possible [Wessling(1995)]. This
property is going to be changed when the flow becomes critical and eventually
turbulent.

Even the Reynolds number in non-Newtonian fluids looses much of its classical
meaning. Traditionally it is defined as the dimensionless ratio of inertial
and viscous forces,

$Re=\frac{U\times L}{\nu}$ (1)

where $U$ [m s-1] is a velocity scale for the mean flow, $L$ [m] a
characteristic length scale (e.g. a distance from a solid wall or body) and
$\nu$ is the kinematic viscosity [m${}^{2}\,$s-1]. This means that $Re$
depends linearly on the state of flow via the velocity, $U$. However, if also
$\nu$ depends on the flow state then things become more difficult.

Relatively simple cases of smart fluids are the so-called shear-thickening
(dilatant) and the shear-thinning (pseudoplastic) fluids. In the former case
viscosity increases with increasing shear rate so that the fluid becomes
‘thicker’. The most prominent example is the fluid in body armor wests of
policemen. Certain water-sand mixtures may behave in similar form. Prominent
examples for pseudoplastic fluids are lava, blood and whipped cream. Also the
notorious quicksand is an example for a shear-thinning material.

In all these cases the turbulent state of those fluids – if at all – is
understood only in an empirical sense.

Figure 1: Examples of shear behavior in a shear-thickening (green, m=0.65), a
shear-thinning (red) and a Newtonian fluid.

The kinematic viscosity of many non-Newtonian fluids may be described as a
power-law function of local shear rate, $S$, by the following slightly
generalized Ostwald-de Waehle ansatz, wherein $\nu_{0}=\nu_{0}(T)$ and
$\gamma=\gamma(T)$ are (generally temperature-dependent) empirical parameters:

$\nu=\nu_{0}+\gamma\times({S}/{S_{*}})^{m}\,.$ (2)

Fig. 1 shows three example fluids with (not completely fictive) parameters
like $\gamma=10^{-6}$ m2 s-1, $S_{*}=0.19$ s-1 (Hz), $\nu_{0}=0.43$ m2 s-1,
and with $m=\\{-0.5;+0.65;+1.5\\}$.

In the following we exclusively deal with shear-thickening fluids ($m>0$) and
use the notion of turbulence in a narrow sense defined later. Nevertheless the
theoretical foundations cover a broader range of applications than only the
shear-thickening case.

### Theoretical foundation

This paper rests essentially on [Baumert(2013)] who demonstrated that the
turbulent viscosity, $\nu_{t}$ [m2 s-1], scaling linearly with turbulent
diffusivity, $\mu_{t}$ [m2 s-1], may be written as follows,

$\nu_{t}=\frac{\cal K}{\pi\,\Omega}$ (3)

where $\cal K$ is turbulent kinetic energy in the narrow sense, TKE [m2 s-2],
$\Omega$ is the r.m.s. turbulent vorticity [s-1] (a relative of enstrophy) and
$\pi$ is the dimensionless circle number. This formula results from a
stochastic-geometric theory of inviscid turbulence wherein turbulent eddies
are taken as singular solutions of the Euler equation and thus as particles
(vortex-tubes dipoles) moving in a stochastic fashion – just like molecules in
Einstein’s theory of Brownian motion – slightly generalized in form of Fokker-
Planck equations wtih generation and dissipation terms. Another result of this
theory is the Karman constant as

${\ae}=1{/}\sqrt{2\pi}=0.40\;.$ (4)

[[, A recent study based on the Princeton Superpipe gave $0.40\pm 0.02$,
see]]baileyetal2014.

##  2 Inviscid and viscous spectra

### The major length scales

Wavenumber spectra of developed viscous turbulence exhibit two distuingished
scales: the larger, so-called energy-containing scale which depends on the
forcing mechanism, $\Lambda$, and the smaller Kolmogorov scale,
$\lambda\ll\Lambda$ , which is related with the kinematic viscosity, $\nu$:

$\lambda=\left(\frac{\nu^{3}}{\varepsilon}\right)^{1/4}\;.$ (5)

Here $\varepsilon$ [m2 s-3] is the total111I.e. the integral over the
turbulent dissipation spectrum. turbulent dissipation rate of TKE, and $\nu$
[m2 s-1] is the kinematic viscosity. Obviously, if $\nu\rightarrow 0$, also
$\lambda\rightarrow 0$.

The energy-containing turbulent scale, $\Lambda$, has been derived on
theoretical grounds as follows [Baumert(2013)],

$\Lambda=\frac{\sqrt{{\cal K}/\pi}}{\Omega}\,.$ (6)

This scale is only indirectly influenced by viscosity.

Figure 2: Inviscid (dashed) and viscous (full) spectra of kinetic-energy
fluctuations, dimensionless. $\kappa_{\Lambda}=2\,\pi$ and
$\kappa_{\lambda}=2\,\pi\,{\Lambda}/\lambda$ denote the dimensionless forms of
the energy-containing and the Kolmogorov wavenumber, respectively.

### Wavenumber spectra

For a better understand of our approach we look at Fig. 2. It exhibits a real-
world viscous spectrum (full) and an idealized spectrum of inviscid turbulence
(dashed) discussed in [Baumert(2013)]. With the dimensionless wave number
$\kappa=k\times{\Lambda}$ [Pope(2000), eqn. 6.246, p. 232] it reads as
follows,

${\cal E}(\kappa)=\frac{E(\kappa)}{\varepsilon^{2/3}{\Lambda}^{5/3}}={C}\times
f_{\Lambda}(\kappa)\times\kappa^{-5/3}\times f_{\lambda}(\kappa,\rho)$ (7)

where

$\rho={\Lambda}/\lambda.$ (8)

$\rho$ is taken below as a rough proxy of the Reynols number.

$C$ in (7) is a so-called universal Kolmogorov constant for the wavenumber
spectrum of turbulence. It has been derived theoretically as follows
[Baumert(2013)]:

$C=\frac{1}{3}(4\pi)^{2/3}=1.80\,.$ (9)

In (7), $f_{\Lambda}(\kappa)$ describes the low-wavenumber, high-energy
interval following Karman’s spectral model [[, loc. cit.]p. 747]pope2000.
$f_{\lambda}(\kappa,\rho)$ denotes the so-called dissipation region wherein
the smallest eddy motions are becoming laminar [[, for details see p. 232, 233
in]]pope2000:

$f_{L}(\kappa)=\left(\frac{\kappa}{\sqrt{\kappa^{2}+c_{L}}}\right)^{17/3}$
(10)

and

$f_{\lambda}(\kappa,\rho)=\exp\left[-\beta\left(\left[\left(\frac{\kappa}{\rho}\right)^{4}+c_{\lambda}^{4}\right]^{1/4}-c_{\lambda}\right)\right]\,.$
(11)

[[, The exponential ansatz for the dissipation region is in a simpler form due
to Kraichnan (1959), loc. cit.]]pope2000.

Figure 3: Spectra of kinetic-energy fluctuations for the Reynolds-number proxy
$Re=L/\lambda$ decreasing from the upper right curve with $Re=\infty$ down to
the lower left with $Re=10$. We see that the area below Kolmogorov’s
$\kappa^{-5/3}$ inertial subrange not only shrinks quantitatively with
decreasing $Re$. It even vanishes completely, indicating the transition from
the turbulent to the laminar state.

##  3 Turbulence in the narrow sense
and other chaotic motions

In the present paper we understand the notion turbulence in a specific narrow
sense and separate it from other forms of chaotic fluid motions. We
distinguish three regions of the kinetic-energy wavenumber spectrum. We namely
interprete Kolmogorov’s inertial subrange in the interval
[$\kappa_{\Lambda},\kappa_{\lambda}$] as turbulence in the narrow sense and
the rest of the spectrum simply as other forms of chaotic fluid motions which
mix only weakly.

We differentiate between the following three wavenumber intervals:

  * (i)

$\kappa\in[0,\kappa_{\Lambda}]$: weakly mixing Karman range,

  * (ii)

$\kappa\in[\kappa_{\lambda},\kappa_{\Lambda}]$: intensely mixing inertial
subrange of Kolmogorov,

  * (iii)

$\kappa\in[\kappa_{\lambda},\infty]$: weakly mixing Kraichnan region.

Arguments for (ii) are given in the next section.

Regarding mixing of momentum in the sense of (3), only the Kolmogorov range
matters which in the inviscid case ($k_{\lambda}\rightarrow\infty$) extends
from the energy-containing scale, $k_{\Lambda}$, until $k_{\lambda}=\infty$ so
that

${\cal
K}=\int_{k_{\Lambda}}^{\infty}E(k)\,dk=C\times\varepsilon^{2/3}\times\int_{k_{\Lambda}}^{\infty}\frac{dk}{k^{5/3}}\;.$
(12)

This last integral gives together with (6) a simple analytical result:

${\cal
K}=C\times\frac{3}{2}\times\left(\frac{\varepsilon\,\Lambda}{2\pi}\right)^{2/3}\;.$
(13)

This result has successfully been tested against high-$Re$ turbulence
observations [[]Fig. 6]baumert2013.

Now we apply this approach to viscous turbulence and write

${\cal K}=\int_{k_{\Lambda}}^{k_{\lambda}}E(k)\,dk\;.$ (14)

Clearly, (12) is a special case of (14).

(14) is analytically integrable and gives the following:

${\cal
K}=C\times\frac{3}{2}\times\left(\frac{\varepsilon}{2\pi}\right)^{2/3}\times\left[{\Lambda}^{2/3}-{\lambda}^{2/3}\right].$
(15)

Due to (5), a kinematic viscosity $\nu$ increasing without bounds gives an
increasing $\lambda$ until $\lambda=\Lambda$ at which point according to (15)
TKE vanishes. This means that turbulence in the narrow sense (Kolmogorov’s
inertial subrange) disappears and chaotic laminar motions of the now joined
Karman-Kraichnan range remain. If ${\cal K}=0$, then eddy viscosity $\nu_{t}$
via (3) and consequently also eddy diffusivity $\mu_{t}$ vanish accordingly.

##  4 Mixing within and outside the Kolmogorov subrange

Our view of Kolmogorov’s inertial subrange has been filled with life through
reports on a series of numerical simulation studies [Herrmann(1990), Herrmann
et al.(1990)Herrmann, Mantica, and Bessis, Manna and Herrmann(1991)] wherein
spectra of space-filling bearing were shown to correspond closely to
Kolmogorov’s 5/3 spectral law. A space-filling bearing is identical with the
densest non-overlapping (Apollonian) circle packing in a plane, with side
condition that the circles are pointwise in contact but able to rotate freely,
without friction or slipping [[, a devil s gear sensu ]]poeppe04. The contact
condition for two different wheels with indices 1 and 2 of the gear reads

$\omega_{1}\times r_{1}=\omega_{2}\times r_{2}.$ (16)

A devil’s gear forms spontaneously when vortex dipoles collide. Those dipoles
move with their inctrinsic velocity governed by their given radius $r$ and
rotation speed $\omega$ [Baumert(2013)]. When they collide they are either
scattered or annihilated. The latter case is realized in form of a dissipative
patch wherein all scales $>0$ are frictionless except the scale zero where all
the dissipation takes place. Because a dissipative patch in form of a devil’s
gear remains at rest, scatter motions scaling with the intrinsic parameters
$r,\omega$ are the only source of place changes and thus of mixing
[Baumert(2013)].

Although in these patches scalar mixing takes also place, the bulk of it is
due to the place changes of dipoles. Hoever, momentum mixing is exclusively
due to place changes of the dipoles. This is the deeper reason why generally
$\mu_{t}>\nu_{t}$. I.e. the turbulent Prandtl number $\sigma=\nu_{t}/\mu_{t}$
is generally less than unity.

##  5 Balance equations for $\cal K$ and $\Omega$
in a homgeneously stirred reactor

### General

The dynamic balances of $\cal K$ and $\Omega$ in a homgeneously stirred
reactor are given as follows [Baumert(2013)]:

$\displaystyle\frac{d\cal K}{dt}$ $\displaystyle=$
$\displaystyle\nu_{t}\,{S}^{2}-\varepsilon_{h}-\varepsilon_{b}\,,$ (17)
$\displaystyle\frac{d\Omega}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\left(\frac{1}{2}{S}^{2}-\Omega^{2}\right).$ (18)

Here $\nu_{t}$ is the eddy viscosity given by (3), $S=\langle
S_{e}^{2}\rangle^{1/2}$ is the (steady) r.m.s. effective shear rate,
$\varepsilon_{h}={\cal K}\,\Omega/\pi$ [Baumert(2005)] is the dissipation rate
within the bulk volume, and $\varepsilon_{b}$ is the bulk dissipation within
the boundary layers (walls etc.) of the reactor.

In the following we denote steady state values ($d{\cal K}/dt=d\Omega/dt=0$)
by an overbar. In particular we have

$\overline{\varepsilon}_{b}=\overline{\varepsilon}_{h}=\frac{\overline{\cal
K}\,\overline{\Omega}}{\pi}\,,$ (19)

fully indepenent on the forcing ${S}$. Whereas the steady-state case of the
dynamic energy balance (17) is useless for the following, the steady-state
r.m.s. vorticity balance finds application further below:

$\overline{\Omega}={S}/\sqrt{2}\,,$ (20)

### Steady state

Until now we have no equation yet which gives us the size of Kolmogorov’s
microscale, $\lambda$. It is derived now, based on the above results and
thoughts. We begin with (15) and replace there the homogeneous volume
dissipation rate via
$\overline{\varepsilon}=\overline{\varepsilon}_{h}=\overline{\cal
K}\;\overline{\Omega}/\pi$ and find with (20) the following relation for the
steady-state turbulent kinetic energy in the reactor:

$\displaystyle\overline{\cal K}$ $\displaystyle=$
$\displaystyle\alpha_{1}\times\overline{S}^{2}\times\left[{\Lambda}^{2/3}-{\overline{\lambda}}^{2/3}\right]^{3}\,,$
(21) $\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle(2\pi)^{-2}\;\;\;=\;\;\;0.0253.$ (22)

This result allows to reformulate the steady-state eddy viscosity (3) in the
reactor as follows:

$\displaystyle\overline{\nu}_{t}$ $\displaystyle=$
$\displaystyle\frac{\overline{\cal
K}}{\pi\,\overline{\Omega}}\;\;=\;\;\alpha_{2}\times{S}\times\left[{\Lambda}^{2/3}-\overline{\lambda}^{2/3}\right]^{3},$
(23) $\displaystyle\alpha_{2}$ $\displaystyle=$ $\displaystyle
2^{-3/2}\times\pi^{-3}\;\;\;=\;\;\;0.0114.$ (24)

Analogously we may derive the steady-state TKE dissipation rate in the reactor
as

$\displaystyle\overline{\varepsilon}$ $\displaystyle=$
$\displaystyle\frac{\overline{\cal
K}\,\overline{\Omega}}{\pi}\;\;=\;\;\alpha_{3}\times{S}^{3}\times\left[{\Lambda}^{2/3}-\overline{\lambda}^{2/3}\right]^{3},$
(25) $\displaystyle\alpha_{3}$ $\displaystyle=$ $\displaystyle
2^{-5/2}\times\pi^{-3}\;\;\;=\;\;\;0.00570.$ (26)

##  6 Microscale bifurcation

### Governing equation

We now consider according to (5) the steady-state microscale,
$\overline{\lambda}$ and replace therein $\varepsilon$ with the reactor’s
steady-state value (25):

$\overline{\lambda}=\left(\frac{\overline{\nu}^{3}}{\overline{\varepsilon}}\right)^{1/4}=\alpha_{4}\times\left(\frac{\overline{\nu}/{S}}{\Lambda^{2/3}-\overline{\lambda}^{2/3}}\right)^{3/4},$
(27) $\alpha_{4}=\alpha_{3}^{-1/4}=1.1509\,.$ (28)

For an easy treatment of (27) we multiply both sides with
$\left[\Lambda^{2/3}-\overline{\lambda}^{2/3}\right]^{3/4}$ and get with some
algebra and using the abbreviation $\xi=\overline{\lambda}^{2}$ the following
working equation:

$(\xi\,\Lambda)^{2/3}=\xi+\alpha_{4}^{4/3}\,\overline{\nu}/{S}\,.$ (29)

Figure 4: Valid steady-state combinations of forcing parameters
$\Lambda,\lambda$. Note that always two states are possible, the upper branch
without turbulence in the narrow sense, and the lower branch with a developed
Kolmogorov spectrum with the property $\lambda\ll\Lambda$. The forcing length
scales are here $\Lambda$ = 10 cm, 2 cm, 1 cm, 5 mm, respectively, from
outside to inside. Figure 5: Similar to Fig. 4 but instead of $\lambda$ the
scale ratio $1/\rho=\lambda/\Lambda$ is shown. $\rho$ proper is a proxy for
the Reynolds number. $\rho=1$ means that the Kolmogorov scale equals the
energy-containing scale so that the highly mixing inertial subrange of the
fluctuation spectrum vanishes and mixing breaks more or less down. The forcing
length scales are here again $\Lambda$ = 10 cm, 2 cm, 1 cm, 5 mm (from outside
to inside). Note that the second horizontal line is situated at
$\lambda/\Lambda\approx 0.53$ and labels the border between the upper and
lower solution branches of (29).

### Physical interpretation

The meaning of (27) is most easily understood by highlighting input and output
variables and assuming mechanical stirring:

  * •

The shear frequency, $\mathbf{S}$, is controlled by the rotation speed of
(e.g.) the propeller(s) used to stir the reactor. I.e. for the solution of
(27), $\mathbf{S}$ is thus a given, prescribed input quantity.

  * •

The turbulent length scale $\mathbf{\Lambda}$ is controlled by geometry (e.g.
radius, length of stirring propeller etc.) and represents also a given,
prescribed input quantity.

  * •

The microscale $\overline{\mathbf{\lambda}}$ however is not given. It is the
dynamic relaxation result of the forcing and thus an output quantity.

As we see in Figures 4 and 5, for a given energy-containing length scale
$\Lambda$ (i.e. along a chosen curve), each admissible shear value $S$ allows
for two values of Kolmogorov’s microscale, $\lambda$:

  * a)

the upper value is associated with the strangling of turbulence due to the
close proximity of $\lambda$ and $\Lambda$. This proximity leaves not enough
space for a developed turbulence spectrum (in the narrow sense from above).

  * b)

The (mostly much) lower value is associated with a good mixing behavior of the
stirred reactor. Unfortunately nature gives no stability guarantee for this
state because there exists the “official” alternative of a strangeled spectrum
– in agreement with the law. Some aspects are discussed in the next Section 7.

##  7 The full dynamic problem

### Smart fluids, smart control?

We have shown in Section 6 that a reactor holding a shear-thickening fluid may
attain two different steady states with respect to scalar mixing and
turbulence. One of them mixes well, the other one does not.

Clearly the system cannot stay in both states at the same time. But non-steady
transitions between them are physically possible and of practical interest. In
Figures 4 and 5 these transitions occur ouside the computed curves.

Questions remain. E.g. how does the system choose, for a given set of initial
conditions $\\{{\cal K}_{0}$, $\Omega_{0}\\}$ and forcing parameters
$\\{\Lambda,S\\}$, between the two admissible steady states, and under which
conditions (including initialization) transitions between them occur. Possibly
smart fluids deserve smart means to keep them under control.

### Relaxation of the “inner fluid”

Due to the inner structure of shear-thickening colloidal systems we have to
expect that a viscosity according to (3) is never instantaneously realized.
Dynamic adjustments to a new micro-mechanical state take certain relaxation
steps. They are governed by the laws of irreversible thermodynamics
[Wessling(1991), Wessling(1993), Wessling(1995)]. Therefore, to get a complete
dynamic picture of the stirred fluid with our focus on learning how one of the
two steady states is selected and how potential transitions between them are
controllable, we actually have at least to augment the equation system (17,
18) by an additional equation for the dynamic relaxation of the kinematic
viscosity (the “inner fluid”) against instantateous changes in the shear rate,
$S=S(t)$.

These results now quantitatively explain the previously only empirically
observed non-Newtonian viscosity due to structure build-up and -breakdown in
colloidal systems under different shear conditions. They are fully in line
with thermodynamical non-equilibrium character of colloidal systems: the
result of chaotic dispersion processes (which are empirically known not to be
easy to control and having quite often unexpected results which require
intensive studies of the process parameters until finally in mass production
reproducible results can be achieved)

### Coda

At least in the present moment we have no detailed knowledge of the inner
relaxation behavior of our test fluid. Therefore we have to apply here the
most simple ansatz – linear or “first order” relaxation (“nudging”). We
augment (17, 18) as follows,

$\displaystyle\frac{d\nu}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}\left(\tilde{\nu}-\nu\right).$ (30)

This equation describes the inner microscopic-dynamic relaxation of
$\nu=\nu(t)$ wherein $\tilde{\nu}$ is the steady-state value of the kinematic
viscosity and known from above as the Ostwald-de Waehle relation,

$\tilde{\nu}=\nu_{0}+\gamma\times({S}/{S_{*}})^{m}\,.$ (31)

Thus finally we have

$\displaystyle\frac{d\nu}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}\left[\nu_{0}+\gamma\times\left(\frac{S}{S_{*}}\right)^{m}-\nu\right].$
(32)

This is our first attempt to describe the effects of changes in nanoscopic
structures of colloidal systems with respect to viscosity.

The value of $\tau$ can be measured in principle but is not yet available. Its
value is important because its relation to intrinsic turbulent time scales
matter. It is also not clear whether the linear ansatz is the correct one
among the many other relaxation models. Further, although (17, 18, 32) form a
closed dynamical problem, (17,18) are possibly not sufficiently representative
for the state of strangled eddy viscosity [[, they have been derived for
inviscid flows,]]baumert2013. These questions deserve specific experimental
studies and further theoretical efforts.

Acknowledgements. HZB thanks Alexander J. Babchin in Tel Aviv for dropping
hints on early fruitful efforts of Yakov Frenkel (1894 – 1952) towards a
kinetic theory of liquids and certain generalized phase transitions.

## References

  * [Bailey et al.(2014)Bailey, Vallikivi, Hultmark, and Smits] Bailey, S. C. C., M. Vallikivi, M. Hultmark, and A. J. Smits, 2014: Estimating the value of von Karman’s constant in turbulent pipe flow. J. Fluid Mech., 749, 79 – 98, doi:dx.doi.org/10.1017/jfm.2014.208.
  * [Baumert(2005)] Baumert, H. Z.: 2005, A novel two-equation turbulence closure for high Reynolds numbers. Part B: Spatially non-uniform conditions. Marine Turbulence: Theories, Observations and Models, H. Z. Baumert, J. H. Simpson, and J. Sündermann, eds., Cambridge University Press, Chapter 4, 31 – 43.
  * [Baumert(2013)] Baumert, H. Z., 2013: Universal equations and constants of turbulent motion. Physica Scripta, T155, 014001 (12pp), doi:10.1088/0031-8949/2013/T155/014001.
  * [Frenkel(1946)] Frenkel, Y., 1946: Kinetic Theory of Liquids. Clarendon Press, 1st engl. ed., 488 pp., russian ed. 1943, Kazan.
  * [Herrmann(1990)] Herrmann, H. J.: 1990, Space-filling bearings. Correlation and Connectivity, H. E. Stanley and N. Ostrowsky, eds., Kluwer Academic Publ., 108 – 120.
  * [Herrmann et al.(1990)Herrmann, Mantica, and Bessis] Herrmann, H. J., G. Mantica, and D. Bessis, 1990: Space-filling bearings. Phys. Rev. Letters, 65, 24 December, 3,223 – 3,226.
  * [Manna and Herrmann(1991)] Manna, S. S. and H. J. Herrmann, 1991: Precise determination of the fractal dimensions of apollonian packing and space-filling bearings. Journal of Physics A, 24.
  * [Pope(2000)] Pope, S. B., 2000: Turbulent flows. Cambridge U. Press, 771 pp.
  * [Pöppe(2004)] Pöppe, C., 2004: Das Getriebe des Teufels. Spektrum d. Wiss., September, 104 – 109.
  * [Wessling(1991)] Wessling, B., 1991: Dispersion hypothesis and non-equilibrium thermodynamics: key elements for a materials science of conductive polymers. A key to understanding polymer blends or other multiphase polymer systems. Synthetic metals, 45, 119 – 149.
  * [Wessling(1993)] — 1993: Dissipative structure formation in collodial systems. Advanced Materials, 5, 300 – 305.
  * [Wessling(1995)] — 1995: Critical shear rate – the instability reason for the creation of dissipative structures in polymers. Z. Phys. Chem., 191, 119 – 135\.