I'm one of the $\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$ value between two (binned) variables as coming from a tilted bivariate normal distribution. The $\chi^2$ needs to pass a certain noise pedestal, else $\phi_k$ is zero.

For low statistics samples $\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$ is affected more than Pearson's correlation constant, because, unlike Pearson, $\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$ simply call:

`df.significance_matrix()`

to get the Z-values. In your example you will see that $\phi_k$ is only evaluated when (roughly) $Z > 0.5$, that the Z scores lie around zero, and also that none of the Z scores are truly significant, say $Z > 5$.

In general, $\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$ will work fine but other correlation constants will be more precise, in particular for very low statistics samples.

Hope this helps!