## 1 Introduction  
In 1924, [PERSON] and [PERSON] introduced the concept of Bose-Einstein statistics and also Bose-Einstein condensation, i.e., the macroscopic occupation of the lowest single-particle state of a system of bosons [1, 2]. In 1938, London suggested that the normal-superfluid phase transition of \({}^{4}\)He is related to the Bose-Einstein condensation and to the existence of a macroscopic wavefunction for the Bose condensate [3, 4]. In 1947, [PERSON] calculated, for a uniform weakly-interacting Bose gas, the quantum depletion, i.e., the fraction of bosons that are not in the Bose condensate at zero temperature due to a repulsive interaction strength [5]. In 1957, [PERSON], [PERSON] and [PERSON] evaluated the first correcting term to the mean-field equation of state of a uniform and weakly-interacting Bose gas [6]. In 1961, [PERSON] and [PERSON] derived the mean-field equation for the space-dependent macroscopic wavefunction of a weakly-interacting Bose gas in the presence of an external trapping potential [7, 8]. The Gross-Pitaevskii equation is the main tool used to describe the properties of the Bose-Einstein condensates, which are now routinely produced with ultracold and dilute alkali-metal atoms [9].  
Some years ago, experiments with atomic gases reported evidence of beyond-mean-field effects on the equation of state of repulsive bosons [10, 11]. These experimental results are quite well reproduced [12] by a modified Gross-Pitaevskii equation, which includes a beyond-mean-field correction that is the local version of the Lee-Huang-Yang term. A few years ago, [PERSON] suggested theoretically the existence of self-bound quantum droplets in an attractive Bose-Bose mixture, where the collapse is suppressed by a beyond-mean-field term [13]. Very recent experiments [14, 15] with two internal states of \({}^{39}\)K atoms in a three-dimensional configuration substantially confirm these theoretical predictions based on a modified Gross-Pitaevskii equation.  
Beyond-mean-field correcting terms into the Gross-Pitaevskii equation [16; 17; 18], or into similar nonlinear Schrodinger equations for superfluids [19; 20], are usually introduced heuristically in the spirit of the density functional theory. Here, we derive the modified Gross-Pitaevskii equation in a self-consistent way, starting from the [PERSON] equation of motion of the bosonic field operator \(\hat{\psi}(\mathbf{r},t)\) and the familiar [PERSON] prescription of writing the quantum field operator as the sum of a classical complex field \(\psi_{0}(\mathbf{r})\), which is the order parameter or macroscopic wavefunction of the Bose-Einstein condensate, and a quantum field \(\hat{\eta}(\mathbf{r},t)\), which takes into account quantum and thermal fluctuations [21; 22]. The presence of a generic external trapping potential \(U(\mathbf{r})\) is circumvented by adopting a semiclassical approximation for the Bogoliubov-de Gennes equations of the fluctuating quantum field [23; 24]. In this way, at zero temperature, we obtain the local density \(\hat{n}(\mathbf{r})\) of the out-of-condensate bosons as a function of the classical field \(\psi_{0}(\mathbf{r})\) and the corresponding equation for \(\psi_{0}(\mathbf{r})\), that is the stationary modified Gross-Pitaevskii equation with beyond-mean-field terms. From these terms, we recover the Lee-Huang-Yang correction [6] in the case of a uniform and real Bose-Einstein order parameter.