I'm one of the phi_k$\phi_k$ authors, so happy to help out.
Indeed there is no closed-form formula for phi_k, but it boils down to interpreting the Pearson chi^2$\chi^2$ value between two (binned) variables as coming from a tilted bivariate normal distribution. The chi^2$\chi^2$ needs to pass a certain noise pedestal, else phi_k$\phi_k$ is zero.
For low statistics samples phi_k$\phi_k$ is indeed affected by statistical fluctuations (like any correlation constant). In case of low statistics the spread in the reported values goes up, but the median is unbiased (see publication). phi_k$\phi_k$ is affected more than Pearson's correlation constant, because, unlike Pearson, phi_k$\phi_k$ does not use exact positional information. For numeric variables it bins the data, and then only uses the number of counts per bin, essentially treating them as categorical variables.
So in case of a low statistics sample, be sure to always check the significance of the correlations found (this holds for any correlation constant!). For phi_k$\phi_k$ simply call:
df.significance_matrix()df.significance_matrix()
to get the Z-values. In your example you will see that phi_k$\phi_k$ is only evaluated when (roughly) Z>0.5$Z > 0.5$, that the Z scores lie around zero, and also that none of the Z scores are truly significant, say Z > 5$Z > 5$.
In general, phi_k$\phi_k$ is useful when you have a set of variables where some are categorical or ordinal. If you only have numeric variables, phi_k$\phi_k$ will work fine but other correlation constants will be more precise, in particular for very low statistics samples.
Hope this helps!