  The random phase approximation (RPA) as formulated as an orbital-dependent,
fifth-rung functional within the density functional theory (DFT) framework
offers a promising approach for calculating the ground-state energies and the
derived properties of real materials. Its widespread use to large-size, complex
materials is however impeded by the significantly increased computational cost,
compared to lower-rung functionals. The standard implementation exhibits an
$\mathcal{O}(N^4)$-scaling behavior with respect to system size $N$. In this
work, we develop a low-scaling RPA algorithm for periodic systems, based on the
numerical atomic orbital (NAO) basis-set framework and a localized variant of
the resolution of identity (RI) approximation. The rate-determining step for
RPA calculations -- the evaluation of non-interacting response function matrix,
is reduced from $\mathcal{O}(N^4)$ to $\mathcal{O}(N^2)$ by just exploiting the
sparsity of the RI expansion coefficients, resultant from localized RI (LRI)
scheme and the strict locality of NAOs. The computational cost of this step can
be further reduced to linear scaling if the decay behavior of the Green's
function in real space can be further taken into account. Benchmark
calculations against existing $\textbf k$-space based implementation confirms
the validity and high numerical precision of the present algorithm and
implementation. The new RPA algorithm allows us to readily handle
three-dimensional, closely-packed solid state materials with over 1000 atoms.
The algorithm and numerical techniques developed in this work also have
implications for developing low-scaling algorithms for other correlated methods
to be applicable to large-scale extended materials.
