The answer to all these questions is Yes. We shall answer this question in the more general case when $X_1,\dots,X_r$ are complex matrices (the inequality for this question looks a little bit different in the complex case since we need to apply conjugations and transposes). This answer shall be in the context of quantum channels. I recommend The Theory of Quantum Information by John Watrous (2018) for more information on quantum information theory. Much of content of this answer is similar to the answers to this related question; fedja remarked that the idea behind the answer to that previous question also applies to this question, so fedja could have answered at least most of this question before I did.
A function $\mathcal{E}:M_d(\mathbb{C})\rightarrow M_d(\mathbb{C})$ is said to be positive if $\mathcal{E}(P)$ is positive semidefinite whenever $P$ is positive semidefinite. If $W$ is a finite dimensional vector space, then let $L(W)$ be the set of all homomorphisms from $W$ to $W$. We say that $\mathcal{E}:M_d(\mathbb{C})\rightarrow M_d(\mathbb{C})$ is completely positive if whenever $V$ is a finite dimensional complex Hilbert space, the mapping $\mathcal{E}\otimes 1_V:L(\mathbb{C}^d\otimes V)\rightarrow L(\mathbb{C}^d\otimes V)$ is positive. We say that a linear mapping $\mathcal{E}:M_d(\mathbb{C})\rightarrow M_d(\mathbb{C})$ is trace preserving if $\text{Tr}(\mathcal{E}(A))=\text{Tr}(A)$ whenever $A\in M_d(\mathbb{C})$. We say that a linear mapping $\mathcal{E}:M_d(\mathbb{C})\rightarrow M_d(\mathbb{C})$ is a quantum channel if it is both completely positive and trace preserving.
If $X_1,\dots,X_r$ are complex matrices, then define a mapping $\Phi(X_1,\dots,X_r):M_{d}(\mathbb{C})\rightarrow M_{d}(\mathbb{C})$ by letting
$$\Phi(X_1,\dots,X_r)(X)=X_1XX_1^*+\dots+X_rXX_r^*.$$ The completely positive mappings from $M_d(\mathbb{C})$ to $M_d(\mathbb{C})$ are precisely the mappings of the form
$\Phi(X_1,\dots,X_r)$, and the quantum channels $\mathcal{E}:M_d(\mathbb{C})\rightarrow M_d(\mathbb{C})$ are precisely the mappings of the form $\Phi(X_1,\dots,X_r)$ where $X_1^*X_1+\dots+X_r^*X_r=1_d$.
If $\mathcal{E}$ is a quantum channel, then $\rho(\mathcal{E})=1$; quantum channels are a lot like stochastic matrices.
I claim that $\rho(X_1\dots X_r)^{2/r}\leq\frac{d}{r}\rho(\Phi(X_1,\dots,X_r))$ whenever $X_1,\dots,X_r$ are $d\times d$-complex matrices.
Let $O_{d,r}$ be the set of all $(X_1,\dots,X_r)\in M_d(\mathbb{C})^r$ such that
$\Phi(X_1,\dots,X_r)$ is not nilpotent. Let $E_{d,r}$ be the set of all $(X_1,\dots,X_r)\in M_d(\mathbb{C})^r$ such that
$\Phi(\lambda BX_1B^{-1},\dots,\lambda BX_rB^{-1})$ is a quantum channel for some complex $\lambda\neq 0$ and invertible matrix $B$. Then by my answer to my other question, $E_{d,r}$ is a dense subset of $O_{d,r}$, and clearly $O_{d,r}$ is dense in $M_d(\mathbb{C})^r$.
Observe that $$\frac{\rho(Y_1\dots Y_r)^{2/r}}{\rho(\Phi(Y_1,\dots,Y_r))}=\frac{\rho(\lambda BY_1B^{-1}\dots \lambda BY_rB^{-1})^{2/r}}{\rho(\Phi(\lambda BY_1B^{-1},\dots,\lambda BY_rB^{-1}))}$$ whenever $\lambda\neq 0$, $B$ is invertible, and $\Phi(Y_1,\dots,Y_r)$ is not nilpotent. Therefore, if $\rho(Z_1\dots Z_r)^{2/r}\leq\frac{d}{r}\rho(\Phi(Z_1,\dots,Z_r))$ whenever $\Phi(Z_1,\dots,Z_r)$ is a quantum channel, then $\rho(Y_1\dots Y_r)^{2/r}\leq\frac{d}{r}\rho(\Phi(Y_1,\dots,Y_r))$ whenever $(Y_1,\dots,Y_r)\in E_{d,r},$
so because $E_{d,r}$ is dense in $M_{d}(\mathbb{C})^{r}$, we will be able to conclude that $\rho(X_1\dots X_r)^{2/r}\leq\frac{d}{r}\rho(\Phi(X_1,\dots,X_r))$ whenever $X_1,\dots,X_r\in M_d(\mathbb{C})^r$.
If $\Phi(Z_1,\dots,Z_r)$ is a quantum channel, then $$d=\text{Tr}(1_d)=\text{Tr}(Z_1^*Z_1+\dots+Z_r^*Z_r)=\|Z_1\|^2_2+\dots+\|Z_r\|_2^2$$ where $\|\cdot\|_2$ denotes the Frobenius norm.
Therefore, by the arithmetic-geometric mean inequality, we have our main inequality, namely
$$\rho(Z_1\dots Z_r)^{2/r}\leq\|Z_1\dots Z_r\|_\infty^{2/r}\leq(\|Z_1\|_\infty^2\dots\|Z_r\|_\infty^2)^{1/r}\leq(\|Z_1\|_2^2\dots\|Z_r\|_2^2)^{1/r}\leq\frac{\|Z_{1}\|_2^2+\dots+\|Z_r\|_2^2}{r}
=\frac{d}{r}=\frac{d}{r}\rho(\Phi(Z_1,\dots,Z_r)).$$
Furthermore, if $\rho(Z_1\dots Z_r)^{2/r}=\frac{d}{r}\rho(\Phi(Z_1,\dots,Z_r))$,
then $\|Z_j\|_\infty^2=\|Z_j\|^2_2=d/r$ for all $j$. However, The only way to get $\|Z_{j}\|_\infty=\|Z_j\|_2$ is if $\text{Rank}(Z_j)=1$ since $\|W\|_{\infty}=\|W\|_2$ precisely when $\text{Rank}(W)=1$.
Therefore, for all $j$, there are column vectors $u_j,v_j$ such that $Z_j=u_jv_j^*$.
In this case, $$\|Z_j\|_2^2=\text{Tr}((u_jv_j^*)(u_jv_j^*)^*)=\text{Tr}(u_jv_j^*v_ju_j^*)=
\text{Tr}(u_j^*u_jv_j^*v_j)=\|u_j\|^2\cdot\|v_j\|^2.$$ Therefore, $\|Z_j\|_2=\|u_j\|\cdot\|v_j\|=\sqrt{d/r}$.
$$\rho(Z_1\dots Z_r)=|\text{Tr}(Z_1\dots Z_r)|=|\text{Tr}(u_1v_1^*\dots u_rv_r^*)|
=|\text{Tr}(v_1^*u_2\dots v_r^*u_1)|$$
$$=|\langle v_1,u_2\rangle|\dots|\langle v_r,u_1\rangle|.$$
Since $\rho(Z_1\dots Z_r)=\|Z_1\|_\infty\dots\|Z_r\|_\infty$, we have
$$|\langle v_1,u_2\rangle|\dots|\langle v_r,u_1\rangle|=\|u_1\|\cdot\|v_1\|\dots \|u_r\|\cdot\|v_r\|.$$
This means that there are constants $\gamma_1,\dots,\gamma_r$ where
$v_j=\overline{\gamma_j}\cdot u_{j+1}$. Therefore,
$Z_j=u_jv_j^*=u_j\gamma_j u_{j+1}^*=\gamma_j u_ju_{j+1}^*$.
Now, set $w_j=\|u_j\|\cdot v_{j}=\|u_j\|\cdot\overline{\gamma_j}u_{j+1}.$
Then $$1_d=\sum_{j=1}^{r}Z_j^*Z_j=\sum_{j=1}^{r}(u_jv_j^*)^*u_jv_j^*=\sum_{j=1}^{r}v_ju_j^*u_jv_j^*=\sum_{j=1}^{r}\|u_{j}\|^2\cdot v_jv_j^*=\sum_{j=1}^{r}w_jw_j^*.$$
Observe that $1_d=\sum_{j=1}^{r}w_jw_j^*$ precisely when $x=\sum_{j=1}^{r}w_{j}\cdot\langle x,w_j\rangle$. We say that a tuple of vectors $(h_1,\dots,h_r)$ is a normalized tight frame if $1_d=\sum_{j=1}^{r}h_jh_j^*$. We say that a normalized tight frame is equal-norm if $\|h_1\|=\dots=\|h_r\|$. Observe that $\|w_j\|^2=d/r$ for all $j$, so $(w_1,\dots,w_r)$ is an equal-norm tight frame.
One can also obtain tuples $(Z_1,\dots,Z_r)$ where $\Phi(Z_1,\dots,Z_r)$ is a quantum channel and where $\rho(Z_1,\dots,Z_r)^{2/r}=d/r$ from equal-norm tight frames. Suppose that $(w_1,\dots,w_r)\in(\mathbb{C}^d)^r$ is an equal-norm tight frame. Then we necessarily have $\|w_j\|^2_2=d/r$ for all $j$. Now, suppose that $|\alpha_j\beta_j|^2=r/d$ for $1\leq j\leq r$. Let $v_j=\alpha_jw_j,u_{j+1}=\beta_{j+1}w_j$, and as always, set $Z_j=u_jv_j^*$ for $1\leq j\leq r$. Then the reader can verify that $\Phi(Z_1,\dots,Z_r)$ is a quantum channel with $\rho(Z_1,\dots,Z_r)^{2/r}=\frac{d}{r}$ and every quantum channel $\Phi(Z_1,\dots,Z_r)$ with $\rho(Z_1,\dots,Z_r)^{2/r}=\frac{d}{r}$ is of this form.
The existence of equal-norm tight frames $(w_1,\dots,w_r)\in(\mathbb{R}^d)^r$ is already a known result. For example, the existence of such frames is Corollary 7.1 in An Introduction to Finite Tight Frames by Shayne F. D. Waldron.
By applying continuity, we can actually show that $\text{Rank}(X_j)=1$ even if we drop the assumption that $\Phi(X_1,\dots,X_r)$ is a quantum channel.
Let $f:O_{d,r}\rightarrow\mathbb{R}$ be the function where
$f(X_1,\dots,X_r)=\frac{\rho(X_1\dots X_r)^{2/r}}{\rho(\Phi(X_1,\dots,X_r))}.$
Let $g^-:M_d(\mathbb{C})\rightarrow\mathbb{R}$ be the function where
$g^-(X)=|\lambda_1\lambda_2|$ where $\lambda_1,\dots,\lambda_d$ are the eigenvalues of $X$ arranged in an order so that $|\lambda_1|\geq|\lambda_2|\geq\dots\geq|\lambda_d|$.
Define $g:O_{d,r}\rightarrow\mathbb{R}$ by letting $$g(X_1,\dots,X_r)=\frac{g^-(X_1)+\dots+g^-(X_r)}{r\cdot\rho(\Phi(X_1,\dots,X_r))}.$$
Suppose that $(X_1,\dots,X_r)\in O_{d,r}$ and $f(X_1,\dots,X_r)=\frac{d}{r}$. Then for each $\delta>0$, there is a neighborhood $U$ of $(X_1,\dots,X_r)$ where if
$(Y_1,\dots,Y_r)\in U$, then $f(Y_1,\dots,Y_r)>\frac{d}{r}-\delta$. Now, let $(Y_1,\dots,Y_r)\in U\cap E_{d,r}$. Then there is some $\lambda$ and invertible $B$ where if we set $Z_j=\lambda BY_jB^{-1}$ for all $j$, then $\Phi(Z_1,\dots,Z_r)$ is a quantum channel. Now, let $\lambda_{j,1},\dots,\lambda_{j,d}$ be the eigenvalues of $Z_j$ ordered in a way so that $|\lambda_{j,1}|\geq|\lambda_{j,2}|\geq\dots\geq|\lambda_{j,d}|$. Let $\sigma_{j,1}\geq\dots\geq\sigma_{j,d}$ denote the singular values of $Z_j$.
Then
$$\frac{d}{r}-\delta<f(Y_1,\dots,Y_r)=f(Z_1,\dots,Z_r)=\rho(Z_1\dots Z_r)^{2/r}
\leq(\sigma_{1,1}^2\dots\sigma_{r,1}^2)^{1/r}$$
$$\leq\frac{\sigma_{1,1}^2+\dots+\sigma_{r,1}^2}{r}
\leq\frac{(\sigma_{1,1}^2+\sigma_{1,2}^2)+\dots+(\sigma_{r,1}^2+\sigma_{r,2}^2)}{r}$$
$$\leq\frac{\|Z_1\|_2^2+\dots+\|Z_r\|_2^2}{r}=\frac{d}{r}.$$
We conclude that $\sigma_{j,1}^2\leq\frac{d}{r}$ and $\sigma_{j,2}^2\leq\delta$. Therefore, by Weyl's inequality, we have
$$|\lambda_{j,1}\lambda_{j,2}|\leq\sigma_{i,1}\sigma_{j,2}\leq\sqrt{\frac{d\delta}{r}}.$$
Therefore, we have
$$g(Y_1,\dots,Y_r)=g(Z_1,\dots,Z_r)\leq\sqrt{\frac{d\delta}{r}}.$$
We therefore conclude that $g(X_1,\dots,X_r)=0$ which means that
$\text{Rank}(X_j)\leq 1$ for $1\leq j\leq r$.
