How to shrink a square with minimal distortion? $\newcommand{\CO}{\text{CO}_2}$
$\newcommand{\euc}{\mathfrak{e}}$
$\newcommand{\SO}{\text{SO}_2}$
$\newcommand{\al}{\alpha}$
$\newcommand{\dist}{\text{dist}}$
$\newcommand{\Lip}{\text{Lip}_{\text{inj}}}$
$\newcommand{\sAverage}[1]{\langle#1\rangle}$
$\newcommand{\IP}[2]{\sAverage{#1,#2}}$
$\newcommand{\pl}{\partial}$
$\newcommand{\Sone}{\mathbb{S}^1}$
Let $0<\sigma_1<\sigma_2$ satisfy $\sigma_1\sigma_2=1$, and let $D=[-1,1]^2$.

Does there exist a Lipschitz injective a.e. map $\phi:D \to D$ such that $d\phi$ has almost everywhere the fixed singular values $\sigma_1,\sigma_2$?

Does there exist such a $\phi$ that maps the boundary onto itself?
By injective a.e. I mean that  $|\phi^{-1}(y)| \le 1$ for a.e. $y \in D$.
Edit
If $[-1,1]^2$ is replaced by a disk, then each $ \phi_t:(r,\theta) \to (r,\theta+t \log r)$ has constant singular values. Unfortunately, an analogous approach doesn't work for the square. Here is why:
Let $\alpha(\theta)$ be an arc length parametrization of $\partial D$, $D=[-1,1]^2$, i.e. $\alpha:\frac{8}{2\pi}\mathbb{S}^1 \to \partial D$, $|\dot \alpha|=1$. Let $h:(0,1) \times \frac{8}{2\pi}\mathbb{S}^1 \to \frac{8}{2\pi}\mathbb{S}^1$, and assume that for every $r \in (0,1)$, $h(r,\cdot)$ is a diffeomorphism of $\frac{8}{2\pi}\mathbb{S}^1$. Define
$$
f\big(r\alpha(\theta)\big)=r\alpha(h(r,\theta)).
$$
$f$ maps each scaled copy $r\text{Image}(\alpha)$ diffeomorphically onto itself. The case of the "logarithmic" map $\phi_t$ above corresponds to choosing $h(r,\theta)=\theta+t \log r$.
Thinking of $f$ as a map
$$
(r,\theta) \mapsto (r,h(r,\theta)),
$$
one gets that
$$
[df]_{E_1,E_2}=\begin{pmatrix} h_{\theta} & \bigg(rh_r(r,\theta)+b(h(r,\theta))-b(\theta) h_{\theta}(r,\theta)\bigg)/s(\theta) \\0 & s(h(r,\theta))/s(\theta) \end{pmatrix}, \tag{1}
$$
where
$
s:=|\al \times \dot \al|, b:=\IP{\al}{\dot \al},
$
and $(E_1,E_2)$ is an orthonormal frame in the domain of parametrization $(0,1] \times \frac{8}{2\pi}\Sone$, given by
$$
(E_1,E_2)_{r,\theta}=\big(\frac{1}{r} \pl_{\theta},  \frac{1}{|\al \times \dot \al|}  (\pl_r-\frac{\IP{\al}{\dot \al}}{r} \pl_{\theta})\big).
$$
More explicitly, we define  $\Phi:(0,1] \times \frac{8}{2\pi}\Sone \to (D,\euc) $, where $\euc$ is the standard Euclidean metric, by
setting $\Phi(r,\theta):=r\al(\theta)$. $(E_1,E_2)$ is an orthonormal frame in $(0,1] \times \frac{8}{2\pi}\Sone$ w.r.t the pullback metric $\Phi^*\euc$.
For the case where the domain $D=[-1,1]^2$ is a square, $s=|\al \times \dot \al|$ is constant, so Equation $(1)$ reduces to
$$
[df]_{E_1,E_2}=\begin{pmatrix} h_{\theta} & \bigg(rh_r(r,\theta)+b(h(r,\theta))-b(\theta) h_{\theta}(r,\theta)\bigg)/s \\0 & 1 \end{pmatrix}, \tag{2}
$$
so $f$ is area-preserving if and only if $h_{\theta}=1$, i.e. $h(r,\theta)=\theta+g(r)$.
Then Equation $(2)$ further reduces to
$$
[df]_{E_1,E_2}=\begin{pmatrix} 1 & \big(rg'(r)+b(\theta+g(r))-b(\theta) \big)/s \\0 & 1 \end{pmatrix}.
$$
Thus $df$ has constant singular values if and only if
$$
A(r,\theta):=rg'(r)+b(\theta+g(r))-b(\theta)
$$
is constant.
$b(\theta)$ is piecewise-affine in $\theta$, with jumps at  the corners, so the locations of jumps in the difference $b(\theta+g(r))-b(\theta)$ depend on both $r,\theta$ and cannot be cancelled by the remaining term $rg'(r)$.
Variational motivation for the question:
Such maps $\phi$ correspond to distortion minimizing maps from the square into a smaller square. Let $0<\lambda \le 1/2$, and let $\Lip(D,\lambda D)$ be the space of injective a.e. Lipschitz maps
$D \to \lambda D$ having nonnegative Jacobian. Set $E:\Lip(D,\lambda D) \to \mathbb{R}$ by
$$E(\phi)=\int_{D} \dist^p( d\phi,\SO )\,\,dx.$$
I proved (Theorem 1.1 here) that for $p>2$, $E(\phi) \ge (1-2\lambda)^{p/2}$, and that equality holds if and only if $\sigma_1(d\phi)+\sigma_2(d\phi)=1, \sigma_1(d\phi)\cdot \sigma_2(d\phi)=\lambda^2$.
$\frac{\phi}{\lambda}:D \to D$ has constant distinct singular values.
The motivation for choosing specifically the square, is that we can tile other shapes in the plane up to an arbitrary accuracy, so if we have a solution on a square, we can glue these solutions to approach an infimal distortion energy for other shapes as well.
 A: Update. I'll show in the end of this post in the proof that the answer to the question is positive for all $\sigma_1\in (0,1)$.
The square is denoted by $\square$.
Set up. Let us fix $x,y\in (0,1)$ with $x^2-2x+y^2<0$ and denote by $\square_{x,y}$ the square with vertices
$$(x,y), (y,-x), (-x,-y), (-y,x)$$
Now let us start do define the map $\varphi_{x,y}:\square\to \square$. The map will be identical on the boundary of the square. First we define it on $\square\setminus \square_{x,y}$. The map will send the boundary of $\square_{x,y}$ to the boundary of $\square_{y,x}$. On vertices it acts as follows:
$$(x,y)\to (x,-y),\, (y,-x)\to (-y-x),\, (-x,-y)\to (-x,y),\, (-y,x)\to (y,x)$$
This is a is depicted on the figure. The set $\square\setminus \square_{x,y}$ is cut into $8$ triangles (with green boundary)  and the map $\varphi_{x,y}$ is linear on them. Now we extend this map to the  $\square_{x,y}\setminus \square_{x',y'}$ and so on.
Easy claim. The constructed map is area preserving and $\sigma_1$ takes only two values - for half of triangles it is one and for the other half it is different.
Finally, we should choose $(x,y)$ so that $\sigma_1$ takes just one value. By simple continuity it is clear that such $(x,y)$ will form a non-empty (algebraic) curve. The equation of the curve is given below.
In the above figure the brown triangle on the left goes to the brown triangle on the right.
Added. So, for this map to have constant $\sigma_1$ we need an equality on $x,y$ to hold. To write down the equality, denote by $A$ the matrix
$$A=\begin{pmatrix}
y+1 & 1-x\\
1-x & 1-y
\end{pmatrix}$$
Then we should have
$$2+\frac{4y^2}{(1-x)^2}=tr\left( A\cdot\begin{pmatrix}
0 & 1\\
1 & 0
\end{pmatrix} \cdot A^{-2}\cdot
 \begin{pmatrix}
0& 1\\
1 & 0
\end{pmatrix} \cdot A\right)
$$
This is the equality between the traces of matrices $BB^T$ and $CC^T$, where  $B$ corresponds to the linear map $\varphi_{x,y}$ induced on a triangle adjacent to an edge of $\square$ and $C$ corresponds to a triangle adjacent to an edge of $\square_{x,y}$.
Proof. So let us denote by $D$ the domain $$y\ge 0,\, x<1,\, x^2-2x+y^2<0$$ where the vertex $(x,y)$ of $\square_{x,y}$ can lie, so that the map $\varphi_{x,y}$ is well defined. On the figure below the domain is in light blue.
Let us introduce two functions on $D$: $\sigma_2^A(x,y)$ and $\sigma_2^B(x,y)$. The first function is equal to $\sigma_2$ of the map $\varphi_{x,y}$ restricted to triangle $A$ and the second is $\sigma_2$ of the map $\varphi_{x,y}$  restricted to triangle $B$. Now, we have the following claims:

*

*$\sigma_2^A(x,y)$ and $\sigma_2^B(x,y)$ are smooth on $D$.


*On the segment $y=0$ both functions are equal to $1$


*When the point $(x,y)$ approaches the circular arc $x^2-2x+y^2=0$ the function $\sigma_2^B(x,y)$ tends to $+\infty$, while $\sigma_2^B(x,y)$ stays bounded.


*When the point $(x,y)$ approaches the segment $x=1$ the function $\sigma_2^A(x,y)$ tends to infinity, while $\sigma_2^A(x,y)$ stays bounded.
Now if one stares a bit more on the figure on trajectories in pink, one sees that the function $\sigma_2^A(x,y)$ goes from $1$ to $+\infty$ when we move along any such trajectory from $(0,0)$ to $x=1$. At the same time  $\sigma_2^B(x,y)$ stays bounded. Moreover, on trajectories that are close to $y=0$ it is close to $1$, while on trajectories close to $x^2-2x+y^2<0$ it is large. This proves that on any trajectory there is a point where $\sigma_2^A(x,y)=\sigma_2^B(x,y)$. By varying further the pink trajectory from the segment $y=0$ to the circular arc we get any $\sigma_2\ge 1$ that we want.

PPS I add here one more picture. It is the previous one rotated by 90 degrees, it is a bit easier to see it this way. In pink I draw the level sets $\sigma_2^B=const$, in green $\sigma_2^A=const$. For $const=1$ these curve coincide with the vertical green-pink segment. Otherwise they are different and they always intersect.

