On diagonal terms in the Coulomb matrix
In terms of how they got the relation for the diagonal elements, I believe it is relatively straightforward. Given a list of atomic energies (energy of atom relative to separated nucleus and electrons), one can try to fit these values to a function of nuclear charge. If you plot a few values, it suggests trying an exponential fit (I don't know why all the papers say polynomial fit, as this could not lead to a noninteger exponent). So, you can try to fit the values to a function of the form $aZ^b$ using, for example, least-squares regression. You should wind up with something pretty close to the equation given in your question, with some variation because I don't know what atomic energies they used to do the fit. (I'm also not sure why the function isn't negative since the atomic energies are negative. My best guess is the sign winds up not being important as a descriptor in the network). 

As to why they choose to use these values along the diagonal, they describe the system in terms of the energies of the individual atoms, along with the energy of interactions between atoms given by the off-diagonal terms. The original Rupp et al. paper compared two Coulomb matrices by looking at the difference between their eigenvalues, with small differences implying similar molecules. These eigenvalues can be thought of as effective atomic energies within the molecule, that is, energies of the individual atoms with the influence of the surrounding atoms included. Similarly, the eigenvectors would describe fragments of the molecule or an effective atom which account for the influence of neighboring atoms.