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.
