For diffeomorphism $f$, if $X$ and $f(X)$ are both Gaussian, then $f$ is affine I am trying to prove the following.

Let $f:\mathbb{R}^{n}\to \mathbb{R}^n$ be a diffeomorphism. If $X$ and
$f(X)$  are both $n$ -dimensional Gaussian variables, then $f$ is
affine. That is, there exists a $n\times n$ matrix $A$ and $b\in \mathbb{R}^n$ such that $f(x)=Ax+b$.

Context
The problem I am trying to solve is the following optimal transport problem in the Monge setting:

Let $\mu_0$ be the probability distribution of $n$-dimensional
Gaussian distribution with mean $0$ and the covariance matrix
$\Sigma_0$. Define $\mu_1$ similarly with mean $0$ and the covariance
matrix $\Sigma_1$.
Find the diffeomorphism $\eta$ on $\mathbb{R}^n$ that minimizes
$\displaystyle \begin{equation}  J(\eta)=\int_{\mathbb{R}^n}d(x,\eta(x))^2d\mu_0\end{equation}\tag*{}$
under the constraint $\eta_*\mu_0=\mu_1$, where $\eta _*\mu_0$ is the
push-forward measure and
$d(x,y)$ is the Euclidean distance.

I have already solved this problem in the case $\eta$ is linear, so I figured I should reduce this problem into a linear one.
The condition $\eta_*\mu_0=\mu_1$ is equivalent to $\mu_0(\eta^{-1}(A))=\mu_1(A)$ for any measurable $A$. In probabilistic notation, we have $P(\eta(X)\in A)=P(Y\in A)$ where $X\sim N(0,\Sigma_0),Y\sim N(0,\Sigma_1)$. This implies $\eta(X)$, and $Y$ have the same distribution. Therefore, both $X$ and $\eta(X)$ are Gaussian (with zero mean).
If $n=1$, then we can easily prove that the diffeomorphism that makes $X$ and $f(X)$ both Gaussian are affine, as we can see here: Gaussian-to-gaussian transformations.
In higher dimensions, the paper Transformations preserving normality and Wishart-ness seems to suggest that we have that $f$ is affine for closely related, general case. (The paper asserts a bijective bimeasurable function that preserves normality for a fixed mean, and any covariance matrix is affine)
My Attempt
We can apply some affine transformation $L$ so that $L\circ f(X)$ has the same distribution as $X$. Thus without loss of generality, we can assume $X$ and $f(X)$ have the same distribution. For simplicity, I will assume $\Sigma_0=\Sigma_1=I$(the identity matrix).
My first guess was that for $f(X)=(Y_1,\cdots, Y_n)$, each $Y_1,\cdots Y_n$ should not depend on more than one $X_i$. This is wrong since $f(X)=\frac{1}{\sqrt{2}}(X_1+X_2,X_1-X_2)$ has the distribution $N(0,I)$ and $f$ is diffeomorphism.
I had a few more failed attempts, but I will omit them to make this post not too long.
 A: $\newcommand\R{\mathbb R}$Ben McKay's idea stated in the comment above is a natural one for a counterexample.
Indeed, for $(x,y)\in\R^2$, let
$$f(x,y):=f((x,y)):=
\left(x \cos \left(r^2\right)-y \sin \left(r^2\right),\ 
x \sin \left(r^2\right)+y \cos
   \left(r^2\right)\right),$$
where $r^2:=x^2+y^2$.
The transformation $f$ is bijective, with
$$f^{-1}((x,y))=
\left(x \cos \left(r^2\right)+y \sin \left(r^2\right),\ 
-x\sin \left(r^2\right)+y \cos
   \left(r^2\right)\right)$$
for all $(x,y)\in\R^2$.
Note also that
$$|f^{-1}((x,y))|^2=r^2 \tag{1}\label{1}$$
for all $(x,y)\in\R^2$, where $|\cdot|$ is the Euclidean norm.
Also,
$$\text{the Jacobian determinant of $f$ is $1\ne0$ everywhere on $\R^2$.} \tag{2}\label{2}$$
So, $f$ is a diffeomorphism.
Also, it follows from \eqref{1}, \eqref{2}, and the formula for the change of variables under the (double) integral sign that $f(X)$ is a standard normal random vector in $\R^2$ provided that $X$ is a standard normal random vector in $\R^2$.
However, it is clear that the transformation $f$ is not affine. (For instance, the partial derivative of the first coordinate of $f(x,y)$ with respect to $x$ equals $1$ at $(x,y)=(0,0)$ and $-1$ at $(x,y)=(0,\sqrt\pi)$.)
