  In most simulations of nonrelativistic nuclear systems, the wave functions
found solving the many-body Schr\"odinger equations describe the
quantum-mechanical amplitudes of the nucleonic degrees of freedom. In those
simulations the pionic contributions are encoded in nuclear potentials and
electroweak currents, and they determine the low-momentum behavior. In this
work we present an alternative quantum Monte Carlo formalism in which both
relativistic pions and nonrelativistic nucleons are explicitly included in the
quantum-mechanical states of the system. We report the renormalization of the
nucleon mass as a function of the momentum cutoff, an Euclidean time density
correlation function that deals with the short-time nucleon diffusion, and the
pion cloud density and momentum distributions. In the two-nucleon sector we
show that the interaction of two static nucleons at large distances reduces to
the one-pion exchange potential, and we fit the low-energy constants of the
contact interactions to reproduce the binding energy of the deuteron and two
neutrons in finite volumes. We show that the method can be readily applied to
light-nuclei.
