The usual definition of a covariance operator on $L_2(D)$ is:
$$ C : L_2(D) \to L_2(D), \qquad (C \psi)(x) = \int_D c(x,y) \psi(y) dy \qquad \forall x\in D, ~~\psi \in L_2(D), $$ where $c(x,y): D \times D \to \mathbb{R}$ is the corresponding covariance function. The operator norm is then $$ \|C\|= \sup_{\|\psi\|_{L_2}=1} \sqrt{\int_D \left(\int_D c(x,y)\psi(y) dy\right)^2dx} $$ I am interested in the covariance operator corresponding to a multivariate covariance function. For example in this paper the author describes the multivariate covariance function $K: D \times D \to \mathbb{R}^p \times \mathbb{R}^p$, so now we have a matrix valued covariance function $K(x,y) = [c_{ij} (x,y)]$ where $c_{ij}$ is a univariate covariance function defined as above.
My question is about the operator corresponding to $K$, and the definition of the operator norm. Is it correct to write this as:
$$ C_K: (L_2(D))^p \to (L_2(D))^p, \qquad (C_K \psi)(x) =\int_{D} (K(x,y))^T \psi(y) dy, \qquad \forall x \in D, ~~ \psi \in (L_2(D))^p $$ and $$ \|C_K\|= \sup_{\| \psi \|_{(L_2(D))^p}=1} \sqrt{\int_{D} \left ( \int_{D} (K(x,y))^T\psi(y) dy \right )^2dx} $$
I am also wondering if there are any good references that discuss this setting in more detail.
I think you are just trying to write out the general operator norm in this setting. (I would say that the choice of $p$ for the dimension is unfortunate, since it too much reminds $L_p(D)$. I will call it by $n$ from now on.)
In your univariate case, your concrete formula is just $$ \|C\| = \sup_{\|\psi\|_{L_2(D)} = 1} \|C\psi\|_{L_2(D)}. $$ Therefore in the multivariate case, the norm should be $$ \|C_K\| = \sup_{\|\psi\|_{L_2(D)^n} = 1} \|C_K\psi\|_{L_2(D)^n}. $$
So it all boils down to how you define your norm on $L_2(D)^n$. The norm on $L_2(D)^n$ comes from $L_2(D)$, just like a norm on ${\mathbb R}^n$ comes from the absolute value norm $|\cdot|$ on ${\mathbb R}$. You can do $L^1, L^\infty$ or $L^2$, and they are all equivalent. So you can define the norm of $\psi = (\psi_1,\dots, \psi_n)\in L_2(D)^n$ as any of the following \begin{gather*} \sum_{i=1}^n \|\psi_i\|_{L_2(D)} = \sum_{i=1}^n \sqrt{\int_D \psi_i^2(x)\,dx},\\ \sqrt{\sum_{i=1}^n \|\psi_i\|^2_{L_2(D)}} = \sqrt{\sum_{i=1}^n \int_D \psi_i^2(x)\,dx},\\ \max_{1\leq i\leq n} \|\psi_i\|_{L_2(D)} = \max_{1\leq i\leq n} \sqrt{\int_D \psi_i^2(x)\,dx}. \end{gather*}
So pick any of these, you have a version of $\|C_K\|$ defined. (Again, they are all equivalent norms.) Say for the $L^2$, i.e. Euclidean norm, we have concretely \begin{align*} \|C_K\| &= \sup_{\|\psi\|_{L_2(D)^n} = 1} \|C_K\psi\|_{L_2(D)^n}\\ &= \sup_{\sum_{i=1}^n \int_D \psi_i^2(x)\,dx = 1} \sqrt{ \sum_{i=1}^n \int_D \Big(\int_D \sum_{j=1}^n c_{ij}(x, y)\psi_j(y) \,dy\Big)^2 \,dx}. \end{align*}
I think the relevant reading is product norm, but more importantly what type of function $c(x, y)$ would define such a map (not discussed in this question).