A non-linear Granger causality test was implemented by Diks and Panchenko (2006). The code can be found here and it is implemented in C. The test work as follows:
Suppose we want to infer about the causality between two variables $X$ and $Y$ using $q$ and $p$ lags of those variables, respectively. Consider the vectors $X_t^q = (X_{t-q+1}, \cdots, X_t)$ and $Y_t^p = (Y_{t-p+1}, \cdots, Y_t)$, with $q, p \geq 1$. The null hypothesis that $X_t^q$ does not contain any additional information about $Y_{t+1}$ is expressed by
$$
H_0 = Y_{t+1}|(X_t^q;Y_t^p) \sim Y_{t+1}|Y_t^p
$$
This null hypothesis is a statement about the invariant distribution of the vector of random variables $W_t = (X_t^q, Y_t^p, Z_t)$, where $Z_t=Y_{t+1}$. If we drop the time indexes, the joint probability density function $f_{X,Y,Z}(x,y,z)$ and its marginals must satisfy the following relationship:
$$
\frac{f_{X,Y,Z}(x,y,z)}{f_Y(y)} = \frac{f_{X,Y}(x,y)}{f_Y(y)} \cdot \frac{f_{Y,Z}(y,z)}{f_Y(y)}
$$
for each vector $(x,y,z)$ in the support of $(X,Y,Z)$. Diks and Panchenko (2006) show that, for a proper choice of weight function, $g(x,y,z)=f_Y^2(y)$, this is equivalent to
\begin{align}
q = E[f_{X,Y,Z}(X,Y,Z)f_Y(Y) - f_{X,Y}(X,Y)f_{Y,Z}(Y,Z)].
\end{align}
They proposed the following estimator for $q$:
\begin{align}
T_n(\varepsilon) = \frac{(n-1)}{n(n-2)} \sum_i (\hat{f}_{X,Y,Z}(X_i,Y_i,Z_i) \hat{f}_Y(Y_i) - \hat{f}_{X,Y}(X_i,Y_i) \hat{f}_{Y,Z}(Y_i,Z_i))
\end{align}
where $n$ is the sample size, and $\hat{f}_W$ is a local density estimator of a $d_W$-variate random vector $W$ at $W_i$ based on indicator functions $I_{ij}^W = I(\|W_i - W_j\| < \varepsilon)$, denoted by
\begin{align}
\hat{f}_W(W_i) = \frac{(2 \varepsilon)^{-d_W}}{n-1} \sum_{j,j \neq i} I_{ij}^W.
\end{align}
In the case of bivariate causality, the test is consistent if the bandwidth $\varepsilon$ is given by $\varepsilon_n = Cn^{-\beta}$, for any positive constant $C$ and $\beta \in (\frac{1}{4}, \frac{1}{3})$. The test statistic is asymptotically normally distributed in the absence of dependence between the vectors $W_i$. For the choice of the bandwidth, Diks and Panchenko (2006) suggest $\varepsilon_n = max(C_n^{-2/7},1.5)$, were $C$ can be calculated based on the ARCH coefficient of the series.
There are other tests of non-linear Granger causality such as in Hiemstra and Jones (1994), but this test in particular suffers from lack of power and over-rejection problems, as stated by Diks and Panchenko here.
As pointed out by @RichardHardy, you should be careful about using local density estimation in small samples. Since Diks and Panchenko showed that in samples smaller than 500 observations their test may under-reject, it would be wise to make further investigations in case the test does not reject the null hypothesis.
References
Diks, C., & Panchenko, V. (2006). A new statistic and practical guidelines for nonparametric Granger causality testing. Journal of Economic Dynamics and Control, 30[9–10], 1647-1669.
Hiemstra, C., & Jones, J. D. (1994). Testing for Linear and Nonlinear Granger Causality in the Stock Price- Volume Relation. The Journal of Finance, 49(5), 1639–1664.
