My favorite reference for these sorts of things that straddle physics and geometry is Frankel's "The geometry of physics".  In the chapter on harmonic forms, you will find what he refers to simply as "Hodge's Theorem".  It's a little more general than you need, because it applies to general $p$-forms, and you only need functions ($0$-forms).  So I'll specialize it to functions.

---

**Hodge's Theorem (for functions):** Let $M^n$ be a closed Riemannian manifold.  Then Poisson's equation
\begin{equation}
  \Delta \alpha = \rho
\end{equation}
(where $\alpha$ and $\rho$ are real-valued functions, and $\Delta$ is the Laplacian) has a solution $\alpha$ if and only if $\rho$ has mean value $0$ on $M^n$:
\begin{equation}
  \int_M \rho\ \mathrm{vol}^n = 0.  \tag{A}
\end{equation}

---

So, to apply this theorem to your question, we have<sup>†</sup>
\begin{align}
  M &= \partial \mathscr{B}, \\
  \alpha &= y, \\
  \rho &= \mathrm{const} - \left( p_{b ; d} \hat{h}^{bd} - R_{ac} Y^{a}_{1} Y^{c}_{2} + R_{adcb} Y^{d}_{1} Y^{c}_{2} Y^{a}_{2} Y^{b}_{1} \right).  \tag{B} \\
\end{align}
To put that another way, there is a solution for $y$ such that
\begin{equation}
  y_{;bd} \hat{h}^{bd} + p_{b ; d} \hat{h}^{bd} - R_{ac} Y^{a}_{1} Y^{c}_{2} + R_{adcb} Y^{d}_{1} Y^{c}_{2} Y^{a}_{2} Y^{b}_{1} = \mathrm{const}  \tag{C}
\end{equation}
if and only if the integral of $\rho$ over $\partial \mathscr{B}$ is zero.  But we get to adjust the value of $\mathrm{const}$, so we can just set it to whatever we need to make that integral zero.

Hawking & Ellis point out that $p_{b ; d} \hat{h}^{bd}$ is a pure divergence.  So you can use [Stokes' theorem](https://en.wikipedia.org/wiki/Stokes%27_theorem) to transform its integral over $\partial \mathscr{B}$ into an integral over the boundary of $\partial \mathscr{B}$.  But the boundary of a boundary is always empty,<sup>††</sup> so that integral has value $0$.  Therefore, this term disappears when you integrate $\rho$.
Now, combining this fact with (A) and (B), the theorem states that $y$ has a solution if and only if
\begin{equation}
  \mathrm{const}\ \int_{\partial \mathscr{B}} \mathrm{vol}^n = \int_{\partial \mathscr{B}} \left( - R_{ac} Y^{a}_{1} Y^{c}_{2} + R_{adcb} Y^{d}_{1} Y^{c}_{2} Y^{a}_{2} Y^{b}_{1} \right)\ \mathrm{vol}^n.  \tag{D}
\end{equation}
The volume (actually, the area of the surface) is assumed to be finite and nonzero, and any volume (area) is non-negative, so the left-hand side is just the constant in question times some positive number.

Now, a math-class way of stating the conclusion would be that given $p_a$, $Y_1^b$, $Y_2^c$, $R_{ijkl}$, and $\hat{h}^{mn}$, one can *choose* a constant [given by rearranging Eq. (D)] such that there exists a function $y$ that solves Eq. (C).  Hawking & Ellis change the emphasis to suit their goals, but the statement is also true: there exists a $y$ such that the first four terms in the original question's Eq. (1) add up to a constant, the sign of which is determined by the integral on the right-hand side of Eq. (D).

---

<sup>†</sup>
Note that I have left out the $p'^a p'_a$ term in (1) from the original question; that term involves derivatives of $y$ other than the Laplacian, so Hodge's theorem doesn't apply to them.  But also note that Hawking & Ellis don't actually claim that it should be included in what is equal to a constant.

<sup>††</sup>
As an interesting side note, the boundary of a boundary is always empty *when dealing with manifolds*.  This is not true of more general topological spaces, because in those settings [the word "boundary" means something different](https://en.wikipedia.org/wiki/Boundary_(topology)#Boundary_of_a_boundary).