Let's simplify a little. Define
$$(U,V) = \frac{1}{\sqrt{\sigma_X^2+\sigma_Y^2}}\left(X+Y,\ \frac{\sigma_Y}{\sigma_X}X - \frac{\sigma_X}{\sigma_Y}Y\right).$$
You can readily check that $U$ and $V$ are uncorrelated standard Normal variables (whence they are independent). In terms of them,
$$X = \frac{\sigma_X}{\sqrt{\sigma_X^2 + \sigma_Y^2}} \left(\sigma_X U + \sigma_Y V\right) = \alpha U + \beta V$$
defines the coefficients of $X$ in terms of $(U,V).$ The question desires a formula for
$$E\left[X^2 \mid |X+Y|\ge k\right] = E\left[\left(\alpha U + \beta V\right)^2 \mid |U| \ge \lambda\right]$$
with $\lambda = k\sqrt{\sigma_X^2 + \sigma_Y^2} \ge 0.$
Expanding the square, we find
$$\begin{aligned}
E\left[\left(\alpha U + \beta V\right)^2 \mid |U| \ge \lambda\right] &= \alpha^2E\left[U^2 \mid |U| \ge \lambda\right] \\&+ 2\alpha\beta E\left[UV \mid |U| \ge \lambda\right] \\&+ \beta^2 E\left[V^2 \mid |U| \ge \lambda\right].
\end{aligned}$$
The second term is zero because $E[V]=0$ and $V$ is independent of $U$. The third term is $\beta^2$ because the independence of $V$ and $U$ gives
$$E\left[V^2\mid |U|\ge \lambda\right] = E\left[V^2\right] = 1.$$
This leaves us to compute the first conditional expectation. The standard (elementary) formula expresses it as the fraction
$$E\left[U^2 \mid |U|\ge \lambda\right] = \frac{\left(2\pi\right)^{-1/2}\int_{|u|\ge \lambda} u^2 e^{-u^2/2}\,\mathrm{d}u}{\left(2\pi\right)^{-1/2}\int_{|u|\ge \lambda} e^{-u^2/2}\,\mathrm{d}u}$$
The denominator is $\Pr(|U|\ge \lambda) = 2\Phi(-\lambda)$ where $\Phi$ is the standard Normal distribution function.To compute the numerator, substitute $x = u^2/2$ to obtain
$$\frac{1}{\sqrt{2\pi}}\int_{|u|\ge \lambda}u^2 e^{-u^2/2}\,\mathrm{d}u = \frac{2^{3/2}}{\sqrt{2\pi}}\int_{\lambda^2/2}^\infty x^{3/2\,-1}\ e^{-x}\,\mathrm{d}x = \frac{1}{\Gamma(3/2)}\int_{\lambda^2/2}^\infty x^{3/2\,-1}\ e^{-x}\,\mathrm{d}x.$$
This equals $\Pr(Z\ge \lambda^2/2)$ where $Z$ has a Gamma$(3/2)$ distribution. It is a regularized incomplete gamma function, $P(3/2, \lambda^2/2).$ Consequently, with $\lambda \ge 0,$
$$E\left[\left(\alpha U + \beta V\right)^2 \mid |U| \ge \lambda\right] =\beta^2 + \frac{\alpha^2 P(3/2, \lambda^2/2)}{2 \Phi(-\lambda)}.$$
To illustrate, this R implementation of the conditional expectation (with a representing $\alpha,$ b representing $\beta,$ and $k$ representing $\lambda$) uses pnorm for $\Phi$ and pgamma for the Gamma distribution:
f <- function(a, b, k) {
b^2 + a^2 * pgamma(k^2/2, 3/2, lower.tail = FALSE) / (2 * pnorm(-k))
}
Best Answer
Exploring
Trying out a few values we see that the values get higher for larger $n$.
(we have drawn a red line based on the bound computed below)
Computing/approximating attempt
Based on the above (tendency to high n and small p) we can try and approximate the sum by replacing the binomial distribution with a Poisson distribution with $\lambda = np$
$$f(k) \approx \frac{\lambda^k e^{-\lambda}}{k!}$$
and let Wolfram alpha compute
$$\sum_{k=1}^\infty \frac{\lambda}{k} \frac{\lambda^k e^{-\lambda}}{ k!} = -e^{-\lambda}\lambda(log(\lambda)-Ei(\lambda)+ \gamma)$$
We arrive at
$$E[np/k] \approx e^{-np}np(1-log(np)+Ei(np)- \gamma) $$
(where we added an extra term $np e^{-np}$ to capture the contribution of $k=0$)
Maybe this can be differentiated, but I had it computed and arrive at a bound of approximately
$$\phi \approx 1.386317$$