As a remark, notice this construction readily generalizes to mixtures with more than two components.

Application

As requested in comments, it would be of interest to choose $\lambda$ to meet a set of criteria, such as

1. Give the components equal weights, which means $$\frac{1}{2}=\pi_\lambda = \int \lambda(z) f(z)\,\mathrm{d}z.$$

2. Since these are intended to model errors in a regression setting (with $\mu=0,$ we would like each of the components also to have zero mean: $0 = E_{F_\lambda}[X].$ In light of (1), that is equivalent to $$0 = \int z\lambda(z) f(z)\,\mathrm{d}z.$$

3. Since regression errors are often assumed to be homoscedastic -- of equal variances -- we would like the variances of $F_\lambda$ and $F_{1-\lambda}$ to be equal. Since they have means of zero, when $f$ is a Normal density, this is achieved when $$\sigma^2 = 2\int z^2\lambda(z) f(z)\,\mathrm{d}z.$$

Although there are many solutions to these equations, one simple (striking) solution is obtained by supposing $\lambda$ and $1-\lambda$ are both simple functions: that is, piecewise constant. By making $\lambda$ symmetric around $0$ we can assure that (2) holds. The simplest of such simple functions is zero except on some positive interval $[a,b]$ and its negative $[-b,-a],$ where it equals $1.$

Without any loss of generality take $\sigma^2=1,$ so that $f = \phi$ is the standard Normal density with the property $\phi^(z) = -z\phi(z).$ Using this fact we may compute

$$\int \lambda(z)\phi(z)\,\mathrm{d}z = 2 \int_a^b \phi(z)\,\mathrm{d}z = 2(\Phi(b)-\Phi(a))$$

(where $\Phi$ is the standard Normal distribution function) and

$$\begin{aligned} \int z^2 \lambda(z)\phi(z)\,\mathrm{d}z &= 2 \int_a^b z^2\phi(z)\,\mathrm{d}z \\ &= 2(\Phi(b) - \Phi(a) + a\phi(b) - b\phi(b)). \end{aligned}$$

This permits numerical solution of (1) and (3). The work is streamlined by noting from (1) that, given $0 \le a\lt \Phi^{-1}(3/4),$

$$b = b(a) = \Phi^{-1}(\Phi(a) + 1/4).$$

That leaves us to solve (3) for $a \ge 0$. Here is an R implementation to illustrate:

f <- function(a) {
  b <- qnorm(1/4 + q <- pnorm(a))
  pnorm(b) - q + a * dnorm(a) - b * dnorm(b) - 1/4
}
uniroot(f, c(0, qnorm(3/4)- 1e-6))$root -> a
qnorm(pnorm(a) + 1/4) -> b

This calculation gives $a \approx 0.508949$ and $b \approx 1.59466.$ Here are plots of the two component densities $f_\lambda$ and $f_{1-\lambda}:$

To illustrate the intended application, here are bivariate data with 150 responses at $X=0$ with errors distributed as $F_\lambda$ and 150 responses at $X=1$ with errors distributed as $F_{1-\lambda}.$ To the right is a quantile plot of the collected residuals.

Although separately neither group of residuals appears Normal, they are both centered at zero, have nearly the same variance, and collectively look perfectly Normal.

Remarks

1. The basic construction readily generalizes to mixtures with more than two components.

2. The example in the application can be extended, by using simple (indicator) functions supported on intervals $[a_i,b_i]$ with $0\le a_1 \lt b_1 \le a_2 \lt b_2 \cdots \lt b_k,$ to create component distributions that match the first $2k$ moments of the Normal distribution their mixture creates. With sufficiently large $k,$ the component distributions will be difficult to discriminate even with largish datasets (at which point one might legitimately wonder whether their non-Normality matters at all).