Self-dual and anti-self dual metrics (and hence locally conformally flat) are also Bach flat. So, for example, the product metric on $S^1 \times S^3$ is Bach flat, even though there is no Einstein metric on $S^1 \times S^3$.

They are equivalent by stereographic projection when you choose a point on $\partial S_+^n$ to go to infinity.

The first case (with the irreducible assumption) is a manifold diffeomorphic to a sphere equipped with a metric whose curvature operator is nonnegative. These do not have to be locally symmetric (i.e. $\nabla R$ need not vanish).

Your conclusion should be diffeomorphic, not isometric, because of Igor's comment (small deformations of the standard metric on the sphere give metric with positive curvature operator which are not locally symmetric).

Correct; the realization of the Bach tensor as the gradient only works in dimension four. As far as I can tell, the generalization of the Bach tensor to other dimensions is motivated by the Fefferman--Graham ambient metric.

I think it is better to call this (a variant of) the Bochner formula. The formula in the original post will (essentially) appear in any textbook under this name. "Reilly formula" usually refers to the result after integrating on a manifold with boundary and performing further algebraic manipulations on the boundary.