A system of harmonic oscillators coupled via nonlinear interaction is a fundamental model in many branches of physics, from biophysics to electronics and condensed matter physics. In quantum optics, weak nonlinear interaction between light modes has enabled, for example, the preparation of squeezed states of light and generation of entangled photon pairs. While strong nonlinear interaction between the modes has been realized in circuit QED systems, achieving significant interaction strength on the level of single quanta in other physical systems remains a challenge. Here we experimentally demonstrate such interaction that is equivalent to photon up- and down-conversion using normal modes of motion in a system of two Yb ions. The nonlinearity is induced by the intrinsic anharmonicity of the Coulomb interaction between the ions and can be used to simulate fully quantum operation of a degenerate optical parametric oscillator. We exploit this interaction to directly measure the parity and Wigner functions of ion motional states. The nonlinear coupling, combined with near perfect control of internal and motional states of trapped ions, can be applied to quantum computing, quantum thermodynamics, and even shed some light on the quantum information aspects of Hawking radiation.