I can't comment yet so I try a short answer:
You could impose a parameter restriction and set $a_0=0$ to gain efficiency.
E.g. in R using the rugarch package you might use the option
fixed.pars = list(omega=0)
in the specification.
The model is still meaningful. The long term volatility in this case would be zero, and as long as $a_1$ (and $a_2$) is significant the model exhibits a mean reversion to the long term volatility of zero. Of course this long-term volatility would never be reached in practice since mean reversion is "slow" and new information will come in and push you away from the long-term vola.
There are two different parametrizations of the GJR-GARCH model in rugarch, and you're applying the formula for the persistence from one parametrization to the other.
For GJR-GARCH(1,1), the first one is the one you've shown, which in the documentation is written like this:
$$\sigma_t^2 = w + (\alpha + \gamma I_{t-1}) \varepsilon_{t-1}^2 + \beta \sigma_{t-1}^2$$
where $I_t = 1$ when $\varepsilon_t \leq 0$, and 0 otherwise.
As you say, when the standardized residuals are drawn from a symmetric distribution, the persistence is given by:
$$\alpha + \beta + \frac{1}{2}\gamma$$
There is an alternate, equivalent parametrization, which arises as a special case of the family-GARCH model, which goes like this:
$$\sigma_t^2 = w + \alpha' \sigma^2_{t-1} (|z_{t-1}|-\eta z_{t-1})^2 + \beta \sigma_{t-1}^2$$
Since the sign of $z_t$ and $\varepsilon_t$ is the same and since $\varepsilon_t = \sigma_t z_t$, you can show that this is the same as:
$$\sigma_t^2 = w + \left(\alpha'(1-\eta)^2 + \alpha'\left[(-1-\eta)^2-(1-\eta)^2\right]I_{t-1} \right)\varepsilon^2_{t-1} + \beta \sigma_{t-1}^2$$
Reading off the coefficients, we have the relationship between the two parametrizations:
$$\alpha = \alpha'(1-\eta)^2$$
$$\gamma = \alpha'\left[(-1-\eta)^2-(1-\eta)^2\right]$$
Therefore, in the second parametrization the persistence is given by:
$$\alpha'(1-\eta)^2 + \beta + \frac{1}{2}\alpha'\left[(-1-\eta)^2-(1-\eta)^2\right]$$
The more general version of this appears in the documentation you've linked to as equations 27-28.
In rugarch, the first parametrization is obtained from rugarch::ugarchspec by passing :
variance.model = list(model = "gjrGARCH", garchOrder = c(1,1))
The second, which is the one you have, by passing:
variance.model = list(model = "fGARCH", garchOrder = c(1,1), submodel = "GJRGARCH")
Plugging in your coefficient estimates (to the precision that you've displayed them), I get that the persistence of your GJR-GARCH model is $0.9575652$, which is very close to what rugarch::persistence gave you.
Best Answer
Interesting. I would not say the external regressor reduces/removes volatility but it is able to fit/explain it better than the previous period's conditional variance $\sigma_{t-1}^2$ can. This reminds me a bit of the realized GARCH model. There, including lagged realized variance into the conditional variance equation tends to make the lagged squared return term not statistically significant; see slide 36 of Peter R. Hansen's "Lecture 3: Realized GARCH Models" (2016). Meanwhile, in your case the external regressor does the same for lagged conditional variance.