Other questions seem to be trying to answer the simpler question of why a
Gaussian mode has a Gouy phase (i.e. compared to plane-waves), rather than the
actual question being asked as to why higher order modes have a faster
accumulation of Gouy phase (i.e. why phase for the $n^{th}$ order mode has a
phase $\Delta\phi=n\psi(z)$ higher than the $0^{th}$ order Gouy phase of
$\psi(z)$), which I will try to answer.
For me, an intuitive understanding comes from understanding what we mean by
"modes." A mode is a configuration of the field that doesn't change upon
propagation, or rather only changes trivially via e.g. a phase or a rescaling.
So upon propagation a "mode" is a field that transforms like
$$U(x) \to U(x')e^{-\phi}.$$
Modes are highly dependent on the physics of wave propagation. The first thing
to note is that propagation from $z=0$ to the far-field $z\to\infty$ is given
by a Fourier transform. Therefore the propagation from $z=-\infty$ to
$z=\infty$ is given by two Fourier Transforms. However, two Fourier transforms
just flip the field, i.e. if $U(x)$ is the field at some transverse point $x$,
then Fourier transforming twice just gives $\mathcal F(\mathcal F(U(x))) =
\mathcal P U(x)= U(-x)$, where $\mathcal F$ is the Fourier operator and
$\mathcal P$ is the parity operator. A ray-picture of propagation from $z=
-\infty\to\infty$ is given in the image below. This can also be thought of by
the picture that at infinity, all sources look like spherical waves.

Now applying the parity operator twice gives the same mode $\mathcal P\mathcal
P U(x) = \mathcal PU(-x) = U(x)$. There are only two types of modes that obey
the parity operator, even and odd functions, i.e.
$$ \mathcal PU_\text{even}(x) = U_\text{even}(x),$$
or
$$\mathcal PU_\text{odd}(x) = -U_\text{odd}(x).$$
This is equivalent to saying that the only phases that a mode can accumulate
are
$$\mathcal P U_n(x) = e^{i\pi n}U(x),$$
where $n$ is an integer representing even or odd modes when $n$ is likewise an
even or odd number.
So physically we see that a "mode," $U(x)$ is either an even or odd function
which simply flips $U(x)\to U(-x)$ when going from $z=-\infty\to\infty$, which
we can represent as a phase $e^{i\pi n}$. The second piece of physics we add is
the fact that propagation is continuous. Therefore, a "mode" acquires a
half the phase from $z=-\infty\to 0$ and half the phase from $z=0\to\infty$.
Therefore, all modes aquire a phase difference from $z=0\to\infty$ of $e^{i\pi
n/2}$, for some integer $n$.
Mathematically, the Fourier transform $\mathcal F$ is simply the square root of
the parity operator $\mathcal P$ (and the inverse Fourier transform can be
thought of as the negative root). If you were to try to mathematically solve
for eigenmodes of $\mathcal F$, you would find modes with exactly the $e^{i\pi
n/2}$ eigenvalues from above. You can actually take this one step further and
consider the Fractional Fourier transform, which is just the continuous
version of $\mathcal F$ where $\pi/2$ is replaced by a continuous parameter
$\alpha$. The Gaussian modes of propagation are also "modes" of this continuous
operator (for the obvious reason that propagation is continuous). For paraxial
propagation $\alpha = \text{arctan}(z/z_R)$, which gives exactly the Guoy phase
as a function of $z$, i.e. $ e^{i\alpha n} =
\text{exp}[in\text{arctan}(z/z_R)]=\psi(z).$
Therefore the intuitive answer for what is the Guoy phase and why is it
different for different modes is the fact that propagation causes fields to
flip (through the origin) and "modes" are defined in terms of fields where this
geometric flipping can be represented by a phase, and different phases are what
actually defines different modes.