  A precise modelling of the dynamics of bubbles nucleated during first-order
phase transitions in the early Universe is pivotal for a quantitative
determination of various cosmic relics, including the stochastic background of
gravitational waves. The equation of motion of the bubble front is affected by
the out-of-equilibrium distributions of particle species in the plasma which,
in turn, are described by the corresponding Boltzmann equations. In this work
we provide a solution to these equations by thoroughly incorporating the
non-linearities arising from the population factors. Moreover, our methodology
relies on a spectral decomposition that leverages the rotational properties of
the collision integral within the Boltzmann equations. This novel approach
allows for an efficient and robust computation of both the bubble speed and
profile. We also refine our analysis by including the contributions from the
electroweak gauge bosons. We find that their impact is dominated by the
infrared modes and proves to be non-negligible, contrary to the naive
expectations.
