I am having some struggle understanding how the change of variable works when using the iso-parametric elements in FEM.
From the weak form of the steady heat equation, I get two integrals
$$ \int_\Omega \nabla \phi_i \nabla \phi_j d\Omega $$ $$ \int_\Omega \phi_i d\Omega $$
In that equations, $\phi$ is defined in the global space (X, Y).
Now, we introduce the natural space. I can define some shape functions in that natural space $\psi$, and for an element of N nodes, I can define
$$ x = \sum_i^N{\psi_i x_i} $$ $$ y = \sum_i^N{\psi_i y_i} $$
Now we need to do the change of variable in the integral, so we need to change three things: integration limits, integration space, and apply the change of variable to the integrands.
Normally, I just see that wherever appears $\phi(x,y)$ in the integral, it is substituted by $\psi(\zeta, \eta)$. Intuitively I can see why $\phi(x,y) = \psi(\zeta, \eta)$. For a point in the global space, the value of $\phi$ must match the value of $\psi$ for its matching point. However, I have never seen a proof of this property.
Furthermore, for the integral limits, I just see that $dxdy$ is subtituted by $|J| d\zeta d\eta$, but I don't fully understand why.
The initial integrals, when transforming to the natural space would convert to
$$ \int_{\Omega_X} \nabla \phi_i \nabla \phi_j dx dy = \int_{\Omega_{\zeta}} (J^{-1} \nabla \psi_i) (J^{-1} \nabla \psi_j) |J| d\zeta d\eta $$
$$ \int_{\Omega_X} \phi_i dx dy = \int_{\Omega_{\zeta}} \psi_i |J| d\zeta d\eta $$
but, how is the change of variables performed in the integral?
Approximation of Function
Suppose you have some differential equations with unknown analytical solution $u(x,y)$. The main idea of FEM is to approximate the analytical solution using functions that are simple and easy to work with, in normal FEM case, polynomials:
$$ u(x,y) \approx \sum_{i=0}^{n} a_{ij} x^{p_{i}}y^{q_{i}} \phantom{x} \forall \phantom{x} (x,y) \in \Omega_{j} $$
You probably notice that this form is not too convenient, because it is quite difficult to see the correlation between the approximation at some points and the coefficients $a_{i}$ without explicitly calculating the approximation. A more convenient form is given below, which is used in FEM:
$$ u(x,y) \approx \sum_{i=0}^{n} \hat{u}(x_{ij},y_{ij})\cdot\phi_{ij}(x,y) \phantom{x} \forall \phantom{x} (x,y)\in\Omega_{j} $$
The polynomials $\phi_{ij}(x,y)$ are defined such that $u(x_{ij},y_{ij})\approx\hat{u}(x_{ij},y_{ij})$. These polynomials are often called shape functions or interpolating functions.
Natural Space
The shape functions $\phi_{ij}(x,y)$ depends on the domain. This is not efficient because you need to formulate different polynomials for different domain.
The solution is mapping to natural space and isoparametric approach. First let's take a look at mapping to natural space:
$$ \mathbf{x}_{j} = \sum_{i=0}^{n} \psi_{i}(\xi,\eta) \cdot \mathbf{x}_{j}(\xi_{i},\eta_{i}) $$
and the approximation in this natural space is given below:
$$ \begin{aligned} u(\xi,\eta) &\approx \sum_{i=0}^{n} \hat{u}\left(x_{j}(\xi_{i},\eta_{i}),y_{j}(\xi_{i},\eta_{i})\right) \cdot \phi_{ij}\left(x_{j}(\xi,\eta),y_{j}(\xi,\eta)\right) \\ &= \sum_{i=0}^{n} \hat{u}(\xi_{ij},\eta_{ij}) \cdot \tilde{\phi}_{ij}(\xi,\eta) \end{aligned} $$
Isoparametric Approach
Up to this point, we have two sets of polynomials for each domain;
$$ \tilde{\phi}_{ij}(\xi,\eta) \phantom{x} \text{and} \phantom{x} \psi_{i}(\xi,\eta) $$
Remember that $\tilde{\phi}_{ij}(\xi,\eta)$ differs from domain to domain and we don't like that. Consequently, some smart guys long time ago proposed that we simply use the same sets of polynomials:
$$ \tilde{\phi}_{ij}(\xi,\eta) = \psi_{i}(\xi,\eta) $$
The term isoparametric approach comes from using the same parametric formulation to map between global & natural space and to approximate the analytical solution. There is no proof for "why the two shape functions must be equal" because it is not necessarily true, it is just a simplification that we decided to make.
Jacobian
I recommend you to revise your linear algebra course or watch the excellent series "Essence of linear algebra" by 3Blue1Brown on YouTube. From there you can understand why determinant of Jacobian gives ratio of surface area (or volume if 3D) in one space and surface area of the same domain when mapped to another space. Also why determinant of Jacobian show up during integration.