  We study the solution $h_\varepsilon$ of the Kardar-Parisi-Zhang (KPZ)
equation for $d \geq 3$: $$ \frac{\partial}{\partial t} h_{\varepsilon} =
\frac12 \Delta h_{\varepsilon} + \bigg[\frac12 |\nabla h_\varepsilon |^2 -
C_\varepsilon\bigg]+
  \beta \varepsilon^{\frac{d-2}2} \xi_{\varepsilon} $$ with
$h_\varepsilon(0,x)=0$. Here $\xi_\varepsilon=\xi\star \phi_\varepsilon$ is a
spatially smoothened (at scale $\varepsilon$) Gaussian space-time white noise
and $C_\varepsilon$ is a divergent constant as $\varepsilon\to 0$. When the
disorder $\beta$ is sufficiently small and $\varepsilon\to 0$,
$h_\varepsilon(t,x)- \mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)\to 0$ in
probability where $\mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)$ is the {\emph
stationary solution} of the KPZ equation - more precisely, $\mathfrak
h^{\mathrm{st}}_{\varepsilon}$solves the above equation with a random initial
condition (that is independent of the driving noise $\xi$) and its law is
constant in $(\varepsilon,t,x)$. In the present article we quantify the rate of
the above convergence in this regime and show that the fluctuation {\emph
about} the stationary solution $$ (\varepsilon^{1-\frac d2} [h_\varepsilon(t,x)
- \mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)])_{x,t} $$ converges pointwise
(with finite dimensional distributions in space and time) to a Gaussian free
field (GFF) evolved by the deterministic heat equation. We also identify the
fluctuations {\it of} the stationary solution itself and show that the rescaled
averages $\int_{\mathbb R^d} {\mathrm d} x \varphi(x) \varepsilon^{1-\frac d2}
[\mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)- \mathbb E(\mathfrak
h^{\mathrm{st}}_{\varepsilon}(t,x))]$ converge to that of the {\emph stationary
solution} of the stochastic heat equation with additive noise, but with
(random) {\emph GFF marginals} (instead of flat initial condition).
