# Weil–Petersson geodesics on the modular surface

Vaibhav Gadre  School of Mathematics and Statistics
University of Glasgow
University Place
Glasgow
G12 8QQ United Kingdom<EMAIL_ADDRESS>
###### Abstract.

We consider the Weil–Petersson (WP) metric on the modular surface. We lift WP
geodesics to the universal cover of the modular surface and analyse geometric
properties of the lifts as paths in the hyperbolic metric on the universal
cover. For any pair of distinct points in the thick part of the universal
cover, we prove that the WP and hyperbolic geodesic segments that connect the
pair, fellow-travel in the thick part and all deviations between these
segments arise during cusp excursions. Furthermore, we give a quantitative
analysis of the deviation during an excursion.

We leverage the fellow traveling to derive a correspondence between recurrent
WP and hyperbolic geodesic rays from a base-point. We show that the
correspondence can be promoted to a homeomorphism on the circle of directions.
By analysing cuspidal winding of a typical WP geodesic ray, we show that the
homeomorphism pushes forward a Lebesgue measure on the circle to a singular
measure. In terms of continued fraction coefficients, the singularity boils
down to a comparison that we prove, namely, the average coefficient is bounded
along a typical WP ray but unbounded along a typical hyperbolic ray.

###### Key words and phrases:

Teichmüller theory, Moduli of Riemann surfaces.

###### 2010 Mathematics Subject Classification:

30F60, 32G15, 37D40, 53D25, 37A25

##  1\. Introduction

We consider the Weil–Petersson metric on the modular surface $S$. Lifting a WP
geodesic segment to the universal cover of $S$, we analyse its geometric
properties in the hyperbolic metric on the universal cover. The thick part of
the universal cover is the complement of lifts of cusp neighbourhoods. For any
pair of distinct points in the thick part, we prove that the WP and hyperbolic
geodesic segments connecting these points fellow travel in the thick part with
all deviations between them arising during excursions in cusp neighbourhoods.
Furthermore, we give a quantitative analysis of these deviations during an
excursion. Our results should extend to mildly constrained general WP-type
metrics. For this reason, we keep our discussions on a more general footing.

We now provide a more precise description. For notation, we denote $S$ with
the WP metric by $X$ and $S$ with the hyperbolic metric by $Y$. We will denote
the universal cover of $S$ by $\mathbb{D}$ and the lift of the metrics by
$d_{\operatorname{WP}}$ and $d_{\operatorname{hyp}}$ respectively.

Let $\delta$ denote the distance to the cusp in $X$. For some $B>0$ small
enough, we may define a cusp neighbourhood by the condition $\delta<B$. We
call the boundary locus $\delta<B$, a WP-horocycle. When $B$ is small enough
the cusp neighbourhood $N=\\{\delta<B\\}$ is convex.

The convexity of $N$ and the compactness of $X\setminus N$ have the following
simple consequence: there exists a smaller neighbourhood $\\{\delta<b\\}$
where $0<b<B$ such that for any pair of points $p,q$ in $X\setminus N$, the
shortest WP geodesic connecting them is disjoint from $\\{\delta<b\\}$.

The consequence stated above implies that a lift of $N$ to the universal cover
$\mathbb{D}$ is also quasi-convex in the hyperbolic metric on $\mathbb{D}$. So
a lift of $N$ to $\mathbb{D}$ is approximated by a horoball in
$d_{\operatorname{hyp}}$. In Section 2, we will clarify in a more precise
technical sense how the hyperbolic horocycle bounding the largest horoball
contained in a lift of $N$ approximates the WP horocycle. The
$\pi_{1}(S)$-orbit of the horoball gives us a collection $\mathcal{H}$ of
disjoint horoballs in $\mathbb{D}$. We shall call the region in $\mathbb{D}$
that is complementary to the collection of these horoballs the _thick part_ of
$\mathbb{D}$. Note that quotient of the thick part by $\pi_{(}S)$ is thick in
either metric, that is, it is disjoint from a cusp neighbourhood in both
$d_{\operatorname{WP}}$ and $d_{\operatorname{hyp}}$.

Let $p,q$ be a pair of distinct points in the thick part. Let
$\gamma_{\operatorname{WP}}[p,q]$ and $\gamma_{\operatorname{hyp}}[p,q]$ be
respectively the WP and hyperbolic geodesic segments between $p$ and $q$.

A _thick-thin decomposition_ of the WP segment
$\gamma_{\operatorname{WP}}[p,q]$ is a partition of
$\gamma_{\operatorname{WP}}[p,q]$ into WP geodesic sub-segments as

(1.1)
$\gamma_{\operatorname{WP}}[p,q]=[y_{0},y_{1}]\cup[y_{1},y_{2}]\cup\dots\cup[y_{2k-1},y_{2k}]\cup[y_{2k},y_{2k+1}],$

where

  * •

$k\geqslant 0$,

  * •

$y_{0}=p$, $y_{2k+1}=q$,

  * •

the segments $[y_{2r},y_{2r+1}]$ for $r=0,k$ are contained in the thick part,
and

  * •

the open segments $(y_{2r-1},y_{2r})$ for $r=1,\cdots,k$ if non-empty, are
contained in a horoball in $\mathcal{H}$.

Analogously, a thick-thin decomposition of $\gamma_{\operatorname{hyp}}[p,q]$
is a partition of $\gamma_{\operatorname{hyp}}[p,q]$ into hyperbolic geodesic
sub-segments as

(1.2)
$\gamma_{\operatorname{hyp}}[p,q]=[z_{0},z_{1}]\cup[z_{1},z_{2}]\cup\dots\cup[z_{2j-1},z_{2j}]\cup[z_{2j},z_{2j+1}],$

where

  * •

$j\geqslant 0$,

  * •

$z_{0}=p$, $z_{2j+1}=q$ and

  * •

the segments $[z_{2r},z_{2r+1}]$ for $r=0,j$ are contained in the thick part,
and

  * •

the open segments $(z_{2r-1},z_{2r})$ for $r=1,\cdots,j$, if non-empty, are
contained in a horoball in $\mathcal{H}$

We prove the following theorem.

######  Theorem 1.3.

There is a constant $R>0$ such that for any pair of distinct points $p,q$ in
the thick part of $\mathbb{D}$, there exist thick-thin decompositions (1.1)
and (1.2) of $\gamma_{\operatorname{WP}}[p,q]$ and
$\gamma_{\operatorname{hyp}}[p,q]$ respectively, such that $j=k$ and the
segments $[z_{2r},z_{2r+1}]$ $R$-fellow travel (in either metric) the segments
$[y_{2r},y_{2r+1}]$ for $r=0,\cdots,k$.

We will also estimate the deviation in $d_{\operatorname{hyp}}$ between the
geodesics during each cusp excursion, that is, we will estimate the hyperbolic
distance between corresponding points on the open segments $(z_{2r-1},z_{2r})$
and $(y_{2r-1},y_{2r})$ when they are both non-empty and large enough.

The coarse geometric viewpoint offers a clarifying perspective of Theorem 1.3.
Recall that the modular surface is the moduli space of hyperbolic once-
punctured tori. The WP metric completion of $\mathbb{D}$ includes noded
surfaces that arise by pinching the systole on a marked once-punctured torus.
It follows that the WP completion is quasi-isometric to the Farey graph with
every edge having length one. In either metric, the thick part of the
universal cover $\mathbb{D}$ is quasi-isometric to the tree dual to the Farey
triangulation. A geodesic segment in the tree can be written as a finite
sequence of alternating rights and lefts. Any pair of triangles that differ
only by rights (or only by lefts) share a vertex in the Farey graph. When the
number of consecutive rights (or lefts) is more than two, the shortest path in
the Farey graph has length two and passes through the shared vertex. In other
words, the WP metric has a shortcut that passes through the noded surface
corresponding to the shared vertex. The actual WP geodesic follows this
shortcut closely by doing a cusp excursion winding around the noded surface.
This feature is shared by the hyperbolic metric except that the cusps are
infinitely far away. The actual hyperbolic geodesic achieves its shortcut by a
cusp excursion of hyperbolic length that is logarithmic in the number of
consecutive rights (or lefts).

###  1.4. The circle map:

Let $p\in\mathbb{D}$ be a base-point. Passing to a smaller WP cusp
neighbourhood if required, we may assume that $p$ is in the thick part.
Consider a WP ray $\gamma$ from $p$ that recurs to the thick part. By using
Theorem 1.3 along recurrence times $\gamma_{t}$, we get hyperbolic geodesic
segments $\gamma_{\operatorname{hyp}}[p,\gamma_{t}]$. We will show that we can
pass to a limit along these hyperbolic geodesic segments to associate a
hyperbolic ray $\gamma^{\prime}$ from $p$ which fellow travels $\gamma$ in the
thick part and all deviations arise during cusp excursions. In particular, the
ray $\gamma^{\prime}$ is also recurrent to the thick part.

We may identify $S^{1}$ with the unit tangent circles
$T_{p}^{1}(\mathbb{D},d_{\operatorname{WP}})$ and
$T_{p}^{1}(\mathbb{D},d_{\operatorname{hyp}})$ in the respective metrics. The
limiting argument on recurrent directions sketched above defines a map $\psi$
from a dense subset of $S^{1}$ to $S^{1}$. In the reverse direction, we can
apply a similar limiting argument to associate a recurrent WP ray to every
recurrent hyperbolic ray. This defines a circle map $\xi$ from a dense subset
of $S^{1}$ to $S^{1}$. It also follows that the compositions $\psi\circ\phi$
and $\phi\circ\psi$ are always defined and $\phi\circ\psi=\text{id}$ and
$\psi\circ\phi=\text{id}$. Thus $\psi$ is a bijection from the dense set of
recurrent WP-directions to the dense set of hyperbolic recurrent directions.
We will further show that $\psi$ and $\phi$ are monotone. From the density of
the domains and the monotonicity, we may conclude that $\psi$ and $\phi$
extend to homeomorphisms of $S^{1}$ that are inverses of each other.

Let $\ell$ be a Lebesgue measure on $S^{1}$. The Liouville measures
$\mu_{\operatorname{WP}}$ or $\mu_{\operatorname{hyp}}$ are absolutely
continuous with respect to each other. After identifying $S^{1}$ with either
of the unit tangent circles $T_{p}^{1}(\mathbb{D},d_{\operatorname{WP}})$ and
$T_{p}^{1}(\mathbb{D},d_{\operatorname{hyp}})$, we may consider the measure
$\ell$ to be the conditional measure on the circle of directions obtained from
either of these Liouville measures. The WP flow is ergodic [9], and in fact,
exponentially mixing [3]. Thus, after identifying
$T_{p}^{1}(\mathbb{D},d_{\operatorname{WP}})$ with $S^{1}$, it follows that
$\ell$-almost every $v\in S^{1}$ determines a recurrent WP ray. Thus the map
$\psi$ is well-defined on a full measure subset of $S^{1}$.

We prove:

######  Theorem 1.5.

The push-forward measure $\psi_{\ast}(\ell)$ is singular with respect to
$\ell$.

Theorem 1.5 boils down to a comparison of the statistics of cuspidal winding
numbers along typical WP and hyperbolic rays. In Section 5, we will
reformulate this comparison in terms of continued fraction coefficients along
typical WP and hyperbolic rays, namely, we will show that the "average"
coefficient is bounded along a typical WP ray but unbounded along a typical
hyperbolic ray.

###  1.6. WP-type metrics on finite type surfaces:

We now mention the constraints under which the results presented here should
extend to general WP-type metrics on surfaces of finite type. The first
requirement is that cusp neighbourhoods can be chosen to be convex and their
horocycles closely approximate the corresponding horocycles for the surfaces
of revolution that model them. The precise nature of the approximation
required is discussed in Section 2.1 and we discuss why it holds for the
modular surface.

For a cusp modelled on the surface of revolution for $y=x^{r}\,:\,r\geqslant
3$, the convexity of a small enough cusp neighbourhood translates to a
constraint on the lower order terms in the expansion of the metric near the
cusp with the metric on a surface of revolution as the leading term. See [3,
Remark on page 242] for the details.

Similarly, the required properties for WP horocycles mentioned in Section 2.1
should also be derivable once there is a reasonable constraint on the lower
order terms in the expansion of the cusp metric. We leave the details of this
derivation to the interested reader.

###  1.7. The WP metric on non-exceptional moduli spaces:

The analogous passage between WP and Teichmüller geodesics on non-exceptional
moduli spaces is likely to be much more involved. For example, in the
exceptional case the existence of the limiting ray relies on the curvatures
being pinched, that is, negative and bounded away from zero. This is no longer
true for non-exceptional moduli.

Brock-Masur-Minsky associated to ending laminations to WP rays. See [1,
Definition 2.5]. It then turns out that there exist recurrent WP rays with
non-uniquely ergodic ending laminations. See [2, Theorem 1]. As such, a naive
limiting argument as we have above cannot hold. It is possible however typical
WP rays do admit a more subtle limiting argument.

###  1.8. Acknowledgements:

I am extremely grateful to Carlos Matheus for the extensive clarifying
discussions that I had with him during the course of this work. I would also
like to thank Scott Wolpert for his helpful comments on an earlier draft of
the paper.

##  2\. Fellow traveling in the thick part

###  2.1. WP and hyperbolic horocycles in model co-ordinates:

To compare the two metrics near the cusp, it will be convenient to work in the
upper half space co-ordinates $\\{z=x+iy)\,:\,y>0\\}$ on $\mathbb{D}$ with the
cusp at infinity. The Poincare metric on the upper half-space descends to the
hyperbolic metric on the modular surface. As such, a hyperbolic cusp
neighbourhood in these co-ordinates can be chosen to be of the form
$\\{y>\text{constant}\\}$.

Viewed as the Teichmüller space of hyperbolic once-punctured tori, the upper
half-space co-ordinates are closely related to the Fenchel-Nielsen co-
ordinates. For a marked hyperbolic once-punctured torus $\tau\in S$, let
$\ell=\ell_{\tau}$ be the hyperbolic length of its systole. Let $t=t_{\tau}$
be the twist parameter about the systole. In [12, Corollary 4.10], Wolpert
showed that the WP distance $\delta(\tau)$ of $\tau$ to the cusp of $S$
satisfies

(2.2) $\delta=\sqrt{2\pi\ell}+O(\ell^{5/2}).$

To keep the model expressions for the WP metric simple, we will allow our
choice of the co-ordinates $(x,y)$ to be flexible up to a homothety of the
form $z\to\alpha z$ for $\alpha\in\mathbb{R}_{>0}$. If $B$ is small enough,
the homothety constants we will need to choose are uniform. This means that
our estimates will change only up to multiplicative constants that depend on
$B$, which is small enough and fixed once and for all.

Up to such a uniform homothety, the co-ordinates $(x,y)$ near $\tau$ written
in terms of the Fenchel–Nielsen co-ordinates $(\ell,t)$ have the following
form:

(2.3) $\ell=\frac{1}{y}+h(x,y),\hskip 15.0ptt=\frac{x}{y}+h^{\prime}(x,y).$

where $h(x,y)$ and $h^{\prime}(x,y)$ comprise of lower order terms. In
particular, we may write $h(x,y)$ as $h(x,y)=g(x,y)/y^{2}$, where we may
assume that both $g$ and its first partials are bounded. We will use
expressions in (2.3) to pass from Fenchel–Nielsen co-ordinates to the upper
half space co-ordinates $(x,y)$.

From equation (2.2) above, a WP horocycle $\\{\delta=B\\}$ corresponds to
$y=2\pi/B^{2}$ to the first order. Up to uniform homothety, this is the same
as $y=1/B^{2}$. As

$\frac{\partial\delta}{\partial
y}=\left(\sqrt{\frac{\pi}{2\ell}}+O(\ell^{3/2})\right)\frac{\partial\ell}{\partial
y}$

if $B$ is small enough then so is $\ell$ and hence $\partial\ell/\partial
y\neq 0$. By the implicit function theorem, the horocycle $\\{\delta=B\\}$ can
be written as $(x,f_{B}(x))$. Now the derivative
$df_{B}/dx=-(\partial\delta/\partial x)/(\partial\delta/\partial
y)=-(\partial\ell/\partial x)/(\partial\ell/\partial y)$. By direct
computation we get

$\frac{\partial\ell/\partial x}{\partial\ell/\partial
y}=\frac{(1/y^{2})\partial g/\partial x}{(-1/y^{2})(1-\partial g/\partial
y+2g/y)}=-\frac{\partial g/\partial x}{1-\partial g/\partial y+2g/y}$

As $B$ and hence $y$ is small, it follows that the right hand side is bounded.
We thus conclude that when the horocycle $\\{\delta=B\\}$ is viewed as a graph
of a function over $x$, the derivative at any point of this graph is bounded.
In particular, this implies that the Euclidean arc-length element along
$\\{\delta=B\\}$ is comparable to $dx$.

The form of the WP metric can be derived by invoking Wolpert’s estimates [12]
in the special case of the modular surface. By converting [12, Page 284] to
the co-ordinates $(x,y)$ up to a uniform homothety, the WP metric is modelled
to the first order by

(2.4) $\frac{dx^{2}+dy^{2}}{y^{3}}.$

Let $I[a,b]$ be the strip $\\{(x,y)\in\mathbb{D}\,:\,a\leqslant x\leqslant
b\\}$. Using the model metric (2.4) and the discussion on the Euclidean arc
length element it follows that if $(b-a)$ is large enough (depending only on
$B$) then the WP length of $\\{\delta=B\\}\cap I[a,b]$ equals $(b-a)/B^{6}$ up
to a multiplicative constant that depends only on $B$. We note that
$(b-a)/B^{6}$ is also the length of the segment $\\{y=1/B^{2}\\}\cap I[a,b]$
in the model metric. On the other hand, the hyperbolic length of
$\\{y=1/B^{2}\\}\cap I[a,b]$ equals $(b-a)/B^{4}$. Thus, up to a
multiplicative constant that depends only on $B$, the WP and hyperbolic
lengths of these related horocycle segments are equal, and furthermore we fix
$B$ once and for all ahead of time. This observation will be crucial in the
discussion that precedes Lemma 2.7.

###  2.5. Projected Paths:

Let $\operatorname{Thick}$ denote the thick part. Let
$\pi_{\operatorname{WP}}:\mathbb{D}\to\operatorname{Thick}$ and
$\pi_{\operatorname{hyp}}:\mathbb{D}\to\operatorname{Thick}$ be the closest
point projections on to the thick part in the corresponding metrics. For any
pair of points $p,q$ in $\operatorname{Thick}$ let
$\gamma_{\operatorname{WP}}[p,q]$ and $\gamma_{\operatorname{hyp}}[p,q]$ be
the WP and hyperbolic geodesics connecting them. We call
$\pi_{\operatorname{hyp}}(\gamma_{\operatorname{hyp}})$ the _hyperbolic
projected path_ of $\gamma_{\operatorname{hyp}}$. Similarly, we call
$\pi_{\operatorname{WP}}(\gamma_{\operatorname{WP}})$ the _WP-projected path_
of $\gamma_{\operatorname{WP}}$.

######  Proposition 2.6.

For any pair of points $p,q$ in $\operatorname{Thick}$ the projected paths
$\pi_{\operatorname{hyp}}(\gamma_{\operatorname{hyp}}[p,q])$ and
$\pi_{\operatorname{WP}}(\gamma_{\operatorname{WP}}[p,q])$ are quasi-geodesics
in $\operatorname{Thick}$ with the corresponding path metrics.

For the hyperbolic metric, this proposition is given in [7, Lemma 2.1]. Here,
we adapt the argument in [7, Lemma 2.1] to the WP metric. The WP metric on
$\mathbb{D}$ is negatively curved and pinched away from 0. Hence, it is
$\text{CAT}(-\kappa)$ for some $\kappa>0$ and thus Gromov hyperbolic. See [10,
Introduction].

Let $H$ be a horoball in the collection $\mathcal{H}$ and let
$\pi_{H}:\mathbb{D}\setminus H\to\partial H$ be the WP closest point
projection. In a Hadamard manifold with pinched negative curvature, the
closest point projection to a horosphere is coarsely contracting. See [5,
Section 4]. Here we give an analogous result for $\pi_{H}$ but focus only on
showing that $\pi_{H}$ is coarsely distance non-increasing as that is
sufficient for our purposes.

Before we state the projection result in Lemma 2.7, we make use of the
observation at the end of Section 2.1. The quotient $H/\pi_{1}(S)$ is a quasi-
convex cusp neighbourhood on the WP surface $X$. Let $0<B^{\prime}<B$ be the
largest constant such that $N^{\prime}=\\{\delta<B^{\prime}\\}$ is contained
in $H/\pi_{1}(S)$. The approximations 2.4 and 2.2 imply that if $B$ is small
enough then $B^{\prime}\geqslant B/2$. Let $H^{\prime}$ be the lift of
$N^{\prime}$ that is contained in $H$. Let
$\pi_{H^{\prime}}:\mathbb{D}\setminus H\to\partial H^{\prime}$ be the WP
closest point projection. By the observation at the end of Section 2.1, it
follows that for any segment $\gamma$ in $\mathbb{D}\setminus H$ the
projections $\pi_{H}(\gamma)$ and $\pi_{H^{\prime}}(\gamma)$ are coarsely the
same, that is, there exists constants $k\geqslant 1$ and $c\geqslant 0$ such
that

$\frac{1}{k}\ell_{\partial H}(\pi_{H}(\gamma))-c<\ell_{\partial
H^{\prime}}\pi_{H^{\prime}}(\gamma)<k\ell_{\partial H}(\pi_{H}(\gamma))+c$

where $\ell_{\partial H}$ and $\ell_{\partial H^{\prime}}$ denote the lengths
in the WP path metrics on $\partial H$ and $\partial H^{\prime}$ respectively.

We now show that $\pi_{H}$ is coarsely distance non-increasing.

######  Lemma 2.7.

The WP closest point projection $\pi_{H}:\mathbb{D}\setminus H\to\partial H$
is coarsely distance non-increasing, that is, there exists constants $K>1,C>0$
such that for any pair of points $p,q\in\mathbb{D}\setminus H$

$d_{\partial H}(\pi_{H}(p),\pi_{H}(q))\leqslant
Kd_{\operatorname{WP}}(p,q)+C.$

###### Proof.

We give a quick sketch. By our preceding observation that the projections
$\pi_{H}$ and $\pi_{H^{\prime}}$ are coarsely the same, it suffices to prove
this property for $\pi_{H^{\prime}}$.

First, suppose $p,q$ are distinct points in $\mathbb{D}\setminus H$ such that
the WP geodesic segment $[p,q]$ is contained in $\mathbb{D}\setminus H$. Let
$z$ be the point that completes the cusp in $H$. Let $[p,z]$ and $[q,z]$ be WP
geodesics from $p$ and $q$ to $z$.

As $B$ was chosen small enough for the metric near the cusp to approximate the
model metric (2.4) closely, the segments $[p,z]$ and $[q,z]$ intersect
$\partial H^{\prime}$ at almost right angles. This implies that there exists a
constant $c^{\prime}\geqslant 0$ that depends only on the metric and the
constant $B$ such that $\pi_{H^{\prime}}(p)$ is $c^{\prime}$-close to the
point $p^{\prime}=[p,z]\cap\partial H^{\prime}$. Similarly
$\pi_{H^{\prime}}(q)$ is $c^{\prime}$-close to the point
$q^{\prime}=[q,z]\cap\partial H^{\prime}$.

The triangle with vertices $z,p,q$ can now be compared to the comparison
triangle in constant negative curvature $-\kappa$ to conclude that the length
of $[p^{\prime}q^{\prime}]$ along $\partial H^{\prime}$ is strictly less than
$d_{\operatorname{WP}}(p,q)$. We deduce that the length of $\pi^{\prime}[p,q]$
along $\partial H^{\prime}$ is coarsely less than
$\ell_{\operatorname{WP}}[p,q]$.

Finally, any path in $\mathbb{D}\setminus H$ can be closely approximated by a
finite concatenation of geodesic segments in $\mathbb{D}\setminus H$. The
lemma follows. ∎

###### Proof of Proposition 2.6 for the WP metric.

Suppose that the number of horoballs in $\mathcal{H}$ that the geodesic
segment $\gamma_{\operatorname{WP}}[p,q]$ passes through is $k$. We decompose
$\pi_{\operatorname{WP}}(\gamma_{\operatorname{WP}})$ into segments
$[y_{0},y_{1}]\cup[y_{1},y_{2}]\cup[y_{2},y_{3}]\cup\dots\cup[y_{2k-1},y_{2k}]\cup[y_{2k},y_{2k+1}]$,
where $y_{0}=$, $y_{2k+1}=q$, all intervals of the form $[y_{2j},y_{2j+1]}$
are WP geodesic segments in $\operatorname{Thick}$ and the remaining intervals
are segments along distinct horocycles.

Let $\beta_{j}$ be the WP geodesics perpendicular to
$\gamma_{\operatorname{WP}}[y_{0},y_{2k+1}]$ at the points $y_{j}$. Each
$\beta_{j}$ separates $y_{0}$ from $y_{2k+1}$ so a thick geodesic from $y_{0}$
to $y_{2k+1}$ must intersect $\beta_{j}$. Moreover, the geodesics $\beta_{j}$
are all disjoint. Otherwise, for some $j$ we would find a triangle bounded by
the geodesic segments of $\beta_{j},\beta_{j+1}$ and
$\gamma_{\operatorname{WP}}[y_{j},y_{j+1}]$. As the WP metric has pinched
negative curvature, this violates the Gauss-Bonnet Theorem as the sides are
geodesic segments and two angles of the triangle are right angles.

In conclusion, the geodesics $\beta_{j}$ divide $\mathbb{D}$ into regions
$R_{j}$ such that each $R_{j}$ contains a single subsegment of
$\gamma_{\operatorname{WP}}[y_{0},y_{2k+1}]$ that is either entirely contained
in $\operatorname{Thick}$ or entirely contained in a horocycle. In each
$R_{j}$, we now show that the subsegment of the projected path is coarsely the
shortest path in the induced path metric, that is any thick geodesic segment
that makes it across $R_{j}$ between the endpoints has length coarsely larger
than the subsegment of the projected path of
$\pi_{\operatorname{WP}}(\gamma_{\operatorname{WP}})$ in $R_{j}$.

First consider a region $R_{j}$ where $[y_{2j},y_{2j+1}]$ is a WP geodesic
segment in $\operatorname{Thick}$. As the metric is $\text{CAT}(-\kappa)$ for
some $\kappa>0$, the nearest point projection to a geodesic is coarsely
contracting. Thus the length of any path that crosses from $\beta_{2j}$ to
$\beta_{2j+1}$ is coarsely at least the length of $[y_{2j},y_{2j+1}]$.

Now consider a region $R_{j}$ where $[y_{2j-1},y_{2j}]$ is a subsegment of a
horocycle $\partial H$. By lemma 2.7, the closest point projection
$\pi_{H}:\operatorname{Thick}\to\partial H$ is also coarsely contracting. Thus
the length of any thick path that crosses $R_{j}$ from $\beta_{2j-1}$ to
$\beta_{2j}$ is coarsely at least the length $\ell_{\partial
H}(y_{2j-1},y_{2j})$.

This finishes the proof of the proposition for the WP metric.

∎

The fundamental group $\pi_{1}(S)$ acts co-compactly on
$\operatorname{Thick}$. Hence, by the Svarć-Milnor lemma,
$\operatorname{Thick}$ with either WP or hyperbolic induced path metric is
quasi-isometric to the fundamental group $\pi_{1}(S)$ with a word metric given
by any finite generating set. In particular,
$(\operatorname{Thick},d_{\operatorname{WP}})$ and
$(\operatorname{Thick},d_{\operatorname{hyp}})$ are quasi-isometric. As a
further corollary of this quasi-isometry and of Proposition 2.6, we get

######  Corollary 2.8.

There exists a constant $R>0$ such that for any pair of points $p,q$ in
$\operatorname{Thick}$ the segments $\pi_{\operatorname{WP}}[p,q]$ and
$\pi_{\operatorname{hyp}}[p,q]$ $R$-fellow travel in $\operatorname{Thick}$.

###### Proof.

By [7, Lemma 2.1], the projected path $\pi_{\operatorname{hyp}}[p,q]$ is a
quasi-geodesic in $(\operatorname{Thick},d_{\operatorname{hyp}})$, and by
Proposition 2.6, the projected path $\pi_{\operatorname{WP}}[p,q]$ is a quasi-
geodesic in $(\operatorname{Thick},d_{\operatorname{WP}})$. As
$(\operatorname{Thick},d_{\operatorname{WP}})$ and
$(\operatorname{Thick},d_{\operatorname{hyp}})$ are quasi-isometric, both
projected paths $\pi_{\operatorname{WP}}[p,q]$ and
$\pi_{\operatorname{hyp}}[p,q]$ being quasi-geodesics, fellow travel.

∎

### Proof of Theorem 1.3:

By Corollary 2.8, the projected paths $\pi_{\operatorname{WP}}[p,q]$ and
$\pi_{\operatorname{hyp}}[p,q]$ $R$-fellow travel. This implies that there is
a constant $R^{\prime}>0$ such that if $\pi_{\operatorname{WP}}(\gamma_{WP})$
has a subsegment of length at least $R^{\prime}$ along a horocycle $\partial
H$ then $\pi_{\operatorname{hyp}}(\gamma_{hyp})$ must also contain a
subsegment of $\partial H$, and vice versa. Moreover, the respective entry and
exit points in $\partial H$ for the two projected paths must be within
distance $R$ of each other in $\operatorname{Thick}$ with the path metric. It
is possible that if $\pi_{\operatorname{WP}}(\gamma_{WP})$ has a subsegment
along $\partial H$ that is less than $R^{\prime}$ in length then
$\gamma_{\operatorname{hyp}}$ does not enter $H$ at all, and vice versa. If
such is the case then we choose the sub-segment along
$\gamma_{\operatorname{hyp}}$ that records the excursion in $H$ to be empty,
and vice versa. This concludes the proof that for any pair of points
$p,q\in\operatorname{Thick}$, the respective geodesics
$\gamma_{\operatorname{WP}}[p,q]$ and $\gamma_{\operatorname{hyp}}[p,q]$
fellow travel in the thick part and all deviations between them arise in cusp
neighbourhoods.

##  3\. Deviations during excursions

We now quantify deviations during a single excursion. Consider a cusp
excursion of a WP geodesic during which the minimum distance to the cusp
satisfies $\delta_{\textrm{min}}=1/D$ for $D>1$.

In upper half space co-ordinates, a WP geodesic during an excursion describes
a curve $\gamma_{\operatorname{WP}}(t)=(x(t),y(t))$ whose tangent vector
$v(t)$ has (up to a uniform multiplicative constant) unit size with respect to
the model metric (2.4). For convenience, we apply a homothety whose constant
depends only on $B$ to assume that we are working with
$\\{z\in\mathbb{C}\,:\,\text{Im}z\geqslant 1\\}$ as the cusp neighbourhood.
This means that all our estimates below are correct up to a multiplicative
constant that depends only on $B$.

An excursion in which $\delta_{\textrm{min}}=1/D$ with $D>1$, corresponds to
$y_{\text{max}}=D^{2}$. By [8, Proposition 5.5], the winding number about the
cusp for such an excursion is $D^{2}$ up to a uniform multiplicative constant.
In particular, this implies that during such an excursion the point that
distance along $\partial H=\\{z\,:\,\text{Im}z=1\\}$ between the entry and
exit points is $D^{2}$ up a uniform multiplicative constant.

We want to compare $\gamma_{\operatorname{WP}}$ with a hyperbolic geodesic
$\gamma_{\operatorname{hyp}}$ that enters and exits
$\\{z\in\mathbb{C}:\,\mid\text{Im}z\mid\geqslant 1\\}$ within bounded
hyperbolic distance of the $\gamma_{\operatorname{WP}}$-entry and exit points.
So we may assume that in the model co-ordinates the entry and exit points of
$\gamma_{\operatorname{WP}}$ are $(0,1)$ and $(2D^{2},1)$ respectively.

We will analyse the two trajectories $\gamma_{\operatorname{WP}}$ and
$\gamma_{\operatorname{hyp}}$ over the first half of their excursions, that is
between the entry point and the point with the largest imaginary part
$y_{\textrm{max}}$. By symmetry, the same analysis holds for the trajectories
on their way out. For the first half, the hyperbolic geodesic is given by the
locus $=\left(x,\sqrt{1+2D^{2}x-x^{2}}\right)$. To prove our deviation bounds,
we will similarly write down the locus for $\gamma_{\operatorname{WP}}$ and
then compare with above locus.

######  Lemma 3.1.

A point on the ingoing half of $\gamma_{\operatorname{WP}}$ is given in co-
ordinates by

$\asymp\left(x,D^{6/5}x^{2/5}\right)$

where $\asymp$ means up the individual co-ordinates may differ up to a uniform
multiplicative constant.

###### Proof.

The co-ordinate $x(\delta)$ is the winding number when
$\gamma_{\operatorname{WP}}$ is distance $\delta$ from the cusp. In the
calculation using dyadic intervals in the proof of [8, Proposition 5.5], let
$k$ be such that

$c2^{k}\asymp\frac{1}{D^{3}\delta^{3}}$

where $c=1/(D^{3}\delta_{0}^{3})$. Recall that $J$ was chosen so that

$c2^{J}\asymp 1$

Then repeating line by line the calculation in [8, Proposition 5.5], we have

$x(\delta)\asymp D^{3}\sum_{j=1}^{k}(c2^{j})^{5/3}\asymp
D^{2}(c2^{k})^{5/3}\asymp\frac{1}{D^{3}\delta^{5}}$

On the other hand,

$y(\delta)\asymp\frac{1}{\delta^{2}}$

By expressing $y$ in terms of $x$, we get that a point on the ingoing
trajectory of $\gamma_{\operatorname{WP}}$ has coordinates

$\asymp\left(x,D^{6/5}x^{2/5}\right)$

finishing the proof. ∎

It is interesting to compare the "shapes" of the hyperbolic and WP
trajectories that we considered above. A hyperbolic geodesic
$\gamma_{\operatorname{hyp}}$ that has coarsely the same entry and exit
points, namely $(0,1)$ and $(2D^{2},1)$, also attains the maximal imaginary
part $D^{2}$, that is, it is about the same height in the upper half space co-
ordinates as $\gamma_{\operatorname{WP}}$. However, from Lemma 3.1, we observe
that quantitatively the WP trajectory has faster initial gradient. So it goes
deep in to the cusp neighbourhood more quickly and winds thereafter. This is
consistent with the standard intuition of how geodesics wind around the cusp
in the WP metric.

##  4\. The circle map

We fix a base-point $p\in\mathbb{D}$. By shrinking the cusp neighbourhoods if
required, we may assume that $p$ lies in the thick part. By the results in
[9], the WP geodesic flow is ergodic. This implies that a dense set of
directions $v\in T^{1}_{p}(\mathbb{D},d_{\operatorname{WP}})$ give WP rays
that are recurrent to the thick part. In fact, identifying
$T^{1}_{p}(\mathbb{D},d_{\operatorname{WP}})$ with $S^{1}$, the set of
recurrent directions has full measure in the Lebesuge measure on $S^{1}$.

For a recurrent direction $v$, let $\gamma^{v}_{\operatorname{WP}}(t)$ denote
the recurrent WP ray that it defines. Let $t_{n}$ be any monotonic sequence of
recurrence times such that $t_{n}\to\infty$. We consider the sequence of
hyperbolic geodesic segments $[x,\gamma^{v}_{\operatorname{WP}}(t_{n})]$. By
Theorem 1.3, the hyperbolic and WP projected paths fellow travel. It follows
that hyperbolic geodesic segments $[x,\gamma^{v}_{\operatorname{WP}}(t_{n})]$
fellow travel. So we can pass to a limit to get a unique limiting hyperbolic
ray $\gamma^{v}_{\operatorname{hyp}}$ that fellow travels
$\gamma^{v}_{\operatorname{WP}}$ in $\operatorname{Thick}$. In particular, the
hyperbolic ray $\gamma^{v}_{\operatorname{hyp}}$ is recurrent to
$\operatorname{Thick}$. Identifying
$T^{1}_{p}(\mathbb{D},d_{\operatorname{hyp}})$ also with $S^{1}$, we get a map
$\psi$ from a dense (full measure) subset of $S^{1}$ to $S^{1}$ by sending $v$
to the unit tangent vector at $p$ to $\gamma^{v}_{\operatorname{hyp}}$.

In similar fashion, we get the reverse map $\xi$ from a dense (full measure)
subset of $S^{1}$ to $S^{1}$. For a recurrent hyperbolic ray $\zeta$
determined by $v\in T^{1}_{p}(\mathbb{D},d_{\operatorname{hyp}})$, we can
consider a monotonic sequence of recurrence times $s_{n}$ such that
$s_{n}\to\infty$ and consider the WP geodesic segments $[x,\zeta(s_{n})]$. As
$(\mathbb{D},d_{\operatorname{WP}})$ is $\text{CAT}(-\kappa)$ for some
$\kappa>0$, the longer and longer WP geodesic segments fellow travel. Thus, we
can pass to a limiting WP ray. By setting up identifications with $S^{1}$ and
sending $v$ to the unit tangent base vector for this WP ray, we define the map
$\xi$ on a dense subset of $S^{1}$ to $S^{1}$.

By Theorem 1.3, the WP ray defined by the direction $\xi\circ\psi(v)$ fellow
travels in $\operatorname{Thick}$ the WP ray defined by $v$. As
$(\mathbb{D},d_{\operatorname{WP}})$ is $\text{CAT}(-\kappa)$ for some
$\kappa>0$, it follows that the ray defined by $\xi\circ\psi(v)$ has to fellow
travel the ray defined by $v$ for all times. Hence, the rays are identical and
$\xi\circ\psi=\text{id}$. By a similar reasoning, the composition
$\psi\circ\xi$ is also well defined and identity. Thus, $\psi$ is a bijection
from the dense set of recurrent WP-directions to the dense set of recurrent
hyperbolic directions.

###  4.1. Monotonicity:

We now derive the monotonicity of the map $\psi$. First, we fix some notation.
Given $\alpha,\beta\in S^{1}$, we denote by $I(\alpha,\beta)$ the (open)
sector of directions in $S^{1}$ going counter-clockwise from $\alpha$ to
$\beta$. We identify the unit tangent spaces
$T^{1}_{p}(\mathbb{D},d_{\operatorname{WP}})$ and
$T^{1}_{p}(\mathbb{D},d_{\operatorname{hyp}})$ with $S^{1}$ in way that
preserves the orientation induced by the surface $S$.

######  Lemma 4.2.

Suppose $v,w\in\text{Domain}(\psi)$ and distinct. For any $\theta\in
I(v,w)\cap\text{Domain}(\psi)$, we have $\psi(\theta)\in I(\psi(v),\psi(w))$.

###### Proof.

The proof follows the same Gauss-Bonnet type arguments. Suppose that there is
some $\theta\in I(v,w)$ such that $\psi(\theta)$ is not contained in
$I(\psi(v),\psi(w))$.

Let $S(\alpha,\beta)$ be the open set in $(\mathbb{D},d_{\operatorname{hyp}})$
swept out by hyperbolic-rays with initial vectors in $I(\alpha,\beta)$. It is
bounded by the hyperbolic rays $\gamma^{\alpha}_{\operatorname{hyp}}$ and
$\gamma^{\beta}_{\operatorname{hyp}}$.

As $\psi(\theta)$ is not in $I(\psi(v),\psi(w))$, the hyperbolic ray
$\gamma^{\psi(\theta)}_{\operatorname{hyp}}$ does not intersect
$S(\psi(v),\psi(w))$. Let us fix a monotonic sequence $t_{n}$ of recurrence
times along $\gamma^{\psi(\theta)}_{\operatorname{hyp}}$ such that
$t_{n}\to\infty$. Suppose that the $R$-ball centred at some
$\gamma^{\psi(\theta)}_{\operatorname{hyp}}(t_{n})$ is disjoint from the rays
$\gamma^{\psi(v)}_{\operatorname{hyp}}$ and
$\gamma^{\psi(w)}_{\operatorname{hyp}}$. As
$\gamma^{\theta}_{\operatorname{WP}}$ intersects this ball it follows that
$\gamma^{\theta}_{\operatorname{WP}}$ must intersect either
$\gamma^{v}_{\operatorname{WP}}$ or $\gamma^{w}_{\operatorname{WP}}$. Breaking
symmetry, let us assume that $\gamma^{\theta}_{\operatorname{WP}}$ intersects
$\gamma^{v}_{\operatorname{WP}}$. In this case the WP rays
$\gamma^{\theta}_{\operatorname{WP}}$ and $\gamma^{v}_{\operatorname{WP}}$
bound a bigon, contradicting Gauss-Bonnet.

So now suppose that for all $t_{n}$, the $R$-balls centred at
$\gamma^{\psi(\theta)}_{\operatorname{hyp}}(t_{n})$ intersect one of
$\gamma^{\psi(v)}_{\operatorname{hyp}}$ or
$\gamma^{\psi(w)}_{\operatorname{hyp}}$. Breaking symmetry, let us assume that
it is $\gamma^{\psi(v)}_{\operatorname{hyp}}$. As this implies the ray
$\gamma^{\psi(\theta)}_{\operatorname{hyp}}$ fellow travels the ray
$\gamma^{\psi(v)}_{\operatorname{hyp}}$ for all time, it forces
$\gamma^{\psi(\theta)}_{\operatorname{hyp}}=\gamma^{\psi(v)}_{\operatorname{hyp}}$.
It follows that $\gamma^{\theta}_{\operatorname{WP}}$ fellow travels
$\gamma^{v}_{\operatorname{WP}}$. As the WP metric is $\text{CAT}(-\kappa)$
for some $\kappa>0$ and $(\mathbb{D},d_{\operatorname{WP}})$ is a geodesic
metric space, it forces
$\gamma^{\theta}_{\operatorname{WP}}=\gamma^{v}_{\operatorname{WP}}$. But this
is a contradiction because $\theta\in I(v,w)$.

∎

Lemme 4.2 implies monotonicity for the map $\psi$. In summary, $\psi$ is a map
defined on a dense subset of $S^{1}$, has dense image, and preserves the
circular order. By a standard exercise, see [4, Lemma 7.3], $\psi$ extends to
a $S^{1}$-homeomorphism.

##  5\. Typical geodesics

Let $\ell$ be a Lebesgue measure on $S^{1}$. After identification with
$T_{p}^{1}(\mathbb{D},d_{\operatorname{WP}})$, we can let $\ell$ be the
conditional measure induced by the WP Liouville measure. Alternatively, $\ell$
could also be the conditional measure induced by the hyperbolic Liouville
measure because the two densities are absolutely continuous.

#### Proof of Theorem 1.5:

We derive Theorem 1.5 by comparing the growth of the total winding number
about the cusp along typical WP and hyperbolic geodesics. For either metric,
we define the total winding number till a particular recurrence time to be the
sum of the winding numbers (without sign) during cusp excursions till that
time.

The statistics of cuspidal winding numbers along typical WP geodesics for
general WP-type metrics were analysed in [8]. The analysis heavily leverages
the precise decay of correlations proved in [3] for the exponential mixing of
the WP flow.

For the total winding number comparison to be valid, we also need the time
changes between the WP and hyperbolic times to work out. To be precise, along
typical geodesics the WP and hyperbolic times at recurrence points need to be
asymptotically comparable. This is also established in [8] for general WP-type
metrics and we shall recall the relevant results at the appropriate steps in
the proof of Theorem 1.5. As the arguments below apply to mildly constrained
WP-type metrics, we keep the discussion on a a general footing as before.

By [8, Theorem 5.9], there exists a constant $K_{1}>1$ such that for almost
every $p$ and $\ell$-almost every $v\in
T_{p}^{1}(\mathbb{D},d_{\operatorname{WP}})$, the total winding number
$W_{v}(T)$ of the Weil–Petersson ray defined by $v$ around the cusp till time
$T$, satisfies

(5.1) $\frac{1}{K_{1}}T<W_{v}(T)<K_{1}T.$

Also by [8, Section 5.10], the hyperbolic distance measured along the
recurrence times of a recurrent Weil–Petersson geodesic ray increases linearly
in time: there exists a constant $K_{2}>1$ such that

(5.2) $\frac{1}{K_{2}}T<d_{hyp}(p,\gamma^{v}_{\operatorname{WP}}(T))<K_{2}T.$

In particular, this implies that when $T$ is a time of recurrence along
$\gamma^{v}_{\operatorname{WP}}$ the hyperbolic time $S$ along
$\gamma^{\psi(v)}_{\operatorname{hyp}}$ is comparable to the Weil–Petersson
time $T$. Let $W_{\psi(v)}(S)$ be the total winding number around the cusp
till time $S$ of the hyperbolic ray $\gamma^{\psi(v)}_{\operatorname{hyp}}$.
From Theorem 1.3 it follows that there exists a constant $K_{3}>1$ such that
when $S$ is the hyperbolic time that corresponds to a time of recurrence $T$
along $\gamma^{v}_{\operatorname{WP}}$, we have

(5.3) $\frac{1}{K_{3}}W_{v}(T)<W_{\psi(v)}(S)<K_{3}W_{v}(T).$

Putting inequalities 5.1, 5.2 and 5.3 together, we can conclude that there
exists a constant $K_{4}>1$ such that

(5.4) $\frac{1}{K_{4}}S<W_{\psi(v)}(S)<K_{3}S$

i.e., the total winding along the corresponding hyperbolic geodesic grows
linearly in hyperbolic time.

On the other hand, by [6, Theorem 1.6], for $\ell$-almost every
$v\in(T_{p}^{1}\mathbb{D},d_{\operatorname{hyp}})$, there is $S\log S$ lower
bound for the growth of total winding number till time $S$ along the
hyperbolic geodesic ray $\gamma^{v}_{\operatorname{hyp}}$. Thus, the
hyperbolic geodesic rays that are the images of typical WP rays under $\psi$
are atypical and this concludes the proof of Theorem 1.5.

###  5.5. Continued fractions:

We can now invoke Theorem 1.5 for the modular surface to give a concrete
description. Series showed that hyperbolic geodesic rays on the modular
surface are coded by the continued fraction expansion of their endpoint at
infinity. See [11]. The continued fraction coefficients track the
combinatorics of how the ray crosses the triangles in the Farey tessellation.
Namely, for each triangle that the geodesic cuts through one records right or
left depending on which of the resulting components is again a triangle. In
this way, one gets an infinite sequence of $R$s and $L$s. Breaking symmetry,
if the sequence has the form $R^{a_{1}}L^{a_{2}}\cdots$ then the number
$a_{j}$ is exactly the $j^{\text{th}}$ coefficient of the endpoint at
infinity. Suppose that $T$ is a recurrence time, $N(T)$ is the number of
excursions till $T$, and $H_{j}\,:\,1\leqslant j\leqslant N(T)$ is the
horoball for the $j^{\text{th}}$-excursion. We may choose upper half-space co-
ordinates $(x,y)$ so that $H_{j}$ is of the form $\\{y>c\\}$ for some constant
$c\geqslant 1$ and the fundamental domain intersected with $H_{j}$ is the
strip $0\leqslant x<1$. Then it readily follows that up to a multiplicative
constant that depends only on $c$, the distance between the entry and exits
points on $\partial H_{j}$ equals the continued fraction coefficient $a_{j}$.
See [6] for the details.

Using the fellow travelling result, namely Theorem 1.3, we can also record
continued fraction coefficients along a recurrent WP ray
$\gamma^{v}_{\operatorname{WP}}$ by tracking the coefficients along the
corresponding hyperbolic ray $\gamma^{\psi(v)}_{\operatorname{WP}}$. We could
also define coefficients using the number of cuspidal fundamental domains
crossed during each excursion of $\gamma^{v}_{\operatorname{WP}}$. Again by
Theorem 1.3, the two definitions produce the same coefficients that differ up
to a bounded additive error. So from now on, we will use the second
definition.

Consider a recurrence time $T$ for $\gamma^{v}_{\operatorname{WP}}$. Then the
segment
$[\gamma^{v}_{\operatorname{WP}}(0),\gamma^{v}_{\operatorname{WP}}(T)]$ has
completed finitely many excursions giving us a finite list of coefficients
$[a_{1},a_{2},\dots,a_{N(T)}]$. Note that there is a constant $K_{5}>1$ such
that for any recurrence time $T$,

$\frac{1}{K_{5}}W_{v}(T)<a_{1}+a_{2}+\dots+a_{N(T)}<K_{5}W_{v}(T)$

where $W_{v}(T)$ is the total cuspidal winding number till time $T$. By the
ergodicity of the Weil–Petersson flow and an argument similar to [6, Lemma
3.8], there exists a constant $K_{6}>1$ such that

$\frac{1}{K_{6}}T<N(T)<K_{6}T,$

that is, the number of excursions till time $T$ along a typical WP ray grow
linearly in $T$. In particular, combining the above facts about winding number
growth and the number of excursions with (5.1), we get that there exists a
constant $K_{7}>1$ such that along a typical WP ray
$\gamma^{v}_{\operatorname{WP}}$, we have

$\frac{1}{K_{7}}<\frac{a_{1}+a_{2}+\dots+a_{N(T)}}{N(T)}<K_{7}$

for all recurrence times $T$ that are large enough (how large depends on $v$).
In other words, the average coefficient along a typical WP ray (and hence
along the corresponding hyperbolic ray) is bounded. On the other hand, the
average coefficient along a typical hyperbolic ray is unbounded and exhibits
an interesting "trim-sum" property. See the introduction in [6].

## References

  * [1] Brock, J., Masur, H. and Minsky, Y. Asymptotics of Weil–Petersson geodesics I: ending laminations, recurrence and flows. Geom. Funct. Anal. 19 (2010), no. 5, 1229-1257.
  * [2] Brock, J. and Modami, B. Recurrent Weil–Petersson geodesic rays with non-uniquely ergodic ending laminations. Geom. Topol. 19 (2015), no. 6, 3565-3601.
  * [3] Burns K., Masur H., Matheus, C. and Wilkinson, A. Rates of mixing for the Weil-Petersson geodesic flow: Exponential mixing in exceptional moduli spaces. Geom. Funct. Anal. 27 (2017), no. 2, 240–288.
  * [4] Dovgoshey, O., Martio, O., Ryazanov, V. , Vuorinen, M. The Cantor function. Expo. Math. 24 (2006), no. 1, 1-37.
  * [5] Farb, B. Relatively hyperbolic groups. Geom. Funct. Anal. 8 (1998), no. 5, 810-840.
  * [6] Gadre, V. Partial sums of excursions along random geodesics and volume asymptotics for thin parts of moduli spaces of quadratic differentials. J. Eur. Math. Soc. 19 (2017), no. 10, 3053-3089.
  * [7] Gadre, V., Maher, J. and Tiozzo, G. Word length statistics and Lyapunov exponents for Fuchsian groups with cusps. New York J. Math 21 (2015), 511-531.
  * [8] Gadre, V. and Matheus, C. Cusp excursions of random geodesics in Weil-Petersson type metrics. to appear in J. Differential Geom.
  * [9] Pollicott, M. and Weiss, H. Ergodicity of the geodesic flow on non-complete negatively curved surfaces. Asian J. Math. 13 (2009), no. 3, 405–419.
  * [10] Pollicott, M., Weiss, H. and Wolpert, S. Topological dynamics of the Weil–Petersson geodesic flow. Adv. Math. 223 (2010), no. 4, 1225-1235.
  * [11] Series, C. The modular surface and continued fractions. J. London Math. Soc. (2) 31 (1985), no. 1, 69-80.
  * [12] Wolpert, S. Behaviour of geodesic length functions on Teichmüller space. J. Differential. Geom. 79 (2008), no. 2, 277-334.


#1 #1

11institutetext: Leibniz-Institut für Astrophysik Potsdam (AIP), An der
Sternwarte 16, D-14482 Potsdam, Germany
11email<EMAIL_ADDRESS>22institutetext: Niels Bohr Institute, University of
Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark

# Does the magneto-thermal instability survive whistler-suppression of thermal
conductivity in galaxy clusters?

Lorenzo Maria Perrone 11 Thomas Berlok 2211 Christoph Pfrommer 11

The hot and dilute intracluster medium (ICM) plays a central role in many key
processes that shape galaxy clusters. Nevertheless, the nature of plasma
turbulence and particle transport in the ICM remain poorly understood and
quantifying the effect of kinetic plasma instabilities on the macroscopic
dynamics represents an outstanding problem. Here we focus on the impact of
whistler-wave suppression of the heat flux on the magneto-thermal instability
(MTI), which is expected to drive significant turbulent motions in the
periphery of galaxy clusters. We perform small-scale Boussinesq simulations
with a sub-grid closure for the thermal diffusivity in the regime of whistler-
wave suppression. Our model is characterized by a single parameter that
quantifies the collisionality of the ICM on the astrophysical scales of
interest that we tune to explore a range appropriate for the periphery of
galaxy clusters. We find that the MTI is qualitatively unchanged for weak
whistler-suppression. Conversely, with strong suppression the magnetic dynamo
is interrupted and MTI-turbulence dies out. In the astrophysically relevant
limit, however, the MTI is likely to be supplemented by additional sources of
turbulence. Investigating this scenario, we show that the inclusion of
external forcing has a beneficial impact and revives even MTI simulations with
strong whistler-suppression. As a result, the plasma remains buoyantly
unstable, with important consequences for turbulent mixing in the ICM.

###### Key Words.:

Galaxies: clusters: intracluster medium – Instabilities – Magnetohydrodynamics
(MHD) – Plasmas – Turbulence – Methods: numerical

##  1 Introduction

Understanding the physics of transport processes in the hot, dilute ICM of
galaxy clusters represents an ongoing theoretical challenge. This has
important implications for the properties of turbulence and the feedback
processes (including from active galactic nuclei, Balbus & Reynolds, 2008;
Sharma et al., 2009; Reynolds et al., 2015; Yang & Reynolds, 2016a). One major
complication arises from the fact that high-$\beta$, weakly-collisional
plasmas are susceptible to kinetic microinstabilities that grow on scales
observationally unresolved (Schekochihin & Cowley, 2006; Kunz et al., 2014).
These instabilities can significantly alter particle motions, affecting
thermal transport (Riquelme et al., 2016; Komarov et al., 2016, 2018), plasma
viscosity (Zhuravleva et al., 2019), particle acceleration (Sironi & Narayan,
2015) and magnetic field amplification (St-Onge et al., 2020). In the presence
of a heat flux (generated, e.g., by a temperature gradient), the whistler
instability grows and saturates by regulating the scattering of electrons via
magnetic fluctuations, reducing their speed along magnetic field lines to
values well below the thermal velocity (Levinson & Eichler, 1992; Pistinner &
Eichler, 1998). As a result, the heat flux can be significantly suppressed
below the collisional (Spitzer) value, depending on the plasma $\beta$
(Riquelme et al., 2016; Roberg-Clark et al., 2016, 2018; Komarov et al., 2018;
Drake et al., 2021). Possible consequences of this suppression are far-
reaching, and include modifications to the Field length (Drake et al., 2021),
and inhibition of plasma instabilities that rely on efficient thermal
conduction along magnetic fields, such as the MTI (Balbus 2000; Kunz 2011).
This instability may be active in the periphery of galaxy clusters and may
contribute to the observed levels of turbulence (Parrish & Stone, 2007;
McCourt et al., 2011).

In previous work, it was shown that the MTI produces a state of turbulence
sustained by the background temperature gradient and composed of density and
velocity fluctuations over a wide range of scales, with rms values at
saturation that follow clear power laws with the thermal diffusivity and the
gradients of gas entropy and temperature (Perrone & Latter, 2022a, b,
hereafter ; , respectively). Through these scaling laws it was estimated that
the MTI is capable of driving turbulent motions of hundreds kilometers per
second, on scales of $100\,$\mathrm{k}\mathrm{p}\mathrm{c}$$, roughly
consistent with existing measurements of turbulence (Hitomi Collaboration et
al. 2018, see Kempf et al. 2023 for a detailed discussion on the observational
constraints on MTI-turbulence).

One of the main uncertainties of these estimates is our still limited
understanding of the impact of microscale instabilities on the macroscopic
heat conduction. A self-consistent model of this impact would require plasma
kinetic simulations spanning these very disparate scales. Because of the
severe computational costs involved, a popular approach is to instead adopt
sub-grid models for the transport coefficient informed by kinetic theory and
simulations (Sharma et al., 2006; Kunz et al., 2012; Foucart et al., 2017;
Kempski et al., 2019; Beckmann et al., 2022; Squire et al., 2023).

Regarding the MTI, studies of its nonlinear saturation with suppression of
thermal conduction by the mirror instability revealed very little difference
compared to the unsuppressed case (Berlok et al., 2021). On the other hand,
the impact of electron scattering at whistler waves could potentially be more
severe (Drake et al., 2021) and remains unexplored. In this work, and in its
companion paper (Perrone et al., 2023, hereafter PBP23), we address this issue
by investigating the survival of the MTI when the plasma is subject to the
expected whistler-suppression of its heat conductivity. We do this using
controlled small-scale simulations. We show that the MTI survives unscathed
for the mild levels of suppression expected in the periphery of galaxy
clusters. With stronger suppression, the MTI instead enters a “spiral-of-
death” scenario. Nevertheless, we show that the presence of external forcing
(the result of, e.g., mergers or substructure accretion) revives the MTI,
rather than further obstructing it. Our findings support the idea that the MTI
could be active in the periphery of galaxy clusters even with suppression of
thermal conduction.

Figure 1: Part (a): 1D simulation illustrating the behaviour of whistler-
suppressed thermal conduction. The temperature (top), plasma $\beta$ (middle
row, kept fixed) and thermal conductivity (bottom row) are shown at different
times (in units of $L^{2}/\chi_{0}$, with $L$ the box size). For this 1D
example the suppression factor has the form
$\chi/\chi_{0}=(1+L\beta\mathrm{d}\ln T/\mathrm{d}x)^{-1}$ such that heat
conduction is suppressed at strong temperature gradients and high $\beta$.
Temperature equilibration thus occurs faster on the right (where $\beta$ is
low) than on the left (where $\beta$ is high). Part (b): snapshots at
saturation of 2D MTI-turbulence with a similar form of whistler-suppression
(Eq. 3). We show temperature fluctuations (top), the Boussinesq equivalent of
plasma beta, $\tilde{\beta}$, and the spatially varying heat conductivity,
$\chi$. The suppression of $\chi$ is seen to spatially correlate with regions
of high $\tilde{\beta}$. The MTI-driven temperature fluctuations reveal that
the MTI continues to work despite the whistler-suppression of the thermal
diffusivity, which amounts to $20\%$ of $\chi_{0}$ (volume-average).

##  2 Model

We perform local simulations of the MTI with the pseudospectral code snoopy
(Lesur, 2015). Similar to PL22a, we adopt the Boussinesq approximation to
study small-scale, subsonic turbulence, assuming a constant negative
background temperature gradient (temperature decreases radially) and a
positive entropy stratification that is stable in the absence of anisotropic
heat conduction. The two gradients define respectively the MTI
($\omega_{\mathrm{T}}$) and Brunt-Väisälä frequency ($N$, see Appendix A.1),
whose values in the periphery of galaxy clusters are roughly comparable and
approximately $(600\,$\mathrm{M}\mathrm{y}\mathrm{r}$)^{-1}$ (PL22a). In our
local Cartesian setup, the radial direction coincides with $\bm{e}_{z}$. The
full MHD Boussinesq equations with anisotropic thermal conduction can be found
in Appendix A.1.

###  2.1 Whistler-suppression of thermal conduction

We employ the following closure for the parallel heat flux in the regime of
electron scattering by whistler waves (Komarov et al., 2018)

$\displaystyle
q_{\parallel}=\frac{q_{\parallel,\mathrm{s}}}{1+\frac{1}{3}\beta_{\mathrm{e}}\lambda_{\text{mfp,e}}/L_{\mathrm{T},\parallel}},$
(1)

where $\beta_{\mathrm{e}}$ is the electron plasma beta (hereafter $\beta$ for
brevity), $\lambda_{\text{mfp,e}}$ is the electron mean-free path, and
$L_{\mathrm{T},\parallel}$ the temperature scale along the magnetic field. In
the periphery of galaxy clusters typical values are $\beta\sim 100$,
$\lambda_{\text{mfp,e}}\sim 10\,$\mathrm{k}\mathrm{p}\mathrm{c}$$, and
$L_{\mathrm{T},\parallel}\sim 100\,$\mathrm{k}\mathrm{p}\mathrm{c}$$. Equation
(1) smoothly interpolates between the Spitzer collisional heat flux
($q_{\parallel,\mathrm{s}}\simeq
0.5n_{\mathrm{e}}m_{\mathrm{e}}\varv_{\text{th,e}}^{3}\lambda_{\text{mfp,e}}/L_{\mathrm{T},\parallel}$;
Spitzer, 1962; Cowie & McKee, 1977, where $n_{\mathrm{e}}$ is the electron
number density, $m_{\mathrm{e}}$ is the electron mass, and
$\varv_{\text{th,e}}$ is the electron thermal speed) at low $\beta$, and the
marginal heat flux controlled by the whistler instability at saturation when
$\beta$ is large ($q_{\parallel,\mathrm{w}}\simeq
1.5\beta^{-1}n_{\mathrm{e}}m_{\mathrm{e}}\varv_{\text{th,e}}^{3}$; Komarov et
al., 2018). Since in the periphery of galaxy clusters
$\lambda_{\text{mfp,e}}/L_{\mathrm{T},\parallel}\ll 1$, in our model we do not
include saturation of the heat flux when the electron mean-free path is larger
than the temperature scale (Cowie & McKee, 1977).

We show the principle of whistler-suppression of thermal conductivity through
a simplified one-dimensional simulation of thermal diffusion (Fig. 1a), where
we evolve the temperature equation with anisotropic diffusion in the presence
of a fixed background magnetic field that has a spatially varying field
strength. As time progresses, sharp temperature gradients survive where
$\beta$ is high but the gradients are gradually erased where the magnetic
field is stronger. On the other hand, steep temperature gradients can be
generated at the interface between regions of high and low $\beta$ where none
was present before (see Fig. 1a).

###  2.2 Whistler-suppression with Boussinesq

To model whistler-suppression of the heat flux in the Boussinesq approximation
we introduce the analogue of the usual plasma beta, $\tilde{\beta}$, as

$\displaystyle\tilde{\beta}=\frac{\varv_{\mathrm{cond}}^{2}}{B^{2}/(8\pi\rho)},$
(2)

($\rho$ is the mass density of the plasma), which we define through the so-
called “conduction speed” $\varv_{\mathrm{cond}}$, rather than the thermal
speed, since the latter is ordered out of the Boussinesq equations. The
conduction speed is given by
$\varv_{\mathrm{cond}}=(\chi_{0}\omega_{\mathrm{T}})^{1/2}$, where $\chi_{0}$
is the Spitzer diffusivity. In typical galaxy clusters,
$\varv_{\mathrm{cond}}\simeq 500\,$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$ and
$\tilde{\beta}\simeq 0.5\beta$. With this substitution, the whistler-
suppressed thermal diffusivity associated with Eq. (1) can be written as (see
Appendix A.2)

$\displaystyle\chi=\frac{\chi_{0}}{1+\frac{1}{3}\tilde{\beta}\alpha\,\Big{|}b_{z}+\omega_{\mathrm{T}}^{-2}\bm{b}\bm{\cdot}\bm{\nabla}\theta\Big{|}},$
(3)

where $\theta$ is the Boussinesq buoyancy variable (proportional to
temperature fluctuation but with opposite sign, PL22a), and $\alpha$ is a
dimensionless parameter that absorbs all the information on the mean-free path
and the typical cluster scales, which is defined as

$\displaystyle\alpha\equiv\frac{1}{\gamma}\left(\frac{\lambda_{\text{mfp,e}}}{H}\right)\frac{H^{2}}{\chi_{0}/\omega_{\mathrm{T}}},$
(4)

where $\gamma$ is the adiabatic index and $H$ is the typical pressure scale-
height of the cluster. By varying $\alpha$ within the range of realistic
values expected in the ICM (see Section 2.3), we can simulate MTI-turbulence
with stronger or weaker whistler-suppression of thermal diffusion. Finally,
the last quantity in the denominator of Eq. (3) represents the projection of
the temperature gradient along the local direction of the magnetic field
$\bm{b}$. Unless otherwise noted, in what follows we take
$\omega_{\mathrm{T}}^{-1}$ as the unit of time and normalize all velocities to
$\varv_{\mathrm{cond}}$. Similarly, we non-dimensionalize $\theta$ as
$\tilde{\theta}\equiv\theta/(N\varv_{\text{cond}})$.

###  2.3 Parameters in the periphery of galaxy clusters

To estimate the value of the $\alpha$ parameter in the periphery of galaxy
clusters we use Eq. (4) together with the definition of the Spitzer
diffusivity and of the electron collisional mean-free path. This yields

$\displaystyle\alpha\simeq
0.01\left(\frac{H}{300\,$\mathrm{k}\mathrm{p}\mathrm{c}$}\right)\left(\frac{\omega_{\mathrm{T}}}{(700\,$\mathrm{M}\mathrm{y}\mathrm{r}$)^{-1}}\right)\left(\frac{T_{\mathrm{e}}}{5\,$\mathrm{k}\mathrm{e}\mathrm{V}$}\right)^{-1/2},$
(5)

where $H$, $\omega_{\mathrm{T}}$ and $T_{\mathrm{e}}$ are measured in units of
$\mathrm{k}\mathrm{p}\mathrm{c}$, $\mathrm{M}\mathrm{y}\mathrm{r}^{-1}$ and
$\mathrm{k}\mathrm{e}\mathrm{V}$, respectively. It is important to note that
these quantities can vary substantially both across different clusters and
within the cluster periphery. Nevertheless, an order-of-magnitude estimate can
be obtained assuming that clusters are in hydrostatic equilibrium, and taking
power-law pressure and temperature profiles (justified in our local model).
For the power-law exponents we take the best-fit profiles for the X-COP sample
reported by Ghirardini et al. (2019). With these assumptions, and taking a
fiducial radius (where these quantities are computed) of the order of several
hundred $\mathrm{k}\mathrm{p}\mathrm{c}$, we find that the range of $\alpha$
allowed by observations is 0.01–0.05, where the lower estimate is obtained
from Eq. (5) with the values given, while the higher estimate is calculated
assuming $H=600\,$\mathrm{k}\mathrm{p}\mathrm{c}$$,
$\omega_{\mathrm{T}}=(600\,$\mathrm{M}\mathrm{y}\mathrm{r}$)^{-1}$ and
$T_{\mathrm{e}}=1\,$\mathrm{k}\mathrm{e}\mathrm{V}$$.

##  3 Suppression of turbulence and demise of the MTI

We run two- and three-dimensional simulations to study the impact of whistler-
suppressed thermal conduction on the MTI (i.e. $\alpha\neq 0)$, and compare
the results against our reference run ($\alpha=0)$.

In Fig. 1b we show a snapshot of a representative 2D MTI run with whistler-
suppression at saturation. The simulation has been initialized with a weak
uniform magnetic field in the horizontal direction and no external forcing
(see Appendix A.3 for further details). We note that despite the significant
suppression of thermal diffusivity, the character of MTI-turbulence remains
qualitatively the same: hot plumes rise across the box, cold plumes sink. The
shape of the plumes traces the morphology of the magnetic field. The volume
average of $\chi/\chi_{0}$ is less than $20\%$, but the spatial distribution
is highly inhomogeneous and correlates strongly with $\tilde{\beta}$.
Conversely, the correlation with $L_{\mathrm{T},\parallel}$ (not shown) is
weaker.

Figure 2: Time series of the volume-averaged specific kinetic energy (top
row), specific magnetic energy (middle row) and thermal diffusivity of 3D
simulations. The left column shows our MTI runs without external forcing:
without suppression (gray starred line), with weak suppression (solid blue),
and with strong suppression (solid gold). In the right panel we introduce
external forcing. For comparison we also show a pure-forcing run (solid red
line) with isotropic thermal diffusion (no MTI). MTI-turbulence decays with
strong whistler-suppression, but the addition of forcing revives it.

A limitation of 2D runs is that the magnetic field strength (and thus the
overall suppression of thermal conduction) is somewhat arbitrary, since they
lack a small-scale dynamo that pins the magnetic energy at the equipartition
value. As we will now show, in 3D the presence of the small-scale dynamo (or
lack thereof) plays a crucial role in the whistler-suppression of MTI-
turbulence. In the left column of Fig. 2 we compare the time evolution of the
volume-averaged specific kinetic ($K\equiv\langle u^{2}\rangle/2$) and
magnetic energy ($M\equiv\langle B^{2}\rangle/(8\pi\rho_{0})$) as well as the
suppression factor $\chi/\chi_{0}$ of a reference 3D MTI run ($\alpha=0$), and
of two 3D runs with weak ($\alpha=0.012$) and strong ($\alpha=0.12$)
suppression. Note that here “weak” and “strong” roughly correspond to the
lower and higher ends of the range of whistler-suppression that can be
expected in the periphery of galaxy clusters (see Section 2.3). As initial
conditions for these latter runs, we use a turbulent snapshot of the reference
run at saturation, and we then turn on whistler-suppression (see Appendix A.3
for the numerical details). The initial spike in the kinetic energy for the
reference MTI run results from the disruption of the fastest-growing MTI mode
(seeded from the initial random perturbations) by parasitic Kelvin-Helmholtz
instabilities (PL22b). Also note that in Fig. 2 magnetic energy reaches
approximate equipartition with the turbulent kinetic energy ($M\sim 0.25K$).
In terms of the modified plasma $\beta$ this approximately corresponds to
$\tilde{\beta}\simeq 25$ (note the equivalence
$M/\varv_{\mathrm{cond}}^{2}=\tilde{\beta}^{-1}$). Converting to the actual
plasma $\beta$ this gives $\beta\sim 50$, depending on the values of
$\varv_{\mathrm{cond}}$ and the sound speed, which is on the lower range of
expected cluster values ($\beta\sim 50-100$ ; Carilli & Taylor, 2002). We
stress, however, that the exact value of the constant of proportionality (and
thus of $\beta$) between kinetic and magnetic energies depends on the
numerical parameters of the simulation.

With weak whistler-suppression, MTI-turbulence carries on, albeit at lower
levels proportional to the average level of suppression ($\simeq 82\%$).
Similarly, the magnetic energy saturates slightly below the reference run.
With stronger suppression, we witness a radically different scenario: here the
kinetic and magnetic energies sharply decrease, with a near-total suppression
of the thermal diffusivity. The physical picture behind the demise of the MTI
can be described as follows: because of the relation of proportionality
between thermal diffusivity and turbulent kinetic energy (; ), a strong
initial suppression of $\chi$ (and therefore a reduction in $K$) means that
magnetic fields cannot be maintained by the small-scale dynamo at the same
equipartition-level as before and that they decay to adjust to the new
turbulence levels. However, this decay leads to further suppression of the
thermal diffusivity through its $\beta$-dependence, which in turns weakens the
turbulence, and so on. The results of this self-reinforcing cycle is a runaway
process that kills the MTI entirely. Note that this runaway process only
occurs in the simulation with strong suppression ($\alpha=0.12$), while in our
weak suppression run ($\alpha=0.012$) the system adjusts to the new
equilibrium without any hint of runaway. We attribute this behavior to the
impairment of the small-scale dynamo which, for MTI-turbulence, requires a
magnetic Reynolds number $\mathrm{Rm}$ (defined as
$\mathrm{Rm}\equiv\mathrm{Re}\times\mathrm{Pm}$, where $\mathrm{Re}$ is the
Reynolds number at the integral scale and $\mathrm{Pm}$ the ratio between the
viscosity and magnetic resistivity) above a critical value of
$\mathrm{Rm}_{\mathrm{c}}\simeq 35$ (PL22b). In our 3D reference run
$\mathrm{Rm}\simeq 77$, well above the critical threshold, with
$\mathrm{Pm}=4$ (see Appendix 2.3 for further details on the numerics). In our
weak suppression run the magnetic Reynolds number remains above the critical
value throughout ($\mathrm{Rm}\simeq 75$), while in our strong suppression run
$\mathrm{Rm}$ quickly drops below $\mathrm{Rm_{c}}$ (around $t=230$ in Fig. 2)
and the small-scale dynamo shuts down. The existence of a critical transition
between the two outcomes is explored in detail in PBP23.

##  4 Revival of the MTI by external turbulence

The results of the previous section may cast a cloud on the contribution of
the MTI to turbulent driving in the periphery of galaxy clusters. However,
other processes are expected to play a role, such as mergers and accretion
events (Simionescu et al., 2019). Even in the event of total suppression of
the MTI, turbulence injected through these other processes can also amplify
magnetic fields and, in so doing, restore thermal conductivity to a
significant fraction of the Spitzer value. This could make it possible for the
MTI to spring back into action.

To assess the feasibility of this scenario, we look at the effect of external
turbulence on the whistler-suppressed MTI by running a number of simulations
with delta-correlated in time, solenoidal, white-noise forcing. We tune energy
injection by the external forcing so as to drive turbulence of comparable
amplitude and at similar scale to that of our reference MTI run. This is
consistent with previous estimates of MTI-turbulence in the periphery of
galaxy clusters (PL22b).

We show in Fig. 2 (right column) the kinetic energy, magnetic energy and the
mean thermal diffusivity of three representative 3D runs before and after
enabling external forcing. These runs have been chosen to illustrate the full
range of possibilities between (i) reference MTI with the addition of external
forcing, and (ii) external forcing only with isotropic thermal conduction (and
thus without the MTI). In between them, the MTI run with strong whistler
suppression and forcing. The addition of forcing results in an increase in the
turbulent kinetic energy of the reference MTI run and an even more significant
increase in the magnetic energy. Similarly, compared to the forcing-only run,
the presence of the MTI clearly shows up as stronger turbulent fluctuations.

It is the run with strong whistler-suppression that exhibits the most dramatic
change: after turning on external forcing, the thermal diffusivity of this
once-decaying run rapidly increases, and the turbulence reaches levels above
those of the forcing-only run, a clear sign that the MTI has been revived.
Eventually, the thermal diffusivity stabilizes at around $64\%\chi_{0}$.

The impact of external forcing on the MTI is particularly evident if we plot
the time evolution of the volume-averaged buoyancy power
$\varepsilon_{\mathrm{b}}\equiv-\langle\theta u_{z}\rangle$ (Fig. 3), i.e. the
rate at which kinetic energy is added or removed by the buoyancy force,
obtained by dotting the momentum equation (Eq. 8) with $\bm{u}$. As $\theta$
is proportional to the density fluctuation, the quantity is positive if there
is an anticorrelation between $\theta$ and $u_{z}$ which is the case when
lighter fluid travels upwards, and denser fluid sinks. A net-positive value of
the buoyancy power thus indicates a buoyantly unstable flow, as in thermal
convection and in the case of reference MTI (where $\varepsilon_{\mathrm{b}}$
is proportional to the fraction of energy drawn from the background
temperature gradient converted into velocity fluctuations by the MTI). Stably-
stratified flows instead have a net-negative buoyancy power, meaning that an
initial perturbation will be damped. Whether the flow is buoyantly stable or
unstable has important implications for the efficiency of plasma and metal
mixing in the cool cores and in the bulk ICM (Sharma et al., 2009; Ruszkowski
& Oh, 2010; Yang & Reynolds, 2016b; Kannan et al., 2017).

The addition of external turbulence clearly interferes with the MTI, which
becomes less efficient at converting thermal energy into turbulence even in
the case of no whistler-suppression (in the reference MTI run
$\varepsilon_{\mathrm{b}}$ is reduced by about $50\%$ after turning on
forcing). Remarkably, the run with strong whistler-suppression goes from
negligible energy injection from buoyancy before forcing (because of the
decaying MTI) to significant energy injection as the external forcing boosts
thermal diffusion and revives the MTI (about $1/2$ of the reference run with
forcing, though the exact figure depends on the value of $\alpha$). This is to
be contrasted with the pure-forcing run which has isotropic thermal diffusion
– and therefore no MTI – where the buoyancy power has an overall stabilizing
effect.

Figure 3: Evolution of the volume-averaged buoyancy power without (left
midplane) and with external forcing (right midplane) in units of
$\varepsilon_{0}\equiv\chi_{0}\omega_{\mathrm{T}}^{4}/N^{2}$, the MTI energy
injection rate (PL22a). The 3D runs shown are the same as in Fig. 2. Despite
strong whistler suppression, the system remains buoyantly unstable thanks to
external forcing.  Figure 4: Shown is the buoyancy power in units of the MTI
energy injection rate for two 3D runs (reference MTI and strongly whistler-
suppressed MTI) before and after turning on external forcing. The volume
average is shown in Fig. 3. External forcing interferes with MTI-turbulence
but the upwelling and downwelling of the plasma is still present (regions in
red).

##  5 Discussion and conclusions

In this work we explored the impact of whistler-suppression of the heat flux
on the MTI through small-scale simulations with a subgrid model for the
thermal diffusivity inspired by kinetic theory and particle-in-cell
simulations. Our main findings are that the impact can be more or less severe
depending on the level of collisionality in the periphery of galaxy clusters.
For weak whistler-suppression, the main features of MTI-driven turbulence are
qualitatively unchanged, with only an overall decrease in the turbulent
fluctuations. For strong whistler-suppression, the MTI loses the ability to
drive a small-scale dynamo and decays. These two extremes span the range which
can be realistically expected in the periphery of galaxy clusters, and it is
thus not possible to rule out a priori whether the MTI would survive on its
own.

The MTI, however, is not the only source of turbulence in the bulk of the ICM.
As we have shown – and contrary to common intuition – these other processes
can have a beneficial effect on the MTI when thermal conduction is suppressed
by the whistler instability. The basic mechanism is summarized in Fig. 4 where
we show snapshots of the buoyancy power for the runs without suppression
($\alpha=0$) and with strong suppression ($\alpha=0.12$) before and after
turning on external forcing. The bottom-left corner shows the standard MTI
without whistler-suppression and external forcing. Turning on strong whistler-
suppression (bottom-right) stalls the MTI-driven small-scale dynamo and as a
result, turbulence dies out. On the other hand, external forcing without
whistler-suppression evidently interferes with the basic workings of the MTI
(top-left), but crucially the convective character of the instability is not
erased. This leads to the main new finding of this work, which is that
external turbulence is able to revive “dead” MTI simulations with strong
whistler-suppression, and that the flow remains overall buoyantly unstable
(top-right). In galaxy clusters, this change in buoyant response leads to an
increased efficiency of turbulent mixing, which could have important
astrophysical implications. For instance, it could explain the increased metal
mixing rates (and more efficient coupling of AGN to the ambient ICM) that was
found in cosmological simulations with anisotropic conduction (Kannan et al.,
2017; Barnes et al., 2019). As we have shown – even with whistler-suppression
– the convectively-unstable nature of the ICM seems to be qualitatively
unchanged. Our work illustrates the need for detailed modelling of the weakly-
collisional outskirts of galaxy clusters and sets the stage for the
development of effective models of turbulence in large-scale cosmological
simulations.

###### Acknowledgements.

The authors acknowledge support by the European Research Council under ERC-AdG
grant PICOGAL-101019746. This project has received support from the European
Union’s Horizon Europe research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 101106080.

## References

  * Balbus (2000) Balbus, S. A. 2000, The Astrophysical Journal, 534, 420
  * Balbus & Reynolds (2008) Balbus, S. A. & Reynolds, C. S. 2008, The Astrophysical Journal, 681, L65
  * Barnes et al. (2019) Barnes, D. J., Kannan, R., Vogelsberger, M., et al. 2019, Monthly Notices of the Royal Astronomical Society, 488, 3003
  * Beckmann et al. (2022) Beckmann, R. S., Dubois, Y., Pellissier, A., Polles, F. L., & Olivares, V. 2022, Astronomy & Astrophysics, 666, A71
  * Berlok et al. (2021) Berlok, T., Quataert, E., Pessah, M. E., & Pfrommer, C. 2021, Monthly Notices of the Royal Astronomical Society, 504, 3435
  * Carilli & Taylor (2002) Carilli, C. L. & Taylor, G. B. 2002, Annual Review of Astronomy and Astrophysics, 40, 319
  * Cowie & McKee (1977) Cowie, L. L. & McKee, C. F. 1977, The Astrophysical Journal, 211, 135
  * Drake et al. (2021) Drake, J. F., Pfrommer, C., Reynolds, C. S., et al. 2021, The Astrophysical Journal, 923, 245
  * Foucart et al. (2017) Foucart, F., Chandra, M., Gammie, C. F., Quataert, E., & Tchekhovskoy, A. 2017, Monthly Notices of the Royal Astronomical Society, 470, 2240
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, Astronomy & Astrophysics, 621, A41
  * Hitomi Collaboration et al. (2018) Hitomi Collaboration, Aharonian, F., Akamatsu, H., et al. 2018, Publications of the Astronomical Society of Japan, 70
  * Kannan et al. (2017) Kannan, R., Vogelsberger, M., Pfrommer, C., et al. 2017, The Astrophysical Journal, 837, L18
  * Kempf et al. (2023) Kempf, J. M., Rincon, F., & Clerc, N. 2023, Dynamical Properties and Detectability of the Magneto-Thermal Instability in the Intracluster Medium
  * Kempski et al. (2019) Kempski, P., Quataert, E., Squire, J., & Kunz, M. W. 2019, Monthly Notices of the Royal Astronomical Society, 486, 4013
  * Komarov et al. (2018) Komarov, S., Schekochihin, A., Churazov, E., & Spitkovsky, A. 2018, Journal of Plasma Physics, 84, 905840305
  * Komarov et al. (2016) Komarov, S. V., Churazov, E. M., Kunz, M. W., & Schekochihin, A. A. 2016, Monthly Notices of the Royal Astronomical Society, 460, 467
  * Kunz (2011) Kunz, M. W. 2011, Monthly Notices of the Royal Astronomical Society, 417, 602
  * Kunz et al. (2012) Kunz, M. W., Bogdanovic, T., Reynolds, C. S., & Stone, J. M. 2012, The Astrophysical Journal, 754, 122
  * Kunz et al. (2014) Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014, Physical Review Letters, 112, 205003
  * Lesur (2015) Lesur, G. 2015, Snoopy: General Purpose Spectral Solver
  * Levinson & Eichler (1992) Levinson, A. & Eichler, D. 1992, The Astrophysical Journal, 387, 212
  * McCourt et al. (2011) McCourt, M., Parrish, I. J., Sharma, P., & Quataert, E. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1295
  * Parrish & Stone (2007) Parrish, I. J. & Stone, J. M. 2007, The Astrophysical Journal, 664, 135
  * Perrone et al. (2023) Perrone, L. M., Berlok, T., & Pfrommer, C. 2023, in preparation
  * Perrone & Latter (2022a) Perrone, L. M. & Latter, H. 2022a, Monthly Notices of the Royal Astronomical Society, 513, 4605
  * Perrone & Latter (2022b) Perrone, L. M. & Latter, H. 2022b, Monthly Notices of the Royal Astronomical Society, 513, 4625
  * Pistinner & Eichler (1998) Pistinner, S. L. & Eichler, D. 1998, Monthly Notices of the Royal Astronomical Society, 301, 49
  * Reynolds et al. (2015) Reynolds, C. S., Balbus, S. A., & Schekochihin, A. A. 2015, The Astrophysical Journal, 815, 41
  * Riquelme et al. (2016) Riquelme, M., Quataert, E., & Verscharen, D. 2016, The Astrophysical Journal, 824, 123
  * Roberg-Clark et al. (2016) Roberg-Clark, G. T., Drake, J. F., Reynolds, C. S., & Swisdak, M. 2016, The Astrophysical Journal, 830, L9
  * Roberg-Clark et al. (2018) Roberg-Clark, G. T., Drake, J. F., Reynolds, C. S., & Swisdak, M. 2018, Physical Review Letters, 120, 035101
  * Ruszkowski & Oh (2010) Ruszkowski, M. & Oh, S. P. 2010, The Astrophysical Journal, 713, 1332
  * Schekochihin & Cowley (2006) Schekochihin, A. A. & Cowley, S. C. 2006, Physics of Plasmas, 13, 056501
  * Sharma et al. (2009) Sharma, P., Chandran, B. D. G., Quataert, E., & Parrish, I. J. 2009, The Astrophysical Journal, 699, 348
  * Sharma et al. (2006) Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. M. 2006, The Astrophysical Journal, 637, 952
  * Simionescu et al. (2019) Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Science Reviews, 215, 24
  * Sironi & Narayan (2015) Sironi, L. & Narayan, R. 2015, The Astrophysical Journal, 800, 88
  * Spitzer (1962) Spitzer, L. 1962, Physics of Fully Ionized Gases
  * Squire et al. (2023) Squire, J., Kunz, M. W., Arzamasskiy, L., et al. 2023, Pressure Anisotropy and Viscous Heating in Weakly Collisional Plasma Turbulence
  * St-Onge et al. (2020) St-Onge, D. A., Kunz, M. W., Squire, J., & Schekochihin, A. A. 2020, Journal of Plasma Physics, 86, 905860503
  * Yang & Reynolds (2016a) Yang, H.-Y. K. & Reynolds, C. S. 2016a, The Astrophysical Journal, 829, 90
  * Yang & Reynolds (2016b) Yang, H.-Y. K. & Reynolds, C. S. 2016b, The Astrophysical Journal, 818, 181
  * Zhuravleva et al. (2019) Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2019, Nature Astronomy, 3, 832

##  Appendix A Methods

###  A.1 Boussinesq approximation to MHD

In the Boussinesq approximations, we evolve small-scale, subsonic fluctuations
on top of a stationary background in hydrostatic equilibrium (denoted by the
subscript “$0$”, i.e., $\rho_{0}$ for the mass density, $T_{0}$ for the
temperature and so on), which represents a local volume of a galaxy cluster at
a reference radius $R_{0}$. The variables evolved by snoopy are the velocity
fluctuation $\bm{u}$, the magnetic field $\bm{B}$ (rescaled by
$\sqrt{4\pi\rho_{0}}$ so that it has the dimensions of a velocity), and the
buoyancy variable

$\displaystyle\theta\equiv g_{0}\frac{\delta\rho}{\rho_{0}},$ (6)

which is given in units of $\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$, where
$g_{0}$ is the gravitational acceleration. The appropriate equations are
derived in PL22a, and are shown below including a spatially-varying
anisotropic thermal diffusivity and an external acceleration $\bm{f}$ in the
momentum equation:

$\displaystyle\bm{\nabla}\bm{\cdot}\bm{u}=\bm{\nabla}\bm{\cdot}\bm{B}$
$\displaystyle=0,$ (7)
$\displaystyle\left(\partial_{t}+\bm{u}\bm{\cdot}\bm{\nabla}\right)\bm{u}$
$\displaystyle=-\frac{\bm{\nabla}p_{\text{tot}}}{\rho_{0}}-\theta\bm{e}_{z}+\left(\bm{B}\bm{\cdot}\bm{\nabla}\right)\bm{B}+\nu\bm{\nabla}^{2}\bm{u}+\bm{f},$
(8) $\displaystyle\left(\partial_{t}+\bm{u}\bm{\cdot}\bm{\nabla}\right)\bm{B}$
$\displaystyle=\left(\bm{B}\bm{\cdot}\bm{\nabla}\right)\bm{u}+\eta\bm{\nabla}^{2}\bm{B},$
(9) $\displaystyle\left(\partial_{t}+\bm{u}\bm{\cdot}\bm{\nabla}\right)\theta$
$\displaystyle=N^{2}u_{z}+\bm{\nabla}\bm{\cdot}\left[\chi\bm{b}\left(\bm{b}\bm{\cdot}\bm{\nabla}\right)\theta\right]+\omega_{\mathrm{T}}^{2}\bm{\nabla}\bm{\cdot}\left(\chi\bm{b}b_{z}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ -(-1)^{n}\nu_{\text{th}}^{n}\bm{\nabla}^{2n}\theta,$ (10)

where $p_{\text{tot}}$ is the sum of thermal and magnetic pressures; $\nu$ and
$\eta$ are the fluid viscosity and magnetic resistivity, while in the buoyancy
equation (Eq. 10) the thermal diffusivity $\chi$ absorbs a factor of
$(\gamma-1)/\gamma$. The (constant in $z$) squared MTI and Brunt-Väisälä
frequencies that appear in Eq. (10) are defined through

$\displaystyle N^{2}$
$\displaystyle=\frac{g_{0}}{\gamma}\left.\frac{\partial(\ln
p\rho^{-\gamma})}{\partial
R}\right|_{0}=\left.-\frac{1}{\gamma\rho}\frac{\partial\ln p}{\partial
R}\frac{\partial\ln(p\rho^{-\gamma})}{\partial R}\right|_{0},$ (11)
$\displaystyle\omega_{\mathrm{T}}^{2}$
$\displaystyle=\left.-g_{0}\frac{\partial\ln T}{\partial
R}\right|_{0}=\left.\frac{1}{\rho}\frac{\partial\ln p}{\partial
R}\frac{\partial\ln T}{\partial R}\right|_{0}$ (12)

and are computed at the reference radius $R_{0}$. Note that in the buoyancy
equation we have also included a higher-order diffusion operator (where
$\nu_{\text{th}}^{n}$ denotes the corresponding hyperdiffusivity, not to be
confused with $\nu$, the fluid viscosity) to regularize the smallest scales in
our simulation in the event of strong spatial variations in the whistler-
suppressed thermal diffusivity. We opted for a third-order ($n=3$) hyper-
diffusion to allow for a wide scale separation between energy injection by the
MTI and dissipation. The introduction of hyperdiffusivity results in a lower
overall energy of the buoyancy fluctuations (for example in our 3D reference
run with hyperdiffusion the thermal energy is $\lesssim 7\%$ less than an
equivalent run without hyperdiffusion, while the kinetic and magnetic energies
are practically unchanged) but – provided this suppression happens near the
smallest scales of the box – the resulting dynamics in spectral space is
mostly unaffected.

###  A.2 Whistler-suppressed thermal conductivity in Boussinesq

The thermal diffusivity associated with Eq. (1) is

$\displaystyle\chi=\frac{\chi_{0}}{1+\frac{1}{3}\beta\lambda_{\text{mfp,e}}/L_{\mathrm{T},\parallel}}.$
(13)

To rewrite it in terms of our Boussinesq variables, we first split the total
temperature gradient parallel to the magnetic field $L_{\mathrm{T},\parallel}$
into the contributions coming from the background profile plus the temperature
perturbation

$\displaystyle L_{\mathrm{T},\parallel}^{-1}$
$\displaystyle=\left\lvert\frac{\bm{b}\bm{\cdot}\bm{\nabla}T}{T}\right\rvert\simeq\left\lvert
b_{z}\left(\frac{\mathrm{d}\ln
T}{\mathrm{d}z}\right)_{0}+(\bm{b}\bm{\cdot}\bm{\nabla})\frac{\delta
T}{T_{0}}\right\rvert$
$\displaystyle=\frac{\omega_{\mathrm{T}}^{2}}{g_{0}}\left\lvert
b_{z}+\omega_{\mathrm{T}}^{-2}\bm{b}\bm{\cdot}\bm{\nabla}\theta\right\rvert.$
(14)

We also assume hydrostatic equilibrium to write the pressure scale-height as
$H=c_{\mathrm{s}}^{2}/g_{0}$, where $c_{\mathrm{s}}^{2}$ is the isothermal
sound speed. Equation (3) is then obtained combining Eq. (14) with the
definitions for $\tilde{\beta}$ (Eq. 2) and $\alpha$ (Eq. 4).

###  A.3 Initial conditions

We run numerical simulations in cubic (or square) boxes with aspect ratio of
unity. We show in this work a subset consisting of one 2D run with weak
whistler-suppression ($\alpha=2\times 10^{-3}$), and four 3D runs as follows:
one reference MTI simulation ($\alpha=0$), one MTI run with weak whistler-
suppression ($\alpha=0.012$), one MTI run with strong whistler-suppression
($\alpha=0.12$), and one forcing-only run with isotropic thermal conduction
and no whistler-suppression ($\alpha=0$). The resolution is $1024^{2}$ for the
2D run, and $288^{3}$ for the 3D simulations. The parameters of the 2D run are
$\chi_{0}/(L^{2}\omega_{\mathrm{T}})=10^{-3}$, $\nu/\chi_{0}=0.01$ and the
physical Prandtl number $\nu/\eta=1$, while in 3D we choose
$\chi_{0}/(L^{2}\omega_{\mathrm{T}})=0.0025$, $\nu/\chi_{0}=0.042$ and
$\nu/\eta=4$ (note that in the forcing-only run, thermal conduction is
isotropic with the same $\chi_{0}$). The hyperdiffusivity is
$\nu_{\text{th}}^{n}/(\Delta h^{4}\chi_{0})=5.2\times 10^{-3}$ for the 2D run,
and $\nu_{\text{th}}^{n}/(\Delta h^{4}\chi_{0})=4.8\times 10^{-3}$ for the 3D
runs (with the exception of the forcing-only run, where it is zero), where
$\Delta h$ is the grid size (uniform in all directions). All our runs have
$N^{2}/\omega_{\mathrm{T}}^{2}=0.1$.

The 2D run, the reference 3D MTI run and the 3D forcing-only run are
initialized with random velocity fluctuations of amplitude $|\delta
u_{i}|/(L\omega_{\mathrm{T}})=10^{-4}$, for $i=1,2,3$. While the 2D run has a
uniform horizontal magnetic field (with initial
$B_{0}/(L\omega_{\mathrm{T}})=10^{-4}$), the reference 3D MTI run and forcing-
only have no-net magnetic flux, and the initial magnetic field is given by
$\bm{B}=B_{0}[\sin(2\pi z/L)\bm{e}_{x}+\cos(2\pi z/L)\bm{e}_{y}]$, with the
same $B_{0}$ as in 2D. These values correspond to $\tilde{\beta}_{0}=2\times
10^{5}$ in 2D and $\tilde{\beta}_{0}=5\times 10^{5}$ in 3D (volume-averaged).

The 3D MTI runs with weak and strong whistler-suppression instead use as
initial conditions a snapshot of the 3D MTI reference run. Similarly, when we
turn on forcing we use as initial condition a snapshot of the corresponding
run without forcing.