The Early Universe
This is more a review material than a text from which beginners can learn about the history of the early universe. Complete knowledge of involved mathematics (primarily differential geometry) and fundamental physics (general relativity, particle physics, statistical mechanics and thermodynamics) is assumed.
Overview

Construction of the FRW metric.
 From observational data it can be assumed that the observed universe is spatially homogeneous and isotropic. From this, a primitive version of FreedmanRobertsonWalker metric, that describes a homogeneous and isotropic universe in general relativity, is constructed,
 From observational data it can be more or less deduced that the universe is expanding, so that the proper time of the cosmic rest frame  the frame in which the CMBR is at rest  needs to be taken as a parameter that scales the metric, leading to the FRW metric.

Standard cosmology. The Einstein equation, with FRW metric and the stressenergy tensor of a perfect fluid, as a tensor equation, has only two independent field equations. The equations usually taken are the socalled Friedmann equation ($0$$0$ component, "involving only time"), and the first law of thermodynamics (the $\mu=0$ component of the conservation of stress energy).
 From the first law of thermodynamics it can be inferred that the early universe was radiation dominated, and then matter dominated, etc.
 From the Friedman equation may be integrated to give the age of the universe in terms of present cosmological parameters.

The $\Lambda$CDM model. A particular model in the Standard cosmology, or the standard model of Big Bang cosmology. Assembled from the observational data, incorporating the accelerating expansion of the Universe, largescale structure in the galaxy distribution, the existence and structure of the CMB, etc.
 Using the machinery in the Standard cosmology, calculate various aspects of the expansion history.
 Nonthermal aspects, such as the time when the deccelerationacceleration and radiation dominationmatter domination transition occured, can be approximated.

From the $\Lambda$CDM model and the Standard Model of particle physics, the Thermal History of the universe can be sketched.
 Plasma of massless thermal particles.
 First spontaneous symmetric breaking (SSB) phase transition: from GUT to standard model ($\sim 10^{16}$GeV), the strong interaction becomes distinct from the electroweak interaction.
 Electroweak symmetry breaking, $\sim 150$GeV, particles acquire mass via Higgs mechanism. Later particles which only interact via massive gauge bosons will decouple from the thermal plasma (but not yet, at around $150$GeV).
 Chiral symmetry breaking and color confinement $\sim 150$MeV. Quarks coalse into hadrons. Local inhomogeneities in the baryon number density
 Neutrino decoupling. $\sim 1$MeV.
 Big Bang nucleosynthesis. $1001$keV.
 Recombination. $\sim 0.4$eV. The origin of CMBR.
Preliminary Thermodynamics and Statistical Mechanics
On the meaning of $\text{eV}$ as temperature
This is a standard practice when dealing with plasma. Via $$ E = k_B T $$ where $k_B$ is the Boltzmann constant, energy is directly linked to temperature. This is because for massless particles (radiations), when $T$ is large, the average energy per particle $E$ is indeed $k_B T$.The probability that a massless particle has energy $E_n = nh\nu$ is given by the Boltzmann factor $$ P(n) = \frac{\exp(E_n/k_B T)}{\sum_{n=0}^\infty \exp(E_n/k_B T)}, $$ the mean energy of a particle with frequency $\nu$ is then $$ \overline{E_\nu} = \frac{h\nu }{e^{h\nu/k_B T} 1}. $$ For large $T$, $e^{h\nu/k_B T} 1 = \tfrac{h\nu}{k_B T}$, and so $$ \overline{E} = k_B T. $$
Even if one is not dealing with massless particles, since for relativistic particle gas $\overline{E} = n k_B T$ where $n$ is the number of degrees of freedom, and $\overline{E} = \frac{n}{2} k_B T$ for nonrelativistic particle gas, the conversion is legitimate. By this we mean, the word 'temperature' really can be intuited to have the usual meaning of temperature.
Numerically $1\text{ GeV} = 1.1605\times 10^{13}\text{ K}$.
Gases of Free Particles in Expanding Universe
Consider the evolution of the effective temperature for free particles. This is the case when the particles decouple from other particles, e.g. after recombination. The effective temperature of massless bosons and fermions decreases in time according to $$ T_{\text{eff}}(t) \propto \frac{1}{a(t)}. $$ Even for massive particles, as long as $T_{\text{eff}} \gg m$ wehere $m$ is the particle mass, the above relation still holds.
While for nonrelativistic particles, $$ T_{\text{eff}}(t) \propto \frac{1}{a^2(t)}. $$
See Gorbunov & Rubakov for details.
Construction of the FRW metric
The materials here are mainly from Wald, General Relativity, 1984.
The observed universe is spatially homogeneous and isotropic,
 A spacetime is spatially homogeneous if there is a 1parameter family of spacelike hypersurfaces $\Sigma_t$ foliating the spacetime, s.t. for each $t$ and any points $p,q\in\Sigma_t$, there is an isometry of the spacetime metric $g_{ab}$ that takes $p$ into $q$.
 A spacetime is said to be spatially isotropic at each point if there is a congruence of timelike curves, whose tangents are denoted $u^a$, filing the space time, s.t. for any tangent vectors $s^a_1, s^a_2\in V_p$ that are orthogonal to $u^a$, there is an isometry of $g_{ab}$ that leaves $p$ and $u^a$ at $p$ fixed but rotates $s^a_1$ to $s^a_2$.
Isotropy implies that the Riemann tensor takes the form of $R_{abc}^{\quad d} = k \delta^c_{\ [a}\delta^d_{\ b]}$, where $k$ is the eigenvalue of the linear map $L:\Omega^2\to \Omega^2$ given by raising the third index of the Riemann tensor (since distinct eigenvalues will destroy the isotropy), and homogeneity implies that the space is of constant curvature, viz. $k$ is invariant in a leaf $\Sigma_t$ of the foliation. From the Bianchi identity of the 4thindexlowered Riemann tensor, it can be seen that actually isotropy forces $k$ to be constant, so the assumption of homogeneity can be dispensed with.
A theorem of Eisenhart says that any two spaces of constant curvature of the same dimension and metric signature which have equal values of $k$ are locally isometric. Thus one needs to classify the possible spatial geometries of $\Sigma_t$ by classifying spaces of constant curvature with different values of $k$. Now $k$ can be positive, negative or zero,
 Geometries with $k>0$ are attained by the 3spheres.
 Geometries with $k=0$ are attained by ordinary 3dimensional flat space.
 Geometries with $k<0$ are attained by hyperboloids, i.e. surfaces in a 4dimensional flat Lorentz signature space with global inertial coordinates $t^2  x^2  y^2  z^2 = R^2$.
The spacetime metric, from the viewpoint of isotropic observables, should then be $$ g_{ab} = u_a u_b  h_{ab}(t) $$ where $h_{ab}(t)$ is the metric of a 3sphere, 3Euclidean space, or a hyperboloid and $t$ is the parameter of the foliation. Label each leaf of the foliation by the proper time $\tau$ of the isotropic observers, the spacetime metric should be $$ ds^2= d\tau^2  a(\tau)^2\left( \frac{dr^2}{1kr^2}+r^2d\theta^2 + r^2\sin^2\theta d\phi^2\right), $$ this is the FLRW metric.
In the following, the parameter $t$ will be dropped, and the proper time $\tau$ will be denoted by $t$ instead. Take $\eta$ to be the conformal time, i.e. $d\eta = dt/a(t)$, then the FRW metric can be written as $$ ds^2= a^2(\eta)\left( d\eta ^2  \frac{dr^2}{1kr^2} r^2d\theta^2  r^2\sin^2\theta d\phi^2\right). $$ Specifically, when $k=0$, $g_{ab} = a^2(\eta)\eta_{ab}$ where $\eta_{ab}$ is the Minkowski metric.
Standard Cosmology
Equations of Cosmology
To solve the Einsein equation $$ G_{ab} = R_{ab} \frac{1}{2}\mathcal{R} g_{ab} = 8\pi T_{ab} $$ one needs, in addition to a metric, the matter content of the universe that is described in terms of its stressenergy tensor $T_{ab}$. This we will take to be the general perfect fluid, $$ T_{ab} = \rho u_a u_b + P(g_{ab}+ u_a u_b). $$
The Einstein equation has only two independent equations, since by isotropy $G^{ab}u_b$ cannot have a spatial component. The independent components are $$ G_{ab}u^a u^b = 8\pi T_{ab}u^a u^b = 8\pi \rho, $$ and $$ G_{ab} s^a s^b = 8\pi T_{ab}s^a s^b = 8\pi P. $$ where $s^a$ is any unit vector tangent to the leaves, leading to
\begin{equation} 3\frac{\ddot{a}}{a} = 4\pi (\rho + 3P);\quad 3\frac{\dot{a}^2}{a^2}=8\pi\rho  \frac{3k}{a^2}. \label{eq} \end{equation}where the second one is the Friedmann equation. The first equation can also be substituted with the covariant conservation of the stressenergy $$ \nabla_a T^{ab} = 0 $$ which yields $$ \dot{\rho} + 3\frac{\dot{a}}{a}(\rho + P) = 0. $$
Supplemented with the equation of state of matter $$ P=P(\rho), $$ the Friedmann equation and the conservation of the stressenergy completely determine dynamics of the cosmological expansion. Summing up,
\begin{gather} \label{eq:cosmology} 3\frac{\dot{a}^2}{a^2}=8\pi\rho  \frac{3k}{a^2}\\\ \dot{\rho} + 3\frac{\dot{a}}{a}(\rho + P) = 0.\\\ P=P(\rho) \end{gather}are the equations of standard cosmology.
Hubble's law and the Redshift
An observations can be made. Provided that $\rho>0$ and $P\geq 0$, or more precisely $\rho+3P > 0$, the universe cannot be static (from the first equation above in $(\ref{eq})$, $\ddot{a}<0$). Either $\dot{a}>0$ or $\dot{a}<0$ with the possible exception of an instant of time when expansion changes over to contraction.
Between two isotropic observers in a homogeneous surface, at proper time $t$, if the distance is $D$, the rate of change of $D$ is $$ \frac{d D}{dt} = \frac{D}{a}\frac{da}{dt} = \frac{\dot{a}}{a} D = DH(t). $$ This is Hubble's law, and $H(t) = \frac{\dot{a}}{a}$. Here $\frac{dD}{da} = \frac{D}{a}$ since $D$ is linear in $a$. $H_0$ would denote the present observed value of $H(t)$.
Since the universe is expanding, $\dot{a}>0$, $\ddot{a}<0$, the universe must have been expanding at a faster rate before. This means that, at a time less than $H_0^{1}$ ago, the universe was in a singular state. The Hubble constant is related to the redshift, $$ z(t) = \frac{a_0}{a(t)}  1. $$ Expanding $a(t)$ around $a_0 = a(t_0)$, to the linear order in $(t_0 t)$, one has $$ z(t) = H_0\cdot(t_0t). $$ and upp to corrections of second order the distance $r=t_0 t$, so $$ z=H_0 r. $$
The $\Lambda$CDM Model, Accelerated Expansion, Age of the Universe
In the Friedmann equation (the first of $\ref{eq:cosmology}$), the curvature term, if any, is small, so we work with spatially flat model, putting$k=0$. Solving the Fiedmann equation requires certain specification of $\rho$, the energy density. As was stated before, the second of the cosmological equations ($\ref{eq:cosmology}$) really needs to be supplemented with the equation of state of the matter, which is the relation between the density and the pressure, $P = P(\rho)$. Now, case by case,
 For nonrelativistic matter (dust), $P = 0$; this leads to $\rho = \operatorname{const}/a^3$. This just says that $\rho$ would be diluted $\propto a^{3}$. Solving the Friedmann equation, $a(t) \propto (tt_s)^{2/3}$, nonsingularity of $\rho$ forces the arbitrary constant $t_s=0$. By $H(t)=\dot{a}/a = 2/(3t)$ and taking the present Hubble parameter $H_0$, $t_0 = 2/(3H_0)$.
 For relativistic matter (radiation), $P = \rho/3$; leading to $\rho = \operatorname{const}/a^4$. The difference between radiation and matter comes, in essense, from the fact that not only the number density is diluted according to $\propto a^{3}$, but the redshift of the energy is $\propto a^{1}$. Solving the Friedmann equation, $a(t) \propto t^{1/2}$, and $H = 1/2t$.
 For vacuum. The vacuum is the same in all inertial frames, thus Lorentz invariance dictates that the stressenergy is of the form $T_{\mu\nu} = \rho \eta_{ij}$. Now $T_{00}=\rho$ and $T_{ij} = P\eta_{ij}$, thus $P=\rho$; the negative sign ensures that the energy density is positive. Thus $\rho = \rho_{vac}$ is constant. Solving the Friedmann equation, $a \propto \exp(\sqrt{8\pi G \rho_{vac}/3}t) := \exp(H_{dS}t)$. $dS$ stands for de Sitter since the de Sitter metric is $$ ds^2 = dt^2 \exp(2H_{dS}t) dx^2. $$
For generic equation of state $P = w\rho$ where $\rho> 1$ is a constant, $\rho = a^{3(1+w)}$, and $$ a \propto t^{\frac{2}{3}\frac{1}{1+w}} := t^{\alpha} $$ The second derivative $$ \ddot{a} \propto \alpha(\alpha1) t^{\alpha 2} $$ implies that for $w>\frac{1}{3}$ the expansion deccelerates, and for $w<\frac{1}{3}$ (dominated by vacuum) the expansion accelerates.
The $\Lambda$CDM Model
Realistic models should, of course, account for nonvanishing spatial curvatures, and combine $\rho_M$ (matter), $\rho_{rad}$ (radiation), $\rho_\Lambda$ (dark energy) and the effective "density" of spatial curvature $\rho_{k}$ given by $\frac{8\pi}{3}G\rho_{k} =  \frac{k}{a^2}$. It should also be noted that $\rho_i$ changes through time: relativistic particles becoming nonrelativistic by cooling, particle creation, etc. In particular, the Universe should have been radiation dominated in its early age, at least in the Hot Big Bang picture, but now it is dark energy dominated.
In the present one can find the numerical value of $\rho_c = \frac{3}{8\pi G} H_0^2$, the critical density, for spatially flat Universe. Should $\rho_c$ equal $\rho_{M,0} + \rho_{rad,0} + \rho_{\Lambda,0}$, then the Universe is actually spatially flat. The current bound on $\rho_{k,0}/\rho_c$ is less than $0.01$, so the Universe is at least approximately flat nowadays. Specifying $\Omega_i := \rho_i / \rho_c$ obtained from observations, one gets the $\Lambda$CDM model.
The present values of $\Omega_i$ are roughly $\Omega_{rad} \lesssim 10^{4}$, $\Omega_{curv} < 0.01$, $\Omega_\Lambda = 0.685$, $\Omega_M = \Omega_B + \Omega_{DM} = 0.05 + 0.265$. $\Omega_{DM}$ is the contribution of dark matter and $\Omega_B$ the baryonic matters. The present universe is dark energy and dark matter dominated. We see that relativistic matter and curvature can all be neglected for the present time in reasonable calculations.
We don't know what dark energy is, and in the framework of $\Lambda$CDM one assumes that $\rho_\Lambda$ is independent of time, since $\Lambda$ is really the cosmological constant that is usually thought of as vacuum energy density. The Friedmann equation in the $\Lambda$CDM model is then
\begin{equation} \label{eq:LCDM} \big(\frac{\dot{a}}{a}\big)^2 = \frac{8\pi}{3}G\rho_c\big[\Omega_M (\frac{a_0}{a})^3 + \Omega_{rad}(\frac{a_0}{a})^4 + \Omega_\Lambda + \Omega_{curv}(\frac{a_0}{a})^2\big]. \end{equation}Age of the Universe
The age of the Universe is really a problem of solving the Friedmann equation under reasonable assumptions about the composition of the energy density. Since relativistic matter is relevant at early times only and spatial curvature is neglectable, the form of the Friedmann equation one employs when approximating the age of the Universe is $$ \big( \frac{\dot{a}}{a} \big)^2 = H^2_0 \big[ \Omega_M (\frac{a_0}{a})^3 + \Omega_\Lambda\big] $$ under the constraint $\Omega_M + \Omega_\Lambda = 1$ with positive $\Omega_\Lambda$. The solution is given by $$ a(t) = a_0 \big(\frac{\Omega_M}{\Omega_\Lambda})^{1/3}\big) \Big[\sinh\big(\frac{3}{2}\sqrt{\Omega_\Lambda}H_0 t\big)\Big]^{2/3}. $$ Then for $a(t_0) = a_0$, $$ t_0 = \frac{2}{3\sqrt{\Omega_\Lambda}}\frac{1}{H_0}\operatorname{Arsh} \sqrt{\frac{\Omega_\Lambda}{\Omega_M}}, $$ and for a redshift $z(t)$ (recall $z(t) = \frac{a_0}{a(t)} 1$) $$ t(z) = \frac{2}{3\sqrt{\Omega_\Lambda}}\frac{1}{H_0}\operatorname{Arsh} \sqrt{\frac{\Omega_\Lambda}{\Omega_M(1+z)^3}}. $$
Transitions in the history of Universe
From decceleration to acceleration
Computing $\ddot{a}$ in the $\Lambda$CDM model via $\ref{eq:LCDM}$, one immediately sees that in the present universe $\ddot{a}>$, the universe is undergoing accelerated expansion. Since the universe is assumed to be raditation/matter dominated in the past, there must have been a particular moment $t_{ac}$ at which $\ddot{a}=0$ occured. Ignoring relativistic matter and curvature, let $a(t_{ac}) = a_{ac}$, it is easy to see that $$ \big(\frac{a_0}{a_{ac}}\big)^3 = \frac{2\Omega_\Lambda}{\Omega_M}, $$ corresponding to the redshift $z_{ac} = (2\Omega_\Lambda/\Omega_M)^{1/3}1$, numerically $z_{ac} \approx 0.63$.
From radiation domination to matter domination
A crude estimation can be given by neglecting dark energy and curvature but use the current value for $\Omega_M$ and $\Omega_{rad}$, repeat the above procedure, one gets the redshift $z_{eq}$ of the moment when radiation and matter contributions are equal $$ z_{eq} = \frac{a_0}{a_{eq}}  1 \sim 10^4, $$ the temperature at that time would be approximately $T_{eq} = T_0(1+z_{eq}) \sim 1 eV \sim 10^4 K$. However at $1 eV$ temperature, all three species of neutrinos are relativistic, and thus must be accounted for. The details are unimportant (see Gorbunov, Rubakov, 4.4), but the estimate given is 52 thousand yrs after Big Bang.
Thermal History (First three minutes)
The reference for this section, until ElectronPositron Annihilation, is Daniel Baumann's lecture notes on Cosmology, URL.
Overview: Equilibrium and Freezeout
The early Universe should be filled with field excitations ("primordial plasma") that interact with each other. If the interaction rate, denoted $\Gamma$, is much larger thant the rate of expansion $H$, i.e. $\Gamma \gg H$, then typical time scaale of interaction is much smaller than the typical expansion time, and the particles can be seen as in thermal equilibrium with each other; they are in local thermal equilibrium.
The interaction rate $\Gamma \equiv n\sigma v$, where $n$ is the number density of particles, $\sigma$ their interaction cross section, and $v$ the average veolicity of the particle (for high temperature, i.e. relativistic particles, $v\sim 1$). For high temperature, dimensional analysis gives $n\sim T^3$, since particle masses can be ignored for ultrarelativistic particles. $\sigma \sim \frac{\alpha^2}{T^2}$, also by dimensional analysis, where $\alpha \equiv g^2_A /2\pi$ for a gauge boson $A$ with coupling $g_A$. Thus for ultrarelativistic particles, $\Gamma \sim \alpha^2 T$.
By the Friedman equation, $H\sim \sqrt{\rho}$, taking dimensions into account $H\sim\sqrt{\rho}/M_{pl}$. $\rho\sim T^4$, and hence $H\sim T^2/M_{pl}$. Thus $$ \frac{\Gamma}{H} \sim \frac{10^{16} \text{GeV}}{T}, $$ where $\alpha\sim 0.01$ is used. When $T$ ranges from $100\text{ GeV}$ to $10^{16}\text{ GeV}$, i.e. for all particles in the Standard Model when they are relativistic the local equilibrium condition is satisfied.
The temperature decreases and relativistic particles become nonrelativistic nevertheless. However, relativistic particles still dominate the density and pressure, and nonrelativistic particles can be ignored from the primordial plasma. The distribution for particles in thermal equilibrium $$ f(p) = \frac{1}{e^{(E(p)\mu)/T}\pm 1} $$ where $p$ is the momentum, implies that when $T\ll m$, the nonrelativistic distribution is supressed by $e^{m/T}$, thus relativistic particles dominate the density and pressure of the primordial plasma. It is sufficient to consider relativistic particles for the approximation of total energy density. Let $g_\ast = g_\ast(T)$ denote the effective number of relativistic degrees of freedom at temperature $T$, then $$ \rho_r = \sum_i \int d^3 p f_i(p) E_i(p) = \frac{\pi^2}{30}g_\ast(T) T^4. $$ In high temperature, all the species of partiles are relativistic, and $$ g_\ast = g_b + \frac{7}{8} g_f = 106.75 $$ the factor $7/8$ comes from integration.
The presence of massive particle species implies that the equilibrium didn't persist, since otherwize massive particles would have been exponentially suppressed and the universe would be mostly photons (at least for nonbaryonic matters, like neutrinos, the suppression would be present). The massive particles that are still abundant today deviated from equilibrium, decoupled from the primordial plasma, and freezed out, leaving relic particles.
From Electroweak Symmetry Breaking to QCD Phase Transition
The suppression of number and mass density of particles when $T \ll m$ can be interpreted as particle antiparticle annihilation; at low temperatures the thermal particle energies aren't sufficient for pair production. Around $80\%$ of the particle antiparticle annihilation takes place in the interval $T = m \to m/6$.
The heaviest particles annihilates first, reducing $g_\ast$.
 The top quarks, with mass aroung $160\text{ GeV}$, annihilates, reducing $g_\ast$ to $106.75\frac{7}{8} \times 12 = 96.25$.
 Then the Higgs boson and the gauge boson of weak interaction annihilates, roughly at the same time. At $T\sim 10\text{ GeV}$ we have $g_\ast = 96.25  (1 + 3\cdot 3) = 86.25$.
 Next annihilate the bottom quarks, leaving $g_\ast = 86.25  \frac{7}{8}\cdot 12 = 75.75$.
 Then the charm quarks and the tau leptons, $g_\ast = 75.75  \frac{7}{8}\times(12+4) = 61.75$.
 Next should have been the strange quarks, if there is no QCD phase transition (or crossover) at around $150\text{ MeV}$, during which the quarks combine into baryons and mesons. These hadrons are nonrelativistic below the temperature of the QCD phase transition, except the pions. Now left are pions, electrons, muons, neutrinos and photons, leaving $g_\ast = 2 + 3 + \frac{7}{8}\times (4+4+6) = 17.25$.
Then should annihilate electrons and positrons, but before that happens, neutrino decouples from the thermal bath.
Neutrino Decoupling
After the Electroweak Symmetry Breaking which occured at around $160\text{ GeV}$, the gauge bosons of weak interactions acquired mass, and the cross section of weak interactions became $\sigma \sim G^2_F T^2$ where the Fermi constant $G_F \sim \alpha/ M_W^2$. Now $$ \Gamma/H \sim \frac{\alpha^2 M_{pl} T^3}{M^4_W} \sim (\frac{T}{1\text{ MeV}})^3, $$ meaning that at around $1\text{ MeV}$, the particles that interact with the primordial plasma only through the weak interaction decouple and freeze out. These are the neutrinos. They are coupled to the primordial plasma via weak interaction processes such as
\begin{gather*} \nu_e + \bar{\nu}_e \leftrightarrow e^+ + e^,\\ e^ + \bar{\nu}_e \leftrightarrow e^ + \bar{\nu}_e. \end{gather*}Electroweak Phase Transition: around 160 GeV
When the temperature is around $160\text{ GeV}$, the electroweak phase transition occurs. This is due to the breaking of the $SU(2)_L \times U(1)_Y$ symmetry by Higgs mechanism to $U(1)_{EM}$ (since the vacuum is not invariant under $SU(2)_L \times U(1)_Y$), from which the Higgs boson acquires vacuum expectation value, the particles, including $W$ and $Z$ bosons, acquire mass, and the electromagnetic $U(1)$ interaction becomes distinct from the weak interaction mediated by $W$ and $Z$ bosons.
There is no general, wellestablished theory of finitetemperature Higgs scalar potential, so the precise dynamics of the electroweak symmetry breaking during the phase transition is still unknown.
Nevertheless, the Higgs field has no vacuum expectation value at high energies since the Higgs potential receives corrections from the thermal bath. When the temperature drops, Higgs field obtains a VEV and the electroweak symmetry breaks.
There are much to be desired in how the theory of finite temperature QFT is formed. It is not clear to me whether one can really add a thermal bath and do perturbation theory as usual, and some claim that perturbation theory is not applicable.
Finite temperature: Vacuum expectation value
The equilibrium state of the interacting medium corresponds to the minimum of the Landau potential (Grand potential). When the temperature $T$ is high ($\sim 1$ GeV), like so in the early Universe, chemical potentials are small, and one can use Helmholtz free energy $F$ instead.
At temperature $T$, let the expectation value of a scalar field $\phi$ be $\langle \phi \rangle_T$. In a system where the average value of $\phi$ is $\langle \phi\rangle_T$ everywhere and in thermal equilibrium otherwise, the free energy $F$ of the system would be $F=\Omega V_{eff}(T,\phi)$ where $\Omega$ is the spatial volume. The effective potential $V_{eff}(T,\phi)$ is the free energy density of the medium at temperature $T$ with the average scalar field $\phi$. In equilibrium, $F$ is at minimum, and hence $\langle \phi \rangle_T$ is the absolute minimum of the effective potential $V_{eff}(T,\phi)$, i.e. the vacuum expectation value (VEV).
At 0 temperature, $F$ is the same as the energy of the system, and $V_{eff}$ coincides with the scalar potential $V(\phi)$ in the Lagrangian of the field theory which $\phi$ takes part in. In the WeinbergSalam model, the absolute minimum of $V(\phi)$, $v = \langle \phi \rangle_{T=0} = \langle \phi \rangle$, is nonzero. The computation of $v$ is given in the next section. At finite temperature $T$, $V_{eff}(T,\phi)\neq V(\phi)$, and hence the absoluate minimum, $\langle \phi \rangle_{T}$ can take other values.
Electroweak Symmetry Breaking (0 temperature)
The WeinbergSalam model is an $SU(2)_L \times U(1)_Y$ gauge theory with tree $SU(2)_L$ gauge bosons $W^i_mu$, $(i=1,2,3)$ and one $U(1)_Y$ gauge boson $B_\mu$. The kinetic Lagrangian density is $$ \mathcal{L}_{K} = \tfrac{1}{4} W^i_{\mu\nu} W^{\mu\nu i} \tfrac{1}{4} B_{\mu\nu}B^{\mu\nu}, $$ where $$ W^i_{\mu\nu} = \partial_{[\nu}W^i_{\mu]} + g\epsilon^{ijk}W^j_\mu W^k_\nu;\quad B_{\mu\nu}=\partial_{[\nu}B_{\mu]}. $$
The most general renoramlizable $SU(2)_L$ invariant potential is $$ V(\Phi) = \mu^2 \Phi^\dagger \Phi+\lambda\Phi^\dagger\Phi^2, $$ with the complex scalar $SU(2)$ doublet $\Phi$ coupled to the gauge fields. If $\mu^2<0$, then the state of minimum energy is not at $\Phi=0$, leading to a vacuum expectation value of the scalar field. The direction of the minimum is not determined since the potential $V$ is ignorant of the direction, thus one can take the VEV to be $$ \langle\Phi\rangle = \frac{1}{\sqrt{2}}\begin{pmatrix} 0 \\\ v\end{pmatrix}, $$ leading to the $U(1)_Y$ charge $Y_\Phi =1$, and the electromagnetic charge $Q=\frac{\tau_3+Y}{2}$ ($\tau_3$ is the third Pauli matrix). Hence $Q\langle \Phi\rangle =0$, and the symmetry is broken from $SU(2)_L\times U(1)_Y$ to $U(1)_{EM}$.
The scalar doublet has the Lagrangian density of $$ \mathcal{L}_s = (D^\mu \Phi)^\dagger (D_\mu \Phi)  V(\Phi) $$ for covariant derivative $$ D_\mu = \partial_\mu + i \frac{g}{2}\tau\cdot W_\mu + i\frac{g'}{2}B_\mu Y $$ which gives physical gauge fields,
\begin{gather*} W^\pm_\mu = \tfrac{1}{\sqrt{2}}(W^1_\mu\mp i W^2_\mu)\\\ Z^\mu = \frac{g' B_\mu + g W^3_\mu}{\sqrt{g^2+g'^2}}\\\ A^\mu = \frac{gB_\mu + g' W^3_\mu}{\sqrt{g^2+g'^2}} \end{gather*}with masses $M^2_W = \tfrac{1}{4}g^2 v^2$, $M^2_Z = \frac{1}{4}(g^2 + g'^2)v^2$ and $M_A = 0$.
The parameter $v$ is found from the charged current of $\mu$decay, $$ \mu \to e\bar{\nu}_e \nu_\mu $$ and the interaction strength for muon decay is measured to be $G_F = 1.16639\times 10^{5}\text{ GeV}^{2}$. The momentum carried by the $W$ boson, of order $m_\mu$, can be neglected in comparison with $M_W$, hence $$ \frac{G_F}{\sqrt{2}}=\frac{g^2}{8 M^2_W} = \frac{1}{2v^2} $$ giving $ v = 246\text{ GeV}$, and the vacuum Higgs VEV $174\text{ GeV}$. It was measured that the Higgs mass $\mu = 125\text{ GeV}$.
The critical temperature
The finite temperature potential, with Higgs field $\phi$, has the form $$ V(\phi,T) = D(T^2  T^2_0)\phi^2  ET\phi^3 + \tfrac{\lambda_T}{4}\phi^4 $$ where $D,E,T^2_0,\lambda_T$ are some factors depending on particle masses, coupling constants, the Higgs VEV in zero temperature, and the temperature. Observe that there is only one minimum of the potential $V(\phi,T)$ sits at $\phi=0$ for large $T$. The critical temperature $T_c$ let $V(\phi,T)$ obtain two minima. Particles obtain masses via the electroweak symmetry breaking mechanism described in the last section. Calculations show that $T_c \approx 178\text{ GeV}$.
A point that needs to be stressed is that some claims that perturbation theory is actually not applicable. What is given here is only a partial and inaccurate but illlustrative picture. The contemporary value for the critical temperature is $T_c \approx 159.5$ GeV (so the finitetemperature perturbation theory actually gives a result that deviates farther w.r.t. the 'actual' value than 0temperature theory), see The Standard Model crossover on the lattice.