This is a continuation of the problem posed here.

Consider a 1-Lipschitz function $f: \mathbb R^n \to \mathbb R$ satisfying the inequality
\begin{align*}
|f(x) - f(y)| \le \|x-y\|_2, \;\forall x,y \in \mathbb R^n.
\end{align*}
For $n \ge 2$, can we find a 1-Lipschitz function that saturates the above inequality *on the average*?

To make the notion of "on the average" precise, let $x$ and $y$ be drawn from distributions with independent sub-Gaussian coordinates with bounded sub-Gaussian norms: $\max_i \|x_i\|_{\psi_2} = O(1)$ and similarly for $\{y_i\}$. We also assume that $\mathbb E x_i = \mathbb E y_i = 1$ and $\mathbb E x_i^2 = \mathbb E y_i^2 = 1$.

One can show that (this follows e.g. from Theorem~3.1.1 here) $$\mathbb E \big|\|x\|_2 - \|y\|_2\big| \asymp 1.$$ while $$\mathbb E\|x-y\|_2 \asymp \sqrt{n}.$$ Is there a 1-Lipschitz function $f : \mathbb R^n \to \mathbb R$ such that $$\mathbb E|f(x) - f(y)| \asymp \sqrt{n}?$$ Here $\asymp$ means inequalities go in both directions up to constants. For the Gaussian case, the answer is negative as was pointed out by Iosif Pinelis in response to the first part of this question. How about general sub-Gaussian distributions?

A related question is determining the order of $$ \sup_{f \in \text{Lip}(1)}\mathbb E|f(x) - f(y)| $$ where $\text{Lip}(1)$ is the set of $1$-Lipschitz functions from $\mathbb R^n$ to $\mathbb R$.