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!