Why is cubic spline interpolation used to analyze light curve of pulsating stars? Is there any scientific justification of such usages? One example of using the cubic spline interpolation is here, section 4.
The main thing a cubic spine lets you do is have a continuous second derivative. If all you need is a continuous curve, connecting the dots with straight lines suffices, with no overshoot, but a very discontinuous derivative. So if you want to think about the derivative of the light curve, you'd have to go at least to a quadratic spline, which allows a continuous first derivative but introduces the potential for overshooting (a curve that goes outside the range of the points themselves). If you actually want to think about the second derivative also, then you'll need that to be continuous, so that will require a cubic spline. But the potential for overshoot keeps getting worse when you add more continuous derivatives. Personally, I find cubic splines to be very unreliable, the overshoot can be terrible unless you have a very nicely resolved curve with not much noise, but that may be just what they do have.
I can't speak for the general case; I'm not aware that cubic splines are always used for modelling light curves. However, a cubic spline is the lowest order of spline that allows for inflection points, thus offering a good combination of flexibility with simplicity. A cubic spline could fit an arbitrary smooth curve with only a few control points. It is also offered in many data analysis and drawing packages.
In the paper, the authors are fitting a light curve by eye. Cubic splines are an "obvious" way of doing this.
Cubic spline interpolations are used because if you have a good pipeline you only need to adjust 2 variables to get a pretty good fit. Granted, a cubic spline have numerous limitations, and offers a ton of inflection points that won't allow for best fitting minimization. That being said, it is fast, dirty, easy, and widely accepted as good enough. Plus, it oscillates just like a pulsating or variable star, so that's another good excuse to go the cubic spline route.
There are literally innumerable scientists and papers that base all of their results on what the cubic spline gives them, in fields from photometry as you've linked, to temperatures derived from Ly$\alpha$ absorption profiles as seen here.
This is a link to a reference on spline interpolation I used in my undergraduate thesis. You might find it helpful.