A comment (well, a few observations), not an answer but too long for that.
Given vector spaces $V,W$ ,a bilinear mapping from $V\times W$ into the reals is one which is linear in each variable separately. Crucial fact: such an operator is uniquely determined by its vales on elements $(x,y)$ where $x$ and $y$ range over bases (yes, I know that this is undergraduate maths but the OP insists on details).
If $X$ and $Y$ are planar vectors and $L(X,Y)$ is the signed area of the square $ABCD$ erected on $A=X$ and $B=Y$, then this is bilinear in $X$ and $Y$ (easy to prove directly in coord8nates using the area formula or by a simple geometrical argument— not my fault—the OP insists on details).
Using these facts, then the case of squares can be reduced to that of determining when two bilinear mappings coincide. This can easily be done with a simple Mathematica programme or by using symmetry and the canonical basis for $2$-space.
This method can be carried over to the more general case of polygons and to gazillions of natural generalisations but I am too terrified of the OP‘s potential comments to even think about posting them here.
Added as an edit—the square case in more detail:
We assume that $A_¡$ is represented by the vector $X_¡$, $B_i$ by $Y_i$—then $C_i$ and $D_i$ are $Y_i+DY_i-DX_i$ and $X_i+DY_i-DX_i$ ($D$ denotes rotation through a right angle).
Now the area of quadrilateral $A_1A_2A_3A_4$ is $$ \frac 12 (X_1\wedge X_2+X_2\wedge X_3 +X_3\wedge X_4+X_4\wedge X_1)$$ with corresponding formulae for the other three quadrilaterals. One can then plug these expressions into the desired equations and complete the proof by routine, if tedious, computations. As I indicated previously, one can, if desired, take advantage of the fact that the expressions are bilinear forms on the corresponding vector spaces.
The general case becomes a battle with indices. If $D$ now denotes rotation through the exterior angle of the regular polygon and $X_i$,$Y_i$ are the vectors corresponding to $A_{i1}$ and $A_{i1}A_{i2}$, then $$ A_{ik}=X_{i}+Y_{i}+DY_i+\dots +D^{k-2}Y_i. $$
We now use the formula for the area of the polygon with vertices $Z_1,\dots,Z_p$,
$$ \frac 12 (Z_1\wedge Z_2+ \dots + Z_p\wedge Z_1),$$ substitute the appropriate expressions and collect terms.
This is just a sketch, I know, but hope that it will be useful to the OP.