One more useful fact that is easy to prove: If $(c,u,\mathrm{d}u)$ is a $D$-parallel section of $E$ on a connected domain $U\subset M$ that is not identically zero, then the zero locus of $u$ is either empty; a nonempty smooth hypersurface, in which case $\tilde\lambda < 0$; or consists of isolated points that are non-degenerate maxima or minima, in which case, $\tilde\lambda = 0$.

With all this understood, let's look at what might reasonably described as a 'conformal connected sum' of compact Einstein manifolds in dimension $n\ge 3$. I'll take it to be this: Three compact Einstein $n$-manifolds $(M_i,g_i)$ for $i=0, 1, 2$, with respective Einstein constants $\lambda_i$, together with smoothly embbeded compact submanifolds $N_i\subset M_i$ where $N_0$ is diffeomorphic to $[0,1]\times S^{n-1}$, and $N_i$ is diffeomorphic to the $n$-ball $B^n$ for $i=1,2$, and, finally, where $M_0{\setminus}N_0$ is conformally diffeomorphic (as Riemannian manfolds) to the disjoint union of $M_1{\setminus}N_1$ and $M_2{\setminus}N_2$. (Without lost of generality, one can take this diffeomorphism to be the identity map, so assume this.)

With all this understood, let's look at what might reasonably described as a 'conformal connected sum' of compact Einstein manifolds in dimension $n\ge 3$. I'll take it to be this: Three compact Einstein $n$-manifolds $(M_i,g_i)$ for $i=0, 1, 2$, with respective Einstein constants $\lambda_i$, together with smoothly embbeded compact submanifolds $N_i\subset M_i$ where $N_0$ is diffeomorphic to $[0,1]\times S^{n-1}$, and $N_i$ is diffeomorphic to the $n$-ball $B^n$ for $i=1,2$, and, finally, where $M_0{\setminus}N_0$ is conformally diffeomorphic (as Riemannian manfolds) to the disjoint union of $M_1{\setminus}N_1$ and $M_2{\setminus}N_2$. (Without lost of generality, one can take this diffeomorphism to be the identity map, so assume this.)

If there are any nontrivial examples with $n\ge3$ beyond what I mentioned in my comment above, I think they are likely to be very few, and there are none that are topologically non-trivial, in the sense that $M_1\# M_2$ will always be homeomorphic to either $M_1$ or $M_2$ (or both).

By the well-known formula for the Ricci tensor of $\tilde g = u^{-2}\,g$, where $u>0$ is a function on $M$, $$ \mathrm{Ric}(\tilde g) = \mathrm{Ric}(g) + (n{-}2)\,u^{-1}\,\nabla(\mathrm{d}u) - \bigl(u^{-1}\,\Delta u + (n{-}1) u^{-2}\,|\mathrm{d}u|^2\bigr)\, g\,, $$ it follows that $\tilde g$ is also Einstein if and only if $u$ is a non-vanishing function on $M$ that satisfies the equation $\nabla(\mathrm{d}u) = v\,g$ for some function $v$. Differentiating this equation, one finds that $v+\lambda u$ must be constant. Thus, by the connectedness of $M$, there exists a constant $c$ such that $\nabla(\mathrm{d}u) = (c{-}\lambda u)\,g$. Consequently, $\Delta u = n(\lambda u{-}c)$, which will be useful below. Moreover, one sees that the Einstein constant of $\tilde g$ is $\tilde\lambda = 2c\,u - \lambda\,u^2 - |\mathrm{d}u|^2$.

Evidently, for a 'generic' Einstein metric, the space of $D$-parallel sections of $E$ consists only of the sections of the form $(c,u,\alpha)=(\lambda u,\,u,\,0)$ where $u$ is a constant. (Such sections are obviously $D$-parallel.)

If $(M^n,g)$ is compact, connected and Einstein with a non-positive Einstein constant, then the only global $D$-parallel sections of $E$ are those with $u$ constant. This is because the defining equation implies $\Delta u = n(\lambda u {-} c)$. When $\lambda=0$, this clearly implies $c=0$ and that $u$ be constant. When $\lambda<0$, this implies that $\lambda u {-} c$ is an eigenvalue of $\Delta$ with negative eigenvalue $n\lambda$ and hence must be $0$.

Meanwhile, if one supposes that $\lambda>0$, the equation $\nabla(\mathrm{d}u) = (c{-}\lambda u)\,g$ implies that, if $p$ is a critical point of $u$, then the Hessian of $u$ at $p$ is $\bigl(c{-}\lambda u(p)\bigr)$ times the metric at $p$. If $c{-}\lambda u(p)=0$, then the $D$-parallel section $(c,u,\mathrm{d}u)$ agrees with the $D$-parallel section $(\lambda\,u(p), u(p),0)$ at $p$ and hence must equal it, i.e., $u$ must be constant. On the other hand, if $c{-}\lambda u(p)\not=0$, then $p$ is a non-degenerate critical point of $u$ that is either a maximum or a minimum. Thus, the only critical points of $u$ on the connected, compact manifold $M^n$ are either strict maxima or minima. It follows from the Mountain Pass Lemma that $u$ can have only one maximum and one minimum. It is a well-known result that this implies that $M^n$ be homeomorphic to the $n$-sphere.

Consequently, on a compact, connected Einstein $n$-manifold $(M^n,g)$ that is not homeomorphic to an $n$-sphere, the only conformal multiples of $g$ that are also Einstein are the constant multiples.

As another example, any connected Einstein $(M^4,g)$ that is not conformally flat also has the property that any $D$-parallel section must be of the form $(c,u,\alpha)=(\lambda u,\,u,\,0)$ for some constant $u$ since the condition $W(\nabla u,X,Y,Z)=0$ would force the Weyl curvature to vanish on the open set where $\nabla u\not=0$, and the real-analyticity of the Einstein metric would then force the Weyl tensor to vanish identically.

At the other extreme, on an $n$-manifold of constant sectional curvature (necessarily conformally flat), the connection $D$ is flat. However, as we have already seen, if the sectional curvature is non-positive and the manifold is compact, then the only global solutions $u$ must be constant. Even when the sectional curvature is positive (and the manifold is compact) the only case in which there are non-constant global solutions is the $n$-sphere. This is because every such space form is a quotient of the round $n$-sphere by a freely acting discrete isometry group $\Gamma$, but such a quotient cannot be homeomorphic to the $n$-sphere when $\Gamma$ is non-trivial.

For any connected Einstein $(M^n,g)$, the sheaf of $D$-parallel sections of $E$ is well-behaved: There is an integer $k$ satisfying $1\le k\le n{+}2$ such that, for any 1-connected open set $U\subset M$, the dimension of the vector space $\mathcal{P}(U)\subset\Gamma(U,E)$ consisting of $D$-parallel sections of $E$ over $U$ is $k$. This follows from the well-known result of DeTurck and Kazdan that $g$, being Einstein, is real-analytic in harmonic coordinates: Since the linear connection $D$ is perforce real-analytic, any germ of a $D$-parallel section defined on an open neighborhood of $p\in M$ can be real-analytically continued as a $D$-parallel section along any real-analytic curve that starts at $p$.

Consequently, $u_i$ for $i=1$ or $2$ must be constant unless $M_i$ is homeomorphic to an $n$-sphere.

If both $u_1$ and $u_2$ are constant, then by scaling $g_1$ and $g_2$ by constants, we can assume that $g_0 = g_i$ on $M_i{\setminus}N_i\subset M_0{\setminus}N_0$ for $i=1$ and $2$, and hence $\lambda_0=\lambda_1=\lambda_2$. Moreover, because $N_i$ is simply-connected and $g_0$ and $g_i$ are real-analytic, it is not hard to show that the isometric inclusion $M_i{\setminus}N_i\subset M_0{\setminus}N_0$ extends to an isometry $\iota_i:(M_i,g_i)\to (M_0,g_0)$, and the compactness of $M_i$ and $M_0$ then implies that $\iota_i$ is a diffeomorphism, mapping $N_i$ diffeomorphically onto $M_0{\setminus}M_i = N_0\cup (M_{3-i}{\setminus}N_{3-i})$. Thus, $M_0$, $M_1$, and $M_2$ are all isometric and homeomorphic to an $n$-sphere. While $(M_0,g_0)$ is indeed a 'conformal connected sum' of $(M_1,g_1)$ and $(M_2,g_2)$ in the above sense, this has to be regarded as a trivial case.

If, say, $u_1$ is constant, we can still reduce to the case that $g_0=g_1$ on $M_1{\setminus}N_1\subset M_0{\setminus}N_0$, and hence $\lambda_0=\lambda_1$. Again, because $N_1$ is simply-connected and $g_0$ and $g_1$ are real-analytic, it follows that the isometric inclusion $M_1{\setminus}N_1\subset M_0{\setminus}N_0$ extends to an isometry $\iota_1:(M_1,g_1)\to (M_0,g_0)$, and the compactness of $M_1$ and $M_0$ then implies that $\iota_1$ is a diffeomorphism, mapping $N_1$ diffeomorphically onto $M_0{\setminus}M_1 = N_0\cup (M_2{\setminus}N_2)$. It follows easily from this that $M_2$ must be homeomorphic to an $n$-sphere. Of course, examples such as this do exist with all of the $(M_i,g_i)$ being conformally flat, and $u_2$ non-constant, in which one attaches a conformally flat bubble to a compact space form, but $(M_0,g_0)$ is isometric (up to a constant multiple) to $(M_1,g_1)$. I don't know of any example of this kind that is not conformally flat, though.

Finally, if neither $u_1$ nor $u_2$ is constant, then, $M_1$ and $M_2$ are both homeomorphic to $n$-spheres, so $M_0$ is as well. Moreover, since $g_i = {u_i}^2 g_0$ for $i=1$ and $2$, and since $M_0$ is simply connected and compact, it follows that $1/u_i$ extends smoothly over the whole $M_0$. From this, it is not hard to argue that $u_i$ cannot vanish on $M_i$ nor can $1/u_i$ vanish on $M_0$. One can then show, using arguments similar to those above, that the compact Einstein manifolds $(M_i,g_i)$ and $(M_0, {u_i}^2 g_0)$ are isometric. Again, this is a case in which the summands of the 'connected sum' are simply global Einstein rescalings of the single compact Einstein manifold $(M_0, g_0)$, which happens to have nontrivial global solutions to the equation $\nabla_0(\mathrm{d}u) = (c-\lambda_0 u) g_0$ on $M_0$ (and hence is homeomorphic to an $n$-sphere). Conversely, such an $(M_0,g_0)$ can always be written as a conformal connect sum of two Einstein manifolds with non-constant 'transition' functions $u_i$.

I do not know whether any examples of this kind exist other than the constant sectional curvature examples on the standard $n$-sphere, but I do not, at the moment, see how to rule them out. A structure equation analysis does not appear to rule out the existence of local Einstein metrics that are not conformally flat but still admit nonconstant Einstein rescalings. Given what was said above, the lowest dimension $n\ge 3$ where this might be possible is $n=5$.

There are no nontrivial examples with $n\ge3$ beyond what I mentioned in my comment above, namely, either a conformal connected sum of a compact space form $(M_1,g_1)$ with the standard round $n$-sphere with the connected sum $M_1\# M_2$ being homothetic to $(M_1,g_1)$ or the case where both $M_1$ and $M_2$ are diffeomorphic to a (possibly exotic) $n$-sphere with homothetic metrics $g_1$ and $g_2$ and the resulting connected sum being again homothetic to $(M_1,g_1)$. (Note that there are many inequivalent Einstein metrics on manifolds homeomorphic to odd-dimensional spheres.)

By the well-known formula for the Ricci tensor of $\tilde g = u^{-2}\,g$, where $u>0$ is a function on $M$, $$ \mathrm{Ric}(\tilde g) = \mathrm{Ric}(g) + (n{-}2)\,u^{-1}\,\nabla(\mathrm{d}u) - \bigl(u^{-1}\,\Delta u + (n{-}1)\, u^{-2}\,|\mathrm{d}u|^2\bigr)\, g\,, $$ it follows that $\tilde g$ is also Einstein if and only if $u$ is a non-vanishing function on $M$ that satisfies the equation $\nabla(\mathrm{d}u) = v\,g$ for some function $v$. Differentiating this equation, one finds that $v+\lambda u$ must be locally constant. Thus, by the connectedness of $M$, there exists a constant $c$ such that $\nabla(\mathrm{d}u) = (c{-}\lambda u)\,g$. Consequently, $\Delta u = n(\lambda u{-}c)$, which will be useful below. Moreover, one sees that the Einstein constant of $\tilde g$ is $\tilde\lambda = 2c\,u - \lambda\,u^2 - |\mathrm{d}u|^2$.

Now, it turns out that, except for the metrics of constant sectional curvature on the $n$-sphere, for a compact, connected Einstein manifold $(M,g)$, the space of global $D$-parallel sections of $E$ consists only of the sections of the form $(\lambda u,\,u,\,0)$ where $u$ is constant. (Such sections are obviously $D$-parallel.)

To see this, first note that if $(M^n,g)$ is compact, connected and Einstein with a non-positive Einstein constant, then the only global $D$-parallel sections of $E$ are those with $u$ constant. This is because the defining equation implies $\Delta u = n(\lambda u {-} c)$. When $\lambda=0$, this clearly implies $c=0$ and that $u$ be constant. When $\lambda<0$, this implies that $\lambda u {-} c$ is an eigenvalue of $\Delta$ with negative eigenvalue $n\lambda$ and hence must be $0$.

Meanwhile, if one supposes that $\lambda>0$, the equation $\nabla(\mathrm{d}u) = (c{-}\lambda u)\,g$ implies that, if $p$ is a critical point of $u$, then the Hessian of $u$ at $p$ is $\bigl(c{-}\lambda u(p)\bigr)$ times the metric at $p$. If $c{-}\lambda u(p)=0$, then the $D$-parallel section $(c,u,\mathrm{d}u)$ agrees with the $D$-parallel section $(\lambda\,u(p), u(p),0)$ at $p$ and hence must equal it, i.e., $u$ must be constant. On the other hand, if $c{-}\lambda u(p)\not=0$, then $p$ is a non-degenerate critical point of $u$ that is either a maximum or a minimum. Thus, on a connected, compact Einstein manifold $(M^n,g)$ the only critical points of a nonconstant $u$ such that $(c,u,\mathrm{d}u)$ is $D$-parallel for some $c$ are either strict maxima or minima. By the Mountain Pass Lemma, it follows that $u$ can only have one maximum and one minimum, which implies that $M$ is homeomorphic to a $n$-sphere.

Suppose that such a non-constant $u$ exists on $(M^n,g)$. By adding a constant to $u$, we can suppose that $u$ is strictly positive on $M$ and by scaling, we can assume that $\lambda=1$. Now by using the structure equations, one can prove that, on a punctured neighborhood of a maximum $p\in M$, $g$ can be written in the 'sine-cone' form $g = \mathrm{d}r^2 + (\sin r)^2\, h$, where $r$ is the distance from $p$, $h$ is an Einstein metric on $S^{n-1}$ with Einstein constant $1$, and $u = a + b \cos r$ for some constants $a$ and $b$. However, unless $h$ is conformally flat, the above sine-cone metric cannot be smooth at $r=0$. Thus, $h$ is conformally flat, which implies that $g$ is also conformally flat. Hence, the only possibility for $(M,g)$ is the round $n$-sphere, as was to be shown.

Now, for any connected Einstein $(M^n,g)$, the sheaf of $D$-parallel sections of $E$ is well-behaved: There is an integer $k$ satisfying $1\le k\le n{+}2$ such that, for any 1-connected open set $U\subset M$, the dimension of the vector space $\mathcal{P}(U)\subset\Gamma(U,E)$ consisting of $D$-parallel sections of $E$ over $U$ is $k$. This follows from the well-known result of DeTurck and Kazdan that $g$, being Einstein, is real-analytic in harmonic coordinates: Since the linear connection $D$ is perforce real-analytic, any germ of a $D$-parallel section defined on an open neighborhood of $p\in M$ can be real-analytically continued as a $D$-parallel section along any real-analytic curve that starts at $p$.

Consequently, for $i=1$ or $2$, $u_i$ must be constant unless $(M_i,g_i)$ is the round $n$-sphere.

If both $u_1$ and $u_2$ are constant, then by scaling $g_1$ and $g_2$ by constants, we can assume that $g_0 = g_i$ on $M_i{\setminus}N_i\subset M_0{\setminus}N_0$ for $i=1$ and $2$, and hence $\lambda_0=\lambda_1=\lambda_2$. Moreover, because $N_i$ is simply-connected and $g_0$ and $g_i$ are real-analytic, it is not hard to show that the isometric inclusion $M_i{\setminus}N_i\subset M_0{\setminus}N_0$ extends to an isometry $\iota_i:(M_i,g_i)\to (M_0,g_0)$, and the compactness of $M_i$ and $M_0$ then implies that $\iota_i$ is a diffeomorphism, mapping $N_i$ diffeomorphically onto $M_0{\setminus}M_i = N_0\cup (M_{3-i}{\setminus}N_{3-i})$. Thus, $M_0$, $M_1$, and $M_2$ are all isometric and homeomorphic to an $n$-sphere. While $(M_0,g_0)$ is indeed a 'conformal connected sum' of $(M_1,g_1)$ and $(M_2,g_2)$ in the above sense, these ought to be regarded as trivial cases. As remarked at the beginning, though, there are many, many such examples that are not conformally flat.

If, say, $u_1$ is constant and $u_2$ is not, we can still reduce to the case that $g_0=g_1$ on $M_1{\setminus}N_1\subset M_0{\setminus}N_0$, and hence $\lambda_0=\lambda_1$. Again, because $N_1$ is simply-connected and $g_0$ and $g_1$ are real-analytic, it follows that the isometric inclusion $M_1{\setminus}N_1\subset M_0{\setminus}N_0$ extends to an isometry $\iota_1:(M_1,g_1)\to (M_0,g_0)$, and the compactness of $M_1$ and $M_0$ then implies that $\iota_1$ is a diffeomorphism, mapping $N_1$ diffeomorphically onto $M_0{\setminus}M_1 = N_0\cup (M_2{\setminus}N_2)$. Meanwhile, $(M_2,g_2)$ must be isometric to the round $n$-sphere, which is conformally flat, implying that $(M_0,g_0)$ is conformally flat. Of course, examples such as this do exist with all of the $(M_i,g_i)$ being conformally flat, and $u_2$ non-constant, in which one attaches a conformally flat bubble to a compact space form, but $(M_0,g_0)$ is isometric (up to a constant multiple) to $(M_1,g_1)$.

Finally, if neither $u_1$ nor $u_2$ is constant, then, $(M_1,g_1)$ and $(M_2,g_2)$ are both round $n$-spheres, so $M_0$ is as well.