  Graphene bilayers display peculiar electronic and mechanical characteristics
associated to their two-dimensional character and relative disposition of the
sheets. Here we study nuclear quantum effects in graphene bilayers by using
path-integral molecular dynamics simulations, which allow us to consider
quantization of vibrational modes and study the effect of anharmonicity on
physical variables. Finite-temperature properties are analyzed in the range
from 12 to 2000~K. Our results for graphene bilayers are compared with those
found for graphene monolayers and graphite. Nuclear quantum effects turn out to
be appreciable in the layer area and interlayer distance at finite
temperatures. Differences in the behavior of in-plane and real areas of the
graphene sheets are discussed. The interlayer spacing has a zero-point
expansion of $1.5 \times 10^{-2}$ \AA\ with respect to the classical minimum.
The compressibility of graphene bilayers in the out-of-plane direction is found
to be similar to that of graphite at low temperature, and increases faster as
temperature is raised. The low-temperature compressibility increases by a 6\%
due to zero-point motion. Especial emphasis is laid upon atomic vibrations in
the out-of-plane direction. Quantum effects are present in these vibrational
modes, but classical thermal motion becomes dominant over quantum
delocalization for large system size. The significance of anharmonicities in
this atomic motion is estimated by comparing with a harmonic approximation for
the vibrational modes in graphene bilayers.
