If you are okay working with the Hitchin map defined on the moduli stack of Higgs bundles rather than the coarse moduli space, then you can find a proof of this statement in the paper The global nilpotent variety is Lagrangian by V. Ginzburg, at least in the case that the genus of the base curve is at least $2$. The precise reference is Corollary 9.
The idea is to show that the global nilpotent cone (the most singular fiber of the Hitchin map) has the same dimension as $Bun_G$. This implies that the Hitchin map is equidimensional and that the stack of Higgs bundles $T^*Bun_G$ is lci. Since the Hitchin base is non-singular, this implies flatness by miracle flatness (Stacks Project Lemma 00R4). Now the stack of semi-stable Higgs bundles is an open substack so the restriction of the Hitchin map to this locus is also flat.
Edit: I think the result also holds for the coarse moduli space of semi-stable Higgs bundles but it seems to be a bit subtle. I think equidimensionality of the coarse Hitchin map follows from that on the stack. Then by miracle flatness the map is flat if and only if the coarse space is Cohen-Macaulay. Locally the coarse space of semi-stable Higgs bundles looks like a quotient $V//H$ where $H$ is reductive and $V$ is affine chart for the stack (in the smooth topology). Thus the question becomes when is $V//H$ Cohen-Macaulay?
For $V$ non-singular this is the Hochster-Roberts theorem. However, it can fail in general even when $V$ itself is Cohen-Macaulay (in fact even when $V$ is a complete intersection). See for example the last paragraph of example $I$ here. In this case we are saved by the fact that moduli space of semi-stable Higgs bundles has symplectic singularities which are in particular Cohen-Macaulay. See for example this paper.
It seems to me then that being in the symplectic setting is used not just for the dimension bounds but also to ensure that the moduli space is Cohen-Macaulay so I'm not sure what to expect for Higgs bundles valued in an arbitrary line bundle $L$.