Perhaps one answer is the appearance of the Gaussian in the
*oscillator* (a.k.a. *Weil*) representation --- which makes sense
even in non-zero characteristic. You're probably
aware of the construction of this infinite dimensional representation
for the double cover of $\mathrm{SL}_2(\Bbb{R})$. The story for $\mathrm{SL}_2 \big( \Bbb{F}_q \big)$
where $q$ a power of a prime strongly parallels the continuous story
except, of course, it's finite:

Choose any non-square $\delta \in \Bbb{F}_q$ and form the quadratic extension
$\Bbb{F}_q(\delta) \cong \Bbb{F}_q^2$. As usual, we identify elements of
$\Bbb{F}_q(\delta)$ as linear combinations $z = a + \sqrt{\delta} b$
with $a, b \in \Bbb{F}_q$ subject the usual formulae for addition
and multiplication; *conjugation* and *norm* are expressed as
$\bar{z} = a - \sqrt{\delta} b$ and $\mathrm{N}(z) = a^2 - \delta b^2$
respectively. We'll need to choose (any) non-trivial additive
character $\psi: \Bbb{F}_q \longrightarrow \Bbb{C}^*$ together with
any multiplicative character
$\chi: \Bbb{F}_q(\delta)^* \longrightarrow \Bbb{C}^*$ which generates the
character group of $\Bbb{F}_q(\delta)^*$. Morally $\psi$ plays the
role of the exponential function $\exp: \Bbb{R} \longrightarrow \Bbb{R}^*$
in this finite context. Now set

\begin{equation}
W_\chi := \ \left\{
\begin{array}{l}
\displaystyle \text{all functions} \
f: \Bbb{F}_q(\delta) \longrightarrow \Bbb{C} \ \ \text{such that} \\
\displaystyle
f(wz) = \overline{\chi(w)} \, f(z) \ \text{whenever $\mathrm{N}(w)=1$}
\end{array} \right\}
\end{equation}

which is $q-1$ dimensional. The *oscillator* representation $\varrho_\chi: \mathrm{SL}_2(\Bbb{F}_q)
\longrightarrow \mathrm{GL}(W_\chi)$ is determined by the action of the
follow elements, which generate $\mathrm{SL}_2(\Bbb{F}_q)$:

\begin{equation}
\begin{array}{ll}
\displaystyle \varrho_\chi { \scriptstyle \begin{pmatrix} x & 0 \\ 0 & x^{-1} \end{pmatrix}}
f(z)
&\displaystyle = \
f(xz) \\
\displaystyle \displaystyle \varrho_\chi { \scriptstyle \begin{pmatrix} 1 & \ \ y \\ 0 & \ \ 1 \end{pmatrix}} f(z)
&\displaystyle = \ \underbrace{\psi \big( y \, \mathrm{N}(z) \big)}_{\text{the Gaussian $G_y(z)$}} \cdot f(z) \\
\displaystyle \displaystyle \varrho_\chi { \scriptstyle \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}} f(z)
&\displaystyle = \ -q \widehat{f}(-\bar{z})
\end{array}
\end{equation}

where $\widehat{f}$ is the Fourier transform of $f$ with respect to the
additive structure of $\Bbb{F}_q(\sqrt{\delta})$. Like the continuous Gaussian,
$G_y(z)$ is an eigenfunction of this Fourier transform.

As far as I know, the construction that I've outlined (and which I learned from Amritanshu Prasad's online notes) can be carried out for
$\mathrm{SL}_2\big( \Bbb{Z}/q\Bbb{Z} \big)$ where $q= p^N$ is still a power of a fixed prime $p$. Furthermore a coherent choice of additive and multiplicative characters
$\psi_N$ and $\chi_N$ can be made for each $N \geq 1$ so that these oscillator representations agree with the inverse system

\begin{equation}
1 \longleftarrow \Bbb{Z}/p\Bbb{Z} \stackrel{\ \text{mod $p$}}{\longleftarrow} \Bbb{Z}/p^2\Bbb{Z} \stackrel{\ \text{mod $p^2$}}{\longleftarrow}
\Bbb{Z}/p^3\Bbb{Z} \longleftarrow \, \cdots
\end{equation}

thus allowing the oscillator representations $\varrho_{\chi_N}$ to be exported to the $p$-adic integers $\Bbb{Z}_p$, which will inherent some kind of Gaussian-like function $\varprojlim G_{\bf y}$ where
${\bf y}= (y_1,y_2,y_3, \dots)$ and $y_N = y_{N+1} \, \text{mod $p^N$}$
for each $N \geq 1$.

3more comments