Trace of a nonlinear matrix equation (cont'd) Let $X_0$ be a trace-one positive definite matrix, i.e. $X_0>0$, $\mathrm{tr}(X_0)=1$. Let $A>0$ and consider the following iteration
$$
X_{k+1} = X_k^{1/2}AX_k^{1/2},\quad k\geq 0,\quad (\star)
$$
where $X_k^{1/2}$ denotes the (principal) square root of $X_k$.

My question: Is it true that if there exists $X_0$ as above such that $(\star)$ is trace-preserving starting from $X_0$, i.e. $\mathrm{tr}(X_{k+1})=\mathrm{tr}(X_{k})=1$ for all $k\geq 0$, then $A=I$?

If $X_0$ and $A$ are scalars this is clearly true. Moreover, it is easy to see that $A$ cannot be such that $A> I$ or $A< I$, but I cannot quite prove that $A$ must be the identity.
Thanks for your help.

Addendum 1. Note that by applying to $(\star)$ an orthogonal change of basis $T$ which diagonalizes $A$, we can rewrite $(\star)$ as
$$
\tilde{X}_{k+1}= \tilde{X}_k^{1/2} D \tilde{X}_k^{1/2} \quad (\star\star)
$$
where $\tilde{X}_k:=T^\top X_k T$ and $D:=T^\top AT>0$ is diagonal.

Addendum 2. A "simpler" version of this question was answered in the affirmative here.

Addendum 3. My attempts so far were based on working with the "simplified" dynamics $(\star\star)$. Namely, consider the partition $D=\left[\begin{smallmatrix}D_1 & 0 \\ 0 & D_2\end{smallmatrix}\right]$ with $D_1>I_{n_1}$ and $D_2<I_{n_2}$ (in the other cases, e.g. $D_1> I_{n_1}$, $D_2=I_{n_2}$, it is easy to show that $\mathrm{tr}(\tilde{X}_1)\neq 1$). Now by partitioning $\tilde X_0^{1/2}$ and $\tilde X_1$ accordingly to the block decomposition of $D$, we get
$$
\tilde{X}_1=\begin{bmatrix}(\tilde X_1)_{11} & (\tilde X_1)_{12} \\ (\tilde X_1)_{12}^\top & (\tilde X_1)_{22}\end{bmatrix} = \begin{bmatrix}(\tilde X_0^{1/2})_{11}D_1(\tilde X_0^{1/2})_{11} + (\tilde X_0^{1/2})_{12}D_2(\tilde X_0^{1/2})_{12}^\top & (\ast) \\ (\ast)^\top & (\tilde X_0^{1/2})_{12}^\top D_1(\tilde X_0^{1/2})_{12} + (\tilde X_0^{1/2})_{22}D_2(\tilde X_0^{1/2})_{22} \end{bmatrix},
$$
with the constraint $\mathrm{tr}(\tilde X_0)=1$ which now reads as 
$$
\mathrm{tr}\left((\tilde X_0^{1/2})_{11}(\tilde X_0^{1/2})_{11} + (\tilde X_0^{1/2})_{12}(\tilde X_0^{1/2})_{12}^\top\right)+\mathrm{tr}\left((\tilde X_0^{1/2})_{12}^\top (\tilde X_0^{1/2})_{12} + (\tilde X_0^{1/2})_{22}(\tilde X_0^{1/2})_{22}\right)=1.
$$
Now my conjecture is that $\mathrm{tr}(\tilde X_2)\neq 1$ for all $\tilde{X_0}>0$, $\mathrm{tr}(\tilde{X_0})=1$, (indeed, it is easy to find examples for which $\mathrm{tr}(\tilde X_1)= 1$). My idea is to use the same block decomposition for $\tilde X_2$ and then exploit some "trace inequalities" applied to the diagonal blocks. However I didn't manage to conclude anything so far.

Addendum 4. A more general version of this question which I suspect hold true (actually, for the $2\times 2$ case, it does hold true) is the following one:

Is it true that if there exists $X_0\ge 0$, $\mathrm{tr}(X_0)=1$, such that $(\star)$ is trace-preserving starting from $X_0$, i.e. $\mathrm{tr}(X_{k+1})=\mathrm{tr}(X_{k})=1$ for all $k\geq 0$, then $X_{k+1}=X_k^{1/2}AX_k^{1/2}=X_k$ for all $k\geq 0$?

 A: Actually, you have completely solved it yourself, just didn't dare to acknowledge it. In my notation, you have $(X\circ X^T)v(A)=(Y\circ Y^T)v(I)$ when $Y^2=XAX$. Similarly, $(Z\circ Z^T)v(I)=(Y\circ Y^T)v(A)$ when $Z^2=YAY$. Taking the trace, we must have 
$$
1=\langle(X\circ X^T)v(I),v(I)\rangle=\langle(Y\circ Y^T)v(I),v(I)\rangle=\langle(Z\circ Z^T)v(I),v(I)\rangle\,.
$$
However,
$$
\langle(Y\circ Y^T)v(I),v(I)\rangle=\langle(X\circ X^T)v(A),v(I)\rangle
$$
and
$$
\langle(Z\circ Z^T)v(I),v(I)\rangle=\langle(Y\circ Y^T)v(A),v(I)\rangle
\\
=
\langle(Y\circ Y^T)v(I),v(A)\rangle=\langle(X\circ X^T)v(A),v(A)\rangle
$$
so 
$$
\langle(X\circ X^T)(v(A)-v(I)),(v(A)-v(I))\rangle=0
$$
whence $X\circ X^T$ and, thereby, $X$ must be degenerate unless $v(A)=v(I)$, i.e., $A=I$.
This story definitely has a few morals but I'll abstain from spelling them out :-).
A: First, let me modify the given problem for the sake of simplicity.  

Given an $n \times n$ matrix $X_0>0$ with unit trace, i.e.,
$$
X_0 = \sum_{1 \le i \le n} \lambda_i v_i v_i^T \;, \quad \sum_{1 \le i \le n} \lambda_i = 1 
$$ where $\{ \lambda_i \}$ and $\{ v_i \}$ are the positive eigenvalues and orthonormal eigenvectors of $X_0$, respectively.  Set $A>0$ to be 
$$
A = \sum_{1 \le i \le n} \alpha_i v_i v_i^T   \tag{$\diamond$}
$$
where $\{ \alpha_i \}$ are positive eigenvalues of $A$, and we stress that $A$ has the same eigenvectors as $X_0$. If the solution to the matrix recurrence relation 
$$
X_{k+1} = X_k^{1/2} A X_k^{1/2}
$$ with initial data $X_0$ satisfies $\operatorname{trace}(X_{k})=1$ for all natural numbers $k \ge 0$, then $A=I_n$, i.e., $\alpha_i = 1$ for all $1 \le i \le n$.

Remark. The difference between this formulation and the OP's formulation is that $A$ has the same eigenvectors as the given seed $X_0$. Admittedly, this form of $A$ is restrictive, but it is nicer to work with because, as we will see, it simplifies the subsequent calculations.  
For the sake of contradiction, suppose that $A \ne I_n$, i.e.,  not all $\alpha_i = 1$ for $1 \le i \le n$. Since
$$
X_0^{1/2} = \sum_{1 \le i \le n} \sqrt{\lambda_i} v_i v_i^T 
$$
after one step of the recurrence relation we have that: 
\begin{align*}
X_1 &= \sum_{1 \le i,j \le n} \sqrt{\lambda_i} \sqrt{\lambda_j} (v_i^T A v_j) v_i v_j^T \\
&= \sum_{1 \le i,j,k \le n}\sqrt{\lambda_i} \sqrt{\lambda_j} \alpha_k (v_i^T v_k) (v_k^T v_j) v_i v_j^T \\
&= \sum_{1 \le i,j \le n}\sqrt{\lambda_i} \sqrt{\lambda_j} \alpha_i \delta_{ij} v_i v_j^T \qquad \text{($\delta_{ij}$ is the Kronecker delta)}\\
&= \sum_{i=1}^n \lambda_i \alpha_i v_i v_i^T
\end{align*} and the unit trace requirement implies that 
$$
\sum_{i=1}^n \lambda_i \alpha_i = 1 \;.
$$ Iterating the above calculation $k$ times, we see that the unit trace requirement  implies that $\sum_{i=1}^n \lambda_i \alpha_i^k = 1 $ which, according to our hypotheses, must hold true for any natural number $k \ge 0$. A more transparent way to write this requirement is as:
$$
\mathbf{V}^T \boldsymbol{\lambda} = \mathbf{1} \tag{$\star$}
$$
where we have introduced an infinite Vandermonde-like matrix and two vectors:
$$
\mathbf{V} = \begin{bmatrix} 1 & \alpha_1 & \cdots & \alpha_1^k & \cdots \\
1 & \alpha_2 & \cdots & \alpha_2^k & \cdots \\
\vdots & \vdots & \ddots & \vdots & \cdots \\
1 & \alpha_n & \cdots & \alpha_n^k & \cdots
\end{bmatrix} \;, \quad \boldsymbol{\lambda} = \begin{bmatrix} \lambda_1 \\
\vdots \\
\lambda_n \end{bmatrix} \;, \quad \mathbf{1} = \begin{bmatrix} 1 \\
1 \\
\vdots \end{bmatrix}
$$ If $\alpha_1 = \cdots =  \alpha_n = 1$, then the solution set of ($\star$) basically contains all unit trace matrices that are positive definite. (This is the trivial case.)   However, if even one of the eigenvalues of $A$ is not equal to one (the case at hand), then the solution set to this infinite (overdetermined) system of equations is empty, and there exists no $X_0$ that satisfies the hypotheses given above, which is a contradiction that is resolved only if $\alpha_1 = \cdots =  \alpha_n = 1$ or $A=I_n$.      
