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.