  Mergers of black hole (BH) and neutron star (NS) binaries are expected to be
observed by gravitational wave observatories in the coming years. Until now,
LIGO has only set an upper limit on this merger rate. BH-NS binaries are
expected to merge in isolation, as their formation is suppressed in star
clusters by the mass segregation and the strong heating by BHs in the cluster
core. Another viable scenario to BH-NS mergers is in triple systems. In this
paper, we carry out a systematic statistical study of the dynamical evolution
of triples comprised of an inner BH-NS binary by means of high-precision
$N$-body simulations, including Post-Newtonian (PN) terms up to 2.5PN order. We
start from the main sequence massive stars and model the supernovae (SN) events
that lead to the formation of BHs and NSs. We adopt different prescriptions for
the natal velocity kicks imparted during the SN processes and illustrate that
large kicks lead to more compact and massive triples that merge on shorter
timescales. We also show that BH-NS merging in triples have a significant
eccentricity in the LIGO band, typically much larger than BH-NS merging in
isolated binaries. Finally, we estimate a rate of $\Gamma_\mathrm{BH-NS}\approx
1.0\times 10^{-3}-3.5\times 10^{-2} \ \mathrm{Gpc}^{-3}\ \mathrm{yr}^{-1}$, for
non-zero velocity kicks, and $\Gamma_\mathrm{BH-NS}=19 \ \mathrm{Gpc}^{-3}\
\mathrm{yr}^{-1}$, for no natal kicks. Our rate estimate overlaps with the
expected BH-NS rate in isolated binaries and within the LIGO upper limit.
