Do roots of these polynomials approach the negative of the Euler-Mascheroni constant? Let $p_n$ be the $n$th degree polynomial that sends $\frac{k(k-1)}{2}$ to $\frac{k(k+1)}{2}$ for $k=1,2,...,n+1$. E.g., $p_2(x) = (6+13x -x^2)/6$ is the unique quadratic polynomial $p(x)$ satisfying $p(0) = 1$, $p(1) = 3$, and $p(3) = 6$. Then it appears that $p_n(x)-x$ always has precisely one negative real root, and moreover this root (as a function of $n$) appears to approach $-0.577$ as $n$ gets large. Do these roots indeed approach a limit, and is this limit indeed the negative of the Euler-Mascheroni constant?
For $n=60$, the root is about $-0.580$; for $n=120$, the root is about $-0.577$. That's as far as I've gone.
(If the question seems unmotivated, I'll give some background, though it might not help anyone find an answer. One way to think about the "values" of certain divergent series is to view them as fixed points of an associated function; e.g., for the series $1+2+4+8+\dots$, the linear map $x \mapsto 2x+1$ sends the sum of the first $n$ terms of the series to the sum of the first $n+1$ terms, and the fixed point of this map is $-1$, which is indeed the natural value to assign to the divergent series. The problem I'm posting arose from trying to analyze $1+2+3+4+\dots$ in a similar way.)
 A: According to Mathematica your polynomials satisfy the recurrence relation
$$
(2 n+1) p(n) \left(n^2+3 n-2 x+2\right)+p(n+1) \left(-4 n^3-18 n^2+4 n x-27 n+2 x-14\right)+(2 n+3) (n+2)^2 p(n+2)=0
$$
with initial conditions $p_1=1+2x,p_2=1-x(x-13)/6$. The solution to this equation is
$$
p_n(x)=\frac{1}{4} \left(\frac{\Gamma \left(n+\frac{1}{2}\right) \cos \left(\frac{1}{2} \pi  \sqrt{8 x+1}\right) \Gamma \left(n+\frac{1}{2} \sqrt{8 x+1}+\frac{3}{2}\right) \Gamma \left(n-\frac{1}{2} \sqrt{8 x+1}+\frac{3}{2}\right) \, _4\tilde{F}_3\left(1,n+\frac{1}{2},n-\frac{1}{2} \sqrt{8 x+1}+\frac{3}{2},n+\frac{1}{2} \sqrt{8 x+1}+\frac{3}{2};n+\frac{3}{2},n+2,n+2;1\right)}{\pi }+2 \, _3F_2\left(-\frac{1}{2},\frac{1}{2}-\frac{1}{2} \sqrt{8 x+1},\frac{1}{2} \sqrt{8 x+1}+\frac{1}{2};\frac{1}{2},1;1\right)+4 x+2\right)
$$
For $n\to\infty$, Mathematica claims that this becomes
$$
p_\infty(x)=\frac{1}{2} \left(\, _3F_2\left(-\frac{1}{2},\frac{1}{2}-\frac{1}{2} \sqrt{8 x+1},\frac{1}{2} \sqrt{8 x+1}+\frac{1}{2};\frac{1}{2},1;1\right)+2 x+1\right)
$$
which solves $p(x)=x$ at around
$$
x=-0.573825523080029241015952733\dots
$$
A: The below answer consists of two parts. At first, we prove that each $p_n(z)-z$ has the unique negative root. Next, we describe the limit as the negative root of a certain function.
Denote $x_k=k(k-1)/2$ for $k=1,2,\ldots$.
Note that $p_n(z)-z$ interpolates the function $\varphi(z):=\sqrt{2z+1/4}+1/2$ at points $x_1,\ldots,x_{n+1}$.

*

*Proving that $p_n(z)-z$ has the unique negative root. Fix $n$ throughout this section, and denote $f(z)=p_n(z)-z$. I claim that $f(z)=1+a_1z+a_2z^2+\ldots+a_n z^n$ where $a_i(-1)^{i-1}>0$. (Then the function $f(-t)$ obviously decreases for $t\geqslant 0$ and has the unique positive root, as needed.) To prove this, assume that on the contrary $a_k(-1)^{k-1}\leqslant 0$ for certain $k\in \{1,2,\ldots,n\}$. Since the function $f-\varphi$ has $n+1$ positive roots, by Rolle's theorem the function $\eta(x):=f^{(k)}(x)-\varphi^{(k)}(x)$ has at least $n+1-k$ distinct positive roots. Note that $(-1)^{(k-1)}\varphi^{(k)}(x)>0$  for all $x>-1/8$. Then, since we assumed $(-1)^{k-1}f^{(k)}(0)\leqslant 0$, we get $(-1)^{(k-1)}\eta(x)<0$. But $(-1)^{(k-1)}\eta(-1/8+0)=+\infty$. Hence $\eta$ has a root on $(-1/8,0)$. Totally, $\eta$ has at least $n+2-k$ roots on $(-1/8,+\infty)$, and $\eta^{(n+1-k)}=-\varphi^{(n+1)}$ has a root on $(-1/8,+\infty)$ again by Rolle. A contradiction.


*How to find the limit of negative roots. Let $\gamma_n$ be any simple closed contour (oriented counterclockwise) on the complex plane containing all $x_1,\ldots,x_{n+1}$ inside but $-1/8$ outside.
Then we may consider $\varphi$ as an analytic function inside $\gamma_n$.
Denoting $H_n(z)=\prod_{k=1}^{n+1} (z-x_k)$, we get by Lagrange interpolation for any $z$ outside $\gamma_n$:
$$
p_n(z)-z=H_n(z)\sum_{k=1}^{n+1}\varphi(x_k) \frac{1}{(z-x_k)H_n'(x_k)}=
H_n(z)\sum_{k=1}^{n+1}{\rm res}_{w=x_k}\frac{\varphi(w)}{(z-w)H_n(w)}\\
=\frac1{2\pi i} H_n(z)\int_{\gamma_n} \frac{\varphi(w)}{(z-w)H_n(w)}dw.
$$
It is more convenient to renormalize now $H_n$: namely, denote $$G_n(z):=z(1-z/x_2)\ldots (1-z/x_{n+1})=(-1)^{n}(x_2\ldots x_{n+1})^{-1}H_n(z).$$ The roots of $p_n(z)-z$ are the same as those of
$$q_n(z):=\int_{\gamma_n}\frac{\varphi(w)}{(z-w)G_n(w)}dw.
$$
For $z<0$ choose $\gamma_n=[-\varepsilon-iR,-\varepsilon+iR]$ plus the corresponding right semicircle for large $R$. The integral over the semicircle tends to 0, thus $$q_n(z)=-\int_{-\varepsilon-i\infty}^{\varepsilon+i\infty}\frac{\varphi(w)}{(z-w)G_n(w)}dw.$$
Denote by $$G(z):=z(1-z/x_2)(1-z/x_3)\ldots $$
the infinite product. It is an entire function, and $G_n$ converge to $G$ uniformly on compact sets. I claim that $q_n$ converge to
$$
q(z):=-\int_{-\varepsilon-i\infty}^{-\varepsilon+i\infty}\frac{\varphi(w)}{(z-w)G(w)}dw=-\frac1{\pi i z}+v.p.
\int_{-i\infty}^{+i\infty}\frac{\varphi(w)}{(z-w)G(w)}dw
$$
uniformly on compact subsets of $(-\infty,-\varepsilon)$. Indeed, this is clear for the integral over any segment $[-\varepsilon-iR,-\varepsilon+iR]$ for a fixed $R$, and for complement of this segment on the line $-\varepsilon+i\mathbb{R}$ we may note that all multiples $1-z/x_k$ have absolute values at least 1, so we may bound from above (in absolute value) the integrand
$\frac{\varphi(w)}{(z-w)G_n(w)}$ by $\frac{|\varphi(w)|}{|(z-w)w}$ which is summable on this vertical line (even uniformly in $z$), so if $R$ is large enough, the integral over $|\Im w|>R$ is small uniformly for all $n$.
I guess that $q(z)$ has a negative root $z_0$ and changes the sign passing through $z_0$ (that would follow from $q(0)$ and $q(-\infty)$ having different sign that looks provable via convergence $q_n\to q$, although some details must be completed), then, the negative roots of $q_n$ must converge to $z_0$ (since for arbitrary $a<z_0<b<0$ we have $q_n(b)q_n(a)<0$ for large $n$, so $q_n$ has a root on $(a,b)$.
