If $a_1=1, a_n=|\cot a_{n-1}|$, then what is $\lim\limits_{n\to\infty}\text{median}\{a_1,a_2,\dots,a_n\}$?

Question

I made up the following sequence: $a_1=1, a_n=|\cot a_{n-1}|$.

What is $\lim\limits_{n\to\infty}\text{median}\{a_1,a_2,\dots,a_n\}$ ?

My thoughts

Here is the graph of $a_n$ against $n$, for $1\le n \le 25$.

Using Excel, we have

$\text{median}\{a_1,a_2,\dots,a_{100}\}\approx 1.107$
$\text{median}\{a_1,a_2,\dots,a_{1000}\}\approx 1.172$
$\text{median}\{a_1,a_2,\dots,a_{10000}\}\approx 1.204$
$\text{median}\{a_1,a_2,\dots,a_{20000}\}\approx 1.201$

It seems that the median converges to something like $1.2$.

Here are the graphs of $y=|\cot x|$ and $y=x$.

We start at $(1, |\cot 1|)$, then go horizontally to the line $y=x$, then go vertically to the curve $y=|\cot x|$, then go horizontally to the line $y=x$, etc. Among the points that we meet on the curve $y=|\cot x|$, apparently (based on Excel) about half of them are below the line $y=1.2$.

Note: If we change $a_n$ from $1$ to any other non-zero integer, the median seems to converge to the same number, approximately $1.2$.

Answer 1

Here is a conditional result, providing a candidate for the value of the limit.

---

Let $x_1 = 1$ and $x_{n+1} = \cot x_n$. Since $\cot$ is an odd function, we have $|x_n| = a_n$ for all $n$. Based on some numerical experiments, we make the following ansatz:

Assumption. The empirical probability measure $\mu_N = \frac{1}{N} \sum_{n=1}^{N} \delta_{x_n}$ converges weakly to the centered Cauchy distribution $\mu$ with shape parameter $\gamma > 0$, that is,

$$ \mu_N(\mathrm{d}x) \to \mu(\mathrm{d}x) = \frac{\gamma}{\pi(x^2 + \gamma^2)} \, \mathrm{d}x. $$

Under this assumption, for $x > 0$ we have

$$ \mathbf{P}_{X\sim \mu_N}(|X| \leq x) \to \mathbf{P}_{X\sim\mu}(|X| \leq x) = \frac{2}{\pi}\arctan\left(\frac{x}{\gamma}\right), $$

and then it can be proved that

\begin{align*} [\text{median of $\{a_1,\ldots,a_n\}$}] &= [\text{median of $|X|$, where $X \sim \mu_N$}] \\ &\to [\text{median of $|X|$, where $X \sim \mu$}] = \gamma. \end{align*}

So it suffices to determine the value of $\gamma$. To this end, consider a random variable $X \sim \mu$. Then $\cot X$ has the same distribution as $X$. By comparing the PDF of $X$ and $\cot X$, we get

\begin{align*} \frac{\gamma}{\pi(x^2 + \gamma^2)} = \frac{1}{1+x^2} \sum_{k\in\mathbb{Z}} \frac{\gamma}{\pi((\operatorname{arccot} x + k \pi)^2 + \gamma^2)} = \frac{\coth \gamma}{\pi(x^2 + \coth^2 \gamma)}. \end{align*}

Therefore $\gamma$ must satisfy $\gamma = \coth \gamma$. This equation has a unique positive solution with the value

$$ \gamma \approx 1.1996786402577338339, $$

which agrees with the numerical observation made by other users.

---

Addendum. Below is the probability histogram of first $10^4$ terms of $(x_n)$ and the graph of PDF of the limiting Cauchy distribution $\mu$:

Caution. Numerical simulation are not totally suited for probing the behavior of the median.

This is because the sequence $(x_n)$, hence $(a_n)$, is very susceptible to small numerical errors, which is significantly amplified when a term gets close to one of the poles of $\cot$. (This is essentially due to loss of significance.)

For instance, simulating $n \mapsto \text{median}(a_1,\ldots,a_n)$ using machine precision and $10000$-digit precision shows fairly different behavior even starting around $n = 20$:

Answer 2

I've expanded my examination of this problem. I have run out the sequence to $n=10^7$ with $a_1=1$, as before. I attempted to find the median for each value of $n$ but it was too time consuming. I stopped the calculation at $n\approx 3.5\cdot10^6$ and then continued the calculation for the last 100,000 points. The results are shown in the two figures below. As you can see in the second figure, the results do not settle down to a constant value at all. It is my opinion that they never will.