{'fraction_non_alphanumeric': 0.042720240157025746, 'fraction_numerical': 0.013768617942500866, 'mean_word_length': 4.894011568560735, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ab-initio study of interacting fermions at finite temperature with neural canonical transformation

We present a variational density matrix approach to the thermal properties of interacting fermions in the continuum. The variational density matrix is parametrized by a permutation equivariant many-body unitary transformation together with a discrete probabilistic model. The unitary transformation is implemented as a quantum counterpart of neural canonical transformation, which incorporates correlation effects via a flow of fermion coordinates. As the first application, we study electrons in a two-dimensional quantum dot with an interaction-induced crossover from Fermi liquid to Wigner molecule. The present approach provides accurate results in the low-temperature regime, where conventional quantum Monte Carlo methods face severe difficulties due to the fermion sign problem. The approach is general and flexible for further extensions, thus holds the promise to deliver new physical results on strongly correlated fermions in the context of ultracold quantum gases, condensed matter, and warm dense matter physics.


I. INTRODUCTION
Consider an interacting quantum system of N fermions in a d-dimensional continuous space with the generic Hamiltonian where (2) (r i − r j ) consists of one-and two-body potentials. We use the shorthand notation x ≡ (r 1 , r 2 , . . . , r N ) ∈ R dN to collectively denote all fermion coordinates. We also assume appropriate natural units so that constants like fermion mass and the Planck constant can be omitted. We would like to study thermodynamic properties of the system encoded in the partition function Z = Tre −βH at inverse temperature β = 1/k B T , which are relevant to a broad range of problems including ultracold Fermi gases [1], condensed matter [2], and warm dense matter [3].
Unfortunately, accurate ab-initio study of Eq. (1) at finite temperatures is generally difficult. As a typical workhorse, quantum Monte Carlo methods suffer from the notorious fermion sign problem at low temperatures [4][5][6][7]. There are extensions of ground state quantum chemistry methods to finitetemperature, e.g. [8]. On the other hand, a fundamental principle to solve quantum system at finite temperature is to minimize the variational free energy F = 1 β Tr(ρ ln ρ) + Tr(ρH) (2) with respect to a variational density matrix ρ. The two terms in Eq. (2) correspond to the entropy and energy of the system, respectively. It can be shown that F − 1 β ln Z is a variational upper bound of the true free energy, where the equality holds only when ρ coincides with the exact density matrix 1 Z e −βH of the system [9].
There are a number of physical constraints on the variational density matrix ρ. Besides basic properties like Hermitian (ρ † = ρ), positive definiteness (ρ 0) and normalization (Trρ = 1), it should also be antisymmetric Px|ρ|x = (−1) P x|ρ|x with respect to permutations Px ≡ (r P1 , r P2 , . . . , r PN ) of the fermion coordinates. The challenge is then to devise a tractable computational scheme to perform optimization of ρ within such a constrained space. In practice, the entropy term in Eq. (2) often turns out to be difficult to compute [10,11]. As a result, most of the previous variational density matrix studies resort to an alternative imaginary-time evolution approach [10,[12][13][14], which is more appropriate at high temperatures.
In the low-temperature regime, the variational density matrix can be reasonably represented by a truncated 1 set of lowenergy many-body basis states |Ψ n : ρ = n µ n (φ)|Ψ n (θ) Ψ n (θ)|, where φ and θ are variational parameters. A discrete probabilistic model µ n , which satisfies 0 < µ n < 1 and n µ n = 1, is used to parametrize the Boltzmann distribution of the basis states |Ψ n . On the other hand, we choose to construct |Ψ n by applying a unitary transformation to a set of reference basis states |Φ n , e.g., the non-interacting Slater determinants: |Ψ n (θ) = U(θ)|Φ n . To make such an ansatz powerful 1 Truncation is necessary for systems in the continuum with infinitedimensional Hilbert space. Given N fermions and M available singleparticle orbitals, the summation in Eq. (3) involves M N terms, which is exponentially large. Nevertheless, such summation appearing in relevant physical quantities can be estimated via Monta Carlo sampling and will not cause big troubles in practice. See also discussions in Sec. IV. enough, the unitary transformation should have many-body nature, so that particle correlations can be effectively introduced into the reference state. In addition, it should also preserve permutation antisymmetry of the reference wavefunction, which we will refer to as the equivariance property. Overall, the modeling of the variational density matrix in the present approach is illustrated in Figure 1.
Parametrizing and optimizing a rich family of equivariant many-body unitary transformations U turn out to be a fairly nontrivial task. In this paper, we present an elegant solution to this problem by constructing U as unitary representation of the canonical transformation of phase space variables in classical mechanics, extending the previous work on neural canonical transformations [15] from classical to quantum domain. The resulting approach naturally generalizes the ground-state variational Monte Carlo (VMC) method [16,17] to finite temperatures and is not hindered by the fermion sign problem. Moreover, based on Born's probabilistic interpretation of wavefunctions, the equivariant unitary transformation turns out to be intimately related to equivariant normalizing flow [18,19], an important class of generative model developed well within the deep-learning community. This way, one can leverage the latest technical advances in probabilistic modeling to efficiently tackle the thermodynamics of strongly correlated fermions in a fully ab-initio way.
It is worth mentioning a variety of related works to put the present contribution into a broader perspective. First, there have been various wavefunction ansatzes for ground-state VMC calculation of fermions, from the traditional Slater-Jastrow [20], backflow [21][22][23][24] to more recent attempts based on neural networks [25][26][27][28][29][30][31][32][33]. However, unitary transformations are not considered in these ansatzes, since only a single wavefunction, instead of a whole basis, is needed in this situation. Second, there have been quantum algorithms for thermal properties of model Hamiltonians [34][35][36], which rely on quantum circuits to construct the unitary transformation. However, they still demand advances in quantum technologies to be practically useful. Third, variational free energy studies of statistical mechanics and field theory problems [15,[37][38][39][40][41][42] can be regarded as the classical counterparts of the present approach. Last but not least, the so-called quantum flow approach [43] also performs a learnable unitary transformation to a single-particle basis. In the many-particle settings considered here, one has to additionally deal with the permutation antisymmetry by imposing the equivariance property into the coordinate transformation carried out by the normalizing flows [44][45][46]. In this way, normalizing flows have also precisely addressed the open problem envisioned in [47]: "The full use of the (coordinate) transform to compute from first principles requires adequate approximation to the Jacobian and the inverse transformation."

II. THEORY
The method of constructing unitary transformations in this work is based on the fact that one can establish a one-toone correspondence between the group of unitary transforma-FIG. 1. Architecture of the variational density matrix representation (3) of the present approach. A discrete probabilistic model µ n parametrizes the Boltzmann distribution of a many-body basis |Ψ n . We construct |Ψ n by applying a parametrized unitary transformation to a reference basis state |Φ n . The unitary transformation corresponds to a permutation equivariant many-body coordinate transformation implemented as a flow of fermion coordinates. tions in quantum mechanics and the group of canonical transformations of phase space variables (x, p) in classical mechanics [48][49][50][51]. One can develop some basic understanding of the idea by inspecting the infinitesimal structure of these two groups. In classical mechanics, one can use an arbitrary generating function G(x, p) to define a continuous family of canonical transformations via the symplectic evolution where λ denotes a continuous parameter. When the canonical transformation is quantized, the generating function is converted to a Hermitian operatorĜ 2 , and the corresponding unitary transformation then takes the form U λ = e −iĜλ . See Appendix A for more details.
An important class of canonical transformation is the socalled point transformations, in which the new generalized coordinates depend solely on the old coordinates, not on the old momenta. The generating function of point transformation reads G = u(x) · p, where u : (r 1 , · · · , r N ) → (u 1 , · · · , u N ) is a function in the dN-dimensional coordinate space. The equation of motion followed by the transformed coordinates then takes the form u can be intuitively viewed as a vector field that guides all the particles to continuously flow in the coordinate space as the parameter λ increases 3 . To see this, one can consider the examples of spatial translation and rotation, which can be generated by the total momentum G = e · N i=1 p i or angu- For systems with such spatial symmetries, these transformations would leave the Hamiltonian unaltered. We also note that the vector fields associated with these two examples are "separable", i.e., the vector field experienced by one particle is independent of the positions of any other particles, thus completely ignores the interactions among them. For practical variational calculations, we need to seek for vector fields u that can effectively introduce particle correlations, as illustrated in the right panel of Fig. 2.
To study the unitary transformation induced by a general point transformation, we simply employ an anticommutator to make the quantized generator Hermitian [52]:Ĝ = 1 2 {u(x),p}. Given a set of basis states Φ n (x), the transformed basis wave- To obtain a physical interpretation, it is instructive to differentiate both sides of Eq. (5) with respect to λ to get which has the form of a continuity equation of probability density. See Appendix A for the derivation details. Pictorially, starting from a family of orthonormal reference states such as Slater determinants, the probability mass of each many-body wavefunction will undergo a continuous evolution guided by the many-body vector field u. During this process, the particles constantly repel or attract each other and correlation effects are gradually cumulated. Crucially, these states remain orthonormal thanks to the unitary nature of the transformation. We thus obtain a correlated many-body basis in the end of the evolution, which can be used to build up the variational density matrix ansatz Eq. (3). Note the particles follow deterministic advection in Eq. (6) rather than random diffusion. In practice, we also integrate the differential equation for a finite amount of time instead of seeking for steady state solutions as diffusion-based approaches [53][54][55].
Transforming probability density continuously in the coordinate space is precisely the idea of continuous normalizing flow [39,54,56,57]. Specifically, the probability p n (z) ≡ |Φ n (z)| 2 associated with the reference state is known as the base distribution, while the more complex model distribution q n (x) ≡ |Ψ n (x)| 2 is obtained from the base by applying a learnable diffeomorphism f : z → x in the dN-dimensional coordinate space. Note we have omitted λ in the notation to avoid cluttering. In the present setting, one builds the diffeomorphism f via the ordinary differential equation (ODE) (4). By making use of the change-of-variable formula, the model probability q n (x) can be written as Taking the square root of both sides yields a more explicit expression for the basis wavefunction 4 : Albeit not so evident at the first sight, Eqs. (8) and (5) are completely equivalent, which can be rigorously proved by more formally establishing the unitary representation of point transformations; see Appendix A. In practice, the diffeomorphism f can be constructed by composing a sequence of point transformations, which is similar to the iterative backflow approach in ground state variational calculations [58,59]. However, an important difference of Eq. (8) from the ground-state backflow wavefunction ansatzes is the presence of a Jacobian determinant factor. This factor is crucial to guarantee orthonormality of the basis states, which is an essential ingredient for the present finite-temperature approach. Finally, as a many-fermion wavefunction, Ψ n (x) should satisfy the permutation antisymmetry property. Since this property holds already for the base Φ n (z), the only requirement is the unitary transformation appearing in Eq. (5) being permutation equivariant: that is, it should commute with the particle permutation operator. This can be achieved simply by requiring the many-body vector field u to be equivariant too, which means that the permutation of particle positions will result in the same permutation of the vector fields they experience: Intuitively, the indistinguishability of the particles is maintained throughout the continuous flow in the coordinate space, since one cannot label them by using the vector fields they experience at any time. The probability density q n (x) = |Ψ n (x)| 2 associated with the transformed wavefunction, on the other hand, is invariant under particle permutations [44][45][46]. A notable feature is that q n (x) inherits nodal lines from the fermionic reference state p n (z) = |Φ n (z)| 2 and, crucially, these nodal lines are deformed by the flow transformation.
To parametrize the permutation equivariant vector field u, one can leverage many recent advances in natural language processing [60], molecular simulation [44,61,62], and point set modeling [45,46]. Moreover, permutation equivariant functions have also been used in various ground-state VMC calculations [24,27,28]. Consequently, one can naturally port these efforts into the present framework almost without any modifications: just use the permutation equivariant layer as the vector field u to drive the flow.

III. IMPLEMENTATION
Substitution of the density matrix ansatz (3) into Eq. (2) yields the following estimator of the variational free energy: Notice the entropy term depends solely on the state occupation probability µ n and can be easily computed, which is a direct consequence of orthonormality of the basis states (8). The second term consists of the local energy associated with each basis state: In Eq. (10) the two-fold expectations correspond to classical thermal average of the Boltzmann distribution and quantum expectation according to the Born rule of wavefunction amplitudes. In the limit β → ∞, only the energy term survives and one naturally restores the ground-state VMC method. The gradients of Eq. (10) with respect to the parameters φ and θ, which appear in the classical and quantum distributions µ n and q n (x) respectively, have the following forms: For both estimators we employ the control variate method [36,38,63] to further reduce their variances. An important observation is that only the non-negative probability density q n (x) = |Ψ n (x)| 2 associated with the wavefunction is involved in the calculation. This is a satisfying feature of working directly in the continuum rather than on a finite basis set or lattice [29,64]: one can deal with the quantum manybody problem completely within the framework of probabilistic modeling. Nevertheless, the sign structure of the fermion wavefunction Ψ n (x) is still important and relevant for the calculation of off-diagonal physical observables such as correlation function and momentum distribution.
In practice, Eqs. (10) and (12) are estimated by sampling a batch of pairs (n, x) from the joint distribution µ n q n (x) following the ancestral sampling strategy. In particular, to sample coordinates x, one can start from samples z from the prior distribution p n (z) (e.g., via Markov chain Monte Carlo) and evolve them according to the ODE (4). The log-likelihood ln q n (x) appearing in the local energy (11) and gradient estimators (12) is evaluated by integrating Eq. (4) jointly with the following ODE [39,57]: To understand this, one can rewrite the continuity equation (6) in the form ∂ ∂λ + u(x) · ∇ ln q n (x, λ) = −∇·u(x) and note that  3. Energy E versus inverse temperature β for 6 spin-polarized electrons in a two-dimensional quantum dot with κ = 0.5. The green and blue points are benchmark data from two different variants of PIMC [68,69], while the red points are results of the present approach, including the zero-temperature limit. d dλ ≡ ∂ ∂λ + u(x) · ∇ is the material derivative associated with the sample x. Furthermore, the gradient and laplacian operations appearing in (11) and (12) can be accurately and efficiently computed by differentiating through the ODE integration using automatic differentiation [65], where the adjoint method with constant memory cost turns out to be useful [57]. Our code implementation based on PyTorch [66] is publicly available [67].

IV. APPLICATION: ELECTRONS IN TWO-DIMENSIONAL QUANTUM DOT
We demonstrate the capability of the present approach by studying electrons in a two-dimensional quantum dot. The one-and two-body potentials take the form of a harmonic trap and repulsive Coulomb interaction, respectively: where κ > 0 is the interaction strength. Despite its simplicity, this model shows rich phenomena due to the interplay of interaction and temperature effects. In particular, as κ increases, the Fermi liquid picture based on the concept of quasiparticles would eventually break down. The resulting phase is usually characterized as a Wigner molecule [2], where the kinetic motion of electrons is largely frozen, and the spatial density distribution would typically exhibit a shell structure. There have been a large number of numerical studies focusing on its ground-state [70][71][72] and finite-temperature properties [73][74][75][76]. However, there have been no reliable method that works for the entire interaction range at low temperatures. Thus, this problem offers an ideal playground for the present method. We consider the spin polarized case. The base wavefunctions Φ n (z) are chosen to be Slater determinants of single-electron orbitals obtained simply by eliminating the two-body term v (2) . Such Slater determinants constitute an exponentially large set of basis of the many-body Hilbert space. Focusing on low-temperature properties of the system, we carry out a truncation of the basis by including only those within an energy cutoff E cut relative to the noninteracting ground state. In the considered parameter region, we found that E cut 4 is sufficient to capture most of the finite-temperature effects.
Since the corresponding number of basis states is no more than 2000, we choose to adopt a simple parametrization of the state probabilities µ n (φ) = e φn m e φm based on the softmax function. Nevertheless, we note that this is not a limiting factor of the present approach because one can capture exponentially large number of basis states by utilizing more sophisticated discrete probabilistic models [38,[77][78][79].
The wavefunction ansatz Ψ n (x) is generated from the base Φ n (z) by the continuous flow guided by a many-body vector field u. We take u to be of the backflow form [21][22][23][24] for simplicity and clear physical interpretation: The many-body nature and permutation equivariance of this vector field can be easily confirmed by inspection. The scalar functions ξ and η can be referred to as the one-and two-body "backflow potential", which capture the "mean field" and electron correlation effects, respectively. Note for a given distance r, η(r) > 0 stands for a repulsive interaction between two electrons, and vice versa; similarly for ξ(r). We parametrize the potentials by two independent neural networks with single hidden layer. Initially, the backflow potentials are set to zero and µ n to Boltzmann distribution of the non-interacting base states Φ n (z). The optimization is performed on a batch of 8000 samples using the Adam stochastic gradient descent algorithm [80] for 3000 iteration steps. As the first benchmark, Figure 3 shows the temperature dependence of the energy for a system of N = 6 spin-polarized electrons with κ = 0.5. For such a weak interaction, the standard path integral Monte Carlo (PIMC) method is severely hindered by the fermion sign problem, since the electrons are largely delocalized and subject to exchange effects. Consequently, it provides reliable results only at relatively high temperatures β 1.5 [6]. Variants of PIMC with alleviated fermion sign problem can access slightly lower temperatures [68,69], where our results agree nicely with the benchmark data as shown in the figure. The slight discrepancy at β = 2 is likely due to insufficiently large E cut in our calculation. On the other hand, the present approach can easily reach even lower temperatures, including the zero-temperature limit β → ∞. We note that alternative Monte Carlo methods based on expansions in the Fock space [81,82] can work more favorably for such weak interactions, but will again suffer from the fermion sign problem in the strong coupling region.
Overall, the present approach serves as a valuable complement of conventional quantum Monte Carlo methods for studying the thermodynamic properties of fermion systems, especially for low temperature, large particle number and in- termediate interaction strength [83]. To demonstrate this, we perform systematic calculations of N = 10 spin-polarized electrons at β = 6 for a wide range of κ from 0.5 to 8. For more benchmark data, see Appendix B. Figure 4 shows the energy dependence on κ, together with electron density profiles at the two limits κ = 0.5 and 8. Notice the density centers around the origin of the trap in the weak coupling regime. On the other hand, stronger repulsive interaction smears out the electron cloud and induces a shell structure, indicating the emergence of the Wigner molecule phase. The observed spatial configuration consisting of two shells for the present parameter settings also agrees with the analysis in the classical limit κ → ∞ [2], where quantum fluctuations arising from the kinetic term in Eq. (1) are ignored. Reaching this result in the strong interaction regime where the density profile is qualitatively different from the weak coupling case is a stringent test to the present method.
To obtain the density profiles as shown in Fig. 4, one starts from the density of non-interacting reference state consisting of a large number of electron coordinate samples z, then evolves them according to the continuous flow specified by the ODE (4) towards the final spatial distribution of x, as described previously in Sec. III. The initial and final values of the continuous parameter λ are conventionally chosen to be 0 and 1, respectively, which are treated as fixed hyperparameters of the model. The many-body vector field u governing such an evolution process is determined by the backflow potentials ξ and η, which are shown in Figure 5. Notice the interactions arising from one-and two-body backflow potentials are both repulsive, which can be viewed as the manifestation of electron repulsion at the level of mean field and two-body correlations, respectively. The overall evolution of the electrons is, nevertheless, jointly determined by the two potentials together. Fig. 5 also shows that in the strong interaction regime, the backflow potentials deviate largely from the values of zero in the non-interacting case. To visualize how such

V. DISCUSSIONS
In essence, this work belongs to the large family of canonical transformation approaches for quantum many-body systems, except that the transformation is directly carried out upon particle coordinates instead of a many-body Fock space within the formulation of second quantization [84][85][86]. Although the framework of neural canonical transformation is general, we have technically restricted to the subgroup of point transformations to make the calculation tractable, which corresponds to the specific choiceĜ = 1 2 {u(x),p} of the generator. The generalization to arbitrary generators, such as the Hamiltonian H in Eq. (1), seems fairly nontrivial in practice. In fact, the tractability of such a general calculation implies that it is possible to accurately simulate the real-time evolution e −iHt |Φ of any many-body systems. In this perspective, the present approach can also be understood as a short-time variational approximation of adiabatic time evolution towards the thermal equilibrium.
One can take a different view of the limit of the present approach by inspecting the basis wavefunction representation Eq. (8). Compared to ground-state variational ansatzes, the coordinate transformation f in this work is implemented by continuous normalizing flow and subject to further limitations due to its invertibility [87][88][89][90]. Moreover, the continuous flow can only deform the nodal surface of reference states without changing its topology [24,91,92]. One example is the number of nodal cells of the ground-state wavefunction, which is conjectured to be always two for spatial dimensions d higher than one [92][93][94][95]. Similar conjecture has also been proposed for thermal density matrices at low temperatures [95]. More thorough characterization of fermion nodes like possible topological obstructions is still lacking and worths further study [96,97]. In practice, one may remedy these issues by increasing the expressibility of the reference states Φ n (z). For example, one can use more physically plausible reference states than the Slater determinants [98,99], or introduce additional parameters into the reference state which are pretrained or trained jointly with the flow transformation.
Although we have made use of continuous normalizing flow in this work, it should be possible to use other class of permutation equivariant normalizing flows [18]. Some examples are the partitioned flow [62] and the invertible residual network [100], which can be more efficient than the present ODE-based implementation. One can also directly carry out Monte Carlo sampling of the electron coordinates other than transforming samples of the reference states. The local energy Eq. (11) resembles the score matching loss function [101] in training generative models, which is known to be expensive to compute. In light of this, advances in efficient score matching training might be beneficial to further reduce the computational efforts when scaling up to larger systems [102].
Building on these technical improvements in the implementation, a promising future direction is to scale up to larger particle number N and study other correlated fermion systems of fundamental importance, such as the homogeneous electron gas [103] and dense hydrogen [10]. Moreover, similar to what was shown in [36], one can also obtain information about the low-lying excited states of these systems as a byproduct of the thermodynamic calculation.
FIG. 6. Several equally-spaced snapshots of the electron density along the continuous flow of coordinates at λ = 0, 1/8, . . . , 1. The system parameters are N = 10, β = 6 and κ = 8. Starting from the non-interacting reference state, the electron correlation effects can be gradually introduced, and one can finally reproduce the shell structure characteristic of a Wigner molecule in the strong interaction regime.
ables (x, p) to a new one (X(x, p), P(x, p)) satisfying the socalled symplectic condition, which is equivalent to saying that all Poisson brackets among the new variables are preserved [104]. To study its implication in the realm of quantum mechanics, we should convert the new (as well as old) variables into Hermitian operators (X, P) → (X,P) following certain quantization procedure. As a result, (X,P) should satisfy the usual commutation relations of coordinate and momenta as (x,p). This largely motivates us to reasonably expect the existence of a unitary transformation U that connects the two set of operators together: U can thus be viewed as the unitary representation of the given canonical transformation.
To further clarify the nature of U, it is instructive to consider the eigenstate |x) of the new coordinate operatorsX defined asX |x) = x|x). (A2) The essential point is that |x) constitutes a different coordinate basis from the old onesx|x = x|x , which is why a slightly different Dirac notation has been used in Eq. (A2). By making use of these two basis, the unitary transformation U can then be formally defined as [48,49] U|x) = |x .
It is then straightforward to verify the operator transformation relations Eq. (A1). The tranformation behavior of U upon the wavefunction Ψ(x) of any given quantum state reads The essential ingredient in evaluating this expression is the state overlap x |x), which can be shown to be uniquely determined (up to a phase factor) by the following set of equations [49,50]:X Unfortunately, for a general canonical transformation, the computation procedure outlined in Eqs. (A4) and (A5) above can be fairly difficult. We would thus restrict ourselves in the subgroup of point transformations, which can be constructed by the type-2 generating function F 2 (x, P) = f (x) · P via the following implicit relations [104]: The resulting transformation formula can be obtained by some simple manipulations: In other words, the new coordinates depend only on the old coordinates through a bijective map f , while the momenta are transformed in a covariant way so as to preserve the Poisson brackets. When such a point transformation is quantized, the corresponding expression for the state overlap x |x) turns out to be simple: The correctness of this result can of course be verified by plugging it into Eq. (A5). One can, however, get the intuitive feeling by noting thatX ≡ f (x) shares the same set of eigenstates asx, i.e., the Dirac delta function, while the additional Jacobian determinant factor is present to make the integration measure in the completeness relation dx|x)(x| = 1 as expected. By substituting (A8) into Eq. (A4), we finally obtain the transformed wavefunction through a point transformation as follows: This result is essentially Eq. (8) in the main text, where it is intuitively obtained from the perspective of normalizing flow.

Infinitesimal unitary transformations
The characteristics of the unitary representation of canonical transformations formulated above can be more clearly revealed by studying the infinitesimal behavior in the vicinity of identity transformation, which is a common practice in physics and of great theoretical importance. An infinitesimal canonical transformation can be constructed by the type-2 generating function F 2 (x, P) = x · P + dλ G(x, P), where the two terms correspond to identity transformation and infinitesimal perturbation, respectively. Substituting this expression into Eq. (A6) and retaining only lowest-order contributions, one could then obtain a continuous family (x(λ), p(λ)) of canonical transformations specified by the Hamilton's equations of motion: The function G(x, p) is usually also called the generating function. When such a canonical transformation is quantized, Eq. (A10) is naturally replaced by the Heisenberg equations of motion. This observation is essential: in light of Eq. (A1), the corresponding unitary transformation can be immediately recognized as where the generatorĜ should be ensured to be Hermitian by the quantization procedure. Equation (A11) is clearly a more explicit and meaningful characterization of the unitary transformation than the formal definition (A3). In particular, the state overlap appearing in Eq. (A4) corresponds precisely to the propagator (x|x = x|e −iĜλ |x . Such a propagator is difficult to evaluate in general cases, so we again concentrate only on the point transformations as in the previous section. Within the formulation presented here, a point transformation corresponds to the choice G(x, p) = u(x) · p, which leads to the type-2 generating function F 2 (x, P) = (x + dλ u(x)) · P. (A12) To obtain the corresponding unitary representation, we employ a simple operator symmetrizationĜ = 1 2 {u(x),p} to make the generator Hermitian, as mentioned in the main text. Comparing Eq. (A12) with the form F 2 (x, P) = f (x) · P discussed in the previous section, one can readily reach the conclusion that the transformed wavefunction (UΨ)(x) = x|e − i 2 {u(x),p}λ |Ψ on a quantum state Ψ(x) can be equivalently written in the form of Eq. (A9), where the coordinate bijection f is specified by the ordinary differential equation We thus rigorously show the equivalence of the two basis wavefunction ansatz expressions Eqs. (5) and (8), which lies at the core of the finite-temperature approach in this work. Finally, we give a few guidelines for the derivation of the continuity equation (6) from Eq. (5) in the main text for read-ers' convenience. Note thatp = −i∇, we have ∂ ∂λ Ψ n (x, λ) = − i 2 {u(x),p} Ψ n (x, λ) = − 1 2 u(x) · ∇Ψ n (x, λ) + ∇ · (u(x)Ψ n (x, λ)) = −u(x) · ∇Ψ n (x, λ) − 1 2 Ψ n (x, λ)∇ · u(x). (A14) To obtain Eq. (6), simply multiply Eq. (A14) by Ψ * n (x, λ) and add the resulting equation to its own complex conjugate.
Appendix B: Some more benchmark data for 2D quantum dot The following table summarizes our results for the energy of a two-dimensional quantum dot at β = 10, for various electron number N and interaction strength κ. PIMC results from [74] are also listed when available. All data correspond to the fully spin-polarized case. Our finite-temperature calculations indicate that the entropy is negligible for a temperature as low as β = 10, so our energy results can be treated as variational. We anticipate these results (as well as those presented in the main text) can be further improved by adopting better model architecture and optimization schemes. We also note the results reported in Figure 4.8 of [83] for β = 10, N = 3, κ = 2 show that the data in [74] (2)