I'll denote by $(S)$ the following system of equations, which is just the system from the statement but with the second equation multiplied by $-1$ for my convenience:
$$(S) :\begin{cases}(f(x) + g(y))^2 = 2 \sqrt{a} \ h(x+y) + (a+b)
\\ (f(x) - g(y))^2 = - 2 \sqrt{a} \ h(x-y) + (a + b)\end{cases}$$
I also won't use qualificators like "$\forall (x,y) \in \mathbb{R}^2, \dots$" at all times but they're implied whenever I use $x$ and $y$ instead of using another name like $y_0$.
We have, by summing the equations of $(S)$:
$$\begin{split}f(x)^2 + g(y)^2 &= \frac{1}{2}\left((f(x) + g(y))^2 + (f(x) - g(y))^2\right)\\
&= \sqrt{a}(h(x+y) - h(x-y)) + (a + b) \end{split}\tag{1}$$
Similarly, by substracting the second line from the first line, it follows that, as observed by OP:
$$2 f(x) g(y) = \sqrt{a} \left( h(x+y) + h(x-y) \right) \tag{2}$$
By substituting $y = 0$ in $(1)$, we gain:
$$f(x)^2 + g(0)^2 = a + b \tag{3}$$
Thus, there are no solutions if $g(0)^2 > a + b$ since $f(x)^2 \geq 0$, hence we can assume that $g(0)^2 \leq a + b$ from now on.
Let $C := a + b - g(0)^2 \geq 0$ such that $f(x)^2 \equiv C$.
Then, we get, after squaring $(2)$ and replacing $f(x)^2$ by $C$:
$$C g(y)^2 = \frac{a}{4}(h(x+y) + h(x-y))^2 \tag{4}$$
By respectively setting $x = 0$ or $y = 0$ in $(3)$, it ensues that:
$$\begin{align*} C g(y)^2 &= \frac{a}{4}(h(2y) + h(0))^2 \tag{5} \\ h(x)^2 &= \frac{Cg(0)^2}{a} =: K\tag{6}\end{align*}$$
Therefore, as a result of $(3)$, $(6)$, and then $(5)$ in light of $(6)$, we know that all three functions $f$, $g$ and $h$ take a finite amount of values, with:
$$\begin{cases}f(x)^2 \equiv C \\ h(x)^2 \equiv K \\ g(y)^2 = \begin{cases} g(0)^2 \quad\text{if } h(2y) = h(0)\\ 0 \quad\quad\,\,\,\,\text{if } h(2y) = - h(0) \end{cases}\end{cases}$$
This is just to sum things up mid-reasoning though. This is not sufficient yet as they're only necessary conditions.
Let's start by looking at the critical case where $K = 0$.
First, assume that $C = 0$. This is equivalent to $g(0)^2 = a+b$, and to $f \equiv 0$.
From the expression of $K$ in $(6)$, $K = 0$, and thus $h \equiv 0$.
Inputting that in $(1)$ gives $g(y)^2 \equiv a + b$. Conversely, one can check that $f \equiv 0$, $g^2 \equiv a + b$, $h \equiv 0$ is solution of $(S)$ (in particular, $g$ doesn't have to be constant).
On the other hand, if $K = 0$ but $C \neq 0$, then, since $aK = C g(0)^2$, we have $g(0)^2 = 0$, thus $g \equiv 0$ due to our past observations. Conversely, $f^2 \equiv a + b$, $g \equiv 0$ and $h \equiv 0$ is solution of $(S)$ (in particular, $f$ didn't have to be constant either), and we are done with the critical case.
Let's assume moving forward that $K \neq 0$, which is equivalent to having $C \neq 0$ and $g(0) \neq 0$.
Take a fixed $y_0$. Going back to $(1)$ with the context of $(3)$ and $(4)$, and using that $C \neq 0$ by assumption, we can see that, if $h(x + y_0) = - h(x-y_0) \in \left\{\pm\sqrt{K}\right\}$, which is equivalent to $h(2y_0) = - h(0) \in \left\{ \pm\sqrt{K}\right\}$ thanks to $(4)$, then:
$$C + g(y_0)^2 = a + b \pm 2\sqrt{aK}$$
hence, since $a + b - C = g(0)^2$:
$$g(y_0)^2 = g(0)^2 \pm 2\sqrt{aK}$$
But we already have $g(y_0)^2 \in \left\{g(0)^2, 0\right\}$, thus this cannot happen unless $g(y_0)^2 = g(0)^2 - 2\sqrt{aK} = 0$, which would impose specifically that $$h(x + y_0) = - h(x - y_0) = h(x' + y_0) = - h(x' - y_0) = -\sqrt{K}$$ for all $(x, x') \in \mathbb{R}^2$.
Yet, if we consider $z \in \mathbb{R}$ and look at $(x, x') = (z - y_0, z + y_0)$, this would imply that: $$h(z) = h(x + y_0) = - h(x' - y_0) = - h(z)$$ absurd since $K \neq 0$ by assumption, hence this subcase cannot happen at all, and we always have:
$$h(x + y) = h(x - y)$$
which consequently provides, because of $(4)$:
$$g(y)^2 = \frac{a}{4C} \cdot \left(2\sqrt{K}\right)^2 = \frac{aK}{C} = g(0)^2$$
and also lets us deduce that:
$$h(x) = h\left(\frac{x}{2} + \frac{x}{2}\right) = h\left(\frac{x}{2} - \frac{x}{2}\right) = h(0)$$
We can now end this by returning to equation $(2)$ which has no squares and nothing is equal to $0$ to see that, by setting $y = 0$ or $x = 0$ respectively we obtain:
$$f(x) = \frac{\sqrt{a}h(0)}{g(0)} \\ g(y) = g(0)$$
where the second line comes from the fact that $f(0) = \frac{\sqrt{a}h(0)}{g(0)}$, hence $g(y) = \frac{\sqrt{a}h(0)}{\frac{\sqrt{a}h(0)}{g(0)}} = g(0)$.
Conversely, $g \equiv \alpha \in (-\sqrt{a+b}, \sqrt{a+b}) \setminus \{0\}, h \equiv \beta := \pm\sqrt{\frac{(a + b - \alpha^2)\alpha^2}{a}}, f \equiv \frac{\sqrt{a}\beta}{\alpha}$ is solution of $(S)$, and we are finished.
Considering this was quite long to write, I would not be surprised if there were mistakes in there, but hopefully this answers your question correctly.
