Does recursively replacing $\frac1n$ by $\frac1n(\frac12+\dots+\frac1{n+1})$ really converge to $\frac1e$?

I was thinking of a problem and have an answer through computer programming, but am unable to prove it mathematically. Start with the following: $$\frac{1}{2}\bigg(\frac{1}{2} + \frac{1}{3}\bigg)\approx 0.41666$$ Replace every unit fraction being summed in the parentheses (in this case, the inner 1/2 and 1/3) with $$\frac{1}{n}\bigg(\frac{1}{2} + \frac{1}{3} + ... + \frac{1}{n+1}\bigg)$$ where $n$ is the denominator of the unit fraction. If we apply this process once, we get $$\frac{1}{2}\bigg(\frac{1}{2}\bigg(\frac{1}{2} + \frac{1}{3}\bigg)+\frac{1}{3}\bigg(\frac{1}{2} + \frac{1}{3}+\frac{1}{4}\bigg)\bigg)\approx 0.3888$$ and if we apply it again, we get $$\frac{1}{2}\bigg(\frac{1}{2}\bigg(\frac{1}{2}\bigg(\frac{1}{2} + \frac{1}{3}\bigg) + \frac{1}{3}\bigg(\frac{1}{2} + \frac{1}{3}+ \frac{1}{4}\bigg)\bigg)+\frac{1}{3}\bigg(\frac{1}{2}\bigg(\frac{1}{2} + \frac{1}{3}\bigg) + \frac{1}{3}\bigg(\frac{1}{2} + \frac{1}{3}+ \frac{1}{4}\bigg)+\frac{1}{4}\bigg(\frac{1}{2} + \frac{1}{3}+ \frac{1}{4}+ \frac{1}{5}\bigg)\bigg)\bigg)\approx 0.37754$$

As we repeat this process, the result seems to approach $e^{-1}$.

If anyone could prove this or direct me to a proof of this, that would be wonderful!

Consider the following random process:

• Start with a natural number $n_0$.
• Choose a new number $n_1$ from $\{2,3,\ldots,n_0+1\}$ uniformly at random.
• Choose a new number $n_2$ from $\{2,3,\ldots,n_1+1\}$ uniformly at random.

and repeat for however many steps. Your expression is ${\mathbb E}[1/n_k]$ - that is, it is the expected value of $1/{n_k}$. This is a Markov process and we can apply our typical analysis tools to it*.

In particular, as $k$ increases, irrespective of our starting position, the distribution of $n_k$ approaches a unique steady state distribution - that is, a distribution invariant under the step.

Let $p_m$ be the probability that, for large $k$ (i.e. in the limit), we have $n_k=m$. It must be a steady state in the sense that $$p_m=\sum_{\ell=m-1}^{\infty}\frac{1}{\ell}\cdot p_{\ell}$$ which just says that, if we started in this distribution, we would stay in it.

Now, we can construct a probability generating function for the steady state; we'll mess up the indexing slightly for convenience. Let: $$f(x)=p_2x+p_3x^2+p_4x^3+\ldots.$$ Note that we can use formal integration to get $$\int_0^x f(t)\,dt=\frac{p_2}2x^2+\frac{p_3}3x^3+\ldots$$ Then, noting $\frac{1}{1-x}=1+x+x^2+\ldots$ and slipping in an extra $x$ for fun, we can define a new function $g(x)$ as follows: $$g(x)=\frac{x\int_0^x f(t)\,dt}{1-x}=\frac{p_2}2x^3+\left(\frac{p_2}2+\frac{p_3}3\right)x^4+\left(\frac{p_2}2+\frac{p_3}3+\frac{p_4}4\right)x^5+\ldots.$$ Why do this? Because if we let $C=\frac{p_2}2+\frac{p_3}3+\ldots$ - that is, $C$ is exactly the quantity we're looking for - we note that $$\frac{Cx}{1-x}-g(x)=\sum_{m=1}^{\infty}\sum_{\ell=m}^{\infty}\frac{1}{\ell}\cdot p_{\ell}\cdot x^m$$ Which we recognize as $f(x)$. Thus, we get an equation: $$(1-x)f(x)=Cx-x\int_0^x f(t)\,dt$$ Now, I'll leave out the details, because I'm not good at calculus, but if you differentiate both sides of the equation twice, you get a degree two differential equation, and if you plug it into Mathematica, you get a solution: $$f(x)=c_1\cdot xe^x+c_2\cdot \frac{e+xe^x\operatorname{Ei}(x)}e$$ where $c_1,c_2$ are as of yet unknown constants. However, we know that $f(0)=0$ from definition of $f$ and $f(1)=1$ since $p_2+p_3+\ldots=1$ since these are probabilities. This fixes $c_2=0$ and $c_1=\frac{1}e$.

Thus $f(x)=\frac{1}exe^x$. This looks promising already! Then, observe that $\int_0^1 f(t)\,dt=C={\mathbb E}[1/n]$ just by looking at the series we derived earlier for the integral! So, then we just integrate that expression and we see $C=1/e$. Hurray! It worked! As a bonus, we know the steady state probabilities are just $p_n=\frac{1}{e(n-2)!}$ for $2\leq n$.

(*I'm ignoring convergence issues, which are relevant since this Markov chain has infinitely many states. The generating function stuff can be justified after the fact, despite involving serious convergence problems - uniqueness of a solution is not so hard to establish - the usual argument of uniqueness and existence of a steady state still gives us uniqueness for infinite Markov processes - and the fact that $f(x)=\frac{1}exe^x$ is a solution is also not too hard to show. Convergence of the chain to its steady state brings more serious issues - I think they can be solved by bounding the probability that $n_k$ is large, which shouldn't be too hard, since the probability that $n_{k+1}>n_k$ given $n_k$ is just $\frac{1}{n_k}$)

We provide a recurrence relation for the numbers converging to $\frac{1}{e}$. We start with some introductory aspects.

Representation via sums

We denote the elements of the sequence which are iteratively generated by OPs rule with $\alpha_n, n\geq 2$ and derive for small $n$

\begin{align*} \alpha_2&=\frac{1}{2}\\ \alpha_3&=\frac{1}{2}\sum_{j_1=2}^3\frac{1}{j_1}\tag{i}\\ &=\frac{1}{2}\left(\frac{1}{2}+\frac{1}{3}\right)=\frac{5}{12}=0,41\dot{6}\\ \alpha_4&=\frac{1}{2}\sum_{j_1=2}^3\frac{1}{j_1}\sum_{j_2}^{j_1+1}\frac{1}{j_2}\tag{ii}\\ &=\frac{1}{2}\left(\frac{1}{2}\sum_{j_2=2}^3\frac{1}{j_2}+\frac{1}{3}\sum_{j_2=2}^4\frac{1}{j_2}\right)\\ &=\frac{1}{2}\left(\frac{1}{2}\left(\frac{1}{2}+\frac{1}{3}\right)+\frac{1}{3}\left(\frac{1}{2}+\frac{1}{3}+\frac{1}{4}\right)\right) =\frac{7}{18}=0,3\dot{8}\\ \alpha_5&=\frac{1}{2}\sum_{j_1=2}^3\frac{1}{j_1}\sum_{j_2}^{j_1+1}\frac{1}{j_2}\sum_{j_3=2}^{j_2+1}\frac{1}{j_3}\tag{iii}\\ &=\frac{1}{2}\left(\frac{1}{2}\sum_{j_2=2}^3\frac{1}{j_2}\sum_{j_3=2}^{j_2+1}\frac{1}{j_3}+\frac{1}{3}\sum_{j_2=2}^4\frac{1}{j_2}\sum_{j_3=2}^{j_2+1}\frac{1}{j_3}\right)\\ &=\frac{1}{2}\left(\frac{1}{2}\left(\frac{1}{2}\sum_{j_3=2}^3\frac{1}{j_3}+\frac{1}{3}\sum_{j_3=2}^4\frac{1}{j_3}\right)\right.\\ &\qquad \left.+\frac{1}{2}\left(\frac{1}{2}\sum_{j_3=2}^3\frac{1}{j_3}+\frac{1}{3}\sum_{j_3=2}^4\frac{1}{j_3}+\frac{1}{4}\sum_{j_3=2}^5\frac{1}{j_3}\right)\right)\\ &=\frac{1}{2}\left(\frac{1}{2}\left(\frac{1}{2}\left(\frac{1}{2}+\frac{1}{3}\right)+\frac{1}{3}\left(\frac{1}{2}+\frac{1}{3}+\frac{1}{4}\right)\right)\right.\\ &\qquad\left.+\frac{1}{3}\left(\frac{1}{2}\left(\frac{1}{2}+\frac{1}{3}\right)+\frac{1}{3}\left(\frac{1}{2}+\frac{1}{3}+\frac{1}{4}\right) +\frac{1}{4}\left(\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+\frac{1}{5}\right)\right)\right)\\ &=\frac{1\,631}{4\,320}=0,377\,54\overline{629} \end{align*} in accordance with OP's calculation.

From (i),(ii),(iii) we obtain for general $n$: \begin{align*} \color{blue}{\alpha_n=\frac{1}{2}\sum_{j_1=2}^3\frac{1}{j_1}\sum_{j_2=2}^{j_1+1}\frac{1}{j_2}\sum_{j_3=2}^{j_2+1}\frac{1}{j_3} \cdots\sum_{j_{n-2}=2}^{j_{n-3}+1}\frac{1}{j_{n-2}}\sum_{j_{n-1}=2}^{j_{n-2}+1}\frac{1}{j_{n-1}}\qquad\qquad n\geq 3}\tag{1} \end{align*}

Note: The ubiquitous Catalan numbers occur here also. If we write $\alpha_n$ as \begin{align*} \alpha_2&=\frac{1}{2},&C_1=1\\ \alpha_3&=\frac{1}{\color{blue}{2}\cdot 2}+\frac{1}{\color{blue}{2}\cdot 3},&C_2=2\tag{2}\\ \alpha_4&=\frac{1}{\color{blue}{2\cdot 2}\cdot 2}+\frac{1}{\color{blue}{2\cdot 2}\cdot 3}+\frac{1}{\color{blue}{2\cdot 3}\cdot 2}+\frac{1}{\color{blue}{2\cdot 3}\cdot 3}+\frac{1}{\color{blue}{2\cdot 3}\cdot 4},&C_3=5\\ &\ \,\vdots \end{align*} we see the number of summands of $\alpha_n$ is $C_{n-1}$. The blue marked factors in $\alpha_n$ come from the factors of $\alpha_{n-1}$.

Recurrence relation

In order to find a recurrence relation we would like to reverse the order of summation in (1). Alas this is somewhat cumbersome since the upper limit of an inner sum does not nicely correspond with the lower limit of the next outer sum. We will instead show the relationship with the help of number triangles which might also provide additional insight $$ \begin{array}{r|cccc} n&w_{\alpha_2}&w_{\alpha_3}&w_{\alpha_4}&w_{\alpha_5}\\ \hline\\ 5&&&&\frac{1}{5}\\\\ 4&&&\frac{1}{4}&\frac{1}{4}\\\\ 3&&\frac{1}{3}&\frac{1}{3}&\frac{1}{3}\\\\ 2&\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&\frac{1}{2}\\ \\ \hline\\ \\ \end{array} \qquad\qquad\qquad \begin{array}{r|rrrr} n&\alpha_2&\alpha_3&\alpha_4&\alpha_5\\ \hline\\ 5&&&&\frac{1}{120}\\\\ 4&&&\frac{1}{4}&\frac{12}{288}\\\\ 3&&\frac{1}{6}&\frac{5}{36}&\frac{7}{54}\\\\ 2&\frac{1}{2}&\frac{1}{4}&\frac{5}{24}&\frac{7}{36}\\ \\ \hline\\ \alpha_n&\frac{1}{2}&\frac{5}{12}&\frac{7}{18}&\frac{1\,631}{4\,320} \end{array} $$

The left table shows a triangular region of weights $\frac{1}{2},\frac{1}{3},\ldots$ with weight $\frac{1}{n}$ in row $n$. These weights are used to build the entries in the right table and $\alpha_n$ which is the sum of the entries in the $n$-th column. We denote rows and columns starting from $2$ so that the bottom left-most element is $a_{22}=\frac{1}{2}$ while the top right-most element is $a_{55}=\frac{1}{5!}=\frac{1}{120}$. We have \begin{align*} \color{blue}{\alpha_2}&=a_{22}\color{blue}{=\frac{1}{2}}\\ \color{blue}{\alpha_3}&=a_{32}+a_{33}\\ &=\frac{1}{2}a_{22}+\frac{1}{3}a_{22}\\ &=\frac{1}{4}+\frac{1}{6}\color{blue}{=\frac{5}{12}}\\ \color{blue}{\alpha_4}&=a_{42}+a_{43}+a_{44}\\ &=\frac{1}{2}\left(a_{32}+a_{33}\right)+\frac{1}{3}\left(a_{32}+a_{33}\right)+\frac{1}{4}a_{33}\\ &=\frac{5}{24}+\frac{5}{36}+\frac{1}{24}\color{blue}{=\frac{7}{18}}\\ \color{blue}{\alpha_5}&=a_{52}+a_{53}+a_{54}+a_{55}\\ &=\frac{1}{2}\left(a_{42}+a_{43}+a_{44}\right)+\frac{1}{3}\left(a_{42}+a_{43}+a_{44}\right) +\frac{1}{4}\left(a_{43}+a_{44}\right)+\frac{1}{5}a_{55}\\ &=\frac{7}{36}+\frac{7}{54}+\frac{13}{288}+\frac{1}{120}\color{blue}{=\frac{1\,631}{4\,320}} \end{align*}

From the calculations above we derive the rule to build the entries $a_{k,l}$:

• Take $\frac{1}{k}$ times the sum of all entries of column $k-1$ and row greater or equal $\max\{l-1,2\}$.

• The recurrence relation is \begin{align*} a_{n,j}&=\frac{1}{j}\sum_{k=\max\{j-1,2\}}^{n-1}a_{n-1,k}\qquad\qquad\qquad\qquad\qquad n\geq 3,j=2,\ldots,n\\ \color{blue}{\alpha_2}&\color{blue}{=\frac{1}{2}}\\ \color{blue}{\alpha_n}&=\sum_{j=2}^na_{n,j}\color{blue}{=\sum_{j=2}^n\frac{1}{j}\sum_{k=\max\{j-1,2\}}^{n-1}a_{n-1,k}\qquad\qquad n\geq 3,j=2,\ldots,n}\tag{3} \end{align*}

We can show the validity of (3) by successively applying the recurrence relation to obtain \begin{align*} \color{blue}{\alpha_n}&=\sum_{j_1=2}^na_{n,j_1}\\ &=\sum_{j_1=2}^n\frac{1}{j_1}\sum_{j_2=\max\{j_1-1,2\}}^{n-1}a_{n-1,j_2}\\ &=\sum_{j_1=2}^n\frac{1}{j_1}\sum_{j_2=\max\{j_1-1,2\}}^{n-1}\frac{1}{j_2}\sum_{j_3=\max\{j_2-1,2\}}^{n-2}a_{n-2,j_3}\\ &\ \,\vdots\\ &=\sum_{j_1=2}^n\frac{1}{j_1}\sum_{j_2=\max\{j_1-1,2\}}^{n-1}\frac{1}{j_2} \cdots\sum_{j_{n-2}=max\{j_{n-3}-1,2\}}^{3}\frac{1}{j_{n-2}}\sum_{j_{n-1}=\max\{j_{n-2}-1,2\}}^{2}a_{2,j_{n-1}}\\ &\,\,\color{blue}{=\frac{1}{2}\sum_{j_1=2}^n\frac{1}{j_1}\sum_{j_2=\max\{j_1-1,2\}}^{n-1}\frac{1}{j_2} \cdots\sum_{j_{n-2}=max\{j_{n-3}-1,2\}}^{3}\frac{1}{j_{n-2}}}\tag{3} \end{align*} and (3) is the same as (1) with the summation in reverse order.

Notes:

• It is interesting, that summing up the columns in the right table above gives a monotonically increasing sequence converging to $\frac{1}{e}$. \begin{align*} \lim_{n\to\infty} \alpha_n=\frac{1}{e} \end{align*} What I find amazing is the usage of the weights $\frac{1}{j}$ used in the expression of $\alpha_n$ (see e.g. (2)) which makes them similar looking to \begin{align*} \sum_{j=0}^n(-1)^j\frac{1}{j!} \end{align*} but without the alternating sign.

• It would be great if someone could show the convergence of $(\alpha_n)$ by solving the recurrence relation (3).

• It might be interesting that the same structure of the recurrence relation (3) but with different weights, namely $n-1$ instead of $\frac{1}{n}$ and different boundary values $a_{2,n}=1$ instead of $\frac{1}{2}$ and $a_{n,n}=n-1$ instead of $\frac{1}{n}$ is archived as A185423 in OEIS denoting the number of Ordered forests of $k$ increasing plane unary-binary trees with $n$ nodes. $$ \begin{array}{r|cccc} n&&&&\\ \hline\\ 5&&&&3\\\\ 4&&&3&3\\\\ 3&&2&2&2\\\\ 2&1&1&1&1\\ \end{array} \qquad\qquad\qquad \begin{array}{r|rrrr} A185423&&&&\\ \hline\\ 5&&&&24\\\\ 4&&&6&36\\\\ 3&&2&6&30\\\\ 2&1&1&3&9\\ \end{array} $$