Energy　flow　calculation　is　a　basic　work　for　the　planning　operation　and　optimal　scheduling　of　integrated　electricity-gas　system.　In　this　paper,　a　holomorphic　embedding　energy　flow　method　for　natural　gas　system　is　proposed　to　overcome　the　issues　of　strong　dependence　on　initial　values　and　difficulties　in　applying　to　looped　networks　in　natural　gas　energy　flow　algorithms.　On　this　basis,　a　holomorphic　embedded　multi-energy　flow　model　is　constructed,　and　the　recursive　relationships　of　power　series　coefficients　of　status　variables　of　power　system　and　gas　system　are　deduced.　The　energy　flow　distribution　of　the　integrated　electricity-gas　system　is　obtained　by　recursively　solution.　Simulation　results　verify　the　advantages　of　the　proposed　method　in　the　aspects　of　initial　value　dependence,　convergence　and　computational　efficiency.　©　2023　IEEE.