The Frehlich polaron model describes a ubiquitous class of problems concerned with understanding the properties of a single mobile particle interacting with a bosonic reservoir. Originally introduced in the context of electrons interacting with phonons in crystals, this model found applications in such diverse areas as strongly correlated electron systems, quantum information, and high energy physics. In the last few years this model has been applied to describe impurity atoms immersed in Bose-Einstein condensates of ultracold atoms. The tunability of microscopic parameters in ensembles of ultracold atoms and the rich experimental toolbox of atomic physics should allow to test many theoretical predictions and give us new insights into equilibrium and dynamical properties of polarons. In these lecture notes we provide an overview of common theoretical approaches that have been used to study BEC polarons, including Rayleigh-Schroedinger and Green's function perturbation theories, self-consistent Born approximation, mean-field approach, Feynman's variational path integral approach, Monte Carlo simulations, renormalization group calculations, and Gaussian variational ansatz. We focus on the renormalization group approach and provide details of analysis that have not been presented in earlier publications. We show that this method helps to resolve the striking discrepancy in polaron energies obtained using mean-field approximation and Monte Carlo simulations. We also discuss applications of this method to the calculation of the effective mass of BEC polarons. As one experimentally relevant example of a non-equililbrium problem we consider Bloch oscillations of Bose polarons and demonstrate that one should find a considerable deviations from the commonly accepted phenomenological Esaki-Tsu model. We review which parameter regimes of Bose polarons can be achieved in various atomic mixtures.