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.
I claim that there is a somewhat abstract notion of a multi-spectral radius and that there is probably an abstract theory behind this abstract notion. I will try to justify this abstract multi-spectral radius by showing that it captures the specific examples of multi-spectral radii that I have mentioned in the question and that the simplest examples of these multi-spectral radii are reasonable mathematical objects. With that being said, there are notions of a multi-spectral radius that I have not shown to fit within this framework, so more research on this topic is needed.
Suppose that $A$ is a complex Banach algebra. We say that a function
$\rho:A^r\rightarrow[0,\infty)$ is a multi-spectral radius if there is an isometric embedding $\iota:A\rightarrow B$ along with a bounded subset $\mathcal{C}\subseteq B^r$ where
$x_j\iota(a)=\iota(a)x_j$ whenever $a\in A,(x_1,\dots,x_r)\in\mathcal{C}$,
if $(x_1,\dots,x_r)\in\mathcal{C}$, then $(\lambda_1x_1,\dots,\lambda_rx_r)\in\mathcal{C}$ whenever $|\lambda_j|=1$ for $1\leq j\leq r$, and
$$\rho(a_1,\dots,a_r)=\rho_{\iota,\mathcal{C}}(a_1,\dots,a_r)=\sup_{(x_1,\dots,x_r)\in\mathcal{C}}\rho(x_1\iota(a_1)+\dots+x_r\iota(a_r))$$ whenever $a_1,\dots,a_r\in A$.
One may also want to require that if $(x_1,\dots,x_r)\in\mathcal{C}$, then $(\lambda_1x_1,\dots,\lambda_rx_r)\in\mathcal{C}$ whenever $|\lambda_j|=1$ for $1\leq j\leq r$, but this condition is not necessary. One can show that if $a,b\in B$, then the mapping from $\mathbb{C}$ to $[-\infty,\infty)$ defined by
$\lambda\mapsto \ln(\rho(a+\lambda b))$ is subharmonic, so by the maximum principle,
$\max\{\rho(\lambda_1x_1\iota(a_1)+\dots+\lambda_rx_r\iota(a_r)):|\lambda_1|=\dots=|\lambda_r|=1\}=\max\{\rho(\lambda_1x_1\iota(a_1)+\dots+\lambda_rx_r\iota(a_r)):|\lambda_1|\leq 1,\dots,|\lambda_r|\leq 1\}.$
The $L_1$-spectral radius can be characterized in terms of our framework.
Theorem: $\rho_1(a_1,\dots,a_r)$ is the maximum value of $\rho(x_1\iota(a_1)+\dots+x_r\iota(a_r))$ where $\iota:A\rightarrow B$ is an isometric embedding of Banach algebras, and $\|x_j\|\leq 1$ for $1\leq j\leq r$.
The proof of the above result is not too hard, and I have given a proof of the above result in this answer.
We say that a multi-spectral radius $\rho$ is unitary invariant if
$\rho(a_1,\dots,a_r)=\rho(b_1,\dots,b_r)$ whenever there is an $n\times n$-unitary matrix $(u_{i,j})_{i,j}$ where $b_j=\sum_{i=1}^ru_{i,j}a_i$ for $1\leq j\leq r$. The following lemma is a standard result from quantum information theory.
Lemma: Suppose that $A_1,\dots,A_r,B_1,\dots,B_r\in M_n(\mathbb{C})$. Then $\Phi(A_1,\dots,A_r)=\Phi(B_1,\dots,B_r)$ if and only if there is an $r\times r$-unitary matrix $(u_{i,j})_{i,j}$ where $B_j=\sum_{i=1}^ru_{i,j}A_i$ for $1\leq j\leq r$.
Therefore, a multi-spectral radius $\rho:M_n(\mathbb{C})^r\rightarrow[0,\infty)$ is unitary invariant if and only if $\rho(A_1,\dots,A_r)=\rho(B_1,\dots,B_r)$ whenever
$\Phi(A_1,\dots,A_r)=\Phi(B_1,\dots,B_r)$. The continuous unitary invariant multi-spectral radii are completely determined by the mapping $\Phi(A_1,\dots,A_r)\mapsto\rho(A_1,\dots,A_r)$ where $\Phi(A_1,\dots,A_r)$ is completely positive and trace preserving (a completely positive trace preserving map is known as a quantum channel).
The following easy lemmas show that how we can always upgrade a multi-spectral radius to a unitary invariant multi-spectral radius.
Lemma: Let $A$ be an algebra over a field $K$. Suppose that $(a_1,\dots,a_r),(b_1,\dots,b_r),(x_1,\dots,x_r),(y_1,\dots,y_r)\in A^r$. Let $(u_{i,j})_{i,j},(v_{i,j})_{i,j}\in M_r(K)$ be inverse matrices. Suppose that
$a_k=\sum_{i=1}^ru_{i,k}b_i$ and $x_k=\sum_{j=1}^rv_{k,j}y_j$ for $1\leq k\leq r$. Then $$\sum_{k=1}^ra_kx_k=\sum_{i=1}^rb_iy_j.$$
Lemma: Suppose that $K$ is a field and $A$ is an algebra over $K$. Let $(u_{i,k})_{i,k}\in M_r(K)$. Suppose furthermore that $a_1,\dots,a_r,b_1,\dots,b_r,x_1,\dots,x_r,y_1,\dots,y_r\in A$ and
$a_k=\sum_{i=1}^ru_{i,k}b_i$ for $1\leq k\leq r$ and $x_i=\sum_{k=1}^ru_{i,k}y_k$ for $1\leq k\leq r$. Then $\sum_{k=1}^ra_ky_k=\sum_{i=1}^rb_ix_i$.
Proposition: Let $A,B$ be Banach algebras. Let $\iota:A\rightarrow B$ be an isometric embedding. Suppose that $\mathcal{C}\subseteq B^r$ is a bounded subset with $x_j\iota(a)=\iota(a)x_j$ for $1\leq j\leq r$. Let $\mathcal{D}$ be the collection of all tuples $(y_1,\dots,y_r)$ where there is some $r\times r$-unitary matrix $(u_{i,j})_{i,j}$ and $(x_1,\dots,x_r)\in\mathcal{C}$ where
$y_j=\sum_{i=1}^ru_{i,j}x_i$ for $1\leq i\leq r$. Then
$$\rho_{\iota,\mathcal{D}}(x_1,\dots,x_r)=\sup\{\rho_{\iota,\mathcal{C}}(\sum_ju_{1,j}x_j,\dots,\sum_ju_{r,j}x_j)\mid (u_{i,j})_{i,j}\in U(r)\}.$$
By using the following version of Holder's inequality that can be proven using the classical Holder's inequality, we can show that the $L_2$-spectral radius is a multi-spectral radius.
Theorem: $\rho(A_1\otimes B_1+\dots+A_r\otimes B_r)\leq
\rho_p(A_1,\dots,A_r)\cdot\rho_q(B_1,\dots,B_r)$ whenever $p,q\in(1,\infty)$ and $\frac{1}{p}+\frac{1}{q}=1$.
As a consequence, if $d\geq n$ and $A_1,\dots,A_r\in M_n(\mathbb{C})$, then
$$\rho_2(A_1,\dots,A_r)=\max_{(X_1,\dots,X_r)\in M_d(\mathbb{C})}\frac{\rho(A_1\otimes X_1+\dots+A_r\otimes X_r)}{\rho_2(X_1,\dots,X_r)}.$$
We have another construction that allows us to show that the $L_p$-spectral radius is a multi-spectral radius for $1\leq p<\infty$. Suppose now that $1\leq p<\infty$.
Now, let $A$ be a Banach algebra. Let $x_1,\dots,x_r$ be non-commutating variables. Let $B$ be the collection of all sums of the form
$\sum_{k=0}^n\sum_{i_1,\dots,i_k\in\{1,\dots,r\}}a_{i_1,\dots,i_k}x_{i_1}\dots x_{i_k}$. We observe that for $p>1$ the Banach space $\ell^p$ indexed with the natural numbers cannot be endowed with a convolution operation since $(1/n)_{n=1}^{\infty}*(1/n)_{n=1}^{\infty}=(+\infty)_{n=1}^\infty$. We can give $B$ a norm that combines the $\ell^p$ and the $\ell^1$ norms that makes the completion of $B$ into a Banach algebra.
Then give $B$ the norm $$\|\sum_{k=0}^n\sum_{i_1,\dots,i_k\in\{1,\dots,r\}}a_{i_1,\dots,i_k}x_{i_1}\dots x_{i_k}\|=\sum_{k=0}^n\|(a_{i_1,\dots,i_k})_{i_1,\dots,i_k}\|_p.$$
Give $B$ the multiplication defined by bilinearity along with the condition that $$(a\cdot x_{i_1}\dots x_{i_m})\cdot (b\cdot x_{j_1}\dots x_{j_n})=
ab\cdot x_{i_1}\dots x_{i_m}x_{j_1}\dots x_{j_n}.$$
In other words, each element in $A$ commutes with each variable $x_j$, but we do not impose any other version of commutativity.
$B$ is submultiplicative: Let $$u=\sum_{j=0}^\infty\sum_{i_1,\dots,i_j\in\{1,\dots,r\}}a_{i_1,\dots,i_j}x_{i_1}\dots x_{i_j}$$ and let $$v=\sum_{j=0}^\infty\sum_{i_1,\dots,i_j\in\{1,\dots,r\}}a_{i_1,\dots,i_j}x_{i_1}\dots x_{i_j}$$ where only finitely many terms of these 'non-commutative polynomials' are non-zero.
Then
$$\|u\cdot v\|$$
$$=\|(\sum_{k=0}^\infty\sum_{i_1,\dots,i_k\in\{1,\dots,r\}}a_{i_1,\dots,i_k}x_{i_1}\dots x_{i_k})\cdot (\sum_{k=0}^\infty\sum_{i_1,\dots,i_k\in\{1,\dots,r\}}b_{i_1,\dots,i_k}x_{i_1}\dots x_{i_k})\|$$
$$=\|\sum_{k=0}^{\infty}\sum_{j=0}^k\sum_{i_1,\dots,i_j\in\{1,\dots,r\}}\sum_{i_{j+1},\dots,i_k}a_{i_1,\dots,i_j}b_{i_{j+1},\dots,i_k}x_{i_1}\dots x_{i_k}\|$$
$$\leq\sum_{k=0}^{\infty}\sum_{j=0}^k\|\sum_{i_1,\dots,i_k\{1,\dots,r\}}a_{i_1,\dots,i_j}b_{i_{j+1},\dots,i_k}x_{i_1}\dots x_{i_k}\|$$
$$=\sum_{k=0}^{\infty}\sum_{j=0}^k\|(a_{i_1,\dots,i_j}\cdot b_{i_{j+1},\dots, i_k})_{i_1,\dots,i_k\in\{1,\dots,r\}}\|_p$$
$$\leq\sum_{k=0}^{\infty}\sum_{j=0}^k\|(a_{i_1,\dots,i_j})_{i_1,\dots,i_j\in\{1,\dots,r\}}\|_p\cdot \|(b_{i_{j+1},\dots,i_k})_{i_{j+1},\dots,i_k\in\{1,\dots,r\}}\|_p$$
$$=\sum_{j=0}^\infty\|(a_{i_1,\dots,i_j})_{i_1,\dots,i_j\in\{1,\dots,r\}}\|_p\cdot\sum_{k=0}^\infty\|(b_{i_1,\dots,i_k})_{i_1,\dots,i_k\in\{1,\dots,r\}}\|_p=\|u\|\cdot\|v\|.$$
Therefore, the completion $\overline{B}$ of $B$ is a Banach algebra, and the original Banach algebra $A$ embeds into $\overline{B}$. In this case, we simply have
$\rho_p(a_1,\dots,a_r)=\rho(a_1x_1+\dots+a_rx_r)$.
One should be able to generalize the above construction to most sensible notions of a multi-spectral radius.
Other examples:
In order for our notion of a multi-spectral radius to be sensible, one would expect that the functions $\rho_{\iota,\mathcal{C}}$ would be coherent and interesting for the simplest possible cases of $\iota,\mathcal{C}$. For example, if $\iota:A\rightarrow A$ is the identity function and $1\leq p\leq\infty$, and $\mathcal{C}$ is the unit ball in $\mathbb{C}^r$ with respect to the $p$-norm, then one should expect for $\rho_{\iota,\mathcal{C}}$ to be about as reasonable of a function as the $L_p$-spectral radii, and experimental computations indicate that this is indeed the case.
Define a mapping $F_{\iota,\mathcal{C},a_1,\dots,a_r}:\mathcal{C}\rightarrow[0,\infty)$ by $F_{\iota,\mathcal{C},a_1,\dots,a_r}(x_1,\dots,x_r)=
\rho(x_1\iota(a_1)+\dots+x_r\iota(a_r))$. Experimental computations suggest that the local maxima $(x_1,\dots,x_r)$ for the function $F_{\iota,\mathcal{C},a_1,\dots,a_r}$ tend to resemble a sort of conjugate of $a_1,\dots,a_r$.
Let $\iota_n:M_n(\mathbb{C})\rightarrow M_n(\mathbb{C})$ be the identity mapping, and let $\mathcal{L}_{r;p}=\{(\lambda_1,\dots,\lambda_r)\in \mathbb{C}^r:\|(\lambda_1,\dots,\lambda_r)\|_p=1\}$.
In some of my experiments with $A_1,\dots,A_r\in M_n(\mathbb{R})$ and in all of my experiments with $A_1,\dots,A_r\in M_n(\mathbb{C})$ that are Hermitian or real symmetric, when I computed $\lambda_1,\dots,\lambda_r$ locally maximizes $F_{\iota_n,S_1^r,A_1,\dots,A_r}$, then one can find a $\lambda\in S_1$ and $e_1,\dots,e_r\in\{-1,1\}$ where
$\lambda_j=\lambda\cdot e_j$ for $1\leq j\leq r$. A similar phenomenon holds when I locally maximized $F_{\iota_n,\mathcal{L}_{r;p},A_1,\dots,A_r}$ for $1\leq p\leq\infty$ even though this phenomenon seems to break down as $p$ gets close to $1$ and it holds better for Hermitian matrices than it does for non-symmetric real matrices.
My computer experiments indicate that if we locally maximize $F_{\iota,\mathcal{C},a_1,\dots,a_r}$, then as $\mathcal{C}$ better approximates $A$, the local maxima $(x_1,\dots,x_r)$ will become more and more similar to a conjugate version of $(a_1,\dots,a_r)$. On the other hand, if $\mathcal{C}$ is too complicated and has too much room to work with, then the local maxima $(x_1,\dots,x_r)$ will again poorly represent the elements in $A$. Therefore, in order to best represent the conjugates of the elements in $A$, it is best if $\mathcal{C}$ is a little bit simpler than $A$.
Let $\iota_{r;n,d}:M_n(\mathbb{C})\rightarrow M_{n\times d}(\mathbb{C})$ be the algebra homomorphism defined by $\iota_{r;n,d}(A)=A\otimes I_d$. Let $\mathcal{C}_{r;n,d}$ be the collection of all tuples $(I_n\otimes X_1,\dots,I_n\otimes X_r)$ where $\rho_2(X_1,\dots,X_r)=1$.
Theorem: Suppose that $A_1,\dots,A_r,B_1,\dots,B_r$ are $n\times n$-complex matrices where $A_1,\dots,A_r$ do not have a common invariant subspace. Suppose furthermore that $\rho_2(A_1,\dots,A_r)>0,\rho_2(B_1,\dots,B_r)>0$. Then
$\rho(A_1\otimes B_1+\dots+A_r\otimes B_r)=\rho_2(A_1,\dots,A_r)\rho_2(B_1,\dots,B_r)$ if and only if there is some $\lambda$ and invertible $C$ where $B_j=\overline{\lambda\cdot C\cdot A_j\cdot C^{-1}}$ for $1\leq j\leq r$.
See this answer or this link for proofs that I gave of the above result.
From the above result, we see that if $I_n\otimes\overline{X_1},\dots,I_n\otimes\overline{X_r}\in M_{n\times n}(\mathbb{C})$ globally maximizes
$F_{\iota_{r;n,n},\mathcal{C}_{r;n,n},A_1,\dots,A_r}$ and $A_1,\dots,A_r$ have no common invariant subspace, then there are $C,\lambda$ where $X_j=\lambda CA_jC^{-1}$ for $1\leq j\leq r.$
If $A_1,\dots,A_r\in M_n(\mathbb{C})$ does not have a common invariant subspace and $I_n\otimes\overline{X_1},\dots,I_n\otimes\overline{X_r}\in M_{n\times d}(\mathbb{C})$ locally maximizes
$F_{\iota_{r;n,d},\mathcal{C}_{r;n,d},A_1,\dots,A_r}$, then the matrices $X_1,\dots,X_r$ will (up-to-similarity and a constant factor) resemble $A_1,\dots,A_r$. For example, if $A_1,\dots,A_r$ are all real, complex symmetric, real symmetric, Hermitian, real positive semidefinite, complex positive semidefinite, quaternionic, rank $\leq k$, etc, and $I_n\otimes\overline{X_1},\dots,I_n\otimes\overline{X_r}$ locally maximizes $F_{\iota_{r;n,d},\mathcal{C}_{r;n,d},A_1,\dots,A_r}$, then one will often be able to find a constant $\lambda$ and invertible matrix $C$ where
$Y_j=\lambda CX_rC^{-1}$ satisfy those properties respectively. Furthermore, one will often be able to find matrices $R,S$ where
$Y_j=RA_jS$ for $1\leq j\leq r$. In this case, $RS=I_d$ and $P=SR$ will be a (non-orthogonal) projection matrix. Define linear operators $F,G:M_n(\mathbb{C})\rightarrow M_n(\mathbb{C})$ by setting $F(X)=\sum_{k=1}^rA_kX(PA_kP)^*$ and $G(X)=\sum_{k=1}^rA_k^*XPA_kP$ (here $F=G^*$). Define $U_0=I_n,V_0=I_n$ and set $U_{n+1}=F(U_n)/\|F(U_n)\|,V_{n+1}=G(V_n)/\|G(V_n)\|$ for $n\geq 0$. Then $U_n,V_n$ experimentally converge to positive semidefinite matrices $U,V$, the dominant eigenvectors of $F$ and $G$. It seems like the strategy that the optimization algorithm chose for locally maximizing the spectral radius was to make the dominant eigenvalues of $F,G$ positive semidefinite matrices of rank $d$, but the best way to retain the positive semidefiniteness of the dominant eigenvectors of $F,G$ is to make the operators $PA_kP$ closely related to the operators $A_k$. It seems like the reason this strategy works is that in order for a spectral radius of a matrix $A$ to be large, the matrix $A$ should be designed to maximize a particular eigenvalue, and by making the operators $PA_kP$ related to $A_k$, we can maximize the spectral radius of $F,G$. Since the local maximum values of $F_{\iota_{r;n,d},\mathcal{C}_{r;n,d},A_1,\dots,A_r}$ are closely related to the tuples $(A_1,\dots,A_r)$ themselves, I would regard the multi-spectral radius $\rho_{\iota_{r;n,d},\mathcal{C}_{r;n,d}}$ as a legitimate generalization of the notion of the spectral radius to multiple operators which I call the $L_{2,d}$-spectral radius $\rho_{2,d}$.
Other multi-spectral radii $\rho_{\iota,\mathcal{C}}$ are probably reasonably well-behaved, but more computer experiments are needed to verify whether other multi-spectral radii $\rho_{\iota,\mathcal{C}}$ behave nearly as well as $\rho_{\iota_{r;n,d},\mathcal{C}_{r;n,d}}$.
One can find more details on $\rho_{\iota_{r;n,d},\mathcal{C}_{r;n,d}}$ at my site here, and here is another page where I apply $\rho_{\iota_{r;n,d},\mathcal{C}_{r;n,d}}$ to evaluate cryptographic algorithms. I also gave some experimental observations of $\rho_{\iota_{r;n,d},\mathcal{C}_{r;n,d}}$ right here.
Multi-spectrum:
There seems to be a somewhat reasonable definition of a multi-spectrum of a collection of operators.
Suppose that $\rho_{\iota,\mathcal{C}}$ is a multi-spectral radius. If
$(x_1,\dots,x_r)\in\mathcal{C}$ and $\rho(x_1\iota(a_1)+\dots+x_r\iota(a_r))=\rho_{\iota,\mathcal{C}}(a_1,\dots,a_r)$, then we say that the spectrum of $x_1\iota(a_1)+\dots+x_r\iota(a_r)$ is a multi-spectrum of $a_1,\dots,a_r$ with respect to the embedding $\iota$ and set $\mathcal{C}$.
Suppose that $(x_1,\dots,x_r)\in\mathcal{C}$ and for every neighborhood $U$ of $(x_1,\dots,x_r)$ with respect to the topology induced by the norm on $A$, then whenever $(y_1,\dots,y_r)\in U\cap\mathcal{C}$, we have
$\rho(x_1\iota(a_1)+\dots+x_r\iota(a_r))\geq\rho(y_1\iota(a_1)+\dots+y_r\iota(a_r))$; then we say that the spectrum of $x_1\iota(a_1)+\dots+x_r\iota(a_r)$ is a local multi-spectrum of $(a_1,\dots,a_r)$ with respect to the embedding $\iota$ and the set $\mathcal{C}$.
This notion of a multi-spectrum depends on the choice of $\iota,\mathcal{C}$ and not only on the multi-spectral radius $\rho_{\iota,\mathcal{C}}$. For example, if $A_1,\dots,A_r$ are complex matrices with no common invariant subspace, then the multi-spectrum of $A_1,\dots,A_r$ with respect to $\iota_{r;n,n}$ and $\mathcal{C}_{r;n,n}$ is simply the spectrum of $A_1\otimes\overline{A_1}+\dots+A_r\otimes\overline{A_r}$. On the other hand, suppose that $x_1,\dots,x_r$ are the non-commuting variables in $\overline{B}$ where $\rho_p(A_1,\dots,A_r)=\rho(x_1\iota(A_1)+\dots+x_r\iota(A_r))$ for all matrices $A_1,\dots,A_r$ and suppose that $\mathcal{C}=\{(\lambda_1 x_1,\dots,\lambda_r x_r)\mid \lambda_1,\dots,\lambda_r\in S_1\}$. Then a multi-spectrum of $(A_1,\dots,A_r)$ with respect to $\iota$ and $\mathcal{C}$ is the spectrum of $x_1A_1+\dots+x_rA_r$. If $\lambda$ is complex number with $|\lambda|=1$, then there is an automorphism $\phi$ of the Banach algebra $B$ with $\phi(x_1A_1+\dots+x_rA_r)=\lambda(x_1A_1+\dots+x_rA_r)$. Therefore, since $x_1A_1+\dots+x_rA_r$ has the same spectrum as $\lambda(x_1A_1+\dots+x_rA_r)$, there is some compact set $C\subseteq[0,\infty)$ with $\sigma(x_1A_1+\dots+x_rA_r)=\{\lambda t:|\lambda|=1,t\in C\}$.
Unresolved properties
If the set $\mathcal{C}$ is compact, then the function $\rho_{\iota,\mathcal{C}}$ is automatically upper-semicontinuous. If $\rho_{\iota,\mathcal{C}}$ is not upper-semicontinuous, then we can just take the upper-semicontinuous regularization of $\rho_{\iota,\mathcal{C}}$, but the possible lack of upper-semicontinuity is a potential problem with the abstract theory that I am proposing. I do not know if we should require $\mathcal{C}$ to always be compact in order to make $\rho_{\iota,\mathcal{C}}$ always upper-semicontinuous.
So far, I have mainly experimental results about multi-spectral radii, but I would like for there to be more theorems about multi-spectral radii.
Best Answer
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})$.