This paper investigates the formation of primordial black holes (PBHs) and gravitational waves (GWs) from slow first-order phase transitions (PTs) in the early universe. The study shows that such transitions generate large inhomogeneities, which can lead to the formation of PBHs and produce a GW spectrum with two components: one from bubble collisions and another from large perturbations. The latter dominates when the parameter β/H₀ is less than 12, significantly impacting the interpretation of recent pulsar timing array (PTA) data. The GW signal associated with a particular PBH population is stronger than in typical scenarios due to a negative non-Gaussianity of the perturbations, and it has a distinguishable shape with two peaks. The paper also explores the formation of inhomogeneities during a slow and strongly supercooled first-order PT. It computes the distribution of perturbations using a new semi-analytic approach and the PBH formation following the standard formalism. The results show that the distribution of perturbations has a negative non-Gaussianity, which suppresses PBH formation. This leads to a stronger GW background associated with a certain abundance of PBHs compared to scenarios with positive non-Gaussianity. The paper further discusses the formation of PBHs from large density perturbations and their mass function. It shows that the fraction of the total energy density that collapses to PBHs depends on the profile of the overdensity and the equation of state of the Universe. The paper also presents the PBH mass function and abundance, showing that models with a PT occurring at T_reh ≈ 10⁶ GeV and β/H₀ ≈ 8 can provide the entirety of the observed dark matter in the form of asteroid mass PBHs and produce a strong GW spectrum within the reach of LISA, ET, AION, and AEDGE. The paper concludes that the secondary GWs dominate the spectrum and have a significant impact on its shape close to the peak, provided the transition is slow enough. The study also shows that the negative non-Gaussianity of the perturbations is important for the interpretation of PTA data, as it disfavors the second-order GWs due to PBH overproduction. The paper highlights the importance of the secondary GWs in the context of current and future gravitational wave observations.This paper investigates the formation of primordial black holes (PBHs) and gravitational waves (GWs) from slow first-order phase transitions (PTs) in the early universe. The study shows that such transitions generate large inhomogeneities, which can lead to the formation of PBHs and produce a GW spectrum with two components: one from bubble collisions and another from large perturbations. The latter dominates when the parameter β/H₀ is less than 12, significantly impacting the interpretation of recent pulsar timing array (PTA) data. The GW signal associated with a particular PBH population is stronger than in typical scenarios due to a negative non-Gaussianity of the perturbations, and it has a distinguishable shape with two peaks. The paper also explores the formation of inhomogeneities during a slow and strongly supercooled first-order PT. It computes the distribution of perturbations using a new semi-analytic approach and the PBH formation following the standard formalism. The results show that the distribution of perturbations has a negative non-Gaussianity, which suppresses PBH formation. This leads to a stronger GW background associated with a certain abundance of PBHs compared to scenarios with positive non-Gaussianity. The paper further discusses the formation of PBHs from large density perturbations and their mass function. It shows that the fraction of the total energy density that collapses to PBHs depends on the profile of the overdensity and the equation of state of the Universe. The paper also presents the PBH mass function and abundance, showing that models with a PT occurring at T_reh ≈ 10⁶ GeV and β/H₀ ≈ 8 can provide the entirety of the observed dark matter in the form of asteroid mass PBHs and produce a strong GW spectrum within the reach of LISA, ET, AION, and AEDGE. The paper concludes that the secondary GWs dominate the spectrum and have a significant impact on its shape close to the peak, provided the transition is slow enough. The study also shows that the negative non-Gaussianity of the perturbations is important for the interpretation of PTA data, as it disfavors the second-order GWs due to PBH overproduction. The paper highlights the importance of the secondary GWs in the context of current and future gravitational wave observations.