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})$.
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
I have a feeling that this specific construction is already known but have no time to do an extensive check. There are a number of ways to start a search.
The book Hadamard Matrices and Their Applications by Horadam is a good reference. There is also an online library of known Hadamard matrices maintained by Neil Sloane.
PS: The idea of group developed Hadamard matrices is defined on page 6 of Padraic O'Cathain's master's thesis. The rest of the thesis may be interesting as well.