Since you are nitpicking, you could well use a sinusoid activation function on the Neural Network, and reach an even smaller loss value.
Splines don't know anything about the nature of a function. They approximate any function with piecewise polynomials.
Maybe you are trying to say that the default activation function (relu) in sklearn is not smooth. No problem, you can add
activation='tanh'
In general you can use a pretty wide set of functions to approximate an arbitrary function. You can do it with polynomials (Taylor expansion), and many others as long as they form a Hilbert space.
Producing a given function from a linear combination of other functions isn't groundbreaking in the least.