This paper explores the formation of primordial black holes (PBHs) and gravitational waves (GWs) from slow first-order phase transitions in the early universe. The authors focus on the scenario where a strongly supercooled first-order phase transition generates large inhomogeneities, leading to the formation of PBHs and GWs. They derive a new contribution to the GW spectrum, the secondary GWs, which arise from large perturbations induced by the transition. For sufficiently slow transitions ($\beta/H_0 < 12$), the secondary GWs dominate the spectrum and exhibit a distinct shape with two peaks, impacting the interpretation of recent pulsar timing array (PTA) data. The authors compute the distribution of perturbations using a semi-analytic approach, finding it to have negative non-Gaussianity. This non-Gaussianity suppresses PBH formation and results in a stronger GW signal compared to typical scenarios. They provide a simple fit to the abundance and mass function of the resulting PBH population, suggesting that models with a transition at $T_{\text{reh}} \approx 10^6 \, \text{GeV}$ and $\beta/H_0 \approx 8$ could explain the observed dark matter and produce a strong GW spectrum detectable by upcoming experiments like LISA, ET, AION, and AEDGE. The paper also discusses the implications for PTA data, showing that the inclusion of the secondary GWs significantly improves the fit to the NANOGrav data, leading to a higher reheating temperature and a different parameter space compared to previous analyses. The authors conclude that the peculiar shape of the GW spectrum, including two peaks, is a smoking gun signature of their scenario.This paper explores the formation of primordial black holes (PBHs) and gravitational waves (GWs) from slow first-order phase transitions in the early universe. The authors focus on the scenario where a strongly supercooled first-order phase transition generates large inhomogeneities, leading to the formation of PBHs and GWs. They derive a new contribution to the GW spectrum, the secondary GWs, which arise from large perturbations induced by the transition. For sufficiently slow transitions ($\beta/H_0 < 12$), the secondary GWs dominate the spectrum and exhibit a distinct shape with two peaks, impacting the interpretation of recent pulsar timing array (PTA) data. The authors compute the distribution of perturbations using a semi-analytic approach, finding it to have negative non-Gaussianity. This non-Gaussianity suppresses PBH formation and results in a stronger GW signal compared to typical scenarios. They provide a simple fit to the abundance and mass function of the resulting PBH population, suggesting that models with a transition at $T_{\text{reh}} \approx 10^6 \, \text{GeV}$ and $\beta/H_0 \approx 8$ could explain the observed dark matter and produce a strong GW spectrum detectable by upcoming experiments like LISA, ET, AION, and AEDGE. The paper also discusses the implications for PTA data, showing that the inclusion of the secondary GWs significantly improves the fit to the NANOGrav data, leading to a higher reheating temperature and a different parameter space compared to previous analyses. The authors conclude that the peculiar shape of the GW spectrum, including two peaks, is a smoking gun signature of their scenario.