Accelerating analysis of Boltzmann equations using Gaussian mixture models: Application to quantum Bose-Fermi mixtures

.


I. INTRODUCTION
The Boltzmann equation is ubiquitous in modern physics as it can model the non-equilibrium dynamics and near-equilibrium collective behavior of a wide range of quantum many-body systems.Phenomena such as hydrodynamics and turbulence in quantum liquids as well as transport properties in metals and semiconductors fall under its purview.Valid for systems where the length scales characterizing quasiparticle interactions are shorter than typical distances quasiparticles travel before colliding [1], the Boltzmann equation simplifies the typically challenging analysis of such interacting manyparticle systems by accurately capturing the many-body dynamics with the evolution of single-particle distribution functions.When studying the near-equilibrium dynamics of a system, this approach gives efficient access to its linear-response properties which characterize the possible quantum phases of matter.Solving the Boltzmann equation, even after it is linearized around equilibrium, however, can often be computationally challenging and precludes analysis of a system.Approaches to the problem must contend with non-Gaussian equilibrium distribution functions encoding quantum statistics as well as quantum correlations in the self-energies such as those arising from particle-exchange processes [1,2].Even with modern computational resources, state-of-the-art techniques grant limited insight into scientifically interesting * Correspondence to: p dolgirev@g.harvard.eduphenomenona such as spin transport in magnetic materials and electron hydrodynamics in 2D van der Waals (vdW) semiconductors.Existing solution methods also struggle to accurately compute transport properties in technologically relevant systems such as lithium-ion batteries [3] and photovoltaic cells [4].
State-of-the-art approaches make different levels of approximations to simplify the numerical task.The test particle method [5,6] assumes that the dynamics of distribution functions of a many-body system composed of around 10 23 particles can be accurately estimated by sampling the trajectories of a much smaller number of 10 4 -10 6 particles.This is typically a safe approximation, but simulation of the system can remain challenging as collisions in the system correlate the trajectories and prevent their parallel computation.The method easily admits details of experimental protocols, but the computational cost often limits analysis of linear response properties as only a few parameter choices can be evaluated.Other methods also solve the Boltzmann equation numerically after simplifying the collision integrals via various assumptions [7][8][9], but are limited in scope to where those assumptions are valid.The method of moments (MoM) [10,11] takes an alternative approach to direct numerical simulation of the dynamics by variationally approximating the response of the system in the frequency domain.The MoM naturally captures collective behavior and, at the cost of some analytical work, reduces the problem to matrix inversion of a semi-analytical solution.It can, however, require additional work to incorporate certain experimental details such as complicated drive protocols.The computational bottleneck of the method  1. Overview of the linearized Boltzmann equation and its applications.In this work, we develop a solution framework based on Gaussian mixture models (GMM) which enables efficient and accurate analysis of the Boltzmann equation.We demonstrate the method in the context of trapped Bose-Fermi mixtures of cold atoms -see Fig. 2.
comes from evaluating the matrix elements of the multidimensional collision integrals.Estimating these matrix elements via Monte Carlo sampling is often feasible for many 2D systems, albeit with difficulty, but fails for 3D systems.
In this work, we present two principal results.Our first result is the development of an efficient approach to analyze the linear response dynamics of many-body systems without resorting to commonly used relaxation time approximations.This method allows one to keep track of the exact collision integrals in Boltzmann equations and treat both 2D and 3D systems.Our second result leverages the above method to present an analysis of collective modes in quantum Bose-Fermi mixtures that are under active experimental investigation in ultracold atom systems [12].We discover a surprising richness in dynamical regimes that this system manifests as the interactions between particles are tuned.
The method we develop allows for efficient computation of collision integrals, which we combine with the MoM framework, thus removing the bottleneck to solving Boltzmann equations.We term our approach GMM-MoM as it centers around using a Gaussian mixture model (GMM) representation of the equilibrium distribution functions around which the Boltzmann equation is linearized.We demonstrate that this representation provides an accurate approximation of the distribution functions and, critical to efficient solution of the Boltzmann equation, enables semi-analytical computation of the multidimensional collision integrals.This latter step is possible as the MoM basis functions are chosen to be polynomials which, when combined with a Gaussian integral kernel, are analytically computable via Wick's theorem [1].The efficiency of our approach results from the insight that both Bose-Einstein and Fermi-Dirac equilibrium distribution functions can be well approximated with a GMM consisting of only a few Gaussians (see Sec. III for details).Our GMM-MoM approach enables efficient analysis of the linear response dynamics of a large range of systems that can be described using Boltzmann equations, thus allowing their characterization through quantities such as properties of collective modes and transport coefficients.We release our fitting codebase online [13] and include tables of Gaussian fits to the equilibrium distribution functions in a broad range of parameter regimes; these fits are generic and can be directly applied to investigate a broad class of quantum many-body systems.

A. Practical advantage and relevance of the GMM framework
Typical approaches to solving the (linearized) Boltzmann equation either resort to approximations that make the analysis tractable at the cost of oversimplification, or employ brute-force numerical methods that are computationally costly and provide only limited information.By rendering multidimensional collision integrals efficiently computable, our GMM approach enables both fast and accurate analysis of the Boltzmann equation, thereby allowing investigation of broad classes of systems that have hereto remained intractable -see Fig. 1.
A ubiquitous approach to simplifying Boltzmann equations is the relaxation time approximation (RTA), often used to analyse transport properties of materials.This approximation assumes that all quantities of interest in the system relax to thermal equilibrium on the same time scale.While the RTA can capture the behavior of relatively simple systems such as bulk crystalline silicon [14], it is often quantitatively inaccurate even in these contexts [15].The approximation also fails to qualitatively describe the physics of 2D semiconductors [16,17], which foster both fundamental [18][19][20][21] and technological interest [22].A variety of other important systems -including high-T c superconductors, resonantly interacting particles in neutron matter, and ultracold atoms at unitarityare also beyond the scope of RTA as these systems are    FIG. 2. Schematic of an atomic quantum Bose-Fermi mixture trapped into a cigar-shaped optical potential.A timedependent (transverse) perturbation of the trap launches collective modes of the system, which characterize the equilibrium state of matter.
expected to exhibit a clear separation of relaxation time scales for quantities such as heat, charge, spin, and momentum.The GMM framework introduced in this work obviates the need for the RTA while preserving computational efficiency.As a concrete demonstration of this capability, we benchmark the method on stronglyinteracting systems in Appendix C. There, we compute the decay rate of a Bose polaron in the unitary limit and the viscosity of a two-component strongly-interacting Fermi fluid.
Two additional simplifications are commonly applied to the Boltzmann equation: the assumption of isotropic scattering and the assumption of position-independent collisions.The former precludes, for instance, analysis of magnetic systems, including simple ferromagnets, as magnon collision processes explicitly depend on the momenta of the scattering particles [23].The latter is strictly applicable only to spatially homogeneous systems, which precludes systems such as trapped gases, including dipolar atoms, that manifest unique transport regimes of fundamental interest [24].The GMM framework is naturally suited to solve problems beyond these simplifications, as evidenced by our detailed analysis of trapped Bose-Fermi mixtures.
In contexts where the typical approximations fail, one must resort to solving the full linearized Boltzmann equation.One such context is material design.Ab initio methods like density functional theory are combined with Boltzmann analysis to design transport properties relevant to a wide range of technologies including thermoelectric materials [15], lithium-ion batteries [3], and photovoltaic cells [4].Accurate computation of transport coefficients often requires going beyond the RTA for most 2D materials [16] as well as 3D materials like diamond where three-phonon or higher processes are relevant [25].Anisotropic and position-dependent scattering must also be accounted for in devices with nontrivial geometries [25][26][27].
Remarkable progress has been made in brute force numerical methods in such systems; recent approaches have, for example, enabled solving the fully coupled electronphonon transport kinetic equations relevant to thermoelectric and thermopower applications [28,29].These state-of-the-art methods, however, are extremely timeconsuming, as typical calculations take several months of CPU time [29].We expect that the GMM framework introduced in this work can significantly ameliorate the cost of such computations.
As a concrete context, we focus on the case of Bose-Fermi mixtures of cold atomic gases trapped in a cigarshaped optical potential [46][47][48][49][50][51][52], illustrated in Fig. 2. We choose this system to demonstrate our GMM-MoM method for three reasons.Firstly, quantum gases of bosonic and fermionic atoms have emerged as a leading quantum simulation platform to study several quantum many-body systems, as cold-atom setups provide a highly controllable setting to investigate different prototypical models [53].Until recently, they have primarily been used to study purely bosonic or fermionic systems due to the experimental difficulty of simultaneously trapping both types of atoms.The collective behavior of such bosonic and fermionic gases has been explored in several experiments, including precise measurements of collective modes [54][55][56], observation of first and second sound modes [57,58], measurements of viscosities in the unitary regime [59,60] and transport properties in strongly-interacting gases [61][62][63][64][65][66][67].Systems in the polaronic regime, where one gas is much more dilute than the other, have been predicted to exhibit novel collective behavior [68] and non-equilibrium dynamics [69,70].Experiments have investigated the emergence of quasiparticles in such systems [71] and, more recently, have started to explore the collective behavior of mixtures of bosonic and fermionic cold atoms [12] beyond the polaronic regime, thus opening the door to studying Bose-Fermi mixtures in a controlled setting leveraging the entire coldatom toolbox [72][73][74].Secondly, on the theoretical side, Bose-Fermi mixtures beyond the zero-temperature impurity regime host potentially rich physics which has been unexplored in prior work due to the difficulty of numer-ical analysis and lack of experiment-guided motivation.Thirdly, the system is described by one of the most challenging types of Boltzmann equations to solve due to its spatial inhomogeneity, resulting from the asymmetric cigar-shaped optical trap, and the particle-conserving nature of collisions in the system.Making experimentally verifiable predictions in this context thus provides confidence in the GMM-MoM analysis method.
We analyze the collective behavior of the trapped Bose-Fermi cold-atom mixture across a variety of relative densities and interaction strengths.Motivated by modern cold-atom experiments, we study monopole and quadrupole collective modes, focusing on a quantum mixture above the Bose-Einstein transition temperature T > T c ; the gases are cold enough for particle statistics to matter and therefore form a quantum Bose-Fermi mixture.We find that this system exhibits a rich phenomenology which we summarize below for various regimes of relative densities of the atomic species.
The system exhibits a collisionless-tohydrodynamic crossover as the boson-boson interaction strength is tuned.Both the monopole and quadrupole mode lifetimes manifest a non-monotonic behavior, being the shortest in the crossover region.We also find frequency mixing between transverse and longitudinal monopole modes for sufficiently strong interactions; these modes are distinct from each other due to the asymmetry of the cigar-shaped trap.
(II) Impurity regime (dilute fermions).When the fermions are dilute, they act as impurities; the bosons remain largely unaffected and behave as the singlecomponent Bose fluid described above.The fermions and bosons have different monopole mode frequencies, and we can therefore understand the system in terms of coupled harmonic oscillators corresponding to the respective collective modes.When the boson-boson interaction is such that the boson fluid is not in its crossover regime and therefore manifests long-lived modes, we find typical Fano interference lineshapes encoding coherent mixing between the bosonic and fermionic oscillators (see Fig. 6 and discussion in Sec.IV C for details).Upon increasing the boson-fermion coupling, initially sharp fermionic resonances become completely featureless in frequency, indicating that the fermions no longer have independent dynamics and instead follow the bosonic atoms.In addition to the above phenomenona, we find that the fermionic spectral functions exhibit Stokes and anti-Stokes sidebands due to the mixing between longitudinal and transverse monopole modes.
(III) Mixture regime (comparable densities of bosons and fermions).The Bose-Fermi mixture exhibits a collisionless-to-hydrodynamic crossover as the bosonfermion interaction strength is tuned.The emergent hydrodynamic behavior of the whole many-body system manifests as mode locking between bosons and fermions.Intriguingly, we find coherent mixing between bosonic and fermionic collective modes even when the bosons are hydrodynamic and fermions are collisionless.Some of the phenomenology summarized above, such as the collisionless-to-hydrodynamic crossover of singlecomponent Bose fluids, has been discussed in previous papers [6,75].Other effects we uncover, however, such as the Fano lineshapes and the emergence of sidebands, are indicative of previously unexplored physics exhibited by Bose-Fermi mixtures.Additionally, many-body mode mixing and locking in the mixture regime has not been confirmed theoretically due to the collision integral computational bottleneck.The above results are therefore both a concrete window into the rich physics of Bose-Fermi cold-atom mixtures, as well as evidence of the power of the GMM-MoM method to analyze previously inaccessible many-body systems.

II. THEORETICAL FRAMEWORK
This section is organized as follows: In Sec.II A, we first introduce the coupled Boltzmann kinetic equations to describe the dynamics of a trapped Bose-Fermi mixture for T > T c .The equilibrium distribution functions are then discussed in Sec.II B. Finally, Sec.II C is devoted to the method of moments for computing the linear response properties of the mixture, including the spectrum of collective modes.Our framework closely follows and substantially extends that in Ref. [24] on twodimensional fermionic dipolar gases.We end this section by writing down expressions (34)- (38) for the 15dimensional collision integrals which represent the main challenge behind the method of moments in three spatial dimensions.These integrals, in general, are not feasible for the existing numerical tools such as direct Monte Carlo sampling that was proven efficient in two dimensions [24].One of our main technical results is that we overcome this challenge and make the method of moments numerically tractable -this will be the subject of the following sections.

A. Kinetic equations for T > Tc
We describe a trapped Bose-Fermi mixture using the following microscopic model (throughout the text, we set Here, the first term represents the bosonic Hamiltonian: The second term is the fermionic one: The two systems interact with each other through the third term, Ĥint , given by: Here ψ B/F (r) and ψ † B/F (r) encode the bosonic/fermionic annihilation and creation operators, respectively; they satisfy the canonical commutation (anti-commutation) relations.We assume that, as in cold-atom experiments, the contact coupling strengths g B and g BF can be tuned over a broad range using the physics of magnetic Feshbach resonances.For future reference, g B = 4πa B /m B and g BF = 2πa BF (m B + m F )/(m B m F ), where a B and a BF are the Bose-Bose and Bose-Fermi scattering lengths, respectively.In what follows, we consider the trapping potentials to have a cigar-like shape: To relate to the MIT experiment [12], we will consider the fermionic trapping potential to be U trap,F = λU trap,B (r), where λ ≃ 1 is an experimental parameter that encodes the fact that the laser coupling to fermions is a bit different from that to bosons.At this stage, we do not need to specify concrete forms of trapping potentials and can let the formalism be generic.We briefly remark that a possible microscopic derivation of Boltzmann equations uses the Keldysh technique of non-equilibrium Green's functions [76].This description is closely related to an appropriate conserving Kadanaff-Baym framework [77].Both these approaches are usually too complicated to work with, either analytically or numerically, and require approximations.In systems such as quantum Bose-Fermi mixtures considered here, a natural simplification occurs because there is a separation of time and length scales between fast microscopic degrees of freedom and macroscopic ones set by the trapping potential.This separation enables one to make a systematic gradient expansion, which dramatically simplifies the kinetic equations.Further progress is possible if one is allowed to employ the quasiparticle approximation, where the complex dynamics of two-time Green's functions is reduced to that of single-particle distribution functions, thus, obtaining the Boltzmann equation [78,79].The quantum statistics of the constituent particles manifests in both incoherent scatterings captured via collision integrals and in self-energies that might encode inherently quantum processes such as particle exchange.For more details on the derivation of Boltzmann equations, we refer the reader to Refs.[76][77][78][79].The primary goal of this paper is to develop a tool to solve Boltzmann equations rather than to provide a microscopic derivation of the equations of motion.To this end, below we consider a Bose-Fermi mixture above the BEC transition temperature T > T c , where the kinetic equations are well-known.
Specifically, they read: where n B (p, r, t) and n F (p, r, t) are the bosonic and fermionic distribution functions, respectively.For the left-hand sides, we employ the Hartree-Fock approximation [80]: where n B/F (r, t) = d 3 p (2π) 3 n B/F (p, r, t) represents the bosonic/fermionic real-space density.We note that the Hartree-Fock self-energies do not depend on p, which in turn implies that the effective masses are not affected by the contact interactions.For the right-hand sides, we use the Born-Markov approximation and write the Bose-Bose, Bose-Fermi, and Fermi-Bose collision integrals as: ))], (12) where ε B/F (p) = p 2 /(2m B/F ).Below we turn to analyse the linear-response properties of the Bose-Fermi mixture within this Boltzmann kinetic theory.Prior to that, we first compute the self-consistent equilibrium distribution functions needed for our subsequent analysis.

B. The equilibrium state
The equilibrium distribution functions are related to the effective potentials U eff,B/F (r) through (β = 1/T ): where z B/F (r) ≡ exp(β(µ B/F − U eff,B/F (r))) is the bosonic (fermionic) local fugacity and µ B/F is the corresponding chemical potential.Since the effective potentials U eff,B/F (r) depend on n B/F,eq , Eq. ( 13) should be solved self-consistently.This latter task is relatively simple because the momentum integrals are evaluated analytically: where ξ 3 2 (x) is the polylogarithm function.The chemical potentials can be fixed via: with N B/F being the total number of Bose/Fermi particles.
For convenience, throughout the rest of the paper we use the following dimensionless variables.As for the units of momentum and length, we choose p ρ = √ m B ω ρ and a ρ = 1/p ρ , respectively; we fix the unit of energy to be ω ρ ; the unit of mass is then naturally set by m B .In the dimensionless units, the trapping potentials read: where κ = ω 2 z /ω 2 ρ is the anisotropy parameter, and bars on top of the symbols indicate that the corresponding variable has been appropriately rescaled.For instance, the dimensionless interactions strengths read: ḡB = g B /(a 3 ρ ω ρ ) and ḡBF = g BF /(a 3 ρ ω ρ ).We finally note that because the two trapping potentials depend only on the combination ρ2 +κz 2 , the equilibrium self-consistency equations can be written as where r2 = β(ρ 2 +κz 2 ).We remark that these equations, which are solved numerically in practice, imply that the two fugacities depend only on the one-dimensional variable r rather than on two independent variables ρ and z.This observation is not crucial (and in principle, the trapping potential for fermions could be very different from the bosonic one) but facilitates our numerical calculations below.
C. Linear response within the method of moments Our primary goal now is to investigate linear response functions, related to collective modes, which represent small amplitude fluctuations on top of the equilibrium state.To this end, we write: (22) where ∆ B/F (r, p) = n B/F,eq (r, p)[1 ± n B/F,eq (r, p)].In what follows, we assume that Φ B (r, p, t) and Φ F (r, p, t) are small and, as such, linearize the equations of motion.In the remainder of the text, we consider only dimensionless variables that have been properly rescaled and omit writing bars on top of the corresponding symbols.
The method of moments is based on the idea of expanding small fluctuations Φ B/F (r, p, t) using a suitable basis set of functions {ϕ α (r, p)} N α=1 : Here N is the total number of basis functions that we take into account, and, as such, this parameter controls the accuracy of our approximations.Motivated by the experiments, below we investigate linear response to perturbations of the bosonic and/or fermionic trapping potentials, which we also write as: where the frequency amplitudes δU B/F,α (ω) are set by the actual experimental driving protocol.Plugging these expansions into the kinetic equations and subsequently projecting the linearized equations onto the basis functions {ϕ α (r, p)}, we arrive at: where we defined Φ = (Φ B,1 , . . ., Φ B,N , Φ F,1 , . . ., Φ F,N ) T and δU = (δU B,1 , . . ., δU B,N , δU F,1 , . . ., δU F,N ) T to be 2N -dimensional vectors.The matrices M and Ĥ are diagonal in the Bose-Fermi space: with the matrix elements given by: where Throughout the derivations, we used the identity {n B/F,eq , f } = −β∆ B/F {H B/F , f }.The matrix Ŝ encodes the changes in self-energies arising from the changes in the bosonic and fermionic distribution functions: where Here we defined Finally, the collision integral matrix Î is given by: where the corresponding matrix elements read: where The matrix I B describes changes in the bosonic distribution function due to boson-boson scatterings, J BF describes changes in the distribution function of bosons due boson-fermion scatterings, etc.
What we have achieved so far is that the initial infinitedimensional (linear) problem is now reduced to a finitedimensional one.To complete the discussion of dynamics, we also need to review an algorithm for evaluating time-dependent expectation values of operators, which can be measured in experiments.As a concrete example, we consider an arbitrary bosonic observable: where we assumed that ⟨O B ⟩ eq = 0, used Eq. ( 21), and expanded O B (r, p) = α O B,α ϕ α (r, p).More generally, we write: To evaluate Φ(ω From this, we finally get: where Equation ( 41) is among the most useful results, as it defines the linear response within the method-of-moments.Physically, the "residue" r a (ω) encodes the relative importance of the corresponding pole Ω a .We remark that the eigenvalues of Ω are not necessarily real (but for stability of the equilibrium state, their imaginary parts must be negative), which implies that the matrix V is generally not unitary.
The main challenge that we encounter in this section is the evaluation of the 15-dimensional collision integral matrix elements in Eqs. ( 34)- (38).We remark that the energy and momentum conservation laws render these integrals to be 11-dimensional.Further simplifications can occur if one uses symmetries of the optical trap.Nevertheless, the resulting expressions are still too complicated for accurate numerical evaluations.So far, this issue has been a major computational challenge, which is the primary reason why the method of moments has not been applied in its full capacity to three-dimensional problems.We note that in two spatial dimensions, such matrix elements encode only 5-dimensional integrals, which could be accurately evaluated using Monte Carlo sampling [24].
In this work, we provide an efficient method for computing these collision integrals.Our approach is based on the observation that the integrals are hard to evaluate because of the complexity of the equilibrium Bose-Einstein and Fermi-Dirac distribution functions.To further appreciate this point, we note that in a classical system, where the distribution functions have Gaussian (Boltzmann-Maxwell) forms, one could compute the collision integrals analytically, by means of Wick's theorem [75].With this in mind, in the following section, we offer an approximation to the equilibrium distribution functions that, on the one hand, can be made arbitrarily accurate and, on the other hand, also acquire a Gaussian-like structure, in turn enabling us to semianalytically evaluate the collision matrix with essentially arbitrary precision.

III. GAUSSIAN MIXTURE MODEL
In this section, we argue that locally, i.e., for a given real-space point r, the equilibrium distribution functions (both bosonic and fermionic) can be accurately approximated using the following ansatz: Here M 0 is the total number of Gaussians that we take into consideration, implying that we have 2 × M 0 fitting parameters a n (r) and γ n (r) that, in general, are allowed to depend on the real-space coordinate r.In the machine learning literature [81], the ansatz in Eq. ( 43) is referred to as Gaussian mixture model (GMM), though in contrast to the most generic GMM, here we set all the means to zero.In the following subsections, we first consider a few situations of increasing complexity, where we explicitly construct such GMM fittings and discuss their limitations.Then, in subsection III D, we briefly summarize and further discuss the presented approach.
A. Dilute fermions, zF ≤ 1 We begin by considering the Fermi-Dirac distribution function, which in the appropriately chosen units, can be written as: The local fugacity z = z(r) is a monotonic function of the density, and, for reasons that will become apparent shortly, we first focus on the case z ≤ 1, corresponding to dilute fermions.It is then suggestive to change variables as z = exp(−q 2 /2) and consider the following one-variable function: In other words, f 1 (x) is nothing but the Fermi-Dirac distribution for z = 1.If this function f 1 (x) can be accurately approximated with a GMM we immediately get a GMM approximation for n F (z, p) for any z ≤ 1: Importantly, the quality of this fit for n F (z, p) is entirely determined by the quality of the approximate form (46) for f 1 (x).Moreover, by comparing Eqs. ( 43) and ( 47), we observe that the amplitudes a n (r) in Eq. ( 43) do explicitly depend on the local fugacity z (and, as such, on the real-space coordinate r), whereas the variances γ n , which are the same for both Eq. ( 43) and Eq. ( 47), are universal, i.e., independent of z.
In our numerics, we optimize the coefficients α n and γ n using standard fitting routines that eventually minimize the absolute error The GMM fits provide excellent approximations to f2(x) in their respective domains of applicability (dashed line for zB = 0.9 and solid dashed line for zB = 0.99): The inset shows that the bosonic error LB, Eq. ( 52), can be as small as 10 −8 with only about 10 ( 17) Gaussians for zB = 0.9 (zB = 0.99).(c) The fermionic function f3(x) (solid line), Eq. ( 55), that captures the Fermi-Dirac distribution function for any 1 ≤ zF ≤ zmax, can be well fitted with the complex version of the GMM, Eq. ( 56).Here, the domain x ∈ [xmin, ∞) is determined by the maximal fermionic fugacity we consider: xmin = − ln(zmax).As shown in the inset, the absolute error of the fit LF decreases with increasing the total number M0 of complex conjugate Gaussian pairs; for M0 ≲ 25, we obtain LF ≃ 10 −8 at zmax = 10 5 .
Specifically, we use the Levenberg-Marquardt optimization algorithm, which we supply with gradients, that is implemented in Matlab via the lsqcurvefit function.The result of such fitting is shown in Fig. 3(a): Remarkably, we find that L F can be made as small as ∼ 10 −10 with only ∼ 10 Gaussians and can be further reduced by increasing the total number of Gaussians M 0 [see the inset of Fig. 3(a)].While f 1 (x) is expected to be well approximated with a GMM [82], what is a bit surprising is how few Gaussians turn out to be needed for such a high quality of fitting.
To further appreciate this point, we note that for z < 1, one can formally Taylor expand the fermionic distribution function: The fact that this expansion exists in the first place confirms the expectation that the GMM in Eq. ( 43) can, in principle, be made arbitrarily accurate.We note, however, that the results in Fig. 3(a) clearly demonstrate that the GMM significantly outperforms the truncated Taylor series, especially for z approaching unity, where more and more Gaussians in the Taylor series are needed to keep the same quality of the approximation.For instance, one needs to consider about 150 Taylor series coefficients for z = 0.9 to get the same precision as one has from only 10 GMM Gaussians.Importantly, the cost of computing the collision matrices scales as ∼ M 4 0 , which implies that the GMM dramatically facilitates otherwise too demanding calculations.We also remark that both the GMM in Eq. ( 47) and the Taylor series in Eq. ( 49) share the property that their variances γ n do not depend on z.
We briefly comment that in our numerics, perhaps the most non-trivial step turns out to be the initial choice of the fitting parameters for the optimization algorithm.Whenever available, one can start from the coefficients of the truncated Taylor series, which we find then lead to excellent GMM fits.In more complicated situations, where Taylor series is no longer an option, one can try to sequentially sample over random Gaussians with subsequent optimization.This approach gives fits of high precision, but it can be further improved by simple postprocessing that then gives a more efficient fit with fewer Gaussians, as we discuss below.

B. Bosons
We switch to the analysis of the Bose-Einstein distribution function: For bosons, the local fugacity has to be smaller than unity z ≤ 1 so that essentially all the discussion from the preceding subsection applies here as well.One notable difference, however, is that the bosonic one-variable function f 2 (x), now defined as diverges for x → 0. In particular, this implies that the function f 2 (x) cannot be fitted with a GMM (43) in the entire range x ∈ [0, ∞).In practice, this is not an issue because one can always choose a maximal bosonic fugacity z max that cuts off the range to x ∈ [x min , ∞), with x min = −2 ln (z max ).Physically, by choosing such a cutoff, one limits the analysis to temperatures T outside a small window near T c since the bosonic fugacity can approach unity only in the immediate vicinity of the BEC transition.For instance, based on mean-field theory, the choice of z max = 0.9 would correspond to temperatures |T −T c |/T c ≳ 5%; and z max = 0.99 implies that |T − T c |/T c ≳ 1%.One can also increase the cutoff z max closer to unity if one wants to investigate the detailed behavior near T c .This will come at the cost that the fitting near T c becomes more challenging and requires more Gaussians [see Fig. 3(b)].Figure 3(b) shows the GMM fits to the bosonic distribution n B (z, p): for z = 0.9 and z = 0.99, we can reduce the bosonic error to 10 −8 with as few as 10 and 17 Gaussians, respectively.If needed, this precision can be further improved by increasing M 0 -see the inset of Fig. 3(b).Once the GMM fit at z max is established we immediately get the same quality GMM approximation for any z ≤ z max : Notably, as Fig. 3(b) illustrates, the fitting at a given z max can be extrapolated (using Eqs. ( 53) and ( 54)) to even get a reasonably good approximation for z > z max .
C. Dense fermions, zF > 1 We turn to discuss the fermionic distribution function n F (z, p) in the regime z > 1, corresponding to dense fermions.Following the preceding analyses, now it is suggestive to change variables as z = exp(q 2 /2) and consider the function so that n F (z, p) = f 3 ((p 2 −q 2 )/2).Since, in principle, one can have arbitrary large fermionic fugacities, it implies that generically x ∈ (−∞, ∞).In practice, however, the range of the variable x is set by the maximal fermionic fugacity z max that we consider, i.e., x ∈ [x min , ∞) with x min = − ln(z max ).Physically, z max is essentially set by the fermionic density and local temperature.To be consistent with the preceding subsections and to allow for negative values of x, a GMM fit to the function f 3 (x) translates as a sum of exponentials (rather then a sum of Gaussians), which we write as: In contrast to what we have developed so far, here we now allow for the coefficients α n and γ n to be complex (and we only require Re γ 2 n > 0), as we find that to fit the function f 3 (x) with purely real parameters, assumed in Eq. ( 43), turns out to be numerically demanding, especially for very large fugacities z ≫ 1.This extended class of functions incorporates exponentially decaying trigonometric functions, which can fit the initial plateau of f 3 (x) substantially more accurately and with fewer resources than real exponentials.Importantly, the collision integrals in Eqs. ( 34)- (38) can still be expressed analytically within this new representation.In what follows, we will refer to the model in Eq. ( 56) as complex GMM (cGMM).We finally remark that to get an approximation for n F (z, p) for any 1 ≤ z ≤ z max , one can use Eq. ( 47) and Eqs. ( 53)-( 54), with the only difference that the corresponding coefficients are now complex.
Since the function f 3 (x) is real, we consider M 0 complex conjugate pairs so that we now have 4 × M 0 fitting parameters.Figure 3(c) summarizes the results of cGMM fits: we find that the absolute fermionic error L F , defined similarly as above, of the fit at z max = 10 5 can be reduced to ∼ 10 −8 with only about M 0 = 21 Gaussian pairs.As we show in the inset of Fig. 3(c), this error L F decreases with increasing M 0 .However, there are regions where L F plateaus at some values, meaning that the addition of new Gaussian pairs does not improve the quality of the fit, which could possibly be improved by using better initialization of the fitting parameters -we leave this issue to future work.Furthermore, this expectation is supported by the fact that a simple post-processing, where we remove as many redundant Gaussians as possible with subsequent re-optimization, gives a notably better result.specifically, we obtain the target precision L F ≤ 10 −7 at z max = 10 5 with only M 0 = 17 complex conjugate Gaussian pairs.

D. Summary of the section
We have numerically shown that one can achieve a highly accurate numerical representation of the equilibrium distribution functions using a mixture of Gaussians.Crucially, this representation makes the evaluation of collision integrals in Eqs. ( 34)- (38) numerically tractable, as we further discuss below.The quality of such fits is controlled solely by the total number of Gaussians (complex Gaussian pairs) M 0 used in the fitting model.In our numerics, we first get self-consistent equilibrium distribution functions and then, for each real-space point r, we obtain a GMM approximation following the procedures outlined in the preceding three subsections; the resulting absolute errors are ensured to be smaller than 10 −7 in all regimes that we investigate in the remainder of the paper.
We have also shown that not only the GMM dramatically outperforms the truncated Taylor series approximation, its complex analogue, cGMM, gives excellent fits in the regimes where the Taylor series expansion is not even available.Since only a few Gaussians are needed to get outstanding-quality fits, it renders calculations of the collision integrals numerically efficient and accurate.
The developed here approach is the central technical result of this paper, as it allows one to study a broad class of quantum many-body systems describable within some Boltzmann equation.In the following section, motivated by the most recent cold-atom experiments, we apply this framework to investigate monopole and quadrupole collective modes of a trapped Bose-Fermi mixture.
We note in passing that a cautious reader might argue that our framework applies specifically to systems with contact s-wave interactions.Indeed, when interactions are explicitly non-local, the resulting Hartree-Fock self-energies, which enter into the equilibrium distribution functions through the self-consistency Eq. ( 13), can become momentum-dependent.In general, this makes the equilibrium forms in Eqs.(44) and (50) no longer applicable.While the full equilibrium distribution functions are generically more complicated, we expect that they still can be accurately approximated by the GMM or cGMM, upon a straightforward extension of our optimization routines.

IV. COLLECTIVE MODES OF A TRAPPED BOSE-FERMI MIXTURE
The formalism developed in the preceding two sections enables us to efficiently evaluate all of the method-ofmoment matrices in Eq. ( 25) for three-dimensional systems.With the GMM and cGMM representations of the equilibrium distribution functions, the multidimensional collision integrals no longer constitute computational obstacles.In practice, however, actual calculations turn out to be quite cumbersome, and, for this reason, detailed derivations, as well as our hybrid numerical-analytical scheme for the collision integrals, are relegated to Appendices A and B. In Appendix C, we present a series of benchmarks and applications of the collision integral calculations, for various physical systems of interest.
Equipped with the method of moments, we turn to investigate the collective modes of a quantum Bose-Fermi mixture trapped in a cigar-shaped optical potential.In Sec.IV A, we first briefly discuss the choice of basis functions and model parameters.Our main findings are then split into three parts: In Sec.IV B, we study a singlecomponent interacting bosonic fluid, with the primary result being the crossover from the collisionless regime to collision-dominated hydrodynamics.The physics of the quantum mixture, when both fluids are present and interact with each other, is rich so that we consider two limiting situations: the case of dilute fermions N F ≪ N B , further refereed to as the impurity regime, is studied in Sec.IV C, and the case with N F ≃ N B , the mixture regime, is discussed in Sec.IV D. There we explore coherent mixing between the bosonic and fermionic monopole modes, their dampening due to incoherent scatterings, and the emergent hydrodynamics characterized by synchronization of the bosonic and fermionic responses.

A. Basis functions and model parameters
So far we have not made any assumptions about the trapping potentials or the basis functions ϕ α (r, p).Below we consider cigar-shaped traps with cylindrical symmetry.In this case, the external trap modulations and the generated by this perturbation responses will have the same z-axis symmetry.This allows us to classify the collective modes and efficiently choose the corresponding basis functions ϕ α (r, p) -see Ref. [24] for a related discussion in two-dimensional pancake-like geometry.In what follows, we investigate two types of collective excitations corresponding to the radial breathing monopole mode [Fig.4a] and to the quadrupole mode [Fig.5a].
The monopole mode represents an excitation with zero z-axis angular momentum, allowing one to restrict the analysis to the following basis functions: and M sets the truncation order.At first order, we have 7 functions {1, ρ 2 , p 2 ρ , ρ • p ρ , z 2 , p 2 z , zp z }, and this simplified choice of basis functions was shown to perfectly capture the monopole mode of a classical gas [75].Specifically, since this gas follows the Boltzmann-Maxwell distribution, one can evaluate all of the method-of-moments matrices analytically, and the result of such an analysis agrees well with more precise molecular-dynamics simulations [75].Besides, the results of Ref. [75] helped us test our collision integral computations -see Appendix C. The situation we consider here, of a quantum system with more complicated distribution functions, is much more challenging and likely not feasible analytically because of the complexity of the collision integral matrix elements.We provide expressions for all of the monopole matrices, for arbitrary M ≥ 1, in Appendix A. Below we investigate the monopole spectral function, which for bosons is defined as: where ) is the response of the transverse bosonic cloud size to a trap modulation of the form: δU B (t) ∝ ρ 2 δ(t).Analogous spectral function can be introduced for fermions as well, but we return to this below.
The quadrupole mode is generated by a perturbation of the type δU B/F ∼ (x 2 − y 2 ), and the corresponding basis functions are given by ξ i=1,2,3 (r, p) × ϕ α (r, p), where We note that 2(ρ • p ρ ) × ξ 2 = p 2 ρ ξ 1 + ρ 2 ξ 3 , implying that a function proportional to (ρ • p ρ ) × ξ 2 should be removed from this basis set since it can be represented as a linear superposition of the rest of the basis functions [24].For instance, this means that for M = 1 we have only 20 functions.In Appendix B, we analytically evaluate all of the corresponding matrices entering the linearized kinetic equations, and the resulting expressions are then used in our numerical simulations.The quadrupole spectral function for bosons is defined as: where χ B x 2 −y 2 (t) ≡ (⟨x 2 ⟩ t −⟨y 2 ⟩ t ) encodes the quadrupole response of the bosonic cloud to a trap modulation of the form: δU B (t) ∝ (x 2 − y 2 )δ(t).
Let us now briefly comment on the validity of the Boltzmann equation to describe the Bose-Fermi mixture.We need to check two conditions: 1) characteristic range of the inter-particle interactions is smaller than the typical travel distance between collisions; 2) it is sufficient to use scattering length to describe the low-energy part of the scattering amplitude, f (k) = −a, and neglect the full momentum dependence of f (k).The lower bound on the travel distance is given by the distance between the atoms.Hence the former condition is guaranteed provided that the distance between bosons (which we always take to be denser than fermions) is larger than both scattering lengths.To discuss the second condition we recall that in three spatial dimensions the s-wave scattering amplitude can be written as [83,84]: where k is the relative wave vector of the scattering particles, a is the scattering length, and r 0 is the effective range of interactions.The use of the Born-Markov approximation for the collision integrals is fully justified if one can neglect the momentum dependence of f (k), i.e., one requires ka ≲ 1.If one uses k = n 1/3 as the typical momentum set by the interparticle separation at low temperatures, then we find that both conditions for the validity of the Boltzmann approach are satisfied provided that a B ≲ 2 × 10 4 a 0 , with a 0 being the Bohr radius.For this estimate, we used the bosonic density at the center of the trap, N B = 10 6 , and parameters of the MIT experiment [12] on a 40 K− 23 Na mixture: ω z = 2π × 12.2 Hz, ω ρ = 2π × 97 Hz, and λ ≈ 2.4.If instead one considers the thermal momentum k ≃ √ mT , then a B ≲ 4 × 10 5 a 0 for T ≃ 500 nK (for comparison, the mean-field transition temperature is T c ≈ 220 nK).The thermal momentum constraint becomes more essential at higher temperatures, away from the quantum regime.In our numerical analyses below, we, therefore, restrict ourselves to |a B |, |a BF | ≤ 10 4 a 0 .In Appendix C, we discuss that our approach can be directly applied to systems with ultrastrong interactions, where one is required to keep the full k-dependence of the scattering amplitude.

B. One-component interacting fluid
We begin by investigating a single-component bosonic fluid, i.e., for now, we assume that there are no fermions.Figure 4 summarizes our results for the monopole mode in this case.We find that when the coupling g B (or equivalently a B ) is weak, the fluid behaves as a non-interacting Bose gas, and, as such, the monopole spectral function A ρ 2 (ω) exhibits a sharp peak in frequency at ω = 2ω ρ .This regime is further referred to as collisionless [Fig.4c].Upon increasing a B (using, for instance, a magnetic Feshbach resonance), we first observe that the resonance in A ρ 2 (ω) becomes dramatically broader [Fig.4d] -the intuition behind this behavior is that the fluid can no longer be considered as non-interacting, since various sound-like excitations of the system mix and, as such, can broaden the signal.Remarkably, however, upon further increase of a B ≳ 5×10 3 a 0 , this peak reforms at a notably different frequency, and it becomes sharper with a B [Fig. 4e].We interpret this result as the system entering the hydrodynamic regime.Interactions are so strong in this regime that the collective response of the fluid to a slow macroscopic perturbation is described via dynamics where the system is locally always in equilibrium.Let us also remark that the width of the monopole mode in the hydrodynamic regime [Fig.4f] is related to the viscosity of the fluid.In other words, by measuring collective modes in optical traps, one can extract information about transport coefficients [59,62,[85][86][87][88][89][90][91][92][93][94][95], and the method we have developed here enables one to accurately evaluate them and compare to such experiments.As a concrete demonstration of this capability of our approach, we analyse in Appendix C the shear viscosity of a two-component strongly-interacting Fermi fluid.
Noteworthy, we find that when the interactions are strong [Fig.4d,e], the spectral function A ρ 2 (ω) exhibits an additional weak resonance at ω ≃ 2ω z .This resonance comes from the longitudinal monopole mode, which is far detuned from the transverse one for strongly anisotropic cigar-shaped traps with κ ≪ 1.The two modes mix with each other for a B ̸ = 0, explaining the appearance of the longitudinal resonance even when the trap perturbation is transverse.
Figure 5 summarizes our results for the quadrupole mode in an interacting bosonic fluid.Here we also find a similar collisionless-to-hydrodynamics crossover in the parameter range consistent with the above discussion.
In the collisionless regime, the quadrupole spectral function A x 2 −y 2 (ω), similarly to A ρ 2 (ω), also exhibits a resonance at ω ≃ 2ω ρ .Interestingly, as the hydrodynamic regime is approached, the quadrupole resonance red-shifts to a frequency notably different from the corresponding monopole peak.Typically, in cold-atom experiments, trapping potentials may be anisotropic, which gives rise to a weak mixing between the monopole and quadrupole modes -this feature might facilitate experiments, as these modes can be seen in a single measurement.Indeed, upon tuning the temperature or interaction strength, one could observe how the two modes, that are almost degenerate in the collisionless regime, split from each other in the hydrodynamic one.Another in- teresting aspect of the quadrupole mode is that it is much broader than the monopole one [see Fig 5f].This has to do with the fact that the phase space, associated with the scattering processes that contribute to the quadrupole linewidth, is significantly larger than the corresponding phase space associated with the monopole mode.
To sum up, we conclude that the overall behavior of the monopole and quadrupole modes is qualitatively similar to each other, except the lifetime of the quadrupole mode can be significantly shorter.For this reason, when considering the Bose-Fermi mixture below, we focus on the monopole mode.

C. Impurity regime NF ≪ NB
Our analysis of the Bose-Fermi mixture starts from the limit of dilute fermions N F ≪ N B .In this case, the bosonic response is barely affected by the presence of fermions, allowing us to focus on the fermionic spectral function.Scenario I,  !% "  (arb.u. ) Scenario II, *  !% "  (arb.u. ) Monopole modes in the impurity regime NF = 10 −2 NB ≪ NB.Effectively, the Bose-Fermi mixture represents a system of coupled harmonic oscillators, where the trap modulation launches both the bosonic and fermionic monopole modes, coupled to each other for aBF ̸ = 0. We consider two types of driving protocols in (a) and (b).The experimental situation would correspond to the scenario I in (a), where the trap modulation perturbs both the bosonic and fermionic traps, in turn directly launching the two modes simultaneously.In the scenario II of (b), the trap modulation couples directly only to the fermionic monopole mode, which then can excite the bosonic one.The bosons remain largely unaffected by the dilute fermions, and, for this reason, we focus here on the fermionic spectral responses.The coherent mixing between the two types of monopole modes competes with the broadening of the fermionic resonance, as the dilute fermions can easily scatter off the bosons -this is illustrated in panels (l) and (k) [corresponding to the scenario II], where the sharp in frequency fermionic peak for aBF = 0 (l) becomes much broader for aBF = −500a0 (k), yet the coupling to the bosonic resonance remains barely notable in (k).Upon further increasing the interaction strength |aBF|, panels (l)-to-(h), the fermionic resonance has largely disappeared, while the bosonic contribution from the back-action of the bosonic cloud becomes notable.In addition to the bosonic resonance, we observe weak sidebands that originate from the longitudinal monopole modes and their mixing with the transverse ones (see text).In the scenario I (a), the optical modulation directly launches the bosonic monopole mode, leading to an interplay between the aforementioned fermionic broadening and coherent driving from the bosonic cloud.This is illustrated in panels (g)-to-(c), as the fermionic resonance becomes weaker and broader until it is completely gone, while the mixing with the bosonic monopole mode quickly becomes strong and completely dominates the fermionic response for large |aBF|.We note characteristic Fano interference profiles that come from the coherent mixing between the two types of monopole modes.For concreteness, here we showed the results for aB = 0, meaning that bosons are collisionless, but we emphasize that similar conclusions hold for other regimes.
For a BF = 0, i.e., when the two fluids are not coupled to each other and the fermionic gas is non-interacting, fermions behave similar to bosons in the collisionless regime studied in the preceding subsection, except now the fermionic response picks up a resonance at Ω F mon = 2ω F ρ [Fig.6(g),(l)].This fermionic monopole frequency is generically different from the bosonic one because of the difference in masses m F /m B ≈ 40/23 and also because we choose λ ≈ 2.4, encoding the fact that in the MIT experiment [12], the trapping potentials for bosons and fermions are different.As such, the Bose-Fermi mixture effectively represents a three-level system shown in Fig. 6(a),(b).In this sense, trap perturbations act as if they excite the two monopole modes, which are coupled to each other for a BF ̸ = 0.
The scenario I in Fig. 6(a) corresponds to an experiment where the trap perturbation simultaneously affects the bosonic and fermionic traps clouds and, as such, directly launches the two monopole modes.The respective fermionic response A F ρ 2 (ω) exhibits two key features.First, as the interaction strength |a BF | is increased, the fermionic resonance at 2ω F quickly becomes broad and weak, until it is no longer observable at large |a BF |.At the same time, A F ρ 2 (ω) acquires an additional contribution from the bosonic resonance, which eventually dominates the fermionic response.Second, for intermediate interactions, which are strong enough so that the bosonic spectral weight is appreciable but not too strong so that the fermionic peak is still sharp, we find coherent mixing between the bosonic and fermionic monopole modes.The spectral lines, even for large |a BF |, acquire characteristic Fano interference profiles typical of a system of coherently coupled harmonic oscillators [96].A prerequisite for observing such profiles is that one of the excited states should be relatively long-lived.This suggests that the bosonic fluid should be either collisionless or hydrodynamic -this aspect will become clearer when we discuss the mixture regime below.
The fermionic behavior at large |a BF | is understood as the interplay between coherent driving from the bosons and incoherent broadening of the fermionic resonance.Indeed, as |a BF | is increased, the dilute fermions can easily scatter off the relatively dense bosons, explaining the eventual disappearance of the fermionic resonance.Simultaneously, the bosonic monopole mode, which has been directly excited by the trap perturbation, acts as a coherent drive to the fermions and, as such, dominates the fermionic response.This physical picture is further supported by the scenario II in Fig. 6(b), where the trap modulation directly launches only the fermionic monopole mode, which then can excite the bosonic one.The corresponding spectral function ÃF ρ 2 (ω) shows that by increasing |a BF | from zero to 500a 0 [Fig.6(l),(k)], the fermionic resonance becomes dramatically broader, whereas the mixing with the bosonic resonance is barely observable.This is because the collision integral, which determines the linewidth of the fermionic resonance, scales as a 2 BF , whereas the (mean-field) coherent coupling to bosons is only linear.Remarkably, upon further increasing |a BF |, the fermionic response becomes so broad in frequency that it has negligible spectral overlap with the initial sharp fermionic resonance.One could even be tempted to argue that as evidenced by this density-density-like response function, the fermion is essentially gone.Such a statement, however, is not entirely correct because the fermion still remains a well-defined quasiparticle, describable within the Landau Fermi-liquid theory that we employ here.
In principle, the fermionic response can become sharp again for larger values of |a BF |, i.e., upon entering the hydrodynamic regime.For dilute fermions, the values required for this are so large that we do not even consider this possibility.Such a situation can occur for reasonable parameters if the fermionic density is appreciable, as we discuss in the following subsection.
Besides this dramatic broadening of the fermionic resonance at large |a BF |, the spectral function ÃF ρ 2 (ω) also acquires a notable contribution from the bosonic resonance.To understand the physical picture, we note that the optical modulation launches the fermionic monopole mode, that is essentially featureless in frequency.Subsequently, this broadband fermionic mode drives the bosonic one, and the back-action from this long-lived excitation in turn gives rise to the relatively sharp bosonic feature in ÃF ρ 2 (ω).Notably, we find that ÃF ρ 2 (ω) also displays weak Stokes and anti-Stokes sidebands near the bosonic resonance.The higher energy sideband has frequency close to Ω B mon +2ω F z , which is due to a process where one quantum of the transverse bosonic mode and one quantum of the longitudinal fermionic mode are both being excited.The lower sideband has frequency close to Ω B mon −2ω F z and corresponds to a process where one quantum of the transverse bosonic mode is being excited and one quantum of the longitudinal fermionic mode is being depopulated.Both types of processes are possible for the thermal equilibrium state, the longitudinal and transverse monopole modes can mix for a BF ̸ = 0, as it was in Fig. 4, and, as such, the system seems to exhibit two additional welldefined collective modes, which manifest in the fermionic response through the same back-action mechanism as for the bosonic monopole mode.These sidebands, together with the Fano interference profiles discussed above, indicate that the Bose-Fermi mixture is expected to give rise to prominent and tunable nonlinearities -this exciting research direction is left for future work.

D. Mixture regime NF ≃ NB
We finally turn to investigate the mixture regime where the fermionic density is comparable to the bosonic one (N F ≃ N B ).
Figure 7 summarizes our results for the case with the bosons being nominally collisionless a B = 50a 0 .We find that a small interaction strength |a BF | ≲ 10 3 a 0 gives rise to both incoherent broadenings of otherwise sharp bosonic and fermionic resonances and coherent mixing between them.For larger interaction strengths, 10 3 a 0 ≲ |a BF | ≲ 2 × 10 3 a 0 , the incoherent part starts to dominate, and the spectral functions appear as rather featureless distinct Lorentzians.At even stronger interactions and in contrast to the impurity regime discussed above, the system becomes hydrodynamic.Both Lorentzians become narrower and they merge into a single, rather than two distinct, resonance -an effect we refer to as mode locking.This synchronization effect represents a hallmark of the hydrodynamic regime of the Bose-Fermi mixture as a whole.
Similar results hold for the case where bosons are nominally hydrodynamic a B = 5 × 10 3 a 0 , as summarized in Fig. 8. Remarkably, when the interaction |a BF | ≲ 600a 0 is not too strong, we now observe that the coherent mixing is between the collisionless fermionic monopole mode and truly many-body hydrodynamic bosonic mode.We also find that the incoherent scatterings play a role and become dominant in the crossover region, where the two response function appear as two broad Lorentzians.Even the bosonic hydrodynamic mode becomes broader with increasing |a BF | in this region.Finally, as |a BF | is further increased, the Bose-Fermi mixture enters the hydro- Monopole modes in the mixture regime NF = NB, with bosons being nominally collisionless aB = 50a0.When the Bose-Fermi coupling is not too strong, |aBF| ≲ 10 3 a0, we find that the bosonic and fermionic monopole modes exhibit both coherent mixing with each other and incoherent broadening of the respective resonances.For stronger interactions, 10 3 a0 ≲ |aBF| ≲ 2 × 10 3 a0, each of the spectral functions displays a single Lorentzian behavior, meaning that the coherent mixing between the two monopole modes is no longer observable -this regime can be associated with the crossover region.Notably, as |aBF| is further increased, the system turns into the hydrodynamic regime, where not only the two Lorentzians become sharper, as it was in Fig. 4, they also merge to have the same frequency -an effect referred to as mode locking.
dynamic regime as a whole, characterized by the aforementioned mode locking and sharpening of the response functions.

V. CONCLUSION AND OUTLOOK
Our work elucidates the collective near-equilibrium dynamics of quantum Bose-Fermi mixtures in trapped cold atom systems.We analyze linear responses of the system to changes in the confining potential and find that they differ dramatically depending on the parameter regimes.We observe several non-trivial phenomena, including coupling and interference between bosonic and fermionic collective modes, strong damping of these modes, and the emergence of hydrodynamics in the strongly interacting regime.For intermediate and strong interactions, we observe appreciable mode mixing, manifesting, for in- FIG. 8.The same analysis for the mixture regime NF = NB as in Fig. 7 but now bosons are nominally hydrodynamic aB = 5000a0.Quite remarkably, for |aBF| ≲ 600a0, we find that the collisionless fermions display coherent mixing with the hydrodynamic bosonic mode.A bit unusual is that this latter mode also exhibits initial broadening as the interaction strength |aBF| is increased.The rest of the results qualitatively follow those in Fig. 7, including eventual mode locking characteristic of the hydrodynamic regime of the mixture as a whole.
stance, as the appearance of a longitudinal monopole resonance even when the optical perturbation is transverse.In the dilute fermion limit, the bare fermionic resonances quickly become featureless in frequency upon increasing the boson-fermion coupling strengths, which is reminiscent of the physics of Bose polarons.When bosonic and fermionic densities are comparable and interactions are strong enough to make the mixture hydrodynamic, we find mode locking between the bosonic and fermionic responses.Remarkably, we also observe coherent mode mixing even when bosons are hydrodynamic and fermions are collisionless.The coherent mode mixing effects we uncover between longitudinal and transverse monopole modes, manifesting as Fano interference profiles and the emergence of the Stokes and anti-Stokes sidebands, indicate that Bose-Fermi mixtures should exhibit strong and tunable nonlinearities.While nonlinear effects are beyond the linearized Boltzmann framework analyzed here, they are a promising direction for future inquiry.Such non-linearities can potentially be used to explore analogues of non-linear optical phenomena using collective excitations rather than Λ-like atomic gases [96,97].Potential applications include generating single and two mode squeezing [98], entangling different collective modes [99], and achieving reflections with phase conjugation [100].Analysis of the dynamics of the cold atomic Bose-Fermi mixture was made possible by the GMM-MoM developed here.This method allows one to efficiently compute the linear-response properties, including the collective modes and transport coefficients, of a large range of quantum many-body systems that are describable within the framework of the Boltzmann equation.While this work focuses on a spatially-inhomogeneous 3D system, our approach will also accelerate investigations of salient 2D systems.Examples of phenomena in such systems that require further elucidation include electron hydrodynamics in vdW semiconductors, spin transport in bernal graphene, and the collisionless-to-hydrodynamic crossover in pancake-like dipolar gases [24,[101][102][103][104].The tables of GMM fits to equilibrium distribution functions we provide are generic and can be directly applied to a broad class of interesting many-body systems [13].
While the trapped cold-atom mixture we analyzed in this work has particle-conserving scattering processes, particle non-conserving collisions are also ubiquitous in solid-state and cold-atom platforms.Such processes are natural for systems containing phonons, magnons, and photons.Particle non-conserving collision integrals, however, are usually much simpler compared to the particleconserving ones considered in this work and therefore are expected to be straightforwardly tractable using the same GMM-MoM approach.An additional step would be required to treat systems that have explicitly non-local (but spherically symmetric) interactions such as those generated by screened Coulomb or dipolar forces.Systems with pure Coulomb interactions may need to be modeled with Boltzmann-Vlasov equations; we expect that our approach facilitates analysis of such models as well.The GMM representaion of the distribution functions may not reduce the collision integrals to fully analytical expressions in these cases, unlike the contact s-wave interaction in the cold atom mixture considered here.This representation would still, however, dramatically simplify the collision integral, reducing it to a low-dimensional one that can easily be computed numerically using limited computational resources.
There are two generalizations of the GMM-MoM approach that may increase its efficiency.For systems that are spatially inhomogeneous, such as optically-trapped cold atomic gases similar to those considered in this work, one may extend the GMM to an ansatz that encodes the full dependence of the equilibrium distribution functions on both real-space and momentum coordinates.This will further facilitate calculations of the collision integrals.Additionally, one may try mixture models composed of functions other than Gaussians that still allow simplification of the collision integral via Wick's theorem.Possible options include poly-Gaussian functions, composed of a sum of polynomials multiplied by a Gaussian, and error functions.Both of these ansatzes still enable semianalytical computation of the collision integrals and may prove advantageous, for example, in the analysis of dense fermionic systems.
Aside from these simple generalizations, there are three non-trivial promising extensions of our GMM-MoM framework.Firstly, one can try to use a GMM with time-dependent Gaussians to investigate the farfrom-equilibrium dynamics of quantum many-body systems.This approach would apply to systems which can be modeled using a Boltzmann equation with timedependent distribution functions.
If such a timedependent GMM-MoM proves accurate, it would enable significant progress in several open topics in modern physics.Examples include the behavior of nonequilibrium many-body steady-states, with applications to spintronics devices and field-effect transistors essential for classical computer hardware, fundamental questions about turbulence and hydrodynamics, and the phenomenology of light-induced phases of matter such as light-induced superconductivity and magnetism.Secondly, the Fano interference profiles and Stokes and anti-Stokes sidebands we find in the Bose-Fermi mixture studied here motivates the extension of the GMM-MoM framework to include nonlinear effects.Such an extension would facilitate analysis of the non-linear electrodynamics of correlated materials, allowing, for example, study of the temperature dependence of photocurrents in photovoltaic semiconductors or analysis of second or third harmonic generation in systems such as striped cuprate semiconductors which exhibit intertwined order parameters [105][106][107][108][109][110].Finally, we can try to extend the GMM-MoM framework to strongly-correlated systems beyond Landau Fermi liquids which are describable using a Boltzmann equation.Such an extension would necessitate a Gaussian mixture ansatz on the level of Green's functions rather than on the single-particle distribution functions at the heart of the Boltzmann model.Generalizing the GMM-MoM in this way would enable the analysis of non-Fermi liquids, which are believed to be the key to understanding strongly-correlated phases of matter.
Our work opens the door to rigorous theoretical investigation of many fundamental and technologicallyimportant quantum many-body systems describable by Boltzmann equations which have previously evaded analysis.Furthermore, extensions of the method may provide an additional tool to analyze quantum many-body systems that are out-of-equilibrium, manifest non-linear behavior, or are strongly-correlated.

Matrix elements of Ĥ
We first note that the relevant Poisson brackets in Eq. ( 28) can be written as: where we defined Plugging this into Eq.( 28) and following the preceding subsection, for the fermionic sector, we arrive at: The expression for the bosonic sector is obtained from this one by substituting n F,eq (1 − n F,eq ) → n B,eq (1 + n B,eq ), putting the mass ratio to one m B m F → 1, and by replacing γ F → γ B .As in the preceding subsection, we evaluate all the momentum integrals analytically and all the real-space integrals numerically.

Matrix elements I αβ
In this subsection, we show how to evaluate the collision matrix I αβ in Eq. ( 34) associated with incoherent bosonboson scatterings.The main idea of the construction below is that the GMM fits allow one to represent the complicated matrix elements as sums of Gaussian integrals that can be evaluated analytically over the 12-dimensional momentum space (the integration over the remaining 3-dimensional real space is then done numerically).Some difficulties still exist, such as taking into account the momentum and energy conservation laws -we discuss this below.
It will be convenient to write S[ϕ α ] = β Nα κ pα+rα/2 ρ 2mα+kα z 2pα+rα S α [P , q, q ′ ], where S α [P , q, q ′ ] is a known polynomial of its arguments.Here, it is implicit that we consider all the in-plane momenta relative to ρ -this is allowed to do since, for the monopole mode, the integrals do not depend on the angle ϕ between ρ and the x-axis.Rescaling the variables in the usual manner, i.e., ρ √ β → ρ, z √ βκ → z, and p √ β → p, we arrive at: We note that the real-space angular integration can be readily done: Here Poly[P , q, q ′ ] = S α [P , q, q ′ ] × S β [P , q, q ′ ] is a polynomial of its arguments.It is worth pointing out that polynomials are particularly suitable for our construction of the GMM fits because, as we mentioned, after all manipulations, Here, the new polynomial Poly[k, k ′ ] is obtained from Poly[q, q ′ ] by simply substituting q = ( Ãk − Bk ′ )/( Ã2 − B2 ) and q ′ = ( Ãk ′ − Bk)/( Ã2 − B2 ), i.e., through the inverse linear transformation.The integral in Eq. (A19) is already feasible for further analytical/numerical evaluations.Indeed, in the spherical coordinates, it becomes: where S αβ [k, ϕ, ϕ ′ , θ, θ ′ ] is obtained from Poly[k, k ′ ] by substituting the spherical coordinates and k = k ′ .Importantly, each function S αβ is some known polynomial of cos θ, sin θ, cos θ ′ , sin θ ′ , cos ϕ, sin ϕ, cos ϕ ′ , and sin ϕ ′ -as such, the angular integration can be done analytically.In our numerics, we write a simple Mathematica code that symbolically evaluates this four-dimensional angular integration.The remaining Gaussian integral over k is then also evaluated analytically (we write a separate Mathematica subroutine for this step as well).This completes our protocol for evaluating the momentum integrals.At this stage, we are left with the integral over r, which is reduced to a one-dimensional integral over r and then evaluated numerically.Explicitly, the first classical contribution to (I B ) αβ reads: and similar expressions hold for the remaining three contributions.

Matrix elements of Ĵ
We turn to discuss the matrix Ĵ associated with collisions between bosons and fermions, Eqs. ( 35)- (38).Now, the corresponding matrix elements are evaluated in exactly the same manner as in the preceding subsection.The only new feature we need to address in this subsection is that the bosonic and fermionic masses might differ.In this case, the transformation to the center-of-mass frame now reads: As above, the momentum conservation enforces P = P ′ , while the energy conservation gives q = q ′ : The rest of the computation follows the preceding subsection step-by-step.

Matrix elements of M
Similar algebra as in the preceding Appendix gives: where we defined: In the derivations here, we used the same rescaling as above, namely: ρ √ β → ρ, z √ βκ → z, and p √ β → p.

Matrix elements of Ĥ
We note that {ξ j ϕ β , H B/F } = ξ j {ϕ β , H B/F } + ϕ β {ξ j , H B/F }.The first term was evaluated above, in Eqs.(A8)-(A9), while the second term can be written as Following the above form for the Poisson brackets, we write (H , and a straightforward algebra in the spirit of preceding calculations gives: We also get: Similar expressions can be obtained for the bosonic sector (or use the same expressions but substitute:

Matrix elements of Ŝ
We note that the self-energies entering the matrix Ŝ in Eq. ( 29) now depend not only on ρ but also on the angle ϕ between ρ and the x-axis: where we rescaled the variables as ρ √ β → ρ, z √ βκ → z, and p √ β → p and introduced: Following the Poisson bracket expansion (see the preceding subsection), we write (S B ) ij αβ = (S B ) ij αβ,1 +(S B ) ij αβ,2 , where: Similar expressions can be obtained for (S BF ) αβ and (S FB ) αβ .

Matrix elements of Î
We compute the matrix Î in the same manner as in the preceding Appendix.The only new element we encounter here is that the polynomial integrands entering collision matrices now depend on the angle ϕ between ρ and the x-axis.To overcome this issue, we write (in the rescaled variables) expressions of the type: S[ϕ α ξ i ] → S α,i [P , q, q ′ , ϕ]ρ 2mα+kα+µi z 2pα+rα , where we also switched to the center-of-mass frame introduced in Appendix A 4. These new functions S α,i [P , q, q ′ , ϕ] are polynomials of cos ϕ and sin ϕ allowing for efficient analytical/numerical computations of In our numerics, we write a Mathematica subroutine that evaluates the functions S ij αβ [P , q, q ′ ], which are now polynomials of P , q, and q ′ .Plugging these polynomials back into, for instance, the matrix I B , we arrive at (compare to Eq. (A14)): × (1 + n B,eq (r, p))(1 + n B,eq (r, p ′ ))n B,eq (r, p 1 )n B,eq (r, p ′ 1 ) × S ij αβ [P , q, q ′ ]. (B11) The rest of the computation, including the computation of the matrix Ĵ which accounts for collisions between bosons and fermions, follows the procedure in Appendix A step-by-step.
(a) (b) Monopole mode of a classical gas, trapped into a harmonic potential with κ = 10 −2 .Solid lines in (a) and (b) are the monopole mode frequency Ωmon and its dampening Γmon, respectively, computed using the approach developed here.Dashed lines were extracted from Ref. [75].Similarly to our results in the main text, cf.Fig. 4, the classical gas also exhibits a collisionless-to-hydrodynamic crossover, as evidenced, for instance, by the non-monotonic behavior of Γmon.
Appendix C: Benchmarking and applications of collision integral calculations

Monopole mode of a classical gas
A simple and straightforward way to benchmark the method developed in this work is to consider a classical gas trapped in a harmonic potential.Such a gas at equilibrium follows the Boltzmann-Maxwell distribution function, which is Gaussian.The Hartree-Fock self-energy of the gas can also be neglected.Correspondingly, in the GMM ansatz in Eq. ( 43), we only need to keep track of a single harmonic (M 0 = 1).In case of a dilute gas, we can further approximate 1 + n eq ≈ 1, which simplifies the calculations below.
We focus on the monopole mode and consider the truncation order to be M = 1, so that in total we have 7 basis functions for the method of moments.All of the matrices entering the kinetic equation ( 25) can be computed analytically: where N tot = e µβ /(β 3 √ κ) is the total number of particles and µ is the chemical potential.Note that following the main text notations, we set m = 1 in this subsection.The collision integral matrix I αβ has only four nonzero matrix elements, namely: I 33 = I 66 = N tot I and I 36 = I 63 = −N tot I, where By switching to the center-of-mass frame following Appendix A 4, we get: Integration over the real-space gives: We note that S[p 2 z ] = 2(q 2 z − q ′2 z ), i.e., it does not depend on P , which can now be easily integrated out: The remaining integration is straightforward: where v T = 8/(πβ) is the thermal velocity, n 0 is the density at the center of the trap, and γ coll = σn 0 v T /2.Our numerical subroutines reproduce the above analytical matrices and reproduce the known results from Ref. [75], as shown in Fig. 9. Specifically, Ref. [75] demonstrated that the method of moments with M = 1 applied to a classical  10.The fermion (Bose polaron) decay rate as a function of temperature for unitary (Ta = 0) and non-unitary Bose-Fermi interactions (Ta = 5Tc).Solid curves correspond to resuming the Taylor series expansion, Eq. (C14) derived in Ref. [95].Dashed curves correspond to the GMM ansatz, Eq. (C12), where the maximal bosonic fugacity is chosen to be zmax = 0.99 (black dashed vertical line) so that we have M0 = 17 Gaussians in total, cf.Fig. 3(b).The domain of applicability of the GMM approach is T ≳ 1.1Tc, where the two types of calculations essentially coincide.Notably, the GMM predictions for T ≲ 1.1Tc (shaded region) are also in remarkable agreement with the more accurate Taylor series resummation.
thermal gas agrees remarkably well with more precise molecular dynamics simulations.From this result, one can also conclude that M = 1 is expected to be an excellent approximation.If one wants to study a quantum gas where the equilibrium distribution function is no longer Gaussian (as done in the main text), one may need to consider M ≥ 1.In the main text, we fixed M = 2 and checked that further increasing M does not result in an appreciable difference in the resulting spectral functions.

Bose polaron decay rate
In this subsection, we demonstrate that the GMM approach can be applied to experimentally relevant quantum many-body systems with strong interactions; in these systems, one cannot disregard the momentum dependence of the scattering amplitude.We follow Ref. [95] and consider the impurity limit where we have a single fermion (Bose polaron) immersed into a homogeneous bosonic bath.The fermion decay rate Γ in this case can be estimated as [95]: where v rel = k/m r = k F /m F − k B /m B , k is the relative momentum between the two scattering particles, and m r = m F m B /(m F + m B ) is the reduced mass.We neglect the Bose-Bose interactions, a B = 0, corresponding to an ideal Bose gas.We will also consider temperature T to be both below and above T c .Below T c , the fermion decay rate comes from scatterings off the thermal bosons.For T ≤ T c , the bosonic chemical potential is zero µ B = 0 (z B = 1).Since the Bose-Fermi coupling can be strong, we take into account the k-dependence of the scattering cross section: where a BF is the Bose-Fermi scattering length.At this stage, we employ the GMM expansion: In what follows, we will avoid writing the subscript B to ease the notations, as it is clear that the GMM ansatz in this subsection refers solely to bosons.The scattering rate Γ is then written as: To evaluate this integral, we introduce the following change of variables from k F and k B to K and k: Here α = m F /m B .This linear transformation (i) has unity Jacobian, (ii) respects that k is the relative momentum, and (iii) allows us to directly integrate out K. After simple algebra, we obtain: where v rel = 8T /(πm r ), T a = 1/(2m r a 2 BF ), and If one were to use a Taylor series expansion instead, then one should substitute: The result in Eq. (C14) exactly reproduces Eq. (32) of Ref. [95].
The physics of the Bose polaron decay rate is quite fascinating, and the interested reader is referred to the thorough discussion in Ref. [95].Here, we focus on the validity of our GMM calculations.Figure 10 shows the comparison between the GMM result in Eq. (C12) and the Taylor series expression (C14) for both unitary (T a = 0) and nonunitary interactions (T a = 5T c ).For the GMM ansatz, we set the maximal bosonic fugacity to be z = 0.99 so that we have M 0 = 17 Gaussians in total, cf.Fig. 3(b).This choice of the maximal fugacity implies that the domain of applicability of the GMM calculations is T ≳ 1.1T c , where the two types of calculations coincide.Even for T ≲ 1.1T c , the two approaches agree remarkably well (see Fig. 10).
To get a sense of how the two calculations differ from each other for T ≲ 1.1T c , we consider the unitary limit with T a = 0.In this case, both approaches predict linear in T behavior for T ≤ T c (where the bosonic fugacity is z = 1), meaning the maximal discrepancy between the two calculations can be estimated to be |Γ(T c ) − Γ TS (T c )|/Γ TS (T c ) ≈ 0.7%.This remarkable agreement, in the most challenging regime at T = T c and with unitary interactions, stems from the fact that the bosonic GMM expansion can be accurately extrapolated outside its domain of applicability, cf.Fig. 3(b), as discussed in Sec.III of the main text.

Shear viscosity of a strongly-interacting two-component Fermi system
In this subsection, we demonstrate that the fermionic version of the GMM, cGMM, readily enables computation of interesting transport coefficients such as the shear viscosity.To this end, we follow Ref. [86] and consider a homogeneous two-component fermionic fluid in its normal state where, within the method of moments, the shear viscosity is given by [86]: √ mTF /(32a 2 √ π)).The solid curve is taken from Ref. [86], while the dashed one is obtained using the cGMM, and the two results reasonably agree.The choice kF |a| = 4.5 implies that Ta ≈ 0.1TF and allows one to see the ∼ T 3/2 scaling (dotted line) at high temperatures, TF , Ta ≪ T .The shaded regions in (a) and (b) are where the chosen cGMM is no longer applicable.
In an above expression, n p = [exp(β(ε p − µ)) + 1] −1 is the Fermi-Dirac distribution function, 1 ⇒ S[p x p y ] = 2(q x q y − q ′ x q ′ y ), and the differential cross section is given by (q = (p − p ′ )/2): where Ω is the solid angle of the relative momentum q ′ = (p 1 − p ′ 1 ) with respect to the direction of q.In the limit where the fermionic distribution function is approximately Gaussian (E F ≪ T ), the shear viscosity can be computed analytically: where T a = 1/(ma 2 ).In the limit T ≪ T a , we have η cl = 5 √ mT /(32a 2 √ π) ∼ √ T .In the opposite limit, T a ≪ T , we get η = 15(mT ) 3/2 /(32 √ π) ∼ T 3/2 (dotted line in Fig. 11(b)).At lower temperatures, where the distribution function is no longer Gaussian, we can employ the cGMM expansion: which is valid for any z = exp(βµ) ≤ z max .We choose z max = 10 5 , cf.Fig. 3(c).This choice limits the lowest temperature in our analysis to be T min ≈ 0.1E F (for this estimate, we used µ = E F ).If one is interested in even lower temperatures, then one should increase z max and redo the cGMM fitting following Sec.III.As a first step, we solve for the temperature dependence of the chemical potential µ(T ), which can be obtained by fixing the total fermionic density: Here, the factor of 2 encodes that we have a two-component Fermi fluid.The momentum integration in Eq. (C19) can also be carried out using Eq.(C18), allowing us to determine the chemical potential within the cGMM approximation.We find that within the applicability of the cGMM expansion (T ≳ 0.1E F ), the two calculations for µ(T ) coincide (see Fig. 11(a)).The primary challenge we encounter when analyzing Eq. (C15) is the evaluation of the denominator which, in the center-of-mass frame, is proportional to: dΩ 4π q q 2 x q 2 y a 2 1 + q 2 a 2 n p n p ′ (1 − n p1 )(1 − n p ′ 1 ).(C20) We normalize this integral using the corresponding classical expression, where one additionally neglects the momentum-dependence of the scattering cross section: As in Appendix A 4, we write I as a sum of four terms I = I 1 − I 2 − I 3 + I 4 and plug in the cGMM expansion (C18).Under the summations over the cGMM indices, we will get the following type of expressions, cf.Eq. (A15): dΩ 4π q q 2 x q 2 y 1 + q 2 T /T a d 3 P (2π) 3 exp − P 2 2γ 2 P − q 2 γ 2 q + 1 γ P (Aq + Bq ′ ) • P e 2βµ d 3 q (2π) 3 dΩ 4π q q 2 x q 2 y d 3 P (2π) 3 exp − Note that we rescale all the momenta as βp 2 /m → p 2 .Integration over P yields, cf.Eq. (A16): Integration over the solid angle Ω yields: Plugging these results in, we rewrite Eq. (C22) as: This one-dimensional integral, together with the summation over the cGMM harmonics, is easy for numerical evaluations on a single laptop.Figure 11 shows the computed temperature dependence of the shear viscosity η(T ): Our result reproduces that of Ref. [86], demonstrating that the cGMM expansion allows for efficient and accurate computation of transport coefficients of strongly-interacting quantum many-body systems.

FIG. 3 .
FIG. 3. Performance of the Gaussian mixture model (GMM).(a)The fermionic function f1(x) (solid line), Eq.(45), that encodes the Fermi-Dirac distribution function for any zF ≤ 1, and its GMM fit (dashed line) are not distinguishable by eye from each other.The inset shows the dependence of the absolute error of the fit LF, Eq. (48), on the total number of Gaussians M0 that are involved in the fitting: for M0 = 10, we obtain LF ≲ 10 −10 .(b) Similar analysis as in (a) but now we consider the bosonic function f2(x) (solid line), Eq. (51), in a range x ∈ [xmin, ∞), where xmin (vertical dashed lines) is set by the maximal fugacity we consider: xmin ≈ 0.46 for zB = 0.9 and xmin ≈ 0.14 for zB = 0.99.The GMM fits provide excellent approximations to f2(x) in their respective domains of applicability (dashed line for zB = 0.9 and solid dashed line for zB = 0.99): The inset shows that the bosonic error LB, Eq. (52), can be as small as 10 −8 with only about 10 (17) Gaussians for zB = 0.9 (zB = 0.99).(c) The fermionic function f3(x) (solid line), Eq. (55), that captures the Fermi-Dirac distribution function for any 1 ≤ zF ≤ zmax, can be well fitted with the complex version of the GMM, Eq. (56).Here, the domain x ∈ [xmin, ∞) is determined by the maximal fermionic fugacity we consider: xmin = − ln(zmax).As shown in the inset, the absolute error of the fit LF decreases with increasing the total number M0 of complex conjugate Gaussian pairs; for M0 ≲ 25, we obtain LF ≃ 10 −8 at zmax = 10 5 .

FIG. 4 .
FIG. 4. Collisionless-to-hydrodynamics crossover in a singlecomponent bosonic fluid, the case of monopole mode schematically shown in (a).(b) The bosonic spectral response A ρ 2 as a function of frequency ω and for different interaction strengths, as encoded in aB.When the Bose-Bose interactions are weak, one gets a sharp response at ω ≃ 2ωρ; as aB is increased, the system starts to exhibit a many-body hydrodynamic response, as evidenced by the fact that first, the spectral function completely broadens but then reforms again into a sharp peak at a frequency notably different from 2ωρ [see also a few cuts shown in (c,d,e)].This is further illustrated in (f) as the behavior of the position Ωmon and width Γmon of the peak in (b), obtained from a single Lorentzian fit.In the hydrodynamic regime, transverse and longitudinal monopole modes mix, resulting in the emergence of an additional weak peak at the longitudinal frequency ω ≃ 2ωz (e).Parameters used: ωz = 2π × 12.2 Hz, ωρ = 2π × 97 Hz, T = 500 nK, NB = 10 6 , and M = 2.

FIG. 5 .
FIG.5.The same analysis as in Fig.4but for the quadrupole mode.(b) The quadrupole spectral function A x 2 −y 2 (ω), similarly to A ρ 2 (ω), exhibits a crossover from the collisionless regime, aB ≲ 10 3 a0, to a hydrodynamic one, aB ≳ 5 × 10 3 a0.Just as it was for the monopole case, for weak interactions, the function A x 2 −y 2 (ω) peaks in frequency at ω ≃ 2ωρ.As the interactions are increased, this peak first significantly broadens [this broadening is much stronger compared to the one seen in the crossover regime for the monopole mode in Fig.4] and then reforms at around ω ≃ √ 2ωρ -this frequency is notably different from the corresponding monopole resonance in the hydrodynamic regime.

3 d 2 .FIG. 11 .
FIG. 11.Performance of the cGMM expansion on a strongly-interacting Fermi fluid.(a) The chemical potential µ as a function of T : the solid line encodes the exact dependence, while the dashed one corresponds to the cGMM.In our simulations, we chose zmax = 10 5 , which implies that the domain of applicability of the cGMM is T ≳ 0.1TF (vertical dashed line), where the two calculation coincide.(b) The shear viscosity η of the fluid as a function of T at kF |a| = 4.5 (η is normalized byη cl (TF ) = 5 √ mTF /(32a 2 √ π)).The solid curve is taken from Ref.[86], while the dashed one is obtained using the cGMM, and the two results reasonably agree.The choice kF |a| = 4.5 implies that Ta ≈ 0.1TF and allows one to see the ∼ T 3/2 scaling (dotted line) at high temperatures, TF , Ta ≪ T .The shaded regions in (a) and (b) are where the chosen cGMM is no longer applicable.