$\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.