  We present the first $N_f=2+1+1$ results for the matrix elements of the
operators describing neutral $K$ and $D$ mixing in the Standard Model and its
extensions. The combination of maximally twisted sea quarks and
Osterwalder-Seiler valence quarks ensures $\mathcal{O}(a)$-improvement and
continuum like renormalization pattern. We have used the $N_f=2+1+1$ dynamical
quark gauge configurations generated by ETMC. Simulations include three lattice
spacings in the interval $[0.06:0.09]$ fm and pseudoscalar meson masses in the
range $[230:500]$ MeV. Our results are extrapolated to the continuum limit and
to the physical quark masses. The calculation of the renormalization constants
has been performed non-perturbatively in the RI-MOM scheme.
