We present a microscopic theory of the second order phase transition in an interacting Bose gas that allows one to describe formation of an ordered condensate phase from a disordered phase across an entire critical region continuously. We derive the exact fundamental equations for a condensate wave function and the Green's functions, which are valid both inside and outside the critical region. They are reduced to the usual Gross-Pitaevskii and Beliaev-Popov equations in a low-temperature limit outside the critical region. The theory is readily extendable to other phase transitions, in particular, in the physics of condensed matter and quantum fields.Published The problem of finding a microscopic theory of the second order phase transitions became one of the major problems in theoretical physics in 1950s after an invention of the powerful Green's functions and quantum field theory methods in the many-body physics. These methods include the Feynman's and Matsubara's diagram techniques, Wick's theorem, and Dyson's equation (see [1][2][3][4][5][6][7][8][9][10][11][12][13][14] and references therein). A major point is a failure of a Landau self-consistent field theory in a critical region, shown by the exact Onsager's solutions. However, all these standard methods turn out to be insufficient to solve the problem: The microscopic theory, which should connect the asymptotics of the ordered and disordered phases across the critical region, has not been found so far even for anyone of the numerous phase transitions.Here we present such a microscopic theory for a BoseEinstein condensation (BEC) in an interacting gas. We derive the fundamental equations for the order parameter and Green's functions, which are valid not only in the fully ordered, low-temperature phase, but also in the entire critical region and in the disordered phase. They resolve a fine structure of the system's statistics and thermodynamics near a critical λ-point. It becomes possible due to the newly developed methods of (a) the nonpolynomial averages and contraction superoperators [15,16], (b) the partial difference (recurrence) equations [17][18][19] (a discrete analog of the partial differential equations) for superoperators, and (c) a characteristic function and cumulant analysis for a joint distribution of the noncommutative observables. They allow us to take into account (I) the constraints in a many-body Hilbert space, which are the integrals of motion prescribed by a broken symmetry in virtue of a Noether's theorem, and constraintcutoff mechanism, responsible for the very existence of a phase transition and its nonanalytical features, [4,20,21] (II) an insufficiency of a grand-canonical-ensemble approximation, which is incorrect in the critical region [2,8] because of averaging over the systems with different numbers of particles, both below and above the critical point, i.e., over the condensed and noncondensed systems at the same time, that implies an error on the order of 100% for any critical function, (III) a necessity to solve the problem ...