# Permanents through probability distributions

Mobolaji Williams<EMAIL_ADDRESS>Jellyfish, Boston, MA 02111

(June 22, 2021)

###### Abstract

We show that the permanent of a matrix can be written as the expectation value
of a function of random variables each with zero mean and unit variance. This
result is used to show that Glynn’s theorem and a simplified MacMahon theorem
extend from a common probabilistic interpretation of the permanent. Combining
the methods in these two proofs, we prove a new result that relates the
permanent of a matrix to the expectation value of a product of hyperbolic
trigonometric functions, or, equivalently, the partition function of a spin
system. We conclude by discussing how the main theorem can be generalized and
how the techniques used to prove it can be applied to more general problems in
combinatorics.

Permanent, Probability, Glynn’s Theorem, MacMahon Master Theorem

##  I Introduction

The permanent is a seemingly odd function in linear algebra. For a matrix $A$
with row $i$ and column $j$ element $a_{i,j}$, it is defined as

$\operatorname{perm}(A)=\sum_{\sigma\in S_{n}}\prod_{i=1}^{n}a_{i,\sigma(i)},$
(1)

where $\sigma$ is an element of the symmetric group $S_{n}$. The permanent
looks similar to the determinant, but it is missing the “sign of the
permutation” factor $\text{sgn}\,(\sigma)$ which enters into a determinant and
provides the latter with much of its computational simplicity. For example,
while algorithms exist to compute the determinant in polynomial time,
calculations (admittedly naive ones) of the permanent require $O(n!)$ time.
Consequently, there is a literature around permanents that seeks alternative
ways of expressing the permanent for the goal of finding either the new ways
to compute it efficiently or new ways to understand it.

Glynn’s theorem [2] falls in the “computational efficiency” category and
provides a way to compute the permanent in $O(2^{n-1}n^{2})$ time. The
MacMahon Master theorem [5] falls in the “new ways of understanding” category
and relates permanents to their linear algebra cousin determinants. And though
both Glynn’s theorem and the MacMahon theorem are related to permanents, they
are derived from quite different starting points.

In this work, we use a probabilistic definition of the permanent to show that
the two results extend from the same foundation, and we then use this
definition to derive a new result combining aspects of both theorems.

We start with the main theorem we seek to prove:

######  Theorem 1.

Let $p_{X}:\Omega_{X}\to\mathbb{R}$ be a probability distribution defined over
the domain $\Omega_{X}$ with zero mean and unit variance. Let $A$ be an
$n\times n$ matrix with elements $a_{i,j}$. Then the permanent of $A$ is

$\operatorname{perm}(A)=\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\prod_{i=1}^{n}p_{X}(x_{i})\,x_{i}\sum_{j=1}^{n}a_{i,j}x_{j},$
(2)

where $\Omega^{n}_{X}=\Omega_{X}\otimes\cdots\otimes\Omega_{X}$ is the
$n$-factor product over the single-variable domain of integration. In
condensed notation, we can write Eq.(10) as the expectation value

$\operatorname{perm}(A)=\left\langle\prod_{i=1}^{n}x_{i}\sum_{j=1}^{n}a_{i,j}x_{j}\right\rangle_{x_{i}\sim
p_{X}},$ (3)

where the average is over $\\{x_{i}\\}$, a set of independent identically
distributed random variables each of which is drawn from $p_{X}$.

###### Proof.

Say we have a single-variable function $p_{X}$, defined over the domain
$\Omega_{X}\subset\mathbb{R}$, which has the following two properties:

$\int_{\Omega_{X}}dx\,p_{X}(x)x=0,\qquad\int_{\Omega_{X}}dx\,p_{X}(x)x^{2}=1.$
(4)

Interpreting $x$ as a random variable, the conditions in Eq.(4) are equivalent
to requiring $p_{X}$ to be a probability distribution with zero mean and unit
variance. We use this interpretation for simplicity in defining the properties
of $p_{X}$. Two such functions are
$p_{X}(x)=\frac{1}{2}\delta(x-1)+\frac{1}{2}\delta(x+1)$ and
$p_{X}(x)=e^{-x^{2}/2}/\sqrt{2\pi}$.

Next we define a quantity called the permutation tensor:

$\Delta_{\textbf{m}}=\begin{dcases}1&\text{if $\,\textbf{m}$ is a permutation
of $(1,2,\ldots,n)$}\\\ 0&\text{otherwise}\end{dcases},$ (5)

where $\textbf{m}=(m_{1},m_{2},\ldots,m_{n})$. That is, $\Delta_{\textbf{m}}$
is a function of $N$ integers $\textbf{m}=(m_{1},m_{2},\ldots,m_{n})$, and it
is non-zero and equal to 1 if and only if the integers are a permutation of
the numbers from $1$ to $n$. The permanent of $A$ Eq.(1) can be written in
terms of the permutation tensor as

$\operatorname{perm}(A)=\sum_{\\{m_{j}\\}}\Delta_{\textbf{m}}\prod_{j=1}^{n}a_{j,m_{j}},$
(6)

where the summation is defined as

$\sum_{\\{m_{j}\\}}\equiv\prod_{i=1}^{n}\sum_{m_{i}=1}^{n}.$ (7)

There are many possible forms for $\Delta_{\textbf{m}}$, but we will use
Eq.(4) to write is as the multiple integral expression

$\Delta_{\textbf{m}}=\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\prod_{i=1}^{n}p_{X}(x_{i})x_{i}x_{m_{i}},$
(8)

where $d^{n}\textbf{x}=dx_{1}dx_{2}\cdots dx_{n}$. Substituting Eq.(8) into
Eq.(6), we obtain

$\displaystyle\operatorname{perm}(A)$
$\displaystyle=\sum_{\\{m_{j}\\}}\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\left[\prod_{i=1}^{n}p_{X}(x_{i})x_{i}x_{m_{i}}\right]\left[\prod_{j=1}^{n}a_{j,m_{j}}\right]$
$\displaystyle=\left[\prod_{i=1}^{n}\sum_{m_{i}=1}^{n}\right]\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\prod_{i=1}^{n}p_{X}(x_{i})x_{i}a_{i,m_{i}}x_{m_{i}}$
$\displaystyle=\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\prod_{i=1}^{n}p_{X}(x_{i})x_{i}\sum_{m_{i}=1}^{n}a_{i,m_{i}}x_{m_{i}},$
(9)

and the theorem is proved.

∎

As a sanity check, we can use Eq.(10) to affirm a basic result in permanent
calculations. We know that if $A=\mathds{1}$, i.e., $A$ is a matrix entirely
composed of $1$s, then the permanent is $n!$. Checking this result, we have

$\displaystyle\operatorname{perm}(\mathds{1})$
$\displaystyle=\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\prod_{i=1}^{n}p_{X}(x_{i})\,x_{i}\sum_{j=1}^{n}x_{j}$
$\displaystyle=\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\left[\prod_{i=1}^{n}p_{X}(x_{i})\,x_{i}\right]\Big{(}\sum_{j=1}^{n}x_{j}\Big{)}^{n}$
$\displaystyle=\sum_{k_{1}+\cdots+k_{n}=n}n!\prod_{i=1}^{n}\frac{1}{k_{i}!}\int_{\Omega_{X}}dx_{i}\,p_{X}(x_{i})x_{i}^{k_{i}+1}$
(10)

This final result is only nonzero for the integrands where $x_{i}$ is
quadratic for each $i$. Given that $k_{1}+\cdots+k_{n}=n$, we can only have
quadratic $x_{i}$ for all $i$ if $k_{i}=1$ for all $i$. Thus, we have

$\displaystyle\operatorname{perm}(\mathds{1})$
$\displaystyle=\sum_{k_{1}+\cdots+k_{n}=n}n!\prod_{i=1}^{n}\frac{1}{k_{i}!}\delta(1,k_{i})=n!,$
(11)

as expected.

Identifying $p_{X}$ as a probability distribution is primarily for convenience
to allow us to reference the integrals over the linear and quadratic variable
as a mean and variance, respectively. The probability-based interpretation of
these conditions is valid, but as we explain in the discussion, there is a
more general way to write this condition.

But first we explore how this result relates to well-known theorems involving
permanents. First, with Theorem 1, it becomes straightforward to prove Glynn’s
theorem for computing permanents[2].

######  Corollary 1.1 (Glynn’s Theorem).

Let $A=(a_{i,j})$ be an $n\times n$ matrix. Then we have

$\operatorname{perm}(A)=\frac{1}{2^{n-1}}\left[\sum_{\\{s_{i}\\};s_{1}=+1}\left(\prod_{k=1}^{n}s_{k}\right)\prod_{j=1}^{n}\sum_{i=1}^{n}a_{i,j}s_{j}\right],$
(12)

where the outer sum is over all $2^{n-1}$ vectors
$s=(s_{1},\ldots,s_{n})\in\\{\pm 1\\}^{n}$ with $s_{1}=1$.

###### Proof.

We start with Theorem 1. We use the probability density

$p_{X}(s)=\frac{1}{2}\delta(s+1)+\frac{1}{2}\delta(s-1),$ (13)

where $\delta(x)$ is the Dirac delta function (pg. 160 [3]). Eq.(13) expresses
$X$ as a Bernoulli distribution with values $X=\pm 1$ having equal
probabilities. It is straightforward to check that Eq.(13) yields zero mean
and unit variance. Inserting Eq.(13) into Eq.(10), and noting that

$\int_{\Omega_{X}}dx\,p_{X}(x)f(x)=\frac{1}{2}\sum_{s\in\\{\pm 1\\}}f(s),$
(14)

we find

$\operatorname{perm}(A)=\frac{1}{2^{n}}\sum_{\\{s_{i}\\}}\prod_{i=1}^{n}\,s_{i}\sum_{j=1}^{n}a_{i,j}s_{j},$
(15)

where

$\sum_{\\{s_{i}\\}}=\prod_{i=1}^{n}\sum_{s_{i}\in\\{-1,+1\\}}$ (16)

Eq.(15) is a preliminary result for the permanent expressed as an an
expectation value over $n$ Bernoulli random variables, and it looks very close
to the standard result Eq.(12).

To derive the standard result, we first recognize that Eq.(15) is unchanged if
we take $s_{i}$ to $-s_{i}$ for all $i$. We can see this most clearly by
noting that the factor $s_{i}\sum_{j=1}^{n}a_{i,j}s_{j}$ does not change when
we take $s_{i},s_{j}\to-s_{i},-s_{j}$. Thus for each term in the summation,
there is another term that is identical to it for which all of its $s_{i}$
values are switched from positive to negative or vice versa. We can separate
these two sets of terms by isolating a particular $s_{i}$ and setting this
$s_{i}=+1$ for one set and $s_{i}=-1$ for the other. Without loss of
generality we take this isolated $s_{i}$ to be at $i=1$. We can then write
Eq.(15) as

$\displaystyle\operatorname{perm}(A)$
$\displaystyle=\frac{1}{2^{n}}\sum_{\\{s_{i}\\}}\prod_{i=1}^{n}s_{i}\sum_{j=1}^{n}a_{i,j}s_{j}$
$\displaystyle=\frac{1}{2^{n}}\sum_{\\{s_{i}\\};s_{1}=+1}\prod_{i=1}^{n}s_{i}\sum_{j=1}^{n}a_{i,j}s_{j}+\frac{1}{2^{n-1}}\sum_{\\{s_{i}\\};s_{1}=-1}\prod_{i=1}^{n}s_{i}\sum_{j=1}^{n}a_{i,j}s_{j},$
$\displaystyle=\frac{1}{2^{n-1}}\sum_{\\{s_{i}\\};s_{1}=+1}\prod_{i=1}^{n}s_{i}\sum_{j=1}^{n}a_{i,j}s_{j}.$
(17)

In the second equality we split the summation into the aforementioned two sets
of terms, and in the third equality we set these two sets equal due to the
noted $\\{s_{i}\\}\to-\\{s_{i}\\}$ symmetry. Thus Glynn’s theorm is proved. ∎

The final result in Eq.(17) involves a summation over $2^{n-1}$ terms each of
which involves an $n$-factor product over an $n$-term sum. Thus the
computation time for the permanent under Eq.(17) is $O(2^{n-1}n^{2})$ which is
much better than the $O(n!)$ computation time obtained from a brute force
calculation according to the definition Eq.(1).

The essential advantage of Theorem 1 is that it is written in terms of the
arbitrary probability distribution $p_{X}$ and thus we can make a convenient
choice for this distribution when we seek to prove any theorem. We make such a
choice in our proof of a simplified version of the MacMahon Master theorem.

######  Corollary 1.2 (Simplified MacMahon Master Theorem).

Let $A=(a_{i,j})$ be an $n\times n$ matrix. Let $X$ be the diagonal matrix
$X=\operatorname{diag}(x_{1},x_{2},\ldots,x_{n})$. Then we have

$\operatorname{perm}(A)=\text{coefficient of {\hbox{x_1x_2\cdots x_n}}
}\;\text{ term in }\;\frac{1}{\det(I-XA)},$ (18)

where $I$ is the $n\times n$ identity matrix.

###### Proof.

We again start with Eq.(10). We choose our probability density to be a
Gaussian $p_{X}(x)=e^{-x^{2}/2}/\sqrt{2\pi}$. Eq.(10) then becomes

$\operatorname{perm}(A)=\frac{1}{(2\pi)^{n/2}}\int_{\mathbb{R}^{n}}d^{n}\textbf{x}\,e^{-\textbf{x}^{2}/2}\prod_{i=1}^{n}\,x_{i}\sum_{j=1}^{n}a_{i,j}x_{j},$
(19)

where $\textbf{x}^{2}=x_{1}^{2}+\cdots+x_{n}^{2}$. Using the contour integral
identity

$\prod_{k=1}^{n}\alpha_{k}=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{k=1}^{n}\frac{dz_{k}}{z_{k}^{2}}\,\right]e^{\sum_{k=1}^{n}\alpha_{k}z_{k}},$
(20)

we write the product in Eq.(19) as

$\prod_{i=1}^{n}\,x_{i}\sum_{j=1}^{n}a_{i,j}x_{j}=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{i=1}^{n}\frac{dz_{i}}{z_{i}^{2}}\right]\exp\left(\sum_{i=1}^{n}z_{i}x_{i}\sum_{j=1}^{n}a_{i,j}x_{j}\right).$
(21)

Inserting this form of the product into Eq.(19), we have

$\displaystyle\operatorname{perm}$ $\displaystyle(A)$
$\displaystyle=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{i=1}^{n}\frac{dz_{i}}{z_{i}^{2}}\right]\frac{1}{(2\pi)^{n/2}}\int_{\mathbb{R}^{n}}d^{n}\textbf{x}\,\exp\left(-\frac{1}{2}\sum_{i=1}^{n}x_{i}^{2}+\sum_{i=1}^{n}z_{i}x_{i}\sum_{j=1}^{n}a_{i,j}x_{j}\right)$
$\displaystyle=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{i=1}^{n}\frac{dz_{i}}{z_{i}^{2}}\right]\frac{1}{(2\pi)^{n/2}}\int_{\mathbb{R}^{n}}d^{n}\textbf{x}\,\exp\left(-\frac{1}{2}\sum_{i,j=1}^{N}x_{i}\left(\delta_{i,j}-2z_{i}a_{i,j}\right)x_{j}\right)$
$\displaystyle=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{i=1}^{n}\frac{dz_{i}}{z_{i}^{2}}\right]\frac{1}{\det{(I-2ZA)^{1/2}}}$
(22)

From Cauchy’s integral theorem [1], we know that the contour integral
$\frac{1}{2\pi i}\oint\frac{dz}{z^{2}}p(z)$ picks out the linear-in-$x$ term
of $p(x)$. Thus we see that the permanent of $A$ is the coefficient of the
$x_{1}x_{2}\cdots x_{n}$ term in ${1}/{\det(I-2XA)^{1/2}}$. Since it is only
the linear-in-$z$ term within the contour integral that contributes to the
permanent, we can add a quadratic-in-$Z$ term to the argument of the
determinant in Eq.(22) with no change to the final result:

$\operatorname{perm}(A)=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{i=1}^{n}\frac{dz_{i}}{z_{i}^{2}}\right]\frac{1}{\det{(I-2ZA+ZAZA)^{1/2}}}.$
(23)

With the fact that $\det(Q^{2})^{1/2}=(\det(Q)\det(Q))^{1/2}=\det(Q)$ for an
arbitrary matrix $Q$, we then have

$\operatorname{perm}(A)=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{i=1}^{n}\frac{dz_{i}}{z_{i}^{2}}\right]\frac{1}{\det{(I-ZA)}}.$
(24)

Using the same reasoning as that presented below Eq.(22), we see that Eq.(24)
allows us to conclude that the permanent of $A$ is the coefficient of the
$x_{1}x_{2}\cdots x_{n}$ term in ${1}/{\det(I-XA)}$.

∎

Figure 1: Theorem map: Map showing how Theorem 1 is related to Glynn’s theorem
and the simplified MacMahon Master theorem. The “*” next to Glynn’s theorem
represents the fact that we are referencing Eq.(15). The “Permanent as a
Gaussian Expectation” refers to Eq.(19). The “Product as a Contour Integral”
result refers to Eq.(20). Although the Gaussian integral analog of Glynn’s
theorem, Eq.(19), is known (e.g., pg 27 of [5]), the Bernoulli distribution
analog of the simplified MacMahon Master theorem is not part of the literature
on permanents. Motivated by this gap we derive the necessary analog, terming
it the Glynn-MacMahon theorem given its placement on the map.

With Glynn’s theorem and a simplified MacMahon Master theorem proved from
Theorem 1, we can create a map (Fig. 1) relating the former two theorems to
their common starting point. From the map, we see that starting from Theorem
1, we followed paths that ultimately led us to Glynn’s theorem and the
MacMahon Master theorem by choosing, respectively, a Bernoulli distribution
and a Gaussian distribution as the particular probability distributions from
which to sample $x_{i}$ in Eq.(10). However, the map suggests there is a gap
in this collection of results. While Glynn’s theorem is the Bernoulli
distribution analog of Eq.(19), there appears to be no Bernoulli distribution
analog to the MacMahon Master theorem.

In the next theorem, we fill in this gap by proving what we term the Glynn-
MacMahon Master theorem given the result’s placement in Fig. 1. The theorem
shows that the permanent of a matrix can be related to the expectation value
of a product over hyperbolic trigonometric functions or, equivalently, the
partition function of a spin system.

######  Corollary 1.3 (Glynn-MacMahon Theorem).

Let $A=(a_{i,j})$ be an $n\times n$ matrix with non-zero determinant, let $X$
be the diagonal matrix $X=\operatorname{diag}(x_{1},x_{2},\ldots,x_{n})$, and
let $(\varphi_{1},\ldots,\varphi_{n})$ to be a vector of $n$ normally
distributed random variables with mean $0$ for all components and covariance
matrix $XA$. Then we have

$\operatorname{perm}(A)=\text{coefficient of {\hbox{x_1x_2\cdots
x_n}}}\;\text{ term in
}\;\left\langle\prod_{i=1}^{n}2\cosh\varphi_{i}\right\rangle_{\varphi_{i}\sim{\cal
N}(0,XA)}.$ (25)

###### Proof.

From an intermediate result in the proof of Glynn’s Theorem, we can write the
permanent of $A$ as an expectation value over Bernoulli random variables:

$\operatorname{perm}(A)=\sum_{\\{s_{i}\\}}\prod_{i=1}^{n}\,\frac{1}{2}\sum_{j=1}^{n}s_{i}a_{i,j}s_{j}\,\qquad\text{where}\quad\sum_{\\{s_{i}\\}}=\prod_{i=1}^{n}\sum_{s_{i}\in\\{-1,+1\\}}.$
(26)

Using Eq.(20), we then obtain

$\displaystyle\operatorname{perm}(A)$ $\displaystyle=\frac{1}{(2\pi
i)^{n}}\oint\left[\prod_{k=1}^{n}\frac{dz_{k}}{z_{k}^{2}}\right]\sum_{\\{s_{i}\\}}\exp\left[\frac{1}{2}\sum_{i,j=1}^{n}s_{i}z_{i}a_{i,j}s_{j}\right].$
(27)

Next, we take $z_{i}a_{i,j}$ to be the component in the $i$th row and $j$th
column of $ZA$ where $Z=\operatorname{diag}(z_{1},\ldots,z_{n})$. For
notational simplicity, we define $W\equiv ZA$ where
$Z=\operatorname{diag}(z_{1},\ldots,z_{n})$. Then, using the evaluation of
Gaussian integrals for multivariate functions, we have

$\displaystyle\exp\left[\frac{1}{2}\sum_{i,j=1}^{n}s_{i}W_{i,j}s_{j}\right]$
$\displaystyle=\left(\frac{\det W^{-1}}{(2\pi)^{n}}\right)^{1/2}$
$\displaystyle\times\int
d^{n}\boldsymbol{\varphi}\,\exp\left[-\frac{1}{2}\sum_{i,j=1}^{n}\varphi_{i}(W^{-1})_{i,j}\varphi_{j}+\sum_{i=1}^{n}\varphi_{i}s_{i}\right].$
(28)

From our definition of $W$, we can show that $\det W=z_{1}\cdots z_{n}\det A$,
and therefore this result is valid only if $\det A\neq 0$. Performing the
summation over spins, and noting that the prefactor in Eq.(28) is a
normalization constant, we have

$\sum_{\\{s_{i}\\}}\exp\left[\frac{1}{2}\sum_{i,j=1}^{n}s_{i}W_{i,j}s_{j}\right]=\left\langle\prod_{i=1}^{n}2\cosh\varphi_{i}\right\rangle_{\varphi_{i}\sim{\cal
N}(0,W)},$ (29)

where $\langle{\cal O}(\\{\varphi_{i}\\})\rangle_{\varphi_{i}\sim{\cal
N}(0,W)}$ indicates we are taking the expectation value of the function ${\cal
O}$ over the normally distributed values $(\varphi_{1},\ldots,\varphi_{n})$
with mean $0$ and covariance matrix $W$. From Eq.(27), and the fact that the
contour integrals pick out the linear parts of their corresponding variable,
the theorem is proved.

∎

An intermediate result in the above proof shows how the permanent is connected
to a canonical system in physics. In statistical physics, if we have a set of
$n$ spin variables denoted $\\{s_{i}\\}$, where each spin $i$ in this set can
have the value $s_{i}=+1$ or $s_{i}=-1$, and these spins interact according to
the matrix $J_{i,j}$, then the thermal partition function of the system is
given by [4]

$Z^{\text{spin}}_{n}(\\{\beta
J\\})=\sum_{\\{s_{i}\\}}\exp\left[-\frac{\beta}{2}\sum_{i,j=1}^{n}s_{i}J_{i,j}s_{j}\right].$
(30)

where $\beta$ is inverse temperature. By comparison with Eq.(27), we find that
the permanent is equal to the $x_{1}\cdots x_{n}$ term of the $n$-spin
partition function with interaction matrix $x_{i}a_{i,j}/\beta$, namely:

$\operatorname{perm}(A)=\text{coefficient of {\hbox{x_1x_2\cdots
x_n}}}\;\text{ term in }\;Z^{\text{spin}}_{n}(\\{-XA\\}).$ (31)

##  II Discussion

In this work we proved a general theorem for permanents, and used the result
to prove Glynn’s theorem and a simplified MacMahon Master theorem. Combining
aspects of each proof, we proved a new result relating permanents to the
expectation value of a product of hyperbolic trigonometric functions. The fact
that we can obtain such disparate results from the same theorem suggests that
the starting theorem can serve as a foundational definition for permanents.

In general, although it is convenient to do so, we do not need to interpret
the $p_{X}$ that appears in Theorem 1 as a probability density. Moreover, it
is possible to write the permanent as any multiple integral where the single-
variable integrals have two key properties: Their results are zero when their
integrands are linear in a function and unity when their integrands are
quadratic in that function. That is, for functions
$f,g:\Omega_{X}\to\mathbb{R}$ and the set of conditions

$\int_{\Omega_{X}}dx\,g(x)f(x)=0\qquad\int_{\Omega_{X}}dx\,g(x)f(x)^{2}=1,$
(32)

we could have just have easily proved

$\operatorname{perm}(A)=\int_{\Omega^{n}_{X}}d^{n}\textbf{x}\,\prod_{i=1}^{n}g(x_{i})\,f(x_{i})\sum_{j=1}^{n}a_{i,j}f(x_{j}),$
(33)

using an approach identical to that for Theorem 1. This more general case
reveals that there is not a coherent probabilistic interpretation for all such
expressions. For example, the permanent of $A$, written as a multiple integral
over a product of sine functions, can be expressed as

$\operatorname{perm}(A)=\frac{1}{\pi^{n}}\int^{2\pi}_{0}d\theta_{1}\ldots
d\theta_{N}\,\prod_{i=1}^{N}\sin(\theta_{i})\sum_{j=1}^{n}a_{i,j}\sin(\theta_{j}),$
(34)

However, in introducing $1/\pi$ to ensure that the single-variable sine
integrals satisfy Eq.(32), we have rendered them incapable of satisfying the
normalization requirement that would allow us to interpret $\theta$ as a
random variable with uniform probability density around the unit circle. Thus
the probabilistic interpretation of $p_{X}$ is useful but not necessary in
Theorem 1 and in its generalizations.

One important aspect of the proof of Theorem 1 is the introduction of the
permutation tensor $\Delta_{\textbf{m}}$, Eq.(8). With Eq.(33) we could
generalize this tensor and define it in terms of products over functions
rather than variables, and similar to how it was used in the main proof, the
tensor could be used to define more general sums, such as those that appear in
combinatorial optimization problems, where the state space consists of
permutations of elements.

##  III Acknowledgements

The author declares no conflicts of interest.

## References

  * [1] Lars V Ahlfors.  Complex analysis: an introduction to the theory of analytic functions of one complex variable.  New York, London, 177, 1953.
  * [2] David G. Glynn.  The permanent of a square matrix.  European Journal of Combinatorics, 31(7):1887–1891, oct 2010.
  * [3] Sadri Hassani.  Mathematical physics: a modern introduction to its foundations.  Springer Science & Business Media, 2013.
  * [4] Hidetoshi Nishimori.  Statistical Physics of Spin Glasses and Information Processing.  Oxford University Press, jul 2001.
  * [5] Jerome K. Percus.  Combinatorial Methods.  Springer New York, 1971.


# Heat Transfer Enhancement in Rotating Rayleigh-Bénard Convection at $Ra\geq
10^{10}$

Mohammad Anas  Department of Mechanical Engineering, Indian Institute of
Technology, Kanpur 208016, India  Pranav Joshi  Department of Mechanical
Engineering, Indian Institute of Technology, Kanpur 208016, India

###### Abstract

Rotation, which stabilizes flow, can enhance the heat transfer in Rayleigh-
Bénard convection (RBC) due to Ekman pumping. In this Letter, we study this
heat transfer enhancement by conducting direct numerical simulations of
rotating RBC for a wide range of parameters: Rayleigh number ($Ra$) from
$2\times 10^{4}$ to $2\times 10^{10}$, Prandtl number ($Pr$) from $1$ to
$1000$, and Taylor number ($Ta$) from 0 to $2\times 10^{12}$. We show that for
each $Ra$, there exists a critical Prandtl number ($Pr_{cr}$) below which no
heat transfer enhancement occurs at any rotation rate, as well as an optimal
Prandtl number ($Pr_{opt}$) at which maximum heat transfer enhancement occurs
at an optimal rotation rate. We demonstrate a significant heat transfer
enhancement for $Ra$ up to $2\times 10^{10}$ and predict that heat transfer
enhancement due to rotation not only would occur at $Ra>10^{10}$ but would
also become more pronounced.

Thermal convection manifests across a broad spectrum of scales, ranging from
the cooling of electronic devices to geophysical and astrophysical flows. To
study the key characteristics of thermal convection, such as heat and momentum
transfer, a simplified model known as Rayleigh-Bénard convection (RBC) is
commonly employed. Within the framework of RBC, fluid motion occurs between a
hot plate (at the bottom) and a cold plate (at the top) as a consequence of
the thermal buoyancy force Chandrasekhar (1961); Drazin and Reid (2004). The
behavior of RBC system is primarily determined by two dimensionless
parameters: the Rayleigh number ($Ra=g\beta\Delta H^{3}/\nu\kappa$), which
represents the strength of the buoyancy forcing, and the Prandtl number
($Pr=\nu/\kappa$), which represents the ratio of the momentum diffusivity to
thermal diffusivity; here, $g$ is the acceleration due to gravity, $\beta$ is
the thermal expansion coefficient, $\Delta$ is the temperature difference
between the hot and cold plates, $H$ is the separation between the plates,
$\nu$ is the kinematic viscosity, and $\kappa$ is the thermal diffusivity. In
RBC research, the primary focus is on the heat transfer and flow intensity,
quantified by the Nusselt number ($Nu$) and the Reynolds number ($Re$),
respectively, and their scaling with $Ra$ and $Pr$ Shraiman and Siggia (1990);
Grossmann and Lohse (2000); Niemela  _et al._ (2000); Chillà and Schumacher
(2012).

Incorporating background rotation into a thermal convection system makes it
relevant to various flows occurring within the Earth’s outer core Glatzmaier
and Roberts (1995), oceans Marshall and Schott (1999), gaseous planets like
Jupiter Ingersoll (1990); Heimpel  _et al._ (2005), and solar interiors
Spiegel (1971). An RBC system with background rotation parallel to the
gravitational force is known as rotating RBC. Rotation, which introduces the
Coriolis force in a convective system, significantly affects its heat and
momentum transfer characteristics Rossby (1969). The strength of the Coriolis
force relative to the viscous force is represented by the dimensionless
parameter Taylor number ($Ta=4\Omega^{2}H^{4}/\nu^{2}$); here, $\Omega$
represents the rotation rate of the system. To characterize the relative
strength of convection over rotation, convective Rossby number
$Ro=\sqrt{Ra/TaPr}$ is commonly used. When $Ro\gg 1$, buoyancy force dominates
over the Coriolis force, and the properties of a rotating RBC system closely
resemble those of non-rotating RBC King  _et al._ (2008, 2012); Stevens  _et
al._ (2009a). On the other hand, when $Ro\ll 1$, rotation becomes dominant and
rotating convection exhibits similarities to geostrophic flow, which is
characterized by a force balance between pressure gradient and the Coriolis
force King  _et al._ (2008); Kunnen (2021); Ecke and Shishkina (2023).

Typically, rotation suppresses the intensity of flow; consequently, the heat
transfer is also reduced Niemela  _et al._ (2010); Ecke and Niemela (2014);
Stellmach  _et al._ (2014). However, it is observed, both in experiments and
simulations, that for a certain range of parameters, rotation can also enhance
the heat transfer Rossby (1969); Zhong  _et al._ (1993); Liu and Ecke (1997);
Kunnen  _et al._ (2006); Stevens  _et al._ (2013); Wei  _et al._ (2015); Zhong
_et al._ (2009); Stevens  _et al._ (2009a); Yang  _et al._ (2020); Vorobieff
and Ecke (2002); Stevens  _et al._ (2009b, 2010a, 2010b, 2011); Joshi  _et
al._ (2017); Weiss  _et al._ (2010, 2016); Chong  _et al._ (2017). This heat
transfer enhancement in rotating RBC as compared to the non-rotating case is
ascribed to Ekman pumping Julien  _et al._ (1996); Kunnen  _et al._ (2006);
Zhong  _et al._ (2009). There are, however, many studies which report no or
negligible enhancement in heat transfer due to rotation at low $Pr$ or high
$Ra$, e.g., Pfotenhauer  _et al._ (1984) for $Pr\approx 0.5$ at low $Ra$, and
Ecke and Niemela (2014) for $Pr=0.7$, $4\times 10^{9}<Ra<4\times 10^{11}$, and
$4\times 10^{9}<Ta<1\times 10^{14}$. In particular, several studies have
reported no or negligible heat transfer enhancement for $Ra\gtrsim 10^{10}$:
Niemela  _et al._ (2010) (for $Pr=0.7-5.9$), Stellmach  _et al._ (2014) (for
$Pr\approx 1-7$), Kunnen  _et al._ (2016) (for $Pr=1$), and Hartmann  _et al._
(2023) (for $Pr=4.38$ and $6.4$). The findings of these studies have led to a
common expectation that there would be no or negligible heat transfer
enhancement due to rotation for $Ra\gtrsim 10^{10}$ Kunnen (2021); Ecke and
Shishkina (2023).

In this Letter, however, we demonstrate a clear evidence of significant heat
transfer enhancement due to rotation even at $Ra\geq 10^{10}$. We show that
for each $Ra$ (at least within the range of $Ra=2\times 10^{4}-2\times
10^{10}$ explored in the present work), significant heat transfer enhancement
will occur only if the Prandtl number is greater than a critical Prandtl
number $Pr_{cr}$ that increases with increasing $Ra$. For $Ra=10^{10}$,
$Pr_{cr}\approx 10$, which is greater than $Pr$ used in most previous studies
that failed to observe any or significant enhancement in the heat transfer for
$Ra\gtrsim 10^{10}$. In this work, we also provide a precise definition of the
optimal Prandtl number $Pr_{opt}$ for obtaining the maximum heat transfer
enhancement at a given $Ra$ and show that $Pr_{opt}$ also increases with $Ra$.
Furthermore, we emphasize that $Ta$ serves as a better parameter compared to
$1/Ro$ in representing the optimal rotation rate at which the maximum
enhancement occurs. Finally, we show that the maximum enhancement in heat
transfer at $Pr_{opt}$ exhibits a monotonic increase with $Ra$. Importantly,
the present findings predict that heat transfer enhancement due to rotation
not only would occur but also would become more pronounced at $Ra>10^{10}$.

For this study, we perform a large number (more than $400$) of direct
numerical simulations (DNS) for a wide range of parameters: $Ra=2\times
10^{4}-2\times 10^{10}$, $Pr=1-1000$, and $Ta=0-2\times 10^{12}$. We perform
simulations in a horizontally periodic rectangular domain of size $L\times
L\times H$ ($L\times L$ in the horizontal directions and $H$ along the
vertical) employing isothermal and no-slip (and impenetrable) boundary
conditions at the top (cold) and bottom (hot) plates. For the simulations of
non-rotating RBC ($Ta=0$) at moderate Ra, we use large aspect ratio
($\Gamma=L/H$) to avoid the effect of confinement on the Nusselt number Huang
_et al._ (2013): $\Gamma=8$ for $Ra=2\times 10^{4}-10^{6}$ and $\Gamma=4$ for
$Ra=10^{7}-10^{8}$. Considering the high computational cost at large $\Gamma$
for high $Ra$, we use $\Gamma=1$ for $Ra=5\times 10^{8}-2.3\times 10^{9}$ and
$\Gamma=0.5$ for $Ra=10^{10}$. We use $Nu\approx 0.12Ra^{0.30}$ (which we
obtain by fitting the $Nu$ data for $Ra=10^{6}-10^{10}$ at $Pr=100$) to
estimate $Nu\approx 148$ for $Ra=2\times 10^{10}$ and $Pr=100$. Since the
horizontal length scale of the flow in rotating convection, $\ell_{c}$,
decreases with $Ta$ as $\ell_{c}=2.4Ta^{-1/6}H$ Chandrasekhar (1961), we use
relatively lower aspect ratios (half or one-fourth of those for the
corresponding non-rotating RBC cases) for some simulations of rotating RBC. In
all simulations, we ensure $\ell_{c}/L\lesssim 1/8$ to mitigate the
nonphysical effects arising due to the periodic boundaries in the horizontal
directions Kunnen  _et al._ (2016); Julien  _et al._ (1996). For more details
about the simulations and the solver used in this study, please refer to the
Supplementary Material SM: .

In Fig. 1, we show the variation of the normalized Nusselt number $Nu/Nu_{0}$
($Nu_{0}$ is the Nusselt number for the non-rotating case) with the Taylor
number, $Ta$, and the inverse of the Rossby number, $1/Ro$, for $Pr=1-1000$ at
$Ra=[10^{7},10^{8},10^{9},10^{10}]$. We observe that when $Pr$ is not too
small, $Nu/Nu_{0}$ first increases and then decreases as the rotation rate is
increased. For each $Ra$, the maximum enhancement in the heat transfer as
compared to the non-rotating case occurs in a certain range of $Pr$ and
rotation rate. This maximum enhancement increases with increasing $Ra$ and can
reach up to approximately $25\%$, $40\%$, and $55\%$ for $Ra=10^{7}$,
$Ra=10^{8}$, and $Ra=10^{9}$, respectively. Note that we observe a significant
heat transfer enhancement (more than $40\%$) even at $Ra=10^{10}$, which will
be discussed later in greater detail. We also observe that the Taylor number
$Ta$ serves as a better parameter than $1/Ro$ in representing the optimal
rotation rate at which the maximum enhancement occurs. Unlike the optimal
rotation rate represented in terms of the inverse of Rossby number,
($1/Ro_{opt}$), the optimal Taylor number, $Ta_{opt}$, is nearly independent
of $Pr$ when significant enhancement is observed at a given $Ra$. Most earlier
studies used $1/Ro_{opt}$ to represent the optimal rotation rate but found
$1/Ro_{opt}$ to be strongly dependent on $Pr$ Zhong  _et al._ (2009); Stevens
_et al._ (2010a); Yang  _et al._ (2020). Since the heat transfer enhancement
due to rotation is largely controlled by the dynamics of the Ekman boundary
layer Stellmach  _et al._ (2014), the thickness of which depends only on the
Taylor number (which represents the ratio of Coriolis to viscous forces), $Ta$
can be expected to represent better the heat transfer enhancement than $1/Ro$
(which represents the ratio of Coriolis to buoyancy forces) at moderate and
high rotation rates. This finding is also in line with the hypothesis of King
_et al._ (2008) that the boundary layer controls the rotation-dominated regime
in rotating RBC, rather than the balance between the buoyancy and Coriolis
forces. Nonetheless, the beginning of the rotation-affected regime, i.e., the
rotation rate at which $Nu/Nu_{0}$ deviates from $1$, is better represented by
$1/Ro$ than by $Ta$, as seen from our results. Specifically, for $Ra=10^{7}$
and $10^{8}$, the rotation-affected regime begins at $1/Ro\approx 0.2$, while
for $Ra=10^{9}$, it begins at $1/Ro\approx 0.4$. These values are consistent
with the findings of Stevens  _et al._ (2009b).

In Fig. 2, we show the variation of $Ta_{opt}$ with $Ra$ for various Prandtl
numbers. As discussed earlier, we observe that at a given $Ra$, the optimal
Taylor number does not vary significantly with $Pr$. Also, $Ta_{opt}$ follows
a power law close to $Ta_{opt}\sim Ra^{1.5}$ up to a certain $Ra$ and this
limiting $Ra$ for the power law seems to increase with increasing $Pr$. Note
that the Taylor number at which convection ceases completely ($Ta_{cs}$) also
follows the scaling $Ta_{cs}\sim Ra^{1.5}$ for $Ta\gg 1$ Chandrasekhar (1961).

In Fig. 3, we show the variation of normalized maximum Nusselt number,
$Nu_{max}/Nu_{0}$, with $Pr$ for $Ra=2\times 10^{4}-2\times 10^{10}$. At any
given $Ra$ and $Pr$, $Nu_{max}$, by definition, corresponds to the Nusselt
number at the optimal Taylor number for that $Ra$ and $Pr$. The maximum heat
transfer enhancement represented by $Nu_{max}/Nu_{0}$ for each $Ra$ increases
with $Pr$ up to a certain $Pr$ and then decreases. Note that for each $Ra$
there exists a Prandtl number (obtained by extrapolating the data for each
$Ra$ to $Nu_{max}/Nu_{0}=1$) below which there will be no (or negligible) heat
transfer enhancement at any rotation rate. We call this $Pr$ the critical
Prandtl number $Pr_{cr}$. In rotating convection, columnar vortical structures
are known to play a critical role in the heat transfer by transporting the
temperature anomaly from one wall to the other Zhong  _et al._ (2009); Julien
_et al._ (1996). However, for $Pr<Pr_{cr}$, it is likely that the lateral
diffusion of the heat/temperature anomaly away from the vortex columns
restricts their ability to transport heat between the top and bottom walls
Stevens  _et al._ (2010a); Hartmann  _et al._ (2023). As $Pr$ increases, this
effect is expected to weaken, and so the heat transfer increases. However, the
heat transfer enhancement decreases again at very large $Pr$.

At any $Ra$, we define the Prandtl number at which $Nu_{max}/Nu_{0}$ reaches
its maximum as the optimal Prandtl number $Pr_{opt}$ for that $Ra$. Some
studies (e.g., Stevens  _et al._ (2010a)), comparing the heat transfer at a
constant $Ro$, have proposed that at large $Pr$, the Ekman boundary layer
($\delta_{u}$) is much thicker than the thermal boundary layer
($\delta_{\theta}$). Consequently, the columnar vortices do not reach the
thermal boundary layer and the fluid entering them is not as hot (or as cold),
leading to a decrease in the heat transfer enhancement at large $Pr$. However,
we observe that the maximum enhancement (which occurs at $Ta_{opt}$) decreases
at large Prandtl number even though
$1.3\lesssim\delta_{u}/\delta_{\theta}\lesssim 1.5$ for $Pr>Pr_{opt}$ (see
Supplementary Material SM: ). We hypothesize that the higher viscous damping
of the flow at very large $Pr$ weakens the vertical advection of heat by the
columnar structures, resulting in maximum heat transfer at an intermediate
(optimal) $Pr$.

Here, we make two important observations. First, we observe a significant heat
transfer enhancement even at $Ra\geq 10^{10}$: approximately $40\%$ at
$Ra=2\times 10^{10}$ for $Pr=100$, and the trends indicate an even higher heat
transfer enhancement for higher $Pr$. This finding challenges a common
expectation that there will be no (or negligible) heat transfer enhancement
due to rotation at $Ra\geq 10^{10}$ Ecke and Shishkina (2023); Kunnen (2021).
The present results predict that heat transfer enhancement due to rotation is
possible even for $Ra>10^{10}$ provided $Pr>Pr_{cr}$. Secondly, we observe a
significant heat transfer enhancement ($>10\%$) even at $Ra=2\times 10^{4}$
(using $\Gamma=8$), in agreement with Rossby (1969)’s experimental results for
$\Gamma\gtrsim 6$. No earlier studies (except Rossby (1969)) have reported any
heat transfer enhancement at such low $Ra$ ($Ra\lesssim 10^{6}$), perhaps due
to the significant effects of lateral confinement on the Nusselt number at low
$\Gamma$ ($\Gamma\approx 1$).

In Fig. 4, we show the variation of $Pr_{cr}$ and $Pr_{opt}$ with $Ra$.
Interestingly, both $Pr_{cr}$ and $Pr_{opt}$ increase monotonically with $Ra$
and approximately follow power-laws: $Pr_{cr}\approx 1.46\times
10^{-3}Ra^{0.375}$ and $Pr_{opt}\approx 8.18\times 10^{-3}Ra^{0.504}$.
Considering the high computational cost at large $Pr$ and large $Ra$, we do
not perform simulations at $Pr>200$ to find $Pr_{opt}$ for $Ra=10^{10}$, which
is estimated to be $Pr_{opt}\approx 900$ by the above power-law fit. As
Rayleigh number increases, the turbulent diffusion of heat is also expected to
become stronger. At low $Pr$, this higher turbulent diffusion will combine
with the large molecular thermal diffusivity to further increase the lateral
diffusion of heat in the bulk, and hence, will decrease the ability of the
vortex columns to transport heat. Thus, a correspondingly larger $Pr$ may be
necessary to counter this effect of the enhanced turbulent thermal diffusivity
to register any enhancement in the heat transfer, i.e., $Pr_{cr}$ will
increase with increasing $Ra$. On the other hand, as the buoyancy forcing
increases with increasing $Ra$, the viscous damping of the flow at large $Pr$
will become weaker and the heat transfer enhancement can be sustained until
larger Prandtl numbers, i.e. $Pr_{opt}$ also increases as $Ra$ is increased.
We again note that $Pr_{cr}\approx 10$ for $Ra=10^{10}$, and it further
increases with $Ra$. Thus, most earlier studies that used $Pr<10$ Niemela  _et
al._ (2010); Stellmach  _et al._ (2014); Kunnen  _et al._ (2016); Ecke and
Niemela (2014); Hartmann  _et al._ (2023) failed to observe any (or
significant) enhancement in the heat transfer due to rotation for $Ra\gtrsim
10^{10}$.

In Fig. 5, we show the variation of $Nu_{max}(Pr_{opt})/Nu_{0}$ with $Ra$.
Here, $Nu_{max}(Pr_{opt})$ is $Nu_{max}$ at $Pr=Pr_{opt}$. Again, we observe a
clear monotonic increase of $Nu_{max}(Pr_{opt})/Nu_{0}$ with $Ra$, and the
trend can be fitted by a power-law $Nu_{max}(Pr_{opt})/Nu_{0}\approx
0.62Ra^{0.044}$. This relationship predicts $Nu_{max}/Nu_{0}\approx 1.7$
(i.e., $70\%$ enhancement in heat transfer) for $Ra=10^{10}$ at
$Pr_{opt}\approx 900$, and even higher enhancement at $Ra>10^{10}$ at higher
$Pr$. Thus, the present trends predict that the maximum heat transfer
enhancement increases with $Ra$, even for $Ra>10^{10}$, provided the Prandtl
number is also increased commensurately. As discussed earlier, at a given
$Ra$, the enhancement increases with $Pr$ up to $Pr_{opt}$ beyond which the
viscous damping of the flow restricts the vertical advection of heat. However,
as $Ra$ and hence buoyancy forcing is increased, the enhancement due to
rotation can increase up to larger $Pr$ before the viscous damping effect
becomes significant. Thus, the maximum enhancement can be expected to increase
with $Ra$.

Note that the effect of finite aspect ratio on $Nu$ in non-rotating
simulations Huang  _et al._ (2013) may lead to some uncertainty in the values
of $Pr_{cr}$ and $Nu_{max}/Nu_{0}$ for $Ra\geq 5\times 10^{8}$. However, this
effect is expected to not alter any of the major findings of this study.

Figure 1: Variation of the normalized Nusselt number $Nu/Nu_{0}$ with Taylor
number $Ta$ (left) and $1/Ro$ (right) for various $Pr$ at (a) $Ra=10^{7}$, (b)
$Ra=10^{8}$, (c) $Ra=10^{9}$, and (d) $Ra=10^{10}$. Solid lines are used to
connect data points, aiding visual interpretation. Figure 2: Variation of the
optimal Taylor number $Ta_{opt}$ with $Ra$ for various $Pr$. Dashed line
represents $Ta_{opt}=0.02Ta_{cs}$, where $Ta_{cs}=(Ra/8.7)^{1.5}$. Figure 3:
Variation of the normalized maximum Nusselt number $Nu_{max}/Nu_{0}$ with $Pr$
for various $Ra$. Solid lines are used to connect data points, aiding visual
interpretation. Figure 4: Variation of the critical Prandtl number $Pr_{cr}$
and the optimal Prandtl number $Pr_{opt}$ with $Ra$. Dot-dashed and dashed
lines represent power-law fits $Pr_{cr}\approx 1.46\times 10^{-3}Ra^{0.375}$
and $Pr_{opt}\approx 8.18\times 10^{-3}Ra^{0.504}$, respectively. Figure 5:
Variation of $Nu_{max}(Pr_{opt})/Nu_{0}$ with $Ra$. Dashed line represents
power-law fit $Nu_{max}(Pr_{opt})/Nu_{0}\approx 0.62Ra^{0.044}$.

The present results show that as $Ra$ is increased in rotating RBC, $Pr_{cr}$,
$Pr_{opt}$, as well as the maximum heat transfer enhancement increase. In
particular, not only enhancement is certainly possible for $Ra>10^{10}$, but
is also expected to be higher than that for lower $Ra$ at the optimal Prandtl
number. However, we do not know up to what $Ra$ these trends will persist.
Simulations and experiments at significantly higher $Ra$ and $Pr$ than
currently feasible may be required to answer this question.

In closing, we briefly discuss the implications of the present results for
identifying the transition from a buoyancy-dominated to a rotation-dominated
regime in rotating RBC. Some previous studies have used the heat transfer
enhancement (the rotation rate for maximum enhancement or for $Nu/Nu_{0}=1$)
by Ekman pumping to define this transition (see Ecke and Shishkina (2023)).
However, according to our findings, these criteria fail when $Pr<Pr_{cr}$, for
which no enhancement occurs, and unconditionally for free-slip boundary
conditions for which there is no Ekman pumping and consequently no enhancement
in the heat transfer. King  _et al._ (2008, 2012) suggested the ratio of the
Ekman and thermal boundary layer thicknesses ($\delta_{u}/\delta_{\theta}$) to
define the transition from the buoyancy-dominated to the rotation-dominated
regime. Although this definition works for $Pr<Pr_{cr}$, it fails for the
free-slip boundary conditions, where no Ekman layer exists. It is evident that
we still need a more general criterion to define the transition from the
buoyancy-dominated regime to the rotation-dominated regime.

We thank Roshan Samuel for his valuable assistance in the development of the
solver used in this study. We also thank Prof. M. K. Verma for inspiring us to
utilize GPUs for scientific computing, and also for providing resources for
the testing of the solver and for running some simulations. Mohammad Anas
thanks Soumyadeep Chatterjee, Shadab Alam, and Manthan Verma for their useful
discussions on the solver and this work. For all the simulations related to
this work, we gratefully acknowledge the support and the resources provided by
Param Sanganak under the National Supercomputing Mission, Government of India
at the Indian Institute of Technology, Kanpur.

## References

  * Chandrasekhar (1961) S. Chandrasekhar,  _Hydrodynamic and Hydromagnetic Stability_ (Oxford University Press, Oxford, 1961).
  * Drazin and Reid (2004) P. G. Drazin and W. H. Reid,  _Hydrodynamic Stability_ (Cambridge University Press, Cambridge, 2004).
  * Shraiman and Siggia (1990) B. I. Shraiman and E. D. Siggia, Phys. Rev. A 42, 3650 (1990).
  * Grossmann and Lohse (2000) S. Grossmann and D. Lohse, J. Fluid Mech. 407, 27 (2000).
  * Niemela  _et al._ (2000) J. J. Niemela, L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly, Nature 404, 837 (2000).
  * Chillà and Schumacher (2012) F. Chillà and J. Schumacher, Eur. Phys. J. E 35, 58 (2012).
  * Glatzmaier and Roberts (1995) G. A. Glatzmaier and P. H. Roberts, Nature 377, 203 (1995).
  * Marshall and Schott (1999) J. Marshall and F. Schott, Rev. Geophy. 37, 1 (1999).
  * Ingersoll (1990) A. P. Ingersoll, Science 248, 308 (1990).
  * Heimpel  _et al._ (2005) M. Heimpel, J. Aurnou, and Wicht, Nature 438, 193–196 (2005).
  * Spiegel (1971) E. A. Spiegel, Annu. Rev. Astron. Astrophys. 9, 323 (1971).
  * Rossby (1969) H. T. Rossby, J. Fluid Mech. 36, 309 (1969).
  * King  _et al._ (2008) E. M. King, S. Stellmach, J. Noir, U. Hansen, and J. M. Aurnou, Nature 457, 301 (2008).
  * King  _et al._ (2012) E. M. King, S. Stellmach, and J. M. Aurnou, J. Fluid Mech. 691, 568 (2012).
  * Stevens  _et al._ (2009a) R. J. A. M. Stevens, J.-Q. Zhong, H. J. H. Clercx, G. Ahlers, and D. Lohse, Phys. Rev. Lett. 103, 024503 (2009a).
  * Kunnen (2021) R. P. J. Kunnen, J. Turbulence 22, 267 (2021).
  * Ecke and Shishkina (2023) R. E. Ecke and O. Shishkina, Annual Review of Fluid Mechanics 55, 603 (2023).
  * Niemela  _et al._ (2010) J. J. Niemela, S. Babuin, and K. R. Sreenivasan, J. Fluid Mech. 649, 509 (2010).
  * Ecke and Niemela (2014) R. E. Ecke and J. J. Niemela, Phys. Rev. Lett. 113, 114301 (2014).
  * Stellmach  _et al._ (2014) S. Stellmach, M. Lischper, K. Julien, G. Vasil, J. S. Cheng, A. Ribeiro, E. M. King, and J. M. Aurnou, Phys. Rev. Lett. 113, 254501 (2014).
  * Zhong  _et al._ (1993) F. Zhong, R. E. Ecke, and V. Steinberg, J. Fluid Mech. 249, 135 (1993).
  * Liu and Ecke (1997) Y. Liu and R. E. Ecke, Phys. Rev. Lett. 79, 2257 (1997).
  * Kunnen  _et al._ (2006) R. P. J. Kunnen, H. J. H. Clercx, and B. J. Geurts, Phys. Rev. E 74, 056306 (2006).
  * Stevens  _et al._ (2013) R. J. Stevens, H. J. Clercx, and D. Lohse, European J. of Mech. - B/Fluids 40, 41 (2013).
  * Wei  _et al._ (2015) P. Wei, S. Weiss, and G. Ahlers, Phys. Rev. Lett. 114, 114506 (2015).
  * Zhong  _et al._ (2009) J.-Q. Zhong, R. J. A. M. Stevens, H. J. H. Clercx, R. Verzicco, D. Lohse, and G. Ahlers, Phys. Rev. Lett. 102, 044502 (2009).
  * Yang  _et al._ (2020) Y. Yang, R. Verzicco, D. Lohse, and R. J. A. M. Stevens, Phys. Rev. Fluids 5, 053501 (2020).
  * Vorobieff and Ecke (2002) P. Vorobieff and R. E. Ecke, J. Fluid Mech. 458, 191–218 (2002).
  * Stevens  _et al._ (2009b) R. J. A. M. Stevens, J.-Q. Zhong, H. J. H. Clercx, G. Ahlers, and D. Lohse, Phys. Rev. Lett. 103, 024503 (2009b).
  * Stevens  _et al._ (2010a) R. J. A. M. Stevens, H. J. H. Clercx, and D. Lohse, New J. Phys. 12, 075005 (2010a).
  * Stevens  _et al._ (2010b) R. J. A. M. Stevens, H. J. H. Clercx, and D. Lohse, Phys. Fluids 22, 085103 (2010b).
  * Stevens  _et al._ (2011) R. J. A. M. Stevens, J. Overkamp, D. Lohse, and H. J. H. Clercx, Phys. Rev. E 84, 056313 (2011).
  * Joshi  _et al._ (2017) P. Joshi, H. Rajaei, R. P. J. Kunnen, and H. J. H. Clercx, J. Fluid Mech. 830, R3 (2017).
  * Weiss  _et al._ (2010) S. Weiss, R. J. A. M. Stevens, J.-Q. Zhong, H. J. H. Clercx, D. Lohse, and G. Ahlers, Phys. Rev. Lett. 105, 224501 (2010).
  * Weiss  _et al._ (2016) S. Weiss, P. Wei, and G. Ahlers, Phys. Rev. E 93, 043102 (2016).
  * Chong  _et al._ (2017) K. L. Chong, Y. Yang, S.-D. Huang, J.-Q. Zhong, R. J. A. M. Stevens, R. Verzicco, D. Lohse, and K.-Q. Xia, Phys. Rev. Lett. 119, 064501 (2017).
  * Julien  _et al._ (1996) K. Julien, S. Legg, J. Mcwilliams, and J. Werne, J. Fluid Mech. 322, 243–273 (1996).
  * Pfotenhauer  _et al._ (1984) J. M. Pfotenhauer, P. G. J. Lucas, and R. J. Donnelly, J. Fluid Mech. 145, 239–252 (1984).
  * Kunnen  _et al._ (2016) R. P. J. Kunnen, R. Ostilla-Mónico, E. P. van der Poel, R. Verzicco, and D. Lohse, J. Fluid Mech. 799, 413–432 (2016).
  * Hartmann  _et al._ (2023) R. Hartmann, G. S. Yerragolam, R. Verzicco, D. Lohse, and R. J. A. M. Stevens, (2023), arXiv:2305.02127 [physics.flu-dyn] .
  * Huang  _et al._ (2013) S.-D. Huang, M. Kaczorowski, R. Ni, and K.-Q. Xia, Phys. Rev. Lett. 111, 104501 (2013).
  * (42) “See Supplemental Material at [URL will be inserted by publisher] for the discussion related to the simulations and the solvers,” .