  In this paper a systematic examination of graphene/hexagonal boron nitride
(g/hBN) bilayers is presented, through a recently developed two-dimensional
phase field crystal model that incorporates out-of-plane deformations. The
system parameters are determined by closely matching the stacking energies and
heights of graphene/hBN bilayers to those obtained from existing
quantum-mechanical density functional theory calculations. Out-of-plane
deformations are shown to reduce the energies of inversion domain boundaries in
hBN, and the coupling between graphene and hBN layers leads to a bilayer defect
configuration consisting of an inversion boundary in hBN and a domain wall in
graphene. Simulations of twisted bilayers reveal the structure, energy, and
elastic properties of the corresponding Moir\'e patterns, and show a crossover,
as the misorientation angle between the layers increases, from a well-defined
hexagonal network of domain boundaries and junctions to smeared-out patterns.
The transition occurs when the thickness of domain walls approaches the size of
the Moir\'e patterns, and coincides with the peaks in the average von Mises and
volumetric stresses of the bilayer.
