This paper presents density evolution (DE) techniques based on reciprocal channel approximation (RCA) for general bit-interleaved coded modulation (BICM) systems. The original RCA-based DE (RCA-DE) technique is implemented under the assumption of binary-input additive white Gaussian noise channels (BI-AWGNCs). In BICM systems, however, $M$ -ary modulation schemes are generally used, and they can be regarded as to be a parallel transmission of $\log _{2} M$ bits. Thus, $\log _{2} M$ separate bit-level channels, which are not Gaussian equivalent, need to be considered. In order to extend the conventional RCA-DE technique to BICM systems, we first establish a model of protograph BICM ensembles. Two methods are developed based on the bit error rate (BER) and the BICM capacity in order to find a corresponding BI-AWGNC that models each of the bit-level channels. Using these methods, we implement the protograph RCA-DE technique for BICM systems and show that it achieves an accurate estimation. As a practical application of the proposed RCA-DE method, we design bit interleavers in order to achieve better BICM performance in the 3GPP New Radio (NR) LDPC coding system. Numerical results show that performance gain of up to 0.3 dB is consistently achieved over a wide range of parameter values.