Aims. We study the properties of slow magneto-acoustic waves that are naturally excited due to turbulent convection and investigate their role in the energy balance of a plage region using three dimensional (3D) radiation-MHD simulations. Methods. To follow slow magneto-acoustic waves traveling along the magnetic field lines, we select 25 seed locations inside a strong magnetic element and track the associated magnetic field lines both in space and time. We calculate the longitudinal component (i.e. parallel to the field) of velocity at each grid point along the field line and compute the temporal power spectra at various heights above the mean solar surface. Additionally, horizontally averaged (over the whole domain) frequency power spectra for both longitudinal and vertical (i.e. the component perpendicular to the surface) components of velocity are calculated using time-series at fixed locations. To compare our results with the observations we degrade the simulation data with Gaussian kernels having FWHM of 100 km and 200 km, and calculate horizontally averaged power spectra for the vertical component of velocity using time-series at fixed locations. Results. The power spectra of the longitudinal component of velocity, averaged over 25 field lines in the core of a kG magnetic flux concentration, reveal that the dominant period of oscillations shifts from ∼ 6.5 minutes in the photosphere to ∼ 4 minutes in the chromosphere. This behavior is consistent with earlier studies restricted to vertically propagating waves. At the same time, the velocity power spectra, averaged horizontally over the whole domain, show that low frequency waves (∼ 6.5 minute period) may reach well into the chromosphere. In addition, the power spectra at high frequencies follow a power law with an exponent close to -5/3, suggestive of turbulent excitation. Importantly, waves with frequencies above 5 mHz propagating along different field lines are found to be out of phase with each other even within a single magnetic concentration. The horizontally averaged power spectra of the vertical component of velocity at various effective resolutions show that the observed acoustic wave energy fluxes are underestimated, by a factor of three even if determined from observations carried out at a high spatial resolution of 200 km. Since the waves propagate along the non-vertical field lines, measuring the velocity component along the lineof-sight, rather than along the field contributes significantly to this underestimate. Moreover, this underestimation of the energy flux indirectly indicates the importance of high-frequency waves that are shown to have a smaller spatial coherence and are thus more strongly influenced by the spatial averaging effect compared to low-frequency waves. Conclusions. Inside a plage region, on average a significant fraction of low frequency waves leak into the chromosphere due to inclined magnetic field lines. Our results show that longitudinal waves carry (just) sufficient energy to heat the chromosphere in solar plage. Howeve...