The Laser Interferometer Space Antenna (LISA) is expected to detect thousands of individually resolved gravitational wave sources, overlapping in time and frequency, on top of unresolved astrophysical and/or primordial backgrounds. Disentangling resolved sources from backgrounds and extracting their parameters in a computationally intensive "global fit" is normally regarded as a necessary step toward reconstructing the properties of the underlying astrophysical populations. Here, we show that it is in principle feasible to infer the population properties of the most numerous of LISA sources-Galactic double white dwarfs- directly from the frequency (or, equivalently, time) strain series by adopting a simulation-based approach, without extracting and estimating the parameters of each single source. By training a normalizing flow on a custom-designed compression of simulated LISA frequency series from the Galactic double white dwarf population, we demonstrate how to infer the posterior distribution of population parameters (e.g., mass function, frequency, and spatial distributions). This allows for extracting information on the population parameters from both resolved and unresolved sources simultaneously and in a computationally efficient manner. This approach can be extended to other source classes (e.g., massive and stellar-mass black holes, extreme mass ratio inspirals) and to scenarios involving non-Gaussian or nonstationary noise (e.g., data gaps), provided that fast and accurate simulations are available.
