diff --git a/exponax/ic/_gaussian_random_field.py b/exponax/ic/_gaussian_random_field.py
index 1af80c4..42d46ee 100644
--- a/exponax/ic/_gaussian_random_field.py
+++ b/exponax/ic/_gaussian_random_field.py
@@ -64,9 +64,12 @@ def __call__(
             self.num_spatial_dims, self.domain_extent, num_points
         )
         wavenumer_norm_grid = jnp.linalg.norm(wavenumber_grid, axis=0, keepdims=True)
+        # Further division by 2.0 in the exponent is because we want to have the
+        # **power-spectrum** follow a **power-law**. See
+        # https://github.com/Ceyron/exponax/issues/9 for more details.
         amplitude = jnp.power(wavenumer_norm_grid, -self.powerlaw_exponent / 2.0)
         amplitude = (
-            amplitude.flatten().at[0].set(0.0).reshape(wavenumer_norm_grid.shape)
+            amplitude.flatten().at[0].set(1.0).reshape(wavenumer_norm_grid.shape)
         )
 
         real_key, imag_key = jr.split(key, 2)