Comparing AIC of Tensor Product Smooths versus Thin Plate Splines I'm comparing the AIC of these two models.
Tensor Product Smooth vs. Thin Plate Spline both fit using REML
library(mgcv)

ThinPlateSpline <- gam(Z ~ s(X1, X2, bs="tp", k = 17), data = DF,  method="REML")
TensorProductSmooth <- gam(Z ~ te(X1, X2, k = 6), data = DF,  method="REML")

AIC(ThinPlateSpline, TensorProductSmooth)

Is the AIC comparison mathematically justifiable, if not, why not?

In theory AIC does not require the models to be from the same family. So you can compare models with different forms.
But in practice some implementations omit constants from the likelihood function as they're not necessary for comparison of models from the same family (for example two linear regressions). This can cause errors when comparing models with different likelihood functions.
So you need to check before using any implementation.
