Machine Learning Universal Bosonic Functionals

The one-body reduced density matrix $\gamma$ plays a fundamental role in describing and predicting quantum features of bosonic systems, such as Bose-Einstein condensation. The recently proposed reduced density matrix functional theory for bosonic ground states establishes the existence of a universal functional $\mathcal{F}[\gamma]$ that recovers quantum correlations exactly. Based on a novel decomposition of $\gamma$, we have developed a method to design reliable approximations for such universal functionals: our results suggest that for translational invariant systems the constrained search approach of functional theories can be transformed into an unconstrained problem through a parametrization of an Euclidian space. This simplification of the search approach allows us to use standard machine-learning methods to perform a quite efficient computation of both $\mathcal{F}[\gamma]$ and its functional derivative. For the Bose-Hubbard model, we present a comparison between our approach and Quantum Monte Carlo.

The one-body reduced density matrix γ plays a fundamental role in describing and predicting quantum features of bosonic systems, such as Bose-Einstein condensation. The recently proposed reduced density matrix functional theory for bosonic ground states establishes the existence of a universal functional F[γ] that recovers quantum correlations exactly. Based on a novel decomposition of γ, we have developed a method to design reliable approximations for such universal functionals: our results suggest that for translational invariant systems the constrained search approach of functional theories can be transformed into an unconstrained problem through a parametrization of an Euclidian space. This simplification of the search approach allows us to use standard machinelearning methods to perform a quite efficient computation of both F[γ] and its functional derivative. For the Bose-Hubbard model, we present a comparison between our approach and Quantum Monte Carlo.
In 1964 Hohenberg and Kohn proved the existence of a universal functional F[ρ] of the particle density ρ, that captures the exact electronic contribution to the groundstate energy of a system of interacting electrons [1]. Due to a remarkable balance of accuracy and computational cost, first principle modeling of electronic systems based on the respective Density Functional Theory (DFT) is nowadays a well established daily practice, with great impact in material science, quantum chemistry or condensed matter [2]. For bosonic systems, however, a fully first-principle description has been elusive. This is due in part to the unsuitability of the particle density to describe fundamental bosonic features as orbital occupations, mode entanglement, or non-diagonal order, which are important for predicting and describing bosonic condensation. As a result, the theoretical treatment of interacting bosonic systems mainly relies on exact diagonalization techniques, which are restricted to few tens of orbitals [3][4][5], or mean-field theories that are particularly suitable for dilute ultracold gases [6][7][8]. Quantum Monte Carlo (QMC) is known to be a powerful family of techniques for computing ground-state properties but is still restricted due to the fermion sign problem. While bosonic systems do not suffer such a sign problem, QMC cannot be applied to, e.g., frustrated quantum spin systems [9]. Nowadays, a renewed interest in the ab-initio description of many-body systems has been motivated by the successful application of artificial neural networks to both fermionic and bosonic problems [10,11].
Since the parameters of correlated bosonic systems (ultracold gases, in particular) can be tuned with a high degree of control, they are powerful platforms to study a wide range of model Hamiltonians, ranging from Hubbard models to bosonic antiferromagnets [12,13]. They * carlosbe@pks.mpg.de have also become an active research field in the context of quantum simulations [14][15][16], and even quantum foundations [17][18][19][20]. Such need to describe quantum correlations of bosonic systems efficiently has motivated quite recently to put forward a novel physical theory for interacting bosonic systems [21,22]. Based on a generalization of the Hohenberg-Kohn theorem [23,24], this reduced density matrix functional theory (RDMFT) for bosons establishes the existence of a universal functional F W [γ] of the one-body reduced density matrix (1RDM): γ ≡ N Tr N −1 [Γ], obtained from the N -boson density operator Γ by integrating out all except one boson, and the two-particle interactionŴ . Since the 1RDM is the natural variable of the theory, RDMFT is particularly wellsuited for the accurate description of Bose-Einstein condensates (BEC), strongly correlated bosonic systems, or fragmented BEC [25,26]. Furthermore, the information contained in the spectra of the 1RDM can also be sufficient to investigate multipartite quantum correlations in those systems [27][28][29][30][31].
Although RDMFT holds the promise of abandoning the complex N -particle wave function as the central object, it does not trivialize the ground-state problem. In fact, the fundamental challenge is to provide reliable approximations to the universal interaction functional F W [γ]. Yet, while the Hohenberg-Kohn-type fundational theorem of RDMFT shows the existence of a universal functional, it does not give any indication of its concrete form. For DFT, the solution to this problem is given in the form of large classes of approximate functionals, hierarchically organized in the so-called Jacob's ladder. In recent years, the number of such approximation has significantly increased thanks to machine learning [32][33][34][35][36][37] and reduced density matrices [38] approaches.
Our work succeeds in providing a strategy on computing approximations for F W [γ]. In this paper we (i) provide an efficient method to capture the essential features arXiv:2104.03208v3 [quant-ph] 25 Aug 2021 of universal functionals for boson lattices, (ii) show how the constrained search approach associated with it can be simplified in the form of an unconstrained problem, and (iii) implement this approach in a standard machinelearning library to compute F W [γ], its derivative, and the ground-state energy. We shall for simplicity describe our method for the Bose-Hubbard model, but the same results apply to any type of interactions for systems with translational symmetry.
Universal bosonic functionals.-In this work we consider Hamiltonians of the form with a one-particle termĥ =t +v, containing the kinetic energy and the external potential terms, and the two-particle interactionŴ . The ground-state energy and 1RDM follow for any choice of the one-particle where Γ → γ indicates that the minimization is carried out over all Γ whose 1RDM is γ. The main challenge of this approach is that the set of Γ such that Γ → γ is in general extremely complex to characterize, and so far only partial results are known for quasi-extremal, twoparticle or translational invariant fermionic systems [40][41][42][43][44]. Even in the extremely popular DFT, the constrained search over many-body wave functions integrating to the same electronic density (i.e., Ψ → ρ) is rarely explicitly carried out [45].
To shed some light on the problem let us represent γ, the 1RDM of a N -boson real wave function |Ψ , with respect to a set of creation and annihilation operators and assume that the dimension of the one-particle Hilbert space is M . Let us also define M (N − 1)-particle wave functions |Φ j ≡b j |Ψ , which satisfy by definition the condition The meaning of these non-normalized wave functions is clear: while their magnitude equals the diagonal entries of γ (i.e., Φ j |Φ j = γ jj ), the angles they form correspond to the non-diagonal entries of γ. Indeed, since Φ i |Φ j = ||Φ i ||||Φ j || cos(θ ij ) = √ γ ii γ jj cos(θ ij ) we have The bound of the non-diagonal entries: |γ ij | 2 ≤ γ ii γ jj , is the Cauchy-Schwarz inequality for operators, and known to be a representability condition for γ [46]. The condition jb † j |Φ j = N |Ψ suggests that the minimizer of the minimization (2) can be written as a set of M vectors in the Hilbert space H N −1 of N − 1 particles, such that their angles and magnitudes are determined by Eq. (4) (see Fig. 1). We now exploit this first insight to explicitly carry out the constrained search approach and find the universal functional of the Bose-Hubbard model, a workhorse in the context of ultracold bosonic atoms [47].
Bose-Hubbard model.-The Hamiltonian of the Bose-Hubbard model reads: where the operatorb † j (b j ) create (annihilate) a boson on site j, andn j is the corresponding number operator. The first term in Eq. (6) describes the hopping between two sites while the second one is the interacting termŴ = U 2 jn j (n j − 1). Since the problem is determined by N spinless bosons and M sites, we write for the functional F N,M [γ]. For a given γ let us take the minimizer of the functional (2) and call it |Ψ γ ∈ H N , the N -particle Hilbert space. Using the prescription discussed above let us define M (N − 1)-particle wave functions |Φ γ,j ≡b j |Ψ γ ∈ H N −1 , which satisfy by definition the condition (4). Due to the translational invariance of the Bose-Hubbard Hamiltonian (6), these wave functions are all normalized to the filling factor, namely, Φ γ,i |Φ γ,i = N/M . The functional is given by Figure 1. Representation of 3 wave functions in the Hilbert space of N − 1 particles giving place to a 1RDM γ. While the magnitude of the vectors is ||Φi|| 2 = γii, the angles they form satisfy Φi|Φj = √ γiiγjj cos(θij).
As shown in the Appendix B, any rotation of the states |Φ γ,i in the subspace spanned by themselves: G γ = span{|Φ γ,1 , . . . , |Φ γ,M }, will give an energy greater or equal than the energy F N,M [γ]. As a consequence, we rewrite the constrained search approach in Eq. (2) as subject to Φ i |Φ i = N/M and Φ i |Φ j = γ ij . This indicates that the constraint in Eq. (2) can be transferred to the subspace G γ . As we will see below, this result leads to a quite efficient optimization problem for the functional. Exact functional of the dimer.-As a first illustration of this novel approach let us take the simple case of the Bose-Hubbard dimer with two particles (N = M = 2). The states of the Hilbert space can be written as two occupations: |n L , n R . For the 1-boson Hilbert space we choose as a basis: . We write these two wave functions as |Φ L = sin(θ L )|1, 0 + cos(θ L )|0, 1 and |Φ R = cos(θ R )|1, 0 + sin(θ R )|0, 1 . As a result of the corresponding minimization (7), the three angles are related: θ L = θ R = (π/2 − θ)/2, and the universal functional equals to: which is one of the few analytical results for a universal functional that can be found in the literature [48,49]. Machine learning.-Despite the spectacular rise of machine learning in the study of quantum many-body systems, no implementation is known so far for the theory of reduced density matrices [50,51]. One of the reasons for this lack of progress is the large amount of constraints swarming in functional theories. We now discuss how our findings will facilitate learning the universal functional of bosonic systems. Notice that a more appealing way of writing the functional in Eq. (7) is the following: Let us choose a basis for the vector space G γ , say: {|m ∈ H N −1 }. A set of wave functions of the sort needed in the minimization (7) can now be written as |Φ j = d jm |m . The condition of Eq. (4) reads: dd † = γ, where we have defined the matrix: Using the singular value decomposition for such a matrix we have d = UΣV † , with U and V M × M unitary matrices. Since ΣΣ † = U † γU, it is now clear that the spectral decomposition of γ equals to ΣΣ T . As a consequence, we obtain where {n α } are the eigenvalues of γ. Collecting these results we obtain that the exact universal functional (7) can be explicitly written in terms of the eigenvalues and the eigenvectors of γ (contained in the matrix U): where ∆ αβ (U, V) = i,mm u * iα u iβ v * mα v m β m |n i |m . The concrete form of the functional presented in Eq. (9) is striking: for fermionic density-matrix-functional theory, the square root of the occupation numbers in Eq. (9) is known to be the optimal choice for Ansätze of the form n α i n 1−α j , compatible with the integral relation between the one-and two-body reduced density matrices [52][53][54][55], even for systems out of equilibirum [56]. As we can see, the only freedom in the functional (9) is the matrix V, which is, unlike U and n α , not fixed by γ. We use this degree of freedom to introduce a standard optimization problem on a connected manifold M: where we have included a sub-index in U γ to remember that such a matrix is defined by γ. Notice that the manifold M is essentially the set of special orthogonal matrices of dimension M , which generates the space G γ . Although the definition of such a space is far from trivial (and we will leave this question open for future research), it is possible to establish some elementary facts. For instance, in the strongly correlation regime U/t 1 with integer filling factor α = N/M , G γ = span{b i |α, ..., α }.
To make further progress on our problem, notice that optimization problems of the form min x∈M f (x) over a connected manifold M can be transformed into an unconstrained one of the form min y∈R n f (φ(y)) by lifting the function f to the current tangent space T x M [57,58]. The map φ : R n → M is called a trivialization map [57]. By letting the minimization in Eq. (10) to run over the set of special orthogonal matrices in dimension M , the relevant minimization space turns out to be an Euclidian space R M . As a result, finding the universal functional of RDMFT presents itself as an unconstrained minimization problem. This is the crucial and last finding of our work, as it finally allows us to compute the universal bosonic functional by solving the problem directly over the set of orthonormal matrices.
Modern machine learning frameworks like pytorch [59] provide fast and rather efficient ways of performing optimizations on connected manifolds of the type we consider here. For the results we will present below, we have implemented the constrained minimization (10) in pytorch with the constrained minimization toolkit GeoTorch [60]. As a first step we implemented an minimization procedure where for each 1RDM the matrix V in Eq. (10) is optimized to produce the universal functional. As a second step we trained a neural network as the universal bosonic functional (see below).
Results.-In Fig. 6 we present the results for the Bose-Hubbard model (6)  been normalized to 0 in the lower point (i.e., η = 0) and to 1 in the upper point (η = 1). The exact known results for M = 2 in Eq. (8) are verified in our calculations. Furthermore, we observe the existence of the Bose-Einstein force discovered in [21], extended in [22] and proved in [61], namely, the divergence of the gradient In functional theories the knowledge of the functional's form is as important as the knowledge of its derivative, as both are needed for a ground state calculation. To perform the derivative of the functional we trained a neural network to output the matrix V using the degrees of freedom of our 1RDM as input. This has multiple advantages. First, once the functional is trained for given particle and site numbers, it can be evaluated for any γ. Secondly, the automatic differentiation allows an exact evaluation of the gradient ∇ γ F N,M [γ] without further work. For the Bose-Hubbard dimer it was sufficient to use the diagonal terms η = γ i(i+1) and its square as inputs. The calculation was structured as follows: Here FCNN N,M,θ is a fully connected network for N par- ticles and M sites with the parameters θ and Remarkably, as shown in Fig. 3, for the Bose-Hubbard dimer (for which we can compare to exact results) the derivatives provided by the neural network only deviate by ∼ 10 −13 from the exact results.
We demonstrate now that our approach allows us to compute the ground-state energy and 1RDM for a large system. Notice first that for a fixed filling factor α = N/M , the energy of the ground state of the Bose-Hubbard model (6) can be computed as the minimum of the energy functional E N, . By performing the minimization ∇ γ E N,M [γ] = 0 on the domain of positive semidefinite matrices, it is then possible to compute the ground-state energy of the system. To generate the functional in that domain, we have optimized the ansatz γ ij = η κ α (0 ≤ η ≤ 1) with 2 ≤ κ ≤ 8, for |j − i| > 1, with η = γ i(i+1) . This is motivated by the fact that when U/t 0, γ ij ≈ 0, ∀i > j, and when U/t 1 γ ij ≈ α, ∀ij. Following this prescription we have computed the ground-state energy for   Fig. (a) the energy is plotted as a function of the relative strength U/t. In blue the results computed using the functional (10) (RDMFT1) and the ansatz |Φi = bi|1, ..., 1 . In black we use as an ansatz for the functional the exact expression of the dimer appropriately rescaled (RDMFT2). In red the results computed with QMC. In the subfigure the functionals FN,N (η) are plotted as a function of η, the nearest-neighbor off-diagonal value of the 1RDM. In Fig. (b) the relative errors ∆/EQMC , where ∆ = ERDMF T 2 − EQMC are plotted as a function of U/t. the 40-site Bose-Hubbard model with 40 bosons. The dimension of the corresponding Hilbert spaces, being of the order of 5.3 × 10 22 , is out of reach for exact diagonalization and prohibits performing the exact constrained search approach. To solve this problem we have (i) ansatzen the space G γ , the subspace generated by the kets |Φ i =b i |Ψ , by choosing |Ψ = |1, ..., 1 (RDMFT1), and (ii) used the exact functional of the dimer (8), appropriately rescaled (RDMFT2). The energy predictions of our machine-learning functionals are quite remarkable, given the subspaces we have chosen. Indeed, the results presented in Fig. 4 indicate that the predicted RDMFT results are in good agreement with the QMC energies: the errors around U/t = 4 are only due to the approximation of the space G γ . In order to check the quality of the approximate functionals more, we have also plotted the relative error in the last panel. We observe that this error is below 8% and practically zero for large and weak interaction. In addition, notice that our implementation is able to approximate the whole range of energies, not only the weakly (the sector easily described by Bogoliubov methods) or the strongly correlation regimens.
Conclusion.-In conclusion, we have demonstrated the viability of approximating universal bosonic functionals in a quite efficient way. The main ingredient of the computation is a simplification of the constrained search approach that we have introduced in this work based on the Schmidt decomposition of the wave function. This formulation of reduced density matrix functional theory (RDMFT) speeds up the design of reliable approximations for the universal functionals for systems with translational symmetry. The quality of the numerical results obtained in this work highlights the potential of RDMFT to become a competitive tool for computing properties of bosonic ground states with large dimensional Hilbert spaces. Strikingly, since RDMFT takes into account the whole range of bosonic correlations, and does not present dimensional or sign problems, it offers a range of new possibilities. For instance, frustrated bosonic systems can be studied in a direct manner [63]. Bosonic systems with impurities, composites of ultra-cold atoms, or even superconducting systems [64][65][66][67] can also potentially be addressed within this framework. As an outlook of this work, we leave open a new line of research based on extending our findings to systems with internal degrees of freedom, finite temperatures, or broken symmetries. We also expect that previous works in the context of twobody reduced density matrix [68,69] will also benefit from our approach.

ACKNOWLEDGMENTS
We thank Jonathan Siegel, Matt Eiles, Adam Sawicki, and Jakob Wolff for helpful discussions. We are most grateful to Peter Karpov for constructive feedback and insight, and for providing us the QMC energies of the Bose-Hubbard models. M. F. was partially supported by the Research Fund of the University of Basel for Excellent Junior Researchers. C. L. B.-R. was supported by the MPI-PKS through a next-step fellowship.
All codes to reproduce, examine and improve our proposed analysis will be made freely available online upon publication.
Appendix A: The γ-representability problem An important problem in the theory of reduced density matrices for indistinguishable particles is the so-called Nrepresentability problem, namely, which conditions should γ, the one-body reduced density matrix (1RDM), satisfy in order to belong to at least one wave function in H N , the Hilbert space of N particles (fermions or bosons). To understand the problem, let us consider for a given wave function |Ψ ∈ H N the corresponding 1RDM as: We have introduced a one-particle basis set |f i = f † i |0 that determines the corresponding set of creation and anhilitation operators. The matrix γ has the following properties: (a) it satisfies the trace condition: i γ Ψ ii = N , (b) it is hermitian, and (c) it is positive semidefinite (i.e., its eigenvalues are non-negative). For fermions the Pauli exclusion principle imposes another constraint: γ ≤ 1 [70], and the generalized Pauli principle imposes even stronger constraints on the eigenvalues [71,72]. Another less explored problem is the following: what is the set of wave functions giving place to the same 1RDM, namely, what is the set of wave functions in the Hilbert space H N ? The characterization of this set is of crucial importance for the ground state problem as seen from the point of view of reduced density matrix functional theory (RDMFT) [24]. Indeed, it allows us to find the universal functional of a two-particle interaction W defined as: As a consequence of this construction, the ground state of a system of indistinguishable particles driven by a Hamiltonian H(h) could be computed in a quite simple way by resorting only to the set of 1RDMs [39]: where h contains all the 1-particle contributions to the full Hamiltonian H(h) = h + W .

Appendix B: The Bose-Hubbard Hamiltonian
For clarity, we will work out the Bose-Hubbard model, but the results can be generalized for systems with translational symmetry. The Hamiltonian of the problem is given by: The two-particle interaction is W = U 2 in i (n i − 1) and the 1-particle Hamiltonian is h = −t i (b † i b i+1 + h.c.). We fix N and M , the number of particles and the number of sites. The filling factor is defined by α = N/M . Notice that the Hamiltonian is translational invariant. We further choose periodic boundary conditions.

Ground-state problem with reduced density matrices
The ground-state energy satisfies by definition: E 0 = min Ψ∈H N Ψ|H|Ψ . For a given wave function |Ψ the expected value of the corresponding Hamiltonian reads: Yet, since the Hamiltonian (B1) is translational invariant, we have in general: γ Ψ ij = γ Ψ (i+m)(j+m) , for m ≥ 0. As a result, the ground state problem is defined by M/2+1 parameters (for an even number of M ), namely, γ ii = α, γ i(i+1) = η 1 , .... The 1RDM of the real ground state can be written as a M × M matrix: While the entries γ i(i+m) = η m for all m > 1 are not relevant for the computation of ground-state energy, since only nearest-neighbour hopping appears in the Hamiltonian (B1), they are crucial for the ground-state minimization of the energy functional. Furthermore, the entries γ ij contain information about the entanglement of the modes b i and b j . Indeed, in the Mott phase γ ij ≈ 0 for i = j. Bose-Einstein condensation appears when γ ij ≈ α for i = j. On this regard, the entries of γ contain relevant physical information of the problem, and impose also a challenge to the theory because the entries should be such that γ is positive semidefinite. As a matter of fact, the ground-state energy can be computed by minimizing the following functional: The universal functional reads as in Eq. (A3). Of course if we knew the expression for the functional F W [γ] the problem would be remarkably simple: the minimum in (B4) could be found by simply computing the derivatives where the subindex α means that it is fixed during the minimization. The minimizers η * i satisfy: and ∂ ∂ηi F W [γ(α)] = 0. The ground-state energy is then given by The main question now is if it is possible to find competitive approximations to the functional. At first sight, it seems an impossible task as it would imply the disregard of an enormous Hilbert space, whose scaling is exponential. We will show now that such an approach is feasible.

Universal functional for the Bose-Hubbard model
We tackle the problem in the following way: Let us take the minimizer of the functional (A3) and call it |Ψ γ ∈ H N . Notice first that |Ψ γ gives place to M wave functions in H N −1 , the Hilbert space of N − 1 particles: where b i is the annihilation operators of the Bose-Hubbard Hamiltonian (B1). These wave functions satisfy i b † i |Φ γ,i = N |Ψ γ , Φ γ,i |Φ γ,i = α, the diagonal of γ, and Φ γ,i |Φ γ,i+1 = η. The two-particle energy of the minimizer is given by where we have usedn i (n i − 1) = b † in i b i . Now we state the following: any rotation of the states |Φ γ,1 , ..., |Φ γ,M in the space spanned by themselves G γ = span{|Φ γ,1 , . . . , |Φ γ,M } will give an energy greater than or equal to Ψ γ |W |Ψ γ . In other words, for any set of states |Φ 1 , ..., |Φ M that result from a rotation of the frame |Φ γ,1 , ..., |Φ γ,M in G γ . To understand the assertion, let us study in detail the Bose-Hubbard dimer with an arbitrary number N of bosons. In such a case we have only two wave functions, namely, |Φ γ,1 and |Φ γ,2 , with Φ γ,1 |Φ γ,1 = Φ γ,2 |Φ γ,2 = α and Φ γ,1 |Φ γ,2 = αη . A rotation of those vectors is defined as where |Φ ⊥ γ,i are orthogonal wave functions to |Φ γ,i on G γ defined by: Hence, taking real wave functions for simplicity, We now develop independently these terms. The second term in the last line of Eq. (B13) can be rewritten as follows: In the third line we have used the fact that |Φ γ,i is an eigenfunction of the operatorn 1 +n 2 . The last term of Eq. (B13) can also be developed: which is zero, due to translational symmetry of the ground state (e.g., Ψ γ |n 2 1 |Ψ γ = Ψ γ |n 2 2 |Ψ γ ). Finally, since the maximum value of the functional F W [γ] is N (N − 1)/2 (see supplemental material of Ref. [21]) and α = N/2, we have > 0, and therefore: for 0 < θ < π, which is what we wanted to prove. A second meaningful example is the Mott phase. For t = 0 and integer filling factor α = N/M the ground state is |α, α, ... . We then have |Φ i = √ α|α, ..., α − 1, ..., α and F W [γ] = M (α − 1)α. A rotation of any pair of those vectors, e.g., |Φ i (θ) = cos(θ)|Φ i + sin(θ)|Φ j , |Φ j (θ) = cos(θ)|Φ j − sin(θ)|Φ i and |Φ k (θ) = |Φ k , for k = i, j, results in the new energy (i = j and j = i): This result allows us to change the constrained search approach in (A3) by the following more appealing unconstrained functional in the Hilbert space H N −1 : While there is still a representability constraint in G γ that cannot be lifted, this construction will facilitate the design and training of a neural network as the universal functional of bosonic RDMFT. Before showing this, we employ our novel approach to explicitly compute the universal functional of the Bose-Hubbard dimer with 2 bosons, which is one of the few (or perhaps the only) analytical result for the universal bosonic functional existing in the literature.
Therefore Φ L |Φ R = cos(θ L ) sin(θ R ) + sin(θ L ) cos(θ R ). By defining γ LR = cos(θ), we have sin(θ L + θ R ) = cos(θ). As a result, θ L + θ R + θ = π/2, and sin(θ R ) = cos(θ L + θ). Our functional (B21) then reads: The minimum is attached when dJ (θ, θ L )/dθ L = 0. An elementary calculation gives as a solution sin 2 (2θ * L ) = cos 2 (θ) for the minimizer θ * L , which results in F(θ) = J (θ, θ * L ) = sin 2 (θ * L ) + cos 2 (θ) cos 2 (θ * L ) + sin 2 (θ * L ) sin 2 (θ) − cos 3 (θ) sin(θ) which is the result found several times in the literature for bosons and fermions [21,48,49,73]. We compare the function in Eq. (B23) w.r.t. the solution of Eq. (B24) in Fig. 6 for θ = 2π/7. There is still a quicker way of computing the same result. In H 1 the number operators can be written asn L = |10 10| andn R = |01 01|. Hence, it is obvious that the minimum of (B23) is reached when θ L = θ R = (π/2 − θ)/2 and the value of the functional is just 2 sin 2 [(π/2 − θ)/2] = 1 − sin(θ). Hence, Therefore DD † = U † γU. As a consequence, U diagonalizes the matrix γ, and DD T is a diagonal matrix with the eigenvalues of γ in the entries. Therefore the diagonal entries of D are √ n 1 , √ n 2 , ..., √ n M , the square root of the eigenvalues of γ. The kets in Eq. (B25) can then be explicitly written as functions of γ: The connection with the original wave function is striking: whereb † α = i b † i u iα and |v α = m v * m,α |m . Notice that the 1RDM γ fixesb † α and √ n α , but not |v α . Quite remarkable, we can recognize in the expression (B29) the well-known Schmidt decomposition of a bipartite system: in this case the system of 1 and (N − 1) indistinguishable particles. Such a decomposition can be used to study fermionic entanglement as quantum resource [74,75] or to describe the correlated electron dynamics in strong laser fields within exact factorization of the wave function [76]. As its is discussed in Section C, the matrix V is used to engineer the bosonic functionals by training a neural network to produce V as a function of γ.