I know that iterating the following incircle construction approaches an equilateral triangle in the limit:
Starting with any triangle $T$, one forms $T'$ by connecting the three points
of tangency of the circle inscribed inside $T$.
Does the analogous process for tetrahedra approach a regular tetrahedron? And the same question may be asked for simplices in $\mathbb{R}^d$. A reference would be appreciated—Thanks!
Here is a partial answer: in $\mathbb R^3$ there are cycles of length $2$. To show this, it is enough to find non-trivial systems of unit vectors $x_i, y_j$ ($i,j=1,\dots,d+1$) such that $(x_i,y_j)=t$ for some $t>0$ and all $i\ne j$. If $d=3$, we can take the normalized copies of $(1,0,4),(1,0,-4),(-1,4,0),(-1,-4,0)$ and $(-2,0,-1),(-2,0,1),(2,-1,0),(2,1,0)$. However, in dimension $d=4$, I do not see how to construct anything like that, so we have a chance for the affirmative answer again.
Edit: OK, let's give it a shot. Unfortunately, the answer is negative in every dimension $d>2$. If we believe dimensional considerations, that is easy to predict: the 2 cycle in $\mathbb R^d$ will require constructing $d+1$ unit $x$ vectors and $d+1$ unit $y$ vectors in a non-trivial decent position so that all scalar products $(x_i,y_j)$, $i\ne j$ are equal. $(d+1)$ vectors on the sphere require $2(d^2-1)$ parameters. We have $d(d+1)-1$ equations plus the action of the orthogonal group of dimension $\frac{d(d-1)}2$, so as soon as $$ 2(d^2-1)>d(d+1)-1+\frac{d(d-1)}2\, $$ which happens for $d>2$, we have a chance. Of course, we need to check that there is no stupid degeneracy anywhere, so the exact analysis is a bit more tiresome.
Let's try to perturb the standard configuration in which $x_i=y_i$ are the vertices of a regular simplex. Let $z_i$ and $w_i$ be the infinitesimal perturbations of $x_i$ and $y_i$ correspondingly. To keep the norms, we need $(z_i,x_i)=(x_i,w_i)=0$. Once we have that, the only linear restrictions are $\sum_j (z_i,x_j)=0$ and $sum_i (x_i,w_j)=0$. The perturbations of the scalar products are $(z_i,x_j)+(x_i,w_j)$, i.e., we are allowed sums of two matrices with zero diagonals, the first of which has row sums $0$ and the second one has column sum $0$. For $d\ge 3$, this allows us to create any linear perturbation of the scalar product matrix with zero diagonal and total sum $0$, and, moreover, to do it with the dimension excess $2(d+1)(d-1)-d(d+1)+1$. Note that the addition of a matrix with total sum $0$ is not enough to achieve arbitrary scalar products, but is enough to balance the scalar products if they are slightly disbalanced.
Now we just play the usual implicit function game: create linear perturbations of order $t$ so that they don't change the scalar products in the first order but the resulting configuration is $ct$ away from any rotation of the initial configuration (this is pure dimension story). Now, the scalar products differ by about $t^2$. By the implicit function theorem (which should rather be called "sufficient linear range" theorem in this case), we can now use another perturbation of size $t^2$ to balance the scalar products perfectly. This will still leave us away from rotations of the regular simplex but create a cycle of length $2$.
The interesting feature here is that the incenter coincides with the circumcenter. I do not have any counterexample to the general conjecture that we always have convergence to this position (which is always a $2$-cycle). Probably, this is how the question should be posed now:
Is it true that no matter what the starting simplex was, the distance between the incenter and the circumcenter of the rescaled to size $1$ iterated simplices tends to $0$?
It looks like you cannot get away with the length $2$ cycle anymore, and to analyze longer options seems quite a bit less trivial task.
Another edit (the penultimate one, I hope :-))
As I said above, if we have a simplex whose circumcenter coincides with its incenter, it generates a cycle of length $2$. Joseph's experiment suggests that we always get attracted to a cycle like this. I think I have an idea how to prove it though one "little detail" is still missing. Here is the sketch.
Let us normalize all iterations so that the circumradius is $1$. Then the problem can be restated as follows:
Given a system of unit vectors $x_i$ such that $0$ can be written as their convex combination with positive coefficients, define the new system $y_j$ of unit vectors by $(x_i+u,y_j)=t$ ($i\ne j$) where $u$ is some vector (whose physical meaning is the vector going from the incenter to the circumcenter) and $t>0$ is some number (the radius of the inscribed sphere). The next step will be given by the equations $(z_i,y_j+v)=s$. What I want to show is that $|v|<|u|$ and, moreover, if we have some quantitative bound on the non-degeneracy of $x_i$, we can improve it to $|v|\le q|u|$ with $q<1$. If we can then show that our configurations are uniformly non-degenerate (the part that is still missing), we'll have geometric convergence of the "correction vectors" $u$ to $0$, which is enough to ensure the geometric convergence of the vector systems themselves to a $2$-cycle.
One key point is that both $x_i+u$ and $z_i$ are orthogonal to $y_j-y_k$ for $j,k\ne i$, which ensures that they are collinear. With a little more thought, we see that $x_i+u=|x_i+u|z_i$ (i.e., the directions are the same).
The other key point is that if we have two systems of vectors $X_i,Y_j$ such that $(X_j,Y_j)=t$ for all $i\ne j$ and $0$ can be written as convex combinations $$ 0=\sum_i a_i X_i=\sum b_j Y_j\,, $$ then $$ a_k(X_k,Y_k)+t\sum_{i:i\ne k}a_i=\left(\sum_i a_iX_i, Y_k\right)=0= \left(X_k,\sum_j b_jY_j\right)=b_k(X_k,Y_k)+t\sum_{j:j\ne k}b_j\,, $$ which is normally enough to conclude that $a_k=b_k$, the "normal" case including our geometric situation.
Now we get $\sum_i a_i(x_i+u)=\sum_i a_i y_i=0$ and $\sum_i b_iz_i=\sum_ib_i(y_i+v)=0$. Since $x_i+u=|x_i+u|z_i$, we get $b_i=\lambda|x_i+u|a_i$ with some $\lambda>0$ chosen so that $\sum_i a_i=\sum_i b_i=1$. Thus, for every $\mu$, we have $$ |v|=\left|\sum_i (b_i-\mu a_i)y_i\right|\le \sum_i|b_i-\mu a_i|\,. $$ Now notice that $1-|u|\le |x_i+u|\le 1+|u|$ and $|u|<1$ (the incenter is inside the circumsphere). Hence the inequality $|v|\le|u|$ will be obtained if we prove the following
Lemma:
Let $a_j>0$, $\sum a_j=1$. Let $\delta\in(0,1)$. Assume that $|u_j|\le\delta$ and $b_j=\frac{a_j(1+u_j)}{1+\sum_i a_i u_i}$. Then $$ \min_{\mu}\sum_j|b_j-\mu a_j|\le \delta $$
Proof
Multiplying both sides by $1+\sum_i a_i u_i$ and adjusting $\mu$, we can rewrite the conclusion of the lemma as $$ \min_{\mu}\sum_j a_j|u_j-\mu|\le \delta\left(1+\sum_i a_i u_i\right) $$ Now the left hand side is convex in $u_j$ and the right hand side is linear in $u_j$, so it is enough to consider the case when each $u_j$ is either $+\delta$ or $-\delta$. Let $A_\pm=\sum_{j:u_j=\pm \delta}a_j$. Then the left hand side is $2\delta\min(A_+,A_-)$ while the right hand side is $\delta[1+\delta(A_+-A_-)]$. Canceling $\delta$, we get the inequality $2\min(A_+,A_-)\le 1+\delta(A_+-A_-)=(1-\delta)A_-+(1+\delta)A_+$, which is obvious.
We need to be a bit more careful to get the strict inequality, of course (the point is that if the system $x_i$ is non-degenerate, the inequalities $1-|u|\le |x_i+u|\le 1+|u|$ are strict most of the time and when we have a quantitative bound for non-degeneracy, the vectors $x_i$ and $u$ cannot align well). However, for now what we've done is good enough.
Last edit (I hope):
OK, let's prove the uniform non-degeneracy claim. We will show that the inradius does not decrease at each step. Since the circumradius is fixed, this will finish the story. In the above notation, we need to show that $s\ge t$. Let's again talk about general systems SX_i,Y_jS satisfying $(X_i,Y_j)=t$ for all $i\ne j$. Let $\sum_i a_iX_i=\sum_i a_i Y_i=0$. For each $j$, we have $0=\sum_i a_i(X_i,Y_j)=a_j(X_j,Y_j)+(1-a_j)t$, whence, summing over $j$, $$ dt=-\sum_j a_j(X_j,Y_j)\,. $$ Note that we can add any fixed vector to each $X_j$ or to each $Y_j$ (but not both) in the scalar products on the right hand side without affecting the identity.
Now, we get $$ ds=-\sum_j b_j(z_j,y_j+v)=-\sum b_j(z_j,y_j)=-\lambda\sum_j a_j(x_j+u,y_j)=\lambda dt\,, $$ so all we need is to show that the normalizing factor $\lambda$ is at least $1$, i.e., that $\sum_j a_j|x_j+u|\le 1$. However, $|x_j+u|=\sqrt{1+2(x_j,u)+|u|^2}\le 1+(x_j,u)+\frac 12|u|^2\le 1+(x_j+u,u)$. Thus $$ \sum_j a_j|x_j+u|\le 1+\sum_ja_j(x_j+u,u)=1 $$ as needed.
The proof may need some minor combing to put everything together in a neat and completely rigorous way, but all that final polishing is totally routine, so, unless somebody requests it explicitly, I'll leave it to the reader.
The upshot is the following
Proposition: Normalized in any meaninful way, the sequence of iterations always converges to a $2$-cycle generated by a simplex whose incenter and circumcenter coincide. Moreover, the convergence speed is geometric, and, if the normalization fixes the circumradius, at each step the distance between the incenter and the circumcenter decreases and the inradius increases.
Note that this also encompasses the classical result for $d=2$ because the only triangle on the plane whose circumcenter coincides with the incenter is the equilateral triangle. I guess we have figured out everything worth figuring out here by now unless Joseph gives this problem some new and unexpected twist again :-).
There is likely a theorem here that, for almost every start tetrahedron (left column below),
very quickly the iterates form an approximate $2$-cycle. Below shows three random
tetrahedra and five iterates each (rescaled):
Fedja's example of a $2$-cycle:
