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?

The original formulation of the AIC requires that the different observations that you use to fit your model are i.i.d.. If this condition is met, you can safely use the AIC to compare both of your models.
I highly recommend the following seminal paper on the AIC and the BIC:
Burnham, K. P., & Anderson, D. R. (2004). Multimodel inference: understanding AIC and BIC in model selection. Sociological methods & research, 33(2), 261-304.
Just a side note: in this very interesting recent paper, the authors extend the use of the BIC to some cases of correlated observations, and especially to the case of time-invariant Hidden Markov Models.
Regarding the "under the hood transformations" you mentioned, and adding on leviva's answer, one should be careful about the AIC function (I remember, when using the spatstats package, that the function used an incorrect number of free parameters). You can easily alleviate this issue by computing the AIC yourself using the likelihood of both your models.
