Since you are assuming $\eta$ is normal what i would do is try to compute the expectation as fast as possible for every $\theta$. So I would compute the expectation using any numerical integration I might know. That is, if $n$ and $N$ are the Gaussian density and distribution respectively,
$$ f( \theta ) = \int_{ -\infty }^\infty R(\theta, \eta ) n( \eta ) d\eta = \int_0^1 R( \theta, N^{-1}(w ) ) dw\\
=\sum_{j=1}^n \int_{w_{j-1} }^{w_j} R( \theta, N^{-1}(w ) ) dw,$$
where $0=w_0 < w_1 < ....< w_n = 1 $. You can use any integration tricks on each of the integrals above. For example, take $w_j = N(z_j)$ for some partition $(z_j)_{j>0 }$ of the real line, and approximate
$$
\int_{w_{j-1} }^{w_j} R( \theta, N^{-1}(w ) ) dw = R( \theta, N^{-1}(w_j ) )( w_j - w_{j-1} ) = R( \theta, z_j )( w_j - w_{j-1} ).
$$
Once you have this step, just use any root finding method like brent's algorithm to find $\theta$. The integration needs to be fast because every step in the root finding is going to perform a new integration. Hence in your implementation you should be able to re-use parameters in the integration step, for example $z$, $w$, etc.