I'll answer #1 only, leaving the rest to you or someone else to figure out. The answer is affirmative.

Dropping the trivial case $m=1$, we may assume by density that we are in the generic position, i.e., that $A_{ij}^T\ne 0$ and together have no non-trivial invariant subspace (this is just to avoid considering degenerate cases, which can also be done using dimension reduction but I prefer not to bother with this extra technicality).

Consider the mapping $X\mapsto F(X)=\sum_{i,j}A_{ij}XA_{ij}^T$ on the cone $K$ of symmetric non-negative definite $d\times d$ matrices $X$. Notice that if $X\ne 0$ and $F(X)=0$, then all $A_{ij}^T$ map the entire space to the kernel of $X$, which contradicts our generic position assumption. Thus the mapping $X\mapsto \frac 1{\operatorname{Tr}F(X)}F(X)$ is well-defined on the convex compact $Q=K\cup\{X:\operatorname{Tr}X=1\}$ and maps $Q$ to itself. So, by the Brouwer fixed point theorem, the linear mapping $F$ has an eigenvector in the cone with positive eigenvalue, which we can normalize to $1$ by scaling $A$, i.e., we have $$ \sum_{i,j}A_{ij}XA_{ij}^T=X\,, $$ which, clearly, implies $\rho(\sum_{i,j}A_{ij}\otimes A_{ij})\ge 1$.

I now claim that $X$ is non-degenerate. Indeed, otherwise the kernel of $X$ would be non-trivial and invariant under all $A_{ij}^T$.

Now for a vector $y=(y_1,\dots,y_m)\in\mathbb R^{md}$ define its norm by $$ \|y\|^2=\sum_j|X^{-1/2}y_j|^2 $$ where $|\cdot|$ stands for the usual Euclidean norm in $\mathbb R^d$. Let $y$ be a unit vector (in this norm) on which the induced operator norm of $A$ is attained.

Then $\|Ay\|^2=\sum_{i} \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2$ for some unit vectors $\xi_i\in\mathbb R^d$. Now use Cauchy-Schwarz: $$ \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2= \left[\sum_j \langle X^{-1/2}A_{ij}y_j,\xi_i\rangle\right]^2 \\ = \left[\sum_j \langle X^{-1/2}y_j,X^{1/2}A_{ij}^TX^{-1/2}\xi_i\rangle\right]^2 \\ \le \left[\sum_j |X^{-1/2}y_j|^2\right]\left[\sum_j |X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2\right] $$ This allows one to estimate the induced operator norm of $A$ and, thereby $\rho(A)$ by $$ \rho(A)^2\le\|A\|^2\le \sum_{i,j}|X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2 \\ =\sum_{i,j}\langle X^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}\xi_i,\xi_i\rangle\,. $$ Now let $L$ be the linear subspace of $\mathbb R^d$ spanned by all vectors $\xi_i$. Obviously, $\operatorname{dim}L\le\min(d,m)$. Let $P$ be the orthogonal projection to $L$. Then for every non-negative definite operator $S$ in $\mathbb R^d$, we have $\langle S\xi_i,\xi_i\rangle\le\operatorname{Tr} (PSP)$, which yields the bound $$ \sum_{i,j} \operatorname{Tr}(PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P) = \operatorname{Tr}\left(\sum_{i,j} PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P\right) \\ = \operatorname{Tr}\left(PX^{-1/2}\left[\sum_{i,j} A_{ij}XA_{ij}^T\right]X^{-1/2}P\right)= \operatorname{Tr}(PX^{-1/2}XX^{-1/2}P) \\ =\operatorname{Tr} P=\operatorname{dim}L\le\min(d,m) $$ as required.

I'll answer #1 only, leaving the rest to you or someone else to figure out. The answer is affirmative.

Dropping the trivial case $m=1$, we may assume by density that we are in the generic position, i.e., that $A_{ij}^T\ne 0$ and together have no non-trivial invariant subspace (this is just to avoid considering degenerate cases, which can also be done using dimension reduction but I prefer not to bother with this extra technicality).

Consider the mapping $X\mapsto F(X)=\sum_{i,j}A_{ij}XA_{ij}^T$ on the cone $K$ of symmetric non-negative definite $d\times d$ matrices $X$. Notice that if $X\ne 0$ and $F(X)=0$, then all $A_{ij}^T$ map the entire space to the kernel of $X$, which contradicts our generic position assumption. Thus the mapping $X\mapsto \frac 1{\operatorname{Tr}F(X)}F(X)$ is well-defined on the convex compact $Q=K\cap\{X:\operatorname{Tr}X=1\}$ and maps $Q$ to itself. So, by the Brouwer fixed point theorem, the linear mapping $F$ has an eigenvector in the cone with positive eigenvalue, which we can normalize to $1$ by scaling $A$, i.e., we have $$ \sum_{i,j}A_{ij}XA_{ij}^T=X\,, $$ which, clearly, implies $\rho(\sum_{i,j}A_{ij}\otimes A_{ij})\ge 1$.

I now claim that $X$ is non-degenerate. Indeed, otherwise the kernel of $X$ would be non-trivial and invariant under all $A_{ij}^T$.

Now for a vector $y=(y_1,\dots,y_m)\in\mathbb R^{md}$ define its norm by $$ \|y\|^2=\sum_j|X^{-1/2}y_j|^2 $$ where $|\cdot|$ stands for the usual Euclidean norm in $\mathbb R^d$. Let $y$ be a unit vector (in this norm) on which the induced operator norm of $A$ is attained.

Then $\|Ay\|^2=\sum_{i} \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2$ for some unit vectors $\xi_i\in\mathbb R^d$. Now use Cauchy-Schwarz: $$ \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2= \left[\sum_j \langle X^{-1/2}A_{ij}y_j,\xi_i\rangle\right]^2 \\ = \left[\sum_j \langle X^{-1/2}y_j,X^{1/2}A_{ij}^TX^{-1/2}\xi_i\rangle\right]^2 \\ \le \left[\sum_j |X^{-1/2}y_j|^2\right]\left[\sum_j |X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2\right] $$ This allows one to estimate the induced operator norm of $A$ and, thereby $\rho(A)$ by $$ \rho(A)^2\le\|A\|^2\le \sum_{i,j}|X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2 \\ =\sum_{i,j}\langle X^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}\xi_i,\xi_i\rangle\,. $$ Now let $L$ be the linear subspace of $\mathbb R^d$ spanned by all vectors $\xi_i$. Obviously, $\operatorname{dim}L\le\min(d,m)$. Let $P$ be the orthogonal projection to $L$. Then for every non-negative definite operator $S$ in $\mathbb R^d$, we have $\langle S\xi_i,\xi_i\rangle\le\operatorname{Tr} (PSP)$, which yields the bound $$ \sum_{i,j} \operatorname{Tr}(PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P) = \operatorname{Tr}\left(\sum_{i,j} PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P\right) \\ = \operatorname{Tr}\left(PX^{-1/2}\left[\sum_{i,j} A_{ij}XA_{ij}^T\right]X^{-1/2}P\right)= \operatorname{Tr}(PX^{-1/2}XX^{-1/2}P) \\ =\operatorname{Tr} P=\operatorname{dim}L\le\min(d,m) $$ as required.

I'll answer #1 only, leaving the rest to you or someone else to figure out. The answer is affirmative.

Dropping the trivial case $m=1$, we may assume by density that we are in the generic position, i.e., that $A_{ij}^T\ne 0$ and together have no non-trivial invariant subspace (this is just to avoid considering degenerate cases, which can also be done using dimension reduction but I prefer not to bother with this extra technicality).

Consider the mapping $X\mapsto F(X)=\sum_{i,j}A_{ij}XA_{ij}^T$ on the cone $K$ of symmetric non-negative definite $d\times d$ matrices $X$. Notice that if $X\ne 0$ and $F(X)=0$, then all $A_{ij}^T$ map the entire space to the kernel of $X$, which contradicts our generic position assumption. Thus the mapping $X\mapsto \frac 1{\operatorname{Tr}F(X)}F(X)$ is well-defined on the convex compact $Q=K\cup\{X:\operatorname{Tr}X=1\}$ and maps $Q$ to itself. So, by the Brouwer fixed point theorem, the linear mapping $F$ has an eigenvector in the cone with positive eigenvalue, which we can normalize to $1$ by scaling $A$, i.e., we have $$ \sum_{i,j}A_{ij}XA_{ij}^T=X\,, $$ which, clearly, implies $\rho(\sum_{i,j}A_{ij}\otimes A_{ij})\ge 1$.

I now claim that $X$ is non-degenerate. Indeed, otherwise the kernel of $X$ would be non-trivial and invariant under all $A_{ij}^T$.

Now for a vector $y=(y_1,\dots,y_m)\in\mathbb R^{md}$ define its norm by $$ \|y\|^2=\sum_j|X^{-1/2}y_j|^2 $$ where $|\cdot|$ stands for the usual Euclidean norm in $\mathbb R^d$.

Then $\|Ay\|^2=\sum_{i} \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2$ for some unit vectors $\xi_i\in\mathbb R^d$. Now use Cauchy-Schwarz: $$ \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2= \left[\sum_j \langle X^{-1/2}A_{ij}y_j,\xi_i\rangle\right]^2 \\ = \left[\sum_j \langle X^{-1/2}y_j,X^{1/2}A_{ij}^TX^{-1/2}\xi_i\rangle\right]^2 \\ \le \left[\sum_j |X^{-1/2}y_j|^2\right]\left[\sum_j |X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2\right] $$ This allows one to estimate thee induced operator norm of $A$ and, thereby $\rho(A)$ by $$ \rho(A)^2\le\|A\|^2\le \sum_{i,j}|X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2 \\ =\sum_{i,j}\langle X^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}\xi_i,\xi_i\rangle\,. $$ Now let $L$ be the linear subspace of $\mathbb R^d$ spanned by all vectors $\xi_i$. Obviously, $\operatorname{dim}L\le\min(d,m)$. Let $P$ be the orthogonal projection to $L$. Then for every non-negative definite operator $S$ in $\mathbb R^d$, we have $\langle S\xi_i,\xi_i\rangle\le\operatorname{Tr} (PSP)$, which yields the bound $$ \sum_{i,j} \operatorname{Tr}(PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P) = \operatorname{Tr}\left(\sum_{i,j} PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P\right) \\ = \operatorname{Tr}\left(PX^{-1/2}\left[\sum_{i,j} A_{ij}XA_{ij}^T\right]X^{-1/2}P\right)= \operatorname{Tr}(PX^{-1/2}XX^{-1/2}P) \\ =\operatorname{Tr} P=\operatorname{dim}L\le\min(d,m) $$ as required.

I'll answer #1 only, leaving the rest to you or someone else to figure out. The answer is affirmative.

Dropping the trivial case $m=1$, we may assume by density that we are in the generic position, i.e., that $A_{ij}^T\ne 0$ and together have no non-trivial invariant subspace (this is just to avoid considering degenerate cases, which can also be done using dimension reduction but I prefer not to bother with this extra technicality).

Consider the mapping $X\mapsto F(X)=\sum_{i,j}A_{ij}XA_{ij}^T$ on the cone $K$ of symmetric non-negative definite $d\times d$ matrices $X$. Notice that if $X\ne 0$ and $F(X)=0$, then all $A_{ij}^T$ map the entire space to the kernel of $X$, which contradicts our generic position assumption. Thus the mapping $X\mapsto \frac 1{\operatorname{Tr}F(X)}F(X)$ is well-defined on the convex compact $Q=K\cup\{X:\operatorname{Tr}X=1\}$ and maps $Q$ to itself. So, by the Brouwer fixed point theorem, the linear mapping $F$ has an eigenvector in the cone with positive eigenvalue, which we can normalize to $1$ by scaling $A$, i.e., we have $$ \sum_{i,j}A_{ij}XA_{ij}^T=X\,, $$ which, clearly, implies $\rho(\sum_{i,j}A_{ij}\otimes A_{ij})\ge 1$.

I now claim that $X$ is non-degenerate. Indeed, otherwise the kernel of $X$ would be non-trivial and invariant under all $A_{ij}^T$.

Now for a vector $y=(y_1,\dots,y_m)\in\mathbb R^{md}$ define its norm by $$ \|y\|^2=\sum_j|X^{-1/2}y_j|^2 $$ where $|\cdot|$ stands for the usual Euclidean norm in $\mathbb R^d$. Let $y$ be a unit vector (in this norm) on which the induced operator norm of $A$ is attained.

Then $\|Ay\|^2=\sum_{i} \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2$ for some unit vectors $\xi_i\in\mathbb R^d$. Now use Cauchy-Schwarz: $$ \langle X^{-1/2}(Ay)_i,\xi_i\rangle^2= \left[\sum_j \langle X^{-1/2}A_{ij}y_j,\xi_i\rangle\right]^2 \\ = \left[\sum_j \langle X^{-1/2}y_j,X^{1/2}A_{ij}^TX^{-1/2}\xi_i\rangle\right]^2 \\ \le \left[\sum_j |X^{-1/2}y_j|^2\right]\left[\sum_j |X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2\right] $$ This allows one to estimate the induced operator norm of $A$ and, thereby $\rho(A)$ by $$ \rho(A)^2\le\|A\|^2\le \sum_{i,j}|X^{1/2}A_{ij}^TX^{-1/2}\xi_i|^2 \\ =\sum_{i,j}\langle X^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}\xi_i,\xi_i\rangle\,. $$ Now let $L$ be the linear subspace of $\mathbb R^d$ spanned by all vectors $\xi_i$. Obviously, $\operatorname{dim}L\le\min(d,m)$. Let $P$ be the orthogonal projection to $L$. Then for every non-negative definite operator $S$ in $\mathbb R^d$, we have $\langle S\xi_i,\xi_i\rangle\le\operatorname{Tr} (PSP)$, which yields the bound $$ \sum_{i,j} \operatorname{Tr}(PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P) = \operatorname{Tr}\left(\sum_{i,j} PX^{-1/2}A_{ij}XA_{ij}^TX^{-1/2}P\right) \\ = \operatorname{Tr}\left(PX^{-1/2}\left[\sum_{i,j} A_{ij}XA_{ij}^T\right]X^{-1/2}P\right)= \operatorname{Tr}(PX^{-1/2}XX^{-1/2}P) \\ =\operatorname{Tr} P=\operatorname{dim}L\le\min(d,m) $$ as required.