The answer to both of your questions is yes, and I think that it works completely analogously to what describe in the first part of your question. More specifically,
The idea of the trace is the same one that you discuss in the
first part of your question. The "$E$-Hessian" $\nabla^2 \varphi$ is a section of $E
\otimes T^\ast M \otimes T^\ast M$. Use the metric (the musical
isomorphism $\sharp$) to identify $T^\ast M$ with $TM$ and obtain a section
of $E \otimes T^\ast M \otimes TM \cong E \otimes \text{End}(TM)$.
Take the trace of the endomorphism piece to obtain a section of $E$.
Note that $\nabla^2 \phi$ may not be symmetric in its entries, so we have a choice
as to which factor of $T^\ast M$ we apply $\sharp$ to, but this choice is irrelevant once
we take the trace.
In terms of a local frame $\{ e_i \}$ for $TM$, we have
$$ \Delta \phi = \text{tr}_g \nabla^2 \varphi = g^{ij} [\nabla^2 \varphi] (e_i, e_j).$$
Of course, if the frame is orthonormal, this reduces to $[\nabla^2 \varphi] (e_i, e_i)$,
as you stated.
We have a connection $\nabla^E$ on $E$ and the Levi-Civita connection $\nabla^{LC}$ on $T^\ast M$. These induce a connection $\tilde{\nabla}^E$ on
$E \otimes T^\ast M$, which is defined by the "product rule":
$$\tilde{\nabla}^E (\varphi \otimes \omega) = (\nabla^E \varphi) \otimes \omega + \varphi \otimes (\nabla^{LC} \omega)$$
where $\varphi$ is a section of $E$ and $\omega $ is a one-form.
This connection gives a map $\tilde{\nabla}^E: \Gamma(E \otimes T^\ast M) \to \Gamma(E \otimes T^\ast M \otimes T^\ast M)$, and as you guessed, $\nabla^2$ as you defined it is precisely the composition of the connections
$$\tilde{\nabla}^E \circ \nabla^E: \Gamma(E) \to \Gamma(E \otimes T^\ast M) \to \Gamma(E \otimes T^\ast M \otimes T^\ast M).$$
It's a good exercise to prove this!
A note, as much for my own understanding as anything: $\tilde{\nabla}^E$ as defined above is an extension of $\nabla^E$ which maps $\Gamma(E \otimes T^\ast M) \to \Gamma(E \otimes T^\ast M \otimes T^\ast M)$. There is another interesting extension of $\nabla^E$ to $E \otimes T^\ast M$, which I will call $d^E$. $d^E$, in contrast to $\tilde{\nabla}^E$, is a map
$$d^E: \Gamma(E \otimes T^\ast M) \to \Gamma(E \otimes \Lambda^2 T^\ast M),$$
i.e., $d^E$ gives a two-form with $E$-coefficients. More generally, one can define $d^E$ as a map on $E$-valued forms of any degree:
$$d^E: \Gamma(E \otimes \Lambda^k T^\ast M) \to \Gamma(E \otimes \Lambda^{k+1} T^\ast M),$$
defined by
$$d^E(\varphi \otimes \omega) = (\nabla^E \varphi) \wedge \omega + \varphi \otimes (d\omega),$$
where $\wedge$ means "wedge the one-form part of $\nabla^E \varphi$ with $\omega$".
I suppose one should think of $d^E$ as a generalization of the de Rham exterior derivative $d$ on forms. (If $E$ is the trivial bundle $M \times \mathbb{R}$ with trivial connection $\nabla^E = d$, we recover $d$.) Note that $d^E$ does not require a connection on $TM$.
The curvature $R^E$ of the connection $\nabla^E$ is the $\text{End}(E)$-valued two-form defined by the composition
$$R^E:=(d^E)^2 \text{ (or }d^E \circ \nabla^E): \Gamma(E) \to \Gamma(E \otimes T^\ast M) \to \Gamma(E \otimes \Lambda^2 T^\ast M) .$$
It's a standard computation to show that $R^E$ really lives in $\text{End}(E)$, i.e., that it's $C^\infty(M)$-linear, and that
$$ R^E(X,Y) = \nabla^E_X \nabla^E_Y - \nabla^E_Y \nabla^E_Y - \nabla^E_{[X,Y]} $$
As a final remark to relate this back to the Hessian, notice that the antisymmetric part of the Hessian is precisely the curvature, i.e.,
$$[\nabla^2 \varphi](X, Y) - [\nabla^2 \varphi](Y, X) = R^E(X,Y) \varphi.$$
The Hessian of a smooth function is symmetric, which is equivalent to the the fact that the de Rham "curvature" $d^2$ is zero.
References: Here are a couple of books I found useful in reminding myself how some of this works:
- Jost, Riemannian Geometry and Geometric Analysis. See chapter 4.
- Taylor, Partial Differential Equations I. See Appendix C (on his website) and chapter 2.
I think the answers to your first two questions are yes, there is always a codifferential; and no, there is not always a Hodge decomposition. I don't know anything about the Einstein manifold case.
Let's say $E$ is a vector bundle over $M$ with metric and compatible connection $\nabla$. ($E$ could be some tensor bundle with $\nabla$ induced from the Levi-Civita connection, for example.) As you say, $\nabla$ gives an exterior derivative $d^\nabla$ on $E$-valued differential forms that obeys the Leibniz rule
$$ d^\nabla (\omega \otimes e) = d\omega \otimes e + (-1)^{\text{deg} \omega} \omega \wedge \nabla e ,$$
where $\omega$ is a homogeneous form and $e$ is a section of $E$, and $\omega \wedge \nabla e$ means $\sum_i (\omega \wedge dx^i) \otimes \nabla_{\partial_i} e$.
For each degree of forms $k$, you can view $d^\nabla$ as a map $d^\nabla: C^\infty( \Lambda^k \otimes E) \to C^\infty( \Lambda^{k+1} \otimes E)$. Using the metrics, each of those spaces of sections are endowed with $L^2$ inner products. Then $d^\nabla$ always has a formal $L^2$-adjoint $\delta^\nabla$ going the other way: $\delta^\nabla: C^\infty( \Lambda^{k+1} \otimes E) \to C^\infty( \Lambda^k \otimes E)$.
(Note the space of smooth sections can be completed to the space of $L^2$ sections, but $d^\nabla$ and $\delta^\nabla$ are unbounded operators and can generally only be defined on some dense subspace of that $L^2$ space.)
So you can always define $d^\nabla$ and its formal adjoint $\delta^\nabla$. The issue is that $d^\nabla$ squares to zero if and only if the connection $\nabla$ on $E$ is flat, meaning its curvature is zero. (Flat vector bundles $E$, i.e., bundles on which there exists a flat connection, are relatively scarce.) We need $d^\nabla$ to square to zero to even attempt to define the de Rham cohomology of $E$-valued forms, and to have nice properties like the Hodge decomposition. (I want to say that the "Dirac operator" $D = d^\nabla + \delta^\nabla$ and the "Laplacian" $\Delta = D^2$ are not elliptic unless $d^\nabla$ squares to zero, but I need to think about that.)
I don't have a good reference for this at hand, although the Wikipedia article on bundle-valued forms might be useful.
Best Answer
The two Laplacians are not the same. One is the other plus a curvature term. The formula is known as the Weitzenböck formula. It is surprisingly difficult to calculate. It and its proof can be found in this note of Petersen.