Does probability independence have anything to do with functional independence (on general case)? On normal-linear case I know it holds for example if $cor(X_1,Y)=0$ then $g_{X_1Y}(X_1,Y)=g_{X_1}(X_1)g_{Y}(Y)$, where $g_{(.)}$ is the distribution of (.), and that imples to $\beta_1=0$ on:
$Y=\beta_1X_1$
, since $cor$ works as measure to evaluate how good a linear approximation between two variables is and those implications hold we have that there is no linear function such that $Y=f(X_1)$ (or there is with $\beta_1=0$ but it is not informative at all) when both $X_1$ and $Y$ are
stochastically independent and normal distribuited.
But I'm not sure if it holds for general case, that is if $g_{X_1Y}(X_1,Y)=g_{X_1}(X_1)g_{Y}(Y)$ so there is no function $f$ such that $Y=f(X_1)$ (I call this functional independence). Am able to go functional dependence as probability depencence or is it pretty wrong? and if I'm not is there any relation between they both?

For random normal data, the sample mean $\bar X$ and the sample variance $S^2$ are independent random variables. But they are not 'functionally' independent because $\bar X$ is used in the definition of $S^2.$
A formal proof of the independence of $\bar X$ and $S^2$ can be done
using a linear transformation in $n$-space or by moment generating functions. [This is a consequence of Basu's Thm; see example. Also, see this Q&A; "Related" page listed in the margin.]
I show a simulation in R below that illustrates the idea (in terms of correlation), using $m$ normal
samples of size $n.$ With $m = 100\,000$ iterations the correlation
should be accurate to one or two decimal places.
    set.seed(2021)
    n = 10;  m = 10^5
    x = rnorm(m*n, 50, 7)
    DTA = matrix(x, ncol=n)  # each row a sample
    a = apply(DTA, 1, mean)  # m sample means
    v = apply(DTA, 1, var)   # m sample variances
    cor(a, v)
    [1] -0.00202886          # consistent w/ 0 correlation
    
    s = sqrt(v)
    plot(a, s, pch=".")


