Yes. The equality can in fact be reached. Our strategy will be to produce a compact set $K_{d,r}\subseteq M_{d}(\mathbb{C})^{r}$ such that if $(X_1,\dots,X_r)\in K_{d,r}$, then $\rho(\Phi(X_1,\dots,X_r))=1$, and where
$$\rho_{2,d}(A_1,\dots,A_r)=\max\{\rho(A_1\otimes X_1+\dots+A_r\otimes X_r)\mid(X_1,\dots,X_r)\in K_{d,r}\}.$$
Let us go over a few definitions and facts to give some context to our construction. These facts can easily be found in John Watrous's 2018 book called The Theory of Quantum Information.
A mapping $\mathcal{E}:L(V)\rightarrow L(W)$ is said to be positive if whenever $X$ is positive semidefinite, $\mathcal{E}(X)$ is also positive semidefinite. A mapping $\mathcal{E}:L(V)\rightarrow L(W)$ is said to be completely positive if the mapping $\mathcal{E}\otimes 1_{U}:L(V\otimes U)\rightarrow L(W\otimes U)$ is positive for each finite dimensional complex Hilbert space $U$.
A mapping $\mathcal{E}:L(V)\rightarrow L(W)$ is said to be trace preserving if $\text{Tr}(\mathcal{E}(X))=\text{Tr}(X)$ for each $X\in L(V)$.
Proposition: Let $\mathcal{E}:L(V)\rightarrow L(W)$ be a linear operator.
$\mathcal{E}$ is completely positive if and only if there are $A_1,\dots,A_r$ where $\mathcal{E}=\Phi(X_1,\dots,X_r)$.
If $\mathcal{E}$ is defined by letting $\mathcal{E}(X)=A_1XB_1^*+\dots+A_rXB_r^*$, then $\mathcal{E}$ is trace preserving if and only if $A_1^*B_1+\dots+A_r^*B_r=1_V$.
A quantum channel is a completely positive trace preserving operator $\mathcal{E}:L(V)\rightarrow L(W)$. In quantum information theory, the quantum channels are the main morphisms between quantum states.
Proposition: If $\mathcal{E}$ is a quantum channel, then $\rho(\mathcal{E})=1$.
Let $K_{d,r}$ be the collection of all tuples $(X_1,\dots,X_r)\in M_{d}(\mathbb{C})^r$ such that $X_1^*X_1+\dots+X_r^*X_r=1_d$. Said differently, $K_{d,r}$ is the collection of all tuples $(X_1,\dots,X_r)\in M_{d}(\mathbb{C})^r$ such that
$\Phi(X_1,\dots,X_r)$ is a quantum channel. The set $K_{d,r}$ is clearly a closed set, and $K_{d,r}$ is bounded since $d=\text{Tr}(X_1^*X_1+\dots+X_r^*X_r)=\|X_1\|_2^2+\dots+\|X_r\|_2^2$ where $\|\cdot\|_2$ denotes the Frobenius norm, so $K_{d,r}$ is compact.
Lemma: Let $A_1,\dots,A_r\in L(V)$. Suppose that there is no $x\in V\setminus\{0\}$ with $A_1x=\dots=A_rx=0$. Furthermore, suppose that there is no subspace $W\subseteq V$ with $W\neq\{0\},W\neq V$, and $W=A_1[W]+\dots+A_r[W]$. Then there is a $\lambda>0$ along with a positive definite $P$ with $\Phi(A_1,\dots,A_r)(P)=\lambda P$.
Proof: Now, let $\mathcal{Q}$ be the collection of all positive semidefinite matrices in $L(V)$ with trace $1$, and let $F:\mathcal{Q}\rightarrow\mathcal{Q}$ be the mapping defined by letting $$F(P)=\frac{\Phi(A_1,\dots,A_r)(P)}{\text{Tr}\big(\Phi(A_1,\dots,A_r)(P)\big)}.$$ Then $\mathcal{Q}$ is convex, and $F$ is a continuous bijection, so by the Brouwer fixed point theorem, there is some $P\in\mathcal{Q}$ with
$F(P)=P$. Therefore, we have $\Phi(A_1,\dots,A_r)(P)=\lambda P$ for some positive $\lambda$.
Now,
$$\text{Im}(P)=\text{Im}(\lambda P)=\text{Im}(A_1PA_1^*+\dots+A_rPA_r^*)$$
$$=\text{Im}(A_1PA_1^*)+\dots+\text{Im}(A_rPA_r^*)=\text{Im}(A_1P)+\dots+\text{Im}(A_rP)$$
$$=A_1[\text{Im}(P)]+\dots+A_r[\text{Im}(P)].$$
This is only possible if $\text{Im}(P)=V$. $\square$
Let $O_{d,r}$ be the collection of all $(X_1,\dots,X_r)\in M_{d}(\mathbb{C})^{r}$ where $\Phi(X_1,\dots,X_r)$ is not nilpotent. Let $E_{d,r}$ be the collection of all $(X_1,\dots,X_r)\in M_d(\mathbb{C})^r$ where there is no subspace $W\subseteq\mathbb{C}^d$ with
$W=X_1^*[W]+\dots+X_r^*[W]$ and $W\neq\{0\},W\neq\mathbb{C}^d$ and where there is no $x$ with $X_j^*x=0$ for $1\leq j\leq r$.
The sets $O_{d,r},E_{d,r}$ are dense subsets of $M_{d}(\mathbb{C})^{r}$ with $E_{d,r}\subseteq O_{d,r}$.
Proposition: Suppose that $(X_1,\dots,X_r)\in E_{d,r}$. Then there is an invertible matrix $B$ and some positive number $\lambda$ with $(\lambda BX_1B^{-1},\dots,\lambda BX_rB^{-1})\in K_{d,r}$.
Proof: Suppose that $B$ is invertible and $\lambda$ is positive. Then $(\lambda BX_1B^{-1},\dots,\lambda BX_rB^{-1})\in K_{d,r}$ if and only if
$$\sum_{k=1}^{r}\lambda^2 (B^{-1})^*X_k^*B^*BX_k^*B^{-1}=I$$ if and only if
$$\sum_{k=1}^{r}\lambda^2 X_k^*B^*BX_k=B^*B$$.
Therefore, there are $\lambda,B$ with $(\lambda BX_1B^{-1},\dots,\lambda BX_rB^{-1})\in K_{d,r}$ if and only if there is a positive $\mu$ and a positive definite $P$ with $$\mu\Phi(X_1^*,\dots,X_k^*)P=P.$$
On the other hand, the existence of such a positive definite $P$ and positive $\mu$ is guaranteed by the above lemma. $\square$
Now, define a mapping $G_{A_1,\dots,A_r}:O_{d,r}\rightarrow\mathbb{R}$ by letting
$$G_{A_1,\dots,A_r}(X_1,\dots,X_r)=\frac{\rho(A_1\otimes X_1+\dots+A_r\otimes X_r)}{\rho(\Phi(X_1,\dots,X_r))^{1/2}}.$$ Since the mapping $G$ is continuous, we have
$$\rho_{2,d}(A_1,\dots,A_r)=\sup\{G_{A_1,\dots,A_r}(X_1,\dots,X_r)\mid (X_1,\dots,X_r)\in E_{d,r}\}.$$
Since $$G_{A_1,\dots,A_r}(X_1,\dots,X_r)=G_{A_1,\dots,A_r}(\lambda BX_1B^{-1},\dots,\lambda BX_rB^{-1})$$ whenever $B$ is invertible and $\lambda$ is a non-zero complex number, by the above proposition, we know that
$$\{G_{A_1,\dots,A_r}(X_1,\dots,X_r)\mid (X_1,\dots,X_r)\in E_{d,r}\}$$
$$=\{G_{A_1,\dots,A_r}(X_1,\dots,X_r)\mid (X_1,\dots,X_r)\in K_{d,r}\}.$$
Therefore, since $K_{d,r}$ is compact, there is some $(Z_1,\dots,Z_r)\in K_{d,r}$ with
$$G_{A_1,\dots,A_r}(Z_1,\dots,Z_r)=\max\{G_{A_1,\dots,A_r}(X_1,\dots,X_r)\mid (X_1,\dots,X_r)\in K_{d,r}\}$$
$$=\rho_{2,d}(A_1,\dots,A_r).$$
I have ran computations that maximize $G_{A_1,\dots,A_r}$, and in these computations the maximum seems to actually be reached.
Yes. We can characterize $\rho_{2,d}(\mathcal{E})$ whenever $\mathcal{E}$ is completely positive without needing to first decompose $\mathcal{E}$ as $\Phi(A)$ or $\Phi(A_1,\dots,A_r).$ As a consequence, we can define $\rho_{2,d}(\mathcal{E})$ for all linear operators $\mathcal{E}:L(V)\rightarrow L(V)$, but if $\mathcal{E}$ is not completely positive, then $\rho_{2,d}(\mathcal{E})$ is usually infinite, so $\rho_{2,d}(\mathcal{E})$ it not well behaved when we do not assume that $\mathcal{E}$ is completely positive.
Let $U_1,U_2,V_1,V_2,U_1^\sharp,U_2^\sharp,V_1^\sharp,V_2^\sharp,U,V$ be finite dimensional complex inner product spaces.
If $A_1,\dots,A_r:U_2\rightarrow V_2,B_1,\dots,B_r:U_1\rightarrow V_1$ are linear, then define a mapping $\Gamma(A_1,\dots,A_r;B_1,\dots,B_r):L(U_1,U_2)\rightarrow L(V_1,V_2)$ by letting $$\Gamma(A_1,\dots,A_r;B_1,\dots,B_r)(X)=\sum_{k=1}^rA_kXB_k^*.$$
Suppose now that $A_1,\dots,A_r:U_2\rightarrow V_2,B_1,\dots,B_r:U_1\rightarrow V_1$. Let $S:L(U_2,V_2)\rightarrow L(U_2^\sharp,V_2^\sharp)$ be linear and let
$T:L(U_1,V_1)\rightarrow L(U_1^\sharp,V_1^\sharp)$ also be linear. Then define
$\mho(\Gamma(A_1,\dots,A_r;B_1,\dots,B_r),S,T):L(U_1^\sharp,U_2^\sharp)\rightarrow L(V_1^\sharp,V_2^\sharp)$ by letting
$$\mho(\Gamma(A_1,\dots,A_r;B_1,\dots,B_r),S,T)=\Gamma(S(A_1),\dots,S(A_r);T(B_1),\dots,T(B_r)).$$ The mapping $\mho$ is well-defined; i.e. it depends on $\Gamma(A_1,\dots,A_r;B_1,\dots,B_r)$ rather than the particular choice of $A_1,\dots,A_r;B_1,\dots,B_r$.
$\rho_{2,d}(\mathcal{E})$ is the maximum value of $$\frac{\rho(\mho(\mathcal{E},1_{L(V)},T))^2}{\rho(\mho(\mathcal{E},T,T))}$$
where $T:L(V)\rightarrow L(U)$ is a linear operator and $\dim(U)=d$, and this definition generalizes to all linear operators $\mathcal{E}:L(V)\rightarrow L(V)$ when we replace the word 'maximum' with 'supremum'.
Suppose now that $\mathcal{E}=\Gamma(A_1,\dots,A_r;B_1,\dots,B_r)$ where $A_1,\dots,A_r,B_1,\dots,B_r$ are linearly independent and $\Phi(A_1,\dots,A_r)$ is not nilpotent. Let $T:L(V)\rightarrow L(V)$ be a linear mapping where $T(A_j)=A_j,T(B_j)=\alpha\cdot A_j$ for $1\leq j\leq r$.
Then
$\mho(\mathcal{E},1_{L(V)},T)=\mho(\mathcal{E},T,T)=\alpha\cdot\Phi(A_1,\dots,A_r)$.
Therefore, $$\rho_{2,d}(\mathcal{E})\geq\frac{\rho(\alpha\cdot\Phi(A_1,\dots,A_r))^2}{\rho(\alpha\cdot\Phi(A_1,\dots,A_r))}=\alpha\cdot\rho(\Phi(A_1,\dots,A_r)).$$
Since $\alpha$ can be made arbitrarily large, we have $\rho_{2,d}(\mathcal{E})=\infty$ in this case.
A modest generalization: added 8/15/2022
There is another way to generalize $\rho_{2,d}(\mathcal{E})$ to some operators that are not completely positive, but this generalization is based on a conjecture.
Conjecture: Suppose that $\mathcal{E}:L(V)\rightarrow L(V)$ is completely positive and $\alpha\geq 0$. Then $\rho_{2,d}(\mathcal{E}+\alpha\cdot 1_{L(V)})=\rho_{2,d}(\mathcal{E})+\alpha$.
From this conjecture, we can define $\rho_{2,d}(\mathcal{E})$ whenever there is some $\alpha\geq 0$ where $\mathcal{E}+\alpha\cdot 1_{L(V)}$ is completely positive by letting $\rho_{2,d}(\mathcal{E})=\rho_{2,d}(\mathcal{E}+\alpha\cdot 1_{L(V)})-\alpha$. Here, we can have negative values of $\rho_{2,d}(\mathcal{E})$, so $\rho_{2,d}$ more closely resembles the maximum value of a Hermitian operator than the spectral radius.
Projective mappings: added 3/21/2023
Define $\rho_{2,d}^P(A_1,\dots,A_r)$ to be the supremum of all values of the form $$\frac{\rho(A_1\otimes\overline{RA_1S}+\dots+A_r\otimes\overline{RA_rS})}{\rho(RA_1S\otimes\overline{RA_1S}+\dots+RA_rS\otimes\overline{RA_rS})^{1/2}}$$
where $R\in M_{d,n}(\mathbb{C}),S\in M_{n,d}(\mathbb{C})$. I conjecture that $\rho_{2,d}^P(A_1,\dots,A_r)=\rho_{2,d}(A_1,\dots,A_r)$ (computer calculations support this conjecture at least sometimes), but it seems like $\rho_{2,d}^P(A_1,\dots,A_r)$ is easier to work with than $\rho_{2,d}(A_1,\dots,A_r)$.
If $S:U\rightarrow V,R:V\rightarrow U$ are linear maps, then define a mapping $T_{R,S}:L(V)\rightarrow L(U)$ by setting $T_{R,S}(A)=RAS$ for $A\in L(V)$.
$$\mho(\Gamma(A_1,\dots,A_r;B_1,\dots,B_r),T_{R,S},T_{R,S})(X)$$
$$=\Gamma(T_{R,S}(A_1),\dots,T_{R,S}(A_r);T_{R,S}(B_1),\dots,T_{R,S}(B_r))(X)$$
$$=\Gamma(RA_1S,\dots,RA_rS;RB_1S,\dots,RB_rS)(X)=\sum_{k=1}^{r}RA_kSX(RB_kS)^*$$
$$=\sum_{k=1}^rRA_kSXS^*B^*_kR^*
=R\big(\Gamma(A_1,\dots,A_r;B_1,\dots,B_r)(SXS^*)\big)R^*.$$
Similarly,
$$\mho(\Gamma(A_1,\dots,A_r;B_1,\dots,B_r),1_{L(V)},T_{R,S})(X)
=\Gamma(A_1,\dots,A_r;RB_1S,\dots,RB_rS)(X)$$
$$=\sum_{k=1}^rA_kX(RB_kS)^*=\sum_{k=1}^rA_kXS^*B_k^*R^*
=\big(\Gamma(A_1,\dots,A_r;B_1,\dots,B_r)(XS^*)\big)R^*.$$
Therefore, $\mho(\mathcal{E},T_{R,S},T_{R,S})(X)=R\big(\mathcal{E}(SXS^*)\big)R^*$ and $\mho(\mathcal{E},1_{L(V)},T_{R,S})(X)=\big(\mathcal{E}(XS^*)\big)R^*$ whenever $\mathcal{E}:L(V)\rightarrow L(V)$ and $R,S,X$ are suitable linear operators. Therefore, if $\mathcal{E}:L(V)\rightarrow L(V)$, then define $\mho_{R,S}^-(\mathcal{E})(X)=\big(\mathcal{E}(XS^*)\big)R^*$ and
$\mho_{R,S}^+(\mathcal{E})(X)=R\big(\mathcal{E}(SXS^*)\big)R^*$.
$\rho_{2,d}^P(A_1,\dots,A_r)$ is the supremum of all values of the form
$$\frac{\rho(\mho_{R,S}^-(\Phi(A_1,\dots,A_r)))}{\rho(\mho_{R,S}^+(\Phi(A_1,\dots,A_r))^{1/2}}.$$ Therefore, we may define $\rho_{2,d}^P(\mathcal{E})$ for all completely positive operators $\mathcal{E}:L(V)\rightarrow L(V)$ by setting
$\rho_{2,d}^P(\mathcal{E})=\rho_{2,d}^P(A_1,\dots,A_r)^2$ whenever $\mathcal{E}=\Phi(A_1,\dots,A_r)$. Then $\rho_{2,d}^P(\mathcal{E})$ is the supremum of all values of the form $\frac{\rho(\mho_{R,S}^-(\mathcal{E}))^2}{\rho(\mho_{R,S}^+(\mathcal{E}))}$ where $R\in L(U,V),S\in L(V,U)$ and $U$ is a complex inner product space $\dim(U)=d$. This definition of $\rho_{2,d}^P(\mathcal{E})$ makes sense even when $\mathcal{E}:L(V)\rightarrow L(V)$ and $V$ is an infinite dimensional complex Hilbert space (we will probably need to use Stinespring's dilation theorem to generalize $\rho_{2,d}$ to infinite dimensional spaces).
Example: $\rho_{2,d}^P(\mathcal{E})$ can be infinite for superoperators which are positive but not completely positive. For example, the transpose map $T:M_n(\mathbb{C})\rightarrow M_n(\mathbb{C})$ defined by $T(X)=X^T$ is positive but not completely positive whenever $n>1$, and we shall prove that $T$ is not completely positive by showing that $\rho_{2,d}^P(T)=\infty$. A straightforward calculation yields $\rho_{2,d}^P(T)=\sup_{R\in M_{d,n}(\mathbb{C}),S\in M_{n,d}(\mathbb{C})}\frac{\rho(SR)^2}{\rho(R\overline{S}\overline{R}S)}.$
If $R,S$ are rank-1 matrices, then we can set $R=uv^*,S=wx^*$ for vectors $u,v,w,x$. In this case,
$\frac{\rho(SR)^2}{\rho(R\overline{S}\overline{R}S)}=|\frac{v^*w}{v^Tw}|^2$. If $v=w=[1,i]^T$, then $\frac{\rho(SR)^2}{\rho(R\overline{S}\overline{R}S)}=(2/0)^2=+\infty.$
It seems like a good way of showing that a superoperator $\mathcal{E}:L(V)\rightarrow L(V)$ is not completely positive is to show that $\rho_{2,d}^P(\mathcal{E})>\rho(\mathcal{E})$ by finding matrices $R,S$ where $\rho(\mho_{R,S}^-(\mathcal{E}))^2>\rho(\mho_{R,S}^+(\mathcal{E}))\cdot\rho(\mathcal{E})$.
Best Answer
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$.