As a start, let us first derive the formula when the unit square is centered at origin.
Let $a = x^2+y^2+\frac12$, the distance to vertex $v$ is:
$$\begin{cases}
r_{++} &= \sqrt{a - x - y}, & v = (+\frac12,+\frac12)\\
r_{+-} &= \sqrt{a - x + y}, & v = (+\frac12,-\frac12)\\
r_{-+} &= \sqrt{a + x - y}, & v = (-\frac12,+\frac12)\\
r_{--} &= \sqrt{a + x + y}, & v = (-\frac12,-\frac12)
\end{cases}$$
If the sum of these 4 distances is $R = \sqrt{z}$, the equation for the $4$-ellipse will be given by
$$r_{++} + r_{+-} + r_{-+} + r_{--} - R = 0\tag{*1}$$
To construct a polynomial curve that contains this $4$-ellipse, the standard trick
is apply all the symmetry operations in the associated Galois group to this expression and then take the product. For the case at hand, this becomes
$$\prod_{(\epsilon_{++},\epsilon_{+-},\epsilon_{-+},\epsilon_{--}) \in \{\pm 1\}^4}
\left( \epsilon_{++}r_{++} + \epsilon_{+-} r_{+-} + \epsilon_{-+} r_{-+} +
\epsilon_{--} r_{--} - \sqrt{z}\right) = 0
$$
If one throw this to an CAS, the product is equal to
$$\begin{array}{rrl}
& z^8\\
-32 & az^7\\
+32 & z^6 & (y^2+x^2+11a^2)\\
-256 & az^5 & (y^2+x^2+7a^2)\\
-256 & z^4 & (7y^4-20x^2y^2-6a^2y^2+7x^4-6a^2x^2-17a^4)\\
+4096 & az^3 & (3y^4-4x^2y^2-2a^2y^2+3x^4-2a^2x^2-a^4)\\
+8192 & z^2 & (2y^6-x^2y^4-4a^2y^4-x^4y^2+3a^2x^2y^2+2a^4y^2+2x^6-4a^2x^4+2a^4x^2)\\
+65536 & az\; & x^2y^2(y^2+x^2-a^2)\\
+65536 & & x^4y^4
\end{array}
\tag{*2}
$$
This is a a big mess even before we substitute $a$ by $x^2+y^2+\frac12$, $z$ by $R^2$ and fully expand it.
To obtain a more manageable formula, consider following two products,
$$\begin{align}
P_{+}(\lambda)
&= \prod_{(\epsilon_{++},\epsilon_{--}) \in \{\pm 1\}^2}
\left( \epsilon_{++}r_{++} + \epsilon_{--} r_{--} - \lambda\right)
= \lambda^4 - 4a\lambda^2 + 4(x + y)^2
\\
P_{-}(\lambda)
&= \prod_{(\epsilon_{+-},\epsilon_{-+}) \in \{\pm 1\}^2}
\left( \epsilon_{+-}r_{+-} + \epsilon_{-+} r_{-+} - \lambda\right)
= \lambda^4 - 4a\lambda^2 + 4(x - y)^2
\end{align}
$$
Notice
- $P_{+}(\lambda)$ contains a factor $r_{++} + r_{--} - \lambda$.
- $P_{-}(\lambda)$ contains a factor $r_{+-} + r_{-+} - \lambda$.
If one compute the resultant
of $P_{+}(\lambda)$ and $P_{-}(\sqrt{z} - \lambda)$
and eliminate $\lambda$, the resulting polynomial in $\sqrt{z}$ will contain
the desire factor $r_{++} + r_{--} + r_{+-} + r_{-+} - \sqrt{z}$ appeared in LHS of $(*1)$.
With help of a CAS again, up to a scaling factor, the resultant equals to
$$
\begin{align}
g(z,A,B) \stackrel{def}{=}
& \;\; z^2 (z^2-4Az+16A-32)(z^2-2Az+A^2-4A+8)^2\\
& + 2Bz (17z^3-20Az^2+11A^2z-28Az+56z-2A^3+8A^2-16A)\\
& + B^2
\end{align}
$$
where
$\displaystyle\;\begin{cases}
A &= 4x^2+4y^2 + 2\\
B &= 256x^2y^2
\end{cases}$.
If you fully expand $g(z,A,B) = 0$ and $(*2)$, you can verify they return the same polynomial.
To obtain the polynomial curve when the vertices are $(0,0)$, $(1,0)$, $(1,0)$ and $(1,1)$, one can replace $x$ by $x - \frac12$ and $y$ by $y - \frac12$ in
above definition of $A,B$ above. The polynomial curve is
$$g(R^2, 4(x^2+y^2-x-y+1), 16(2x-1)^2(2y-1)^2) = 0\tag{*3}$$
If I count correctly, when you expand this out, you will obtain a polynomial in $R,x,y$ with $289$ terms! This is way too big to include it directly in this answer.
If you really want the coefficients, I will recommend you either input $(*3)$ into an CAS and ask the CAS to expand for you or even better, reproduce the derivations here by an CAS.
For the case where the vertices are located at $(0,0), (n,0), (0,n), (n,n)$,
the formula can be obtained from $(*3)$ by substituting $(x,y,R)$ with $(\frac{x}{n},\frac{y}{n},\frac{R}{n})$.
eliminate([a^2=(x+1)^2+(y+1)^2,b^2=(x+1)^2+(y-1)^2,c^2=(x-1)^2+(y+1)^2,d^2=(x-1)^2+(y-1)^2,a+b+c+d=R],[a,b,c,d]);returns a polynomial $f(x,y,R)$ raised to some irrelevant power that I stripped off. The remaining $f(x,y,R)$ (implicitly set to zero) is an irreducible polynomial of degrees $(10,10,16)$ in $(x,y,R)$. Quite a huge expression; you wouldn't like it. – ccorn May 16 '16 at 00:35