This is not an answer, just some sketchy thoughts that are too long for the comment box. I and some other HHS enthusiasts are very interested in this question being answered; we've tried a fair bit and have set it aside, so I don't think they'll mind me trying to recall what some of the strategies and issues are.

It's indeed open for cocompact special groups, as far as I know. (Part of the motivation for the question comes from Coxeter groups, which are known to be hyperbolic relative to thick subgroups [Theorem VII here], and many of which are virtually compact special.)

There's a very similar question: for the classes of groups in question (HHG, cocompactly cubulated, virtually compact special), is it true that any group $G$ in the class is hyperbolic relative to a collection of subgroups, each of which has at most polynomial divergence (we're allowing $G$ to be hyperbolic relative to itself)?

In particular, does super-polynomial divergence imply exponential divergence for such $G$? (A nontrivially relatively hyperbolic group has at least exponential divergence [Theorem 6.13 here].)

It's also interesting to ask whether $G$ is hyperbolic relative to NRH subgroups (whether because the peripheral subgroups are thick, or NRH for some other reason).

Here is a naive approach that has been tried a couple of times (once for cubical groups, once for HHG) unsuccessfully. I think it's worth sketching because it's probably the first thing one might try, so it might save someone some work to see what the problems with it are. (Maybe they are surmountable.)

The idea is to build candidate thick/NRH subspaces inductively by "brute force", show that this construction terminates, and then show that the resulting subspaces give a valid peripheral structure.

More precisely: let $G$ act geometrically on a proper CAT(0) cube complex $X$, or a proper hierarchically hyperbolic space $X$. In either case, either $X$ is already hyperbolic, or there is a nontrivial product region $P$. (In the cubical case, $P$ is a convex subcomplex decomposing as the product of two unbounded CAT(0) cube complexes. In the HHS case, $P$ is a hierarchically quasiconvex subspace admitting a coarse-median-preserving quasi-isometry to the product of two unbounded hierarchically hyperbolic spaces. The HHS fact follows from Corollary 2.16 here.)

Let $\mathcal P_0$ be the ($G$--invariant) set of such product regions. Call $P,P'\in\mathcal P_0$ "elementary equivalent" if they have unbounded coarse intersection. Taking the transitive closure gives an equivalence relation on $\mathcal P_0$. For each equivalence class, one can take the union $U$ of the subspaces in the class. These $U$ are your candidate thick-of-order-at-most-$1$ pieces.

Now iterate this process: take the set of such $U$, impose the same equivalence relation, and make candidate thick-of-order-at-most-$2$ pieces, etc.

(If I remember right, the equivalence relation used at the $n^{th}$ stage, $n\geq 1$, is that $P,P'\in\mathcal P_{n-1}$ are elementary equivalent if their *hierarchically quasiconvex hulls* have unbounded coarse intersection. These are described nicely by Russell-Spriano-Tran in Section 5 of this paper; in the cubical case, you use cubical convex hulls.)

Suppose the process terminates, i.e. at some point you have a $G$--invariant collection $\mathcal P_n$ of subspaces, no two of which have hierarchically quasiconvex hulls with unbounded coarse intersection.

At this point, you should:

(1) Verify that each $P$ is "quasiconvex" in the sense appropriate to the category you're working in. In the cubical case, you should check that $P$ is at finite Hausdorff distance from its cubical convex hull; in the HHS case, you should check that $P$ is \emph{hierarchically quasiconvex} --- Proposition 5.11 in this paper by Russell-Spriano-Tran is probably the best way to go about this.

(2) Hope that the stabiliser of each $P\in\mathcal P_n$ acts coboundedly. ($G$ acts on $X$ coboundedly, $\mathcal P_n$ is $G$--invariant; you want to use properness of $X$ and boundedness of the coarse intersections between the elements of $\mathcal P_n$ to show that only finitely many elements of $\mathcal P_n$ intersect any given ball.)

Already there is something tricky here.

(3) Verify that adding combinatorial horoballs over the elements of $\mathcal P_n$ in $X$ gives you something hyperbolic.

I think one approach was to verify hyperbolicity by finding a hierarchically hyperbolic structure with no interesting "product regions".

(4) Hope that, more or less by construction, each $P\in\mathcal P_n$ is thick of order at most $n$.

At that point, you'd have what you want.

Difficulties include:

(a) It's not remotely clear that this process terminates. (But I think it's worth nailing down whether or not it does).

(b) If it doesn't terminate, one can still hope to be able to pass to some sort of limit, and find a collection $\mathcal P_\infty$ of subspaces, all with bounded coarse intersection, to get candidate peripherals in a relatively hyperbolic structure. If this works, one doesn't get that the peripherals are thick, but the hope is that they are at least not relatively hyperbolic.

(c) If I remember correctly, one big difficulty is that, at each (finite) stage of the induction, one might have to pass to larger and larger neighbourhoods to get something convex/hierarchically quasiconvex, and this creates problems if the process doesn't terminate.

More vaguely, at each step, there are various constants to control, and they blow up if the process doesn't terminate.

There was a quite different approach (for CAT(0) groups) discussed by many people at the AIM in 2016. There was some subsequent work on it and some useful ideas, but I'll have to check with those involved to see if it's okay to mention the idea/subsequent developments here --- it's possible that one or more of them are still actively working on it, although I'm not.

So, the short answer, is that it's a very interesting question for which the naive approach is a bit of a mess.

My feeling is that if there's a counterexample, then there's a counterexample where $G$ is $\pi_1$ of a CAT(0) square complex.