Recovering correlations

I haven't looked at your whole simulation, but I can make a few quick
comments that hopefully help

I have thought about the following models

m2 <- lmer(y ~ 1 + (0+R1|f) + (0+R2|f), data=dat)

m3 <- lmer(y ~ 1 + (1|f) + (0+R2|f), data=dat)

but I've been struggling to figure out a way to recover the correlations r1
and r2 with the variances from the random effects based on the models m2
and m3.

Splitting these out into separate |f terms suppresses the computation of
the correlations. If you want the associated correlations, then you need
to have them in one group, e.g. (0 + R1 + R2 | f).

lme4 uses an unstructured covariance structure, so there's no real
middle ground between estimating all correlations (one grouping term) or
estimating no correlations (split things into separate grouping terms)
without having repeated/overlapping grouping terms.

There's nothing wrong per se with a correlation you assume is zero --
maybe you'll even get zero back. However, if you want to make that
assumption stronger and selectively constrain certain elements in the
covariance matrix to be zero, you'll probably need to use nlme (I think
you could hack a solution there?), handroll something in TMB or use a
Bayesian method where you're expressing exactly the model structure you
want (either through the likelihood or the prior). If you're willing to
go over to the Julia side and use MixedModels.jl, it wouldn't be too
difficult to construct a model, then selectively constrain certain
elements of the covariance matrix to zero. We don't have a interface for
this at the moment, but the code for ZeroCorr (roughly equal to || in
lme4) just does this for every off-diagonal element, and it wouldn't be
too hard to adapt that code for a single model to only zero-out certain
off-diagonal elements.

Also, for the construction of dummy variables in the random effects,
make sure to check out lme4::dummy.

Hope this helps

Phillip