Boundary conditions for star matter and other periodic fermionic systems
Francesca Gulminelli, Takuya Furuta, Olivier Juillet, Christian Leclercq
aa r X i v : . [ nu c l - t h ] O c t November 21, 2018
Boundary conditions for star matter and other periodic fermionic systems
F. Gulminelli , , , T.Furuta , , , , O.Juillet , , , C.Leclercq , , ENSICAEN, UMR6534, LPC, F-14050 Caen, France Universit´e de Caen-Basse Normandie, UMR6534, LPC, F-14032 Caen, France CNRS, UMR6534, LPC, F-14050 Caen, France present address: RIKEN Nishina Center, RIKEN, Wako, Japan Bulk fermionic matter, as it can be notably found in supernova matter and neutrons stars, issubject to correlations of infinite range due to the antisymmetrisation of the N-body wave function,which cannot be explicitly accounted for in a practical simulation. This problem is usually addressedin condensed matter physics by means of the so-called Twist Averaged Boundary Condition method.A different ansatz based on the localized Wannier representation has been proposed in the contextof antisymmetrized molecular dynamics. In this paper we work out the formal relation between thetwo approaches. We show that, while the two coincide when working with exact eigenstates of the N-body Hamiltonian, differences appear in the case of variational approaches, which are currently usedfor the description of stellar matter. Some model applications with Fermionic Molecular Dynamicsare shown.
PACS numbers: 21.60.-n 26.60.-c 71.10.Ca
I. INTRODUCTION
Interacting fermionic systems in the bulk limit are a standard object of theoretical study in condensed matterphysics. Electrons are subject to an external periodic potential in the presence of a cristalline ionic structure andthe bulk limit can be seen as an infinite number of spatial replicas of a finite system within a specific geometry[1].In that case, the observables of the bulk system can be obtained from the modelization of one single elementary cell,provided adequate boundary conditions are applied to the many-body wave function. This amounts to introduce aBloch phase or twist to each wave function in the single-particle basis, and average the twisted observables over thedifferent phases within the first Brillouin zone[2, 3]. In practical applications, this technique has been applied toQuantum Monte-Carlo simulations of electron systems also in the absence of any external periodic potential[2–4]. Inthis case, the introduction of Bloch phases has to be understood as a technique to accelerate the convergence towardsthe thermodynamic limit of these very expensive numerical calculations, which would otherwise become prohibitivein computation time. In the absence of an external potential, the periodic cell is just a computation cell with nophysical meaning, and independence of the results respect to its size has to be checked.Coming to the strongly interacting fermionic systems studied in nuclear physics, the bulk limit has not attractedmuch interest in the community since nuclei are finite. However there are physical situations where the propertisof fermionic matter composed of protons and neutrons in the thermodynamic limit are essential such as for core-collapsing supernova, and for the crust of the neutron stars which are left over by the explosion. This nucleonic stellarmatter covers a very wide domain of densities ranging from ρ ≈ g · cm − to a few times the normal saturationnuclear density ρ ≈ g · cm − , temperatures between less than 1 and more than 20 MeV, and proton fractionsvarying between 0. and 0.5 In the sub-saturation density regime, it is well established that matter is charge neutraland mainly composed of neutrons, protons, electrons, positrons and photons in thermal and typically also chemicalequilibrium [5, 6]. Depending on the thermodynamic condition, neutrinos and anti-neutrinos can also participate tothe equilibrium.A very large amount of literature exists on the microscopic modelizations of the neutron star outer and innercrust[7–12]. In this regime at zero temperature and relatively low density, matter consists of a lattice of Wigner-Seitzcells, each cell containing a spherical neutron-rich nucleus immersed in a sea of dilute gas of neutrons and relativisticelectrons uniformly distributed inside the cell[13, 14]. The linear size of the cell is of several hundreds of fermis inthe outer crust. The size decreases going towards the center of the neutron star and the central nucleus becomesheavier and increasingly neutron rich. Typical values[13] in the inner crust range from about 200 particles in a cellof linear size of about 50 fm for ρ ≈ − f m − to about 1500 particles in a cell of linear size of about 15 fm for ρ ≈ · − f m − . Most of the existing calculations are based on variational approaches (Hartree-Fock or Hartree-Fock-Bogoliubov) and employ mixed Dirichlet-Neumann boundary conditions at the edge of the cell[13], chosen toproduce a flat density close to the cell border and thus to simulate a uniform neutron gas. The specific way of fixingthese mixed boundary conditions is not completely clear, and discrepancies due to the choice of fixing these conditionsincrease with density[15]. A more conceptual problem with Dirichlet or Neumann boundary conditions, is that theyneglect antisymmetrization correlations which extend beyond the cell size. This has been pointed out in refs.[16–19],where Bloch boundary conditions, similar to the ones used for the QMC modelizations of bulk electron systems, havebeen employed. In these works it was shown that both the neutron specific heat and the motion of unbound neutronsare affected by these effects.In the intermediate region between crust and core, complex phases are predicted that can break the sphericalsymmetry of the Wigner-Seitz cell, and can violate the associated translational invariance. This is even more true atfinite temperature, where the periodicity of the Wigner-Seitz cell is broken by thermal agitation. These conditions aremet in stellar matter in the pre- and post-bounce supernova dynamics, as well as in the cooling process of proto-neutronstars. Even when periodicity cannot rigourously be assumed in these thermodynamic conditions, it can be kept as apractical working hypothesis allowing to address the thermodynamic limit. As already mentioned, the drawback ofthat is that convergence with respect to the cell size has to be systematically checked. Time-dependent variationalmicroscopic calculations have been proposed to address this region, where statistical averages are calculated fromtime averages assuming ergodicity[20–23]. These works have shown that structures, though they can be degeneratein energy[23], are nonetheless approximately periodic in space even at finite temperature.In these calculations simple periodic boundary conditions are employed, completely neglecting the antisymmetriza-tion correlations beyond the calculation grid. The importance of properly accounting for these correlations wasrecently stressed in ref.[24]. In this work, a specific ansatz for the boundary conditions based on the localized Wannierrepresentation has been proposed in the context of antisymmetrized molecular dynamics. It was shown that suchboundary conditions allow obtaining a distribution similar to the one of a free Fermi gas, if gaussian wave packets offixed width are periodically disposed on a two-dimensional grid, while an artificial Pauli potential is needed in orderto obtain the same result with classical molecular dynamics. This result implies that, in the inner crust region wherestellar matter contains an important component of quasi-free neutrons, properly accounting for the periodic characterof the system may be of importance.To conclude these introductory remarks, it appears that simple Dirichlet, Neumann, or periodic boundary conditionsare not adapted to the description of bulk fermionic matter. Different solutions are proposed in the literature fordifferent applications, namely the Twist Averaged Boundary Conditions (TABC) in QMC[2, 3], the Bloch methodin mean-field calculations[16–19], and the Wannier replica method for molecular dynamics approaches[24], but theequivalence of the different techniques and/or their domain of validity is not completely clear.In this paper, we will formally develop the link between the different methods and show some model applicationsfor simple non-interacting nuclear systems described through the variational Fermionic molecular Dynamic (FMD)method. We will show that Bloch (or TABC) and the replica method are equivalent when they are applied to the exacteigenstates of the many-body Hamiltonian. When this is not the case, as for variational mean-field theories appliedto interacting systems, the equivalence is broken. In this case, the Bloch or TABC technique can produce solutionswhich differ from the exact results more than if simple periodic boundary conditions are applied. Conversely, thereplica method appears more powerful and can easily be applied to any variational based mean-field treatment withan affordable extra computational cost. Our numerical applications will concern one-dimensional systems, for whichit has been recently argued[25] that additional drawbacks appear in the TABC method. However it is important toremark that one-dimensional modelizations can often be used in the stellar matter case where (except at very highdensity close to saturation) spherical symmetry is often a good approximation. II. THE BLOCH THEOREM AND TWIST AVERAGED BOUNDARY CONDITIONS
We want to address the physical problem of an infinite system (specifically: the neutron star crust) constituted ofan infinite number of spatial replicas of finite systems of linear size L (the Wigner-Seitz cells). In the absence of anyperiodic external potential, for the replicated system to be equivalent to the infinite system, each Wigner-Seitz cellis supposed to contain an integer number of structures (spherical nuclei or more exotic ’pasta’ structures[20–23]) ofsize ~l , L i = n i l i , i = x, y, z . The minimal choice is to take a box containing exactly one structure, and in this casethe symmetry of the box will follow the symmety of the physical system (a cylindrical box for cylindrical structures,etc).To simplify the notations we will work in one dimension only and write L = nl . The extension to three dimensionsis straightforward. A. The Bloch theorem at the N-body level
We are interested in the translational invariance properties of the Hamiltonian induced by the imposed periodicity.The Hamiltonian of the global system is obviously invariant respect to a simultaneous translation of all particlecoordinates of the same arbitrary length r . This invariance physically corresponds to the general statement that thecenter of mass momentum is a good quantum number. This symmetry has no influence in an infinite system, and wewill not consider it further. The fact of working in a finite box of linear size L which is replicated induces additionallyan extra non trivial invariance, which we now discuss. The total Hamiltonian of the infinitely replicated system readsˆ H = ∞ X m =1 H L (cid:16) ˆ ~x ( m ) (cid:17) (1)where ~x ( m ) ≡ ( x m , . . . , x mN ) denotes coordinates of particles belonging to the m -th replica and the cell Hamiltonian isˆ H L = N X i =1 ˆ t i + ∞ X m = −∞ N X i>j v (ˆ x i − ˆ x j − mL ) , (2)Note that the cell Hamiltonian depends on a finite number N of particles even if these particles are not necessarilyconfined into the specific volume of the cell. This Hamiltonian is invariant under the translation of any particlecoordinate x k of a length mL , with m integer. This invariance can be expressed as h ˆ H L , ˆ T k ( m ) i = 0 (3)where ˆ T k ( m ) = exp (cid:2) − i ~ mL · ˆ p k (cid:3) is the L -translational operator of particle k . Translational invariance implies thatthe eigenfunctions Ψ ( x , . . . , x N ) of the cell Hamiltonian eq.(2) can be written asˆ T k ( m )Ψ θ k,m = exp ( − iθ k,m ) Ψ θ k,m (4)where θ k,m is the eigenvalue associated to the translation x k → x k − mL . Because of the property of translationaloperators ˆ T k ( m ) ˆ T k ( m ) = ˆ T k ( m + m ), the eigenvalues satisfy θ k,m + θ k,m = θ k,m + m . This means that we canassociate the translational invariance of periodicity L for particle k with an eigenvalue θ k such that θ k,m = mθ k (5)Let us define an auxiliary wave function asΦ ( x , . . . , x N ) = exp (cid:18) − i θ k L x k (cid:19) Ψ ( x , . . . , x N ) (6)Using the fact that Ψ is an eigenfunction of the translation operator eq.(4), we getΦ ( x , . . . , x N ) = exp (cid:18) − i θ k L ( x k − mL ) (cid:19) exp ( − imθ k ) Ψ ( x , . . . , x N )= Φ ( x , . . . , x k − mL, . . . , x N ) (7)We have shown that, if Ψ is an eigenfunction, then the function Φ defined by eq.(6) is a periodic function of period L . The same reasoning can be done for any particle k = 1 , . . . , N . This means that the eigenfunctions of thetranslationally invariant Hamiltonian eq.(2) can be written asΨ ( x , . . . , x N ) = exp i L N X k =1 θ k x k ! Φ ( x , . . . , x N ) (8)where Φ is a periodic function, that is invariant under the translation mL of any particle coordinate.Since the Hamiltonian eq.(2) is invariant under the translation of any particle coordinate separately, one couldconsider in principle a different phase θ k for each particle. However the indistinguishability of nucleons imposes thatall the phases must be equal, as we now show[4]. Let us consider a phase θ for the L -translation of particle 1, and aphase θ for the L -translation of particle 2:Ψ ( x + L, x , . . . , x N ) = exp ( iθ ) Ψ ( x , x , . . . , x N ) (9)Ψ ( x , x + L, . . . , x N ) = exp ( iθ ) Ψ ( x , x , . . . , x N ) (10)Applying the permutation symmetry to eq.(9) givesΨ ( x , x , . . . , x N ) = − exp ( − iθ ) Ψ ( x , x + L, . . . , x N ) (11)Applying a translation − L to the second coordinate givesΨ ( x , x , . . . , x N ) = − exp ( − i ( θ − θ )) Ψ ( x , x , . . . , x N ) (12)Applying the permutation symmetry once again we getΨ ( x , x , . . . , x N ) = exp ( − i ( θ − θ )) Ψ ( x , x , . . . , x N ) (13)which shows that the two phases must be equal, θ = θ = θ L . To be precise we could have θ − θ = 2 nπ with anyinteger n , but the phase θ L can be taken without any loss of generality in the interval (0 , π ] or ( − π, π ], that is in thefirst Brillouin zone. Indeed if we consider a very large number N r of replicas of the WS cell, N r → ∞ , the effect ofthe Bloch phase is negligible and we can write global periodic boundary conditions for the wave functionΨ ( x + N r L, . . . , x N + N r L ) = Ψ ( x , . . . , x N ) (14)This implies exp( iθ L N r ) = 1 or θ L = (2 n − N r ) π/N r with n integer. Let us write n = n ′ N r + n ” with 1 ≤ n ” ≤ N r and n ′ integer, then θ L = (2 n ′ − π + θ (15)with − π < θ ≤ π . Finally the wave function isΨ ( x , . . . , x N ) = exp i θL N X i =1 x i ! Φ ( x , . . . , x N ) (16)This Bloch formulation is a possible way to represent the eigenfunction of the exact N-body Hamiltonian in the cell,accounting for the infinite range antisymmetrization correlations with the replicated system. It amounts to introducingan extra quantum number θ , Ψ ( x , . . . , x N ) = h x , . . . , x N | ~n, θ i (17)where ~n denotes the ensemble of the other good quantum numbers. The system observables will thus explicitlydepend on the phase θ .The form of the periodic Hamiltonian eq.(1) implies that the wave function of the global system can be expressedas an antisymmetrized product of N -body Bloch wave functions (16) according toΨ tot ( x , . . . , x ∞ ) = ˆ A ∞ Y m =1 Ψ θ m ( x m , . . . , x mN ) (18)The Twist Averaged Boundary Condition method consists in assuming that the different angles θ m in the productare all different. Then the physical value of any observable ˆ O can be computed from eq.(18) with the extra restriction θ m = θ n giving h ˆ O i totn = ∞ X m =1 h ~n, θ m | ˆ O | ~n, θ m i (19)We can see that within this decoherence hypothesis among the different phases, possible non-diagonal terms disappearand the observables per unit cell are obtained as simple averages over phases h ˆ O i Ln = 12 π Z π − π dθ h ~n, θ | ˆ O | ~n, θ i (20)We will explicitly show in the following that the parts of the Fock space corresponding to different phases areindeed disjoint in the case of Slater determinants. For more complex eigenstates, eq.(20) has to be considered as anapproximation which however appears very well verified in practical applications[3]. B. Application to Slater determinants
In the following we look for an ansatz for Φ ( x , . . . , x N ) to be introduced as a variational approximation to the fullN-body problem. In particular mean field approaches, including Antisymmetrized Molecular Dynamics (AMD) andFermionic Molecular Dynamics (FMD), are based on a Slater approximation. This means that the variational ansatzcan be written as Φ ( x , . . . , x N ) = ˆ A N Y k =1 u k ( x k ) (21)where ˆ A is the antisymmetrization operator among the N particles and u k ( x + mL ) = u k ( x ) ∀ k . To fulfill thisperiodicity condition we can write u k ( x ) = N ∞ X m = −∞ g k ( x − mL ) (22)where g is an arbitrary function and N is a normalization factor. In the following we will show applications with thenon-orthogonal single-particle FMD/AMD basis set[26, 27] given by non-normalized gaussian wave packets g Z k ( x ) ≡ exp (cid:18) − a k ( Z k − x ) (cid:19) (23)where Z k , a k are complex variational parameters. For this specific choice the translation of the argument is equivalentto a translation of the gaussian centroid Z k u k ( x ) = lim N r →∞ √ N r N r / X m = − N r / g Z k + mL ( x ) (24)where the dependence on a k is implicit and omitted to simplify the notations. It is important to remark that thechoice of the normalization in eq.(24) guarantees that the norm of the wave function is finite, which will allow thenumerical evaluations below. In the special case where the gaussian width is sufficiently small respect to the cell size,such that the overlaps between gaussians can be neglected this norm is readily evaluated as < u k | u k > = p π | a k | .Let us consider a generic one-body operator ˆ O = P Nk =1 ˆ o k . The matrix element in r-space representation reads o jk = < u j | ˆ o | u k > = Z ∞−∞ dx Z ∞−∞ dx ′ u ∗ j ( x ′ ) o ( x ′ , x ) u k ( x ) = lim N r →∞ Z Nr L − Nr L dx Z ∞−∞ dx ′ u ∗ j ( x ′ ) o ( x ′ , x ) u k ( x ) (25)Because of the periodicity of the u k this simplifies to o jk = lim N r →∞ N r Z L/ − L/ dx Z ∞−∞ dx ′ u ∗ j ( x ′ ) o ( x ′ , x ) u k ( x ) (26)where we have used the fact that the operator is translationally invariant. For a local one-body operator ˆ A = P Nk =1 ˆ a k the matrix element simplifies to a jk = lim N r →∞ N r Z L/ − L/ dx u ∗ j ( x ) a ( x ) u k ( x ) (27)This shows that, because of the periodicity, only integrals over one single box are needed. The expectation value ofthe observable ˆ O is given by (for the moment we are ignoring the Bloch phase, and considering expectations only overthe auxiliary function Φ, which we indicate <> Φ ): < ˆ O > Φ = < Φ | ˆ O | Φ > = N X k =1 o kk (28)if the u k constitute an orthonormal basis, and < ˆ O > Φ = < Φ | ˆ O | Φ >< Φ | Φ > = N X j,k =1 o jk B − kj (29)if not. Here B jk is the overlap matrix B jk = < u j | u k > . (30)Implementing the ansatz (24), the matrix element of a generic local operator is calculated by a jk = lim N r ,N ′ r ,N ′′ r →∞ N ′′ r p N r N ′ r N r / X m = − N r / N ′ r / X m ′ = − N ′ r / Z L/ − L/ dxg ∗ Z j + mL ( x ) a ( x ) g Z k + m ′ L ( x ) (31)Because of the finite norm of u k , the double sum is a convergent quantity with increasing number of replicaslim N r ,N ′ r →∞ N r / X m = − N r / N ′ r / X m ′ = − N ′ r / Z L/ − L/ dxg ∗ Z j + mL ( x ) a ( x ) g Z k + m ′ L ( x ) = (32)lim N r ,N ′ r →∞ p N r N ′ r N r / X m = − N r / N ′ r / X m ′ = − N ′ r / Z ∞−∞ dxg ∗ Z j + mL ( x ) a ( x ) g Z k + m ′ L ( x ) (33)= < u k | ˆ a | u k > (34)This means that the limits in eq.(31) can be performed separately giving: a jk = Z L/ − L/ dx ∞ X m,m ′ = −∞ g ∗ Z j + mL ( x ) a ( x ) g Z k + m ′ L ( x ) (35)In the hypothesis that the spatial extension of the wave function is negligible beyond a linear size L max = M L ,lim x →± L max h x | g Z i = 0, we can write a jk = Z L/ − L/ dx M X m,m ′ = − M g ∗ Z j + mL ( x ) a ( x ) g Z k + m ′ L ( x ) (36)which is a feasible integral. Eq.(36) is readily generalized to the case of a non-local operator.Because of the translational invariance the wave function has to be multiplied by the Bloch phase factor eq.(16).This means that the expectation values of observables have to be corrected respect to eq.(28),(29) < ˆ O > Ψ = N X j,k =1 o jkθ B − kjθ = < Ψ θ | ˆ O | Ψ θ >< Ψ θ | Ψ θ > (37)where the notation <> Ψ corresponds to a specific choice for the (arbitrary) Bloch phase − π < θ ≤ π ,Ψ θ ( x , . . . , x N ) = ˆ A N Y k =1 ψ kθ ( x k ) , (38) ψ kθ ( x ) = exp ( iθx/L ) u k ( x ) (39)We can finally write the expression of the different matrix elements for the infinite periodic system, for a given choiceof the Bloch phase θ . The overlap matrix is not influenced by the Bloch phase, and the same is true for the matrixelement of any local operator: a jkθ = M X m,m ′ = − M Z L/ − L/ dxg ∗ Z k + mL ( x ) a ( x ) g Z k + m ′ L ( x ) = a jk (40)In particular the local density ρ θ ( x ) = A X j,k =1 < x | ψ k >< ψ j | x > B − kj = A X j,k =1 < x | u k >< u j | x > B − kj = ρ ( x )is independent of θ .The situation is different for non-local operators, because in this case the Bloch phase does not cancel any more o jkθ = ∞ X m = −∞ M X m ′ = − M Z L/ − L/ dx Z ∞−∞ dx ′ g ∗ Z j + mL ( x ′ ) o ( x, x ′ ) g Z k + m ′ L ( x ) exp ( iθ ( x − x ′ ) /L ) (41)Let us take the exemple of the momentum operatorˆ kψ j ( x ) = − i M X m = − M exp ( iθx/L ) (cid:18) ∂∂x + iθL (cid:19) g Z j + mL ( x ) (42)The expectation value of the total momentum in the cell ˆ K = P Nj =1 ˆ k j is given by < ˆ K > Ψ = N X ji =1 < ψ j | ˆ k | ψ i > B − ij = N X ji =1 < u j | ˆ k | u i > B − ij + N θL = < ˆ K > Φ + N θL (43)We can see that the Bloch phase physically represents a boost to the center of mass motion. It will cancel whenperforming the average between the different phases − π < θ ≤ π (TABC), but this will give a finite contribution tothe kinetic energy. Indeed if we consider the square momentumˆ k ψ j ( x ) = − exp ( iθx/L ) M X m = − M (cid:18) ∂ ∂x + 2 iθL ∂∂x − θ L (cid:19) g Z j + mL ( x ) (44)This gives rise to a total kinetic energy expectation value < ˆ E K > Ψ = ~ m N X jk =1 < ψ j | ˆ k | ψ k > B − kj ++ ~ m θL < ˆ K > Φ + ~ m N θ L = < ˆ E K > Φ + ~ θmL < ˆ K > Φ + N m (cid:18) ~ θL (cid:19) (45)This extra kinetic energy term linked to the Bloch phase explicitly enters in the energy variation (or in the equationsof motion, in the case of dynamical models[22, 23, 26, 27]). If we consider the typical size of a Wigner-Seitz cellas calculated by Negele and Vautherin[13], in the inner crust ( L ≈
15 fm) the phase effect gives rise to an energycontribution ∆ E ≈
11 MeV, which is far from being negligible if we compare to the corresponding Fermi energy E F ≈
24 MeV [16–18].
C. One-dimensional model cases
As a first model case, we consider one-dimensional free particles. Even if this system is clearly very far fromthe correlated dishomogeneous solutions relevant for stellar matter, it has the advantage of being exactly solvable.Moreover, it is an especially interesting test case for molecular dynamics models [24, 26, 27] and more general modelsemploying localized wave functions[23]: these models are optimized to treat density fluctuations but not necessarilyadapted to treat the coupling to the continuum which is needed for the extreme neutron-proton ratio which is foundin the inner crust of neutron stars.In one dimension, the energy per particle of a degenerate ideal Fermi gas at density ρ = N/L = k F π is given by e F G = ~ π ρ m . In the case of a finite system of N particles within a length L with periodic boundary conditions, theenergy levels are given by E i = ~ m (cid:18) πn i L (cid:19) , (46)where n i is an integer within the interval [ − N , N ]. Due to this level degeneracy, the total cell momentum h ˆ K i L isequal to ± πNL for an even number of particles, while it is zero for an odd number.Because of the periodicity of the system, it is easier to work in Fourier space. All one-body observables can beexpressed as a function of the one-body density in momentum space ρ Ψ ( k, k ′ ) = P Ni,j =1 h k | ψ i ih ψ j | k ′ i B − ij , where thesingle-particle states are given in the momentum representation. These states can be obtained by taking the Fouriertransform of the wave function. Considering that a periodic function can always be expressed as a Fourier series, u n ( x ) = ∞ X l = −∞ c nl e ( − πilx/L ) (47)we get ψ nθ ( k ) = √ π ∞ X l = −∞ c nl δ ( θ/L − k − πl/L )= lim N r →∞ ∞ X l = −∞ Nr − X m = − Nr − r πN r L Z L − L dx g Z n + mL ( x ) e iπlx/L δ ( θ/L − k − πl/L )= lim N r →∞ ∞ X l = −∞ r πN r L Z NrL − NrL dx g Z n ( x ) e iπlx/L δ ( θ/L − k − πl/L )= 2 πL r a n N r ∞ X l = −∞ e − an ( πlL ) e i πlL Z n δ ( θ/L − k − πl/L ) (48)At variance with the spatial density, the momentum density thus explicitly depends on the value of the phase θ = Lk : ρ θ ( k ) = lim N r →∞ π N r L N X n,p =1 B − n,p p a n a ∗ p e − an + a ∗ p ( πlL ) e i πlL ( Z n − Z ∗ p ) δ ( θ/L − k − πl/L ) (49)The momentum density is obtained by averaging over the different phase values, giving ρ Ψ ( k ) = lim N r →∞ πN r L Nr X n = − Nr − N X n,p =1 B − n,p p a n a ∗ p e − an + a ∗ p ( πlL ) e i πlL ( Z n − Z ∗ p ) (50)We can see that the value of ρ in k is solely determined by the wave function Ψ corresponding to the phase value θ = k − πm , which corresponds to a unique value in the first Brillouin zone. This means that the number of pointschosen for the discretization of the integral over the Bloch phase defines the resolution of the momentum density: theBloch phase generates new plane waves which ’fill up’ the distribution as needed to obtain a continuous Fermi seastarting from the discrete levels of a finite system.The total energy and momentum easily follow D ˆ E E Ψ = ~ m Z ∞−∞ dkρ Ψ ( k ) k , D ˆ K E Ψ = Z ∞−∞ dkρ Ψ ( k ) k. (51)The result for the FMD model is shown in Fig.1. We have chosen a particle density ρ = . π . The size L of theelementary cell is directly proportional to the number of particles in the cell according to L = Nρ . The FMD groundstate is obtained minimizing the total energy with respect to the variational parameters ~z = { Z n , a n , n = 1 , . . . , N } .This can be obtained using the gradient method, or alternatively by the numerical solution of the FMD equations ofmotion [26] d h ˆ H i Ψ dz i = X j C ∗ ij dz ∗ i dt , (52) C ij = i ~ ∂ ∂z ∗ i ∂z j ln h Ψ | Ψ i (53)It is possible to show[26] that if a friction term is added as a multiplicative factor on the r.h.s. of eq.(52) d h ˆ H i Ψ dz i = A ∗ X j C ∗ ij dz ∗ i dt , (54)it gives a dynamics of the variational parameters leading to a decrease in time of the total energy. The resultingdissipative dynamics can be written as: d h ˆ H i Ψ dt = 2 i ℑ ( A ) X ij C ij dz ∗ i dt . dz j dt (55)and this equation is numerically followed until convergence. ρ [ f m − ] x [ fm ] 00.20.40.60.811.2 -1 0 1 < k | ˆ ρ | k > k F k Nk F -101 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 -6-4-20246 h ˆ K i Ψ k F h ˆ K i Φ k F k Nk F h ˆ H i Ψ N h ˆ H i Φ N k Nk F FIG. 1: FMD calculation of a system of three free particles in a one dimensional box of size L with Bloch boundary conditions.Upper part: behavior of the total linear momentum (left) and total energy (right) as a function of the Bloch wave number k = θ/L Full lines: total expectation values h ˆ K i Ψ , h ˆ H i Ψ (see text); dashed lines: averages over the variational part of thewave function h ˆ K i Φ , h ˆ H i Φ without the contribution of the Bloch phase (see text); dash-dotted line: total energy average overthe Bloch phases. Lower part: spatial (left) and momentum (right) density. As expected, Figure 1 shows that the periodicity in k -space is reproduced by the calculation. This is due to theimplicit dependence of the expection values h ˆ P i Φ , h ˆ H i Φ on the choice of the Bloch quantum number k = θ/L obtainedthrough the variational procedure.The exact results of a free Fermi gas are recovered in the model. This means that the FMD ansatz for the wavefunction appears adapted to describe the free particle system. This was not a-priori obvious because the choice ofsingle-particle wave packets represents a restriction of the complete Slater determinant variational space, and the0plane wave is obtained only as a mathematical limit a → ∞ of the wave function, which is not exactly accessible innumerical simulations. This implies that the FMD model, which has been introduced essentially to describe clusterdegrees of freedom, can also address continuum states.The other interesting point is that the exact Fermi gas result is obtained with any arbitrary number of particles inthe simulation. For the 3-particles system shown in Figure 1, the triangles in the lower right part of the figure showthe momentum values allowed for this finite system. These values are eigenvalues of the periodic part of the wavefunction Φ, and their spacing decreases with increasing number of particles. As shown by eq.(50), the effect of the θ quantum number is to produce all the other missing values for the k quantum number which would be accessible tothe infinite system, thus reproducing the continuous Fermi distribution.In conclusion, the application of TABC to Bloch single particle wave functions appears as a very powerful method torecover the correct kinetic energies and wave functions at the thermodynamic limit with variational methods appliedto small systems even in one dimension. A word of caution is however necessary: to demonstrate the Bloch theoremwe have explicitly used the fact that the states are eigenvectors of the many-body Hamiltonian, which is in generalnot true when using variational methods. For this reason we will develop an alternative scheme in the next section. III. THE REPLICA METHOD
The twist averaged boundary condition are currently applied in QMC calculations of correlated electron systemsat the thermodynamic limit. In the specific framework of fermionic molecular dynamics, an alternative methodto deal with the infinite range antisymmetrisation correlations has been proposed, based on the localized Wannierrepresentation[24]. Taking advantage of the nested structure of the overlap matrix in a periodic system, an analyticalmethod of inversion of this matrix is proposed in ref.[24]. We present here a slightly different derivation of the sameequations, for an application to any arbitrary single particle basis.
A. General formulation
Let us consider that our WS cell contains already many different replicas of the physical system, N tot = N r N wherefinally we will let N r → ∞ . Similarly, we now interpret the simulation cell L tot = nL as a huge length, n ≫
1, suchthat working with a finite L tot can be considered as equivalent to the thermodynamic limit, meaning that the scalinghas to be fulfilled for all observables ˆ O < ˆ O > = N r < ˆ O > L (56)with < ˆ O > L = 1 N r N tot X j,k =1 o jk B − kj (57)independent of the periodicity properties of the Hamiltonian. If we consider the N tot -body wave function, the invari-ance property of the system under the translation mL tot of any particle coordinate is still verified. The differencerespect to the previous section is that the system composed of N tot particles has an additional symmetry which isnot directly linked with the use of a finite box L but is rather due to the physical periodicity of the pasta structuresover a lenght scale L = L tot /n . We have already observed that, because of the absence of an external field, the cellHamiltonian eq.(2), which now represents the full Hamiltonian eq.(1),ˆ H = N tot X i =1 ˆ t i + ∞ X m = −∞ N tot X i>j v (ˆ x i − ˆ x j − mL tot ) (58)is invariant under the simultaneous translation of every particle coordinate of any arbitrary length. The translationalinvariance which will be relevant to us is the invariance under the simultaneous translation of every particle coordinateof the physical length L , defined as the periodic length of the one body density ρ ( x + mL ) = ρ ( x ) (59)If the Hartree-Fock field as an external field, this translational invariance would be the only one verified by the mean-field hamiltonian, and it would be verified for each particle coordinate separately. The situation is different in the1case of an interacting system as described by the self-consistent Hartree-Fock approach. In this case, translationalinvariance with respect to any arbirary length is respected because of the self-consistency of the mean-field approach, U ij = δE HF /δρ ji where E HF is the interaction part of the Hartree-Fock energy. However, this translational invarianceapplies only to the simultaneous translation of all particle coordinates. The only special feature of the length L isthat eq.(59) is fulfilled, which will allow us to impose a simplified expression for the one-body wave functions, as wewill see in a moment.Following the derivation of the previous section, these two invariance properties lead to the definition of a Blochphase for the N tot -body wave function according toΨ ( x , . . . , x N tot ) = exp i θ l LN tot N tot X k =1 x k ! Φ ( x , . . . , x N tot ) (60)where Φ is a periodic function, that is invariant both under the translation mL of all N tot particle coordinatessimultaneously, and under the translation mL tot of any particle coordinate.The difference respect to the derivation of the previous section is that now we have a factor 1 /N tot multiplying theBloch phase. Since we are at the thermodynamic limit N tot → ∞ , the correction on kinetic observables induced bythis Bloch phase can be neglected and we can consider that the total N tot -body wave function Ψ shares the periodicityproperties of Φ.Let us now look for an ansatz for this L − and L tot − periodic Ψ ( x , . . . , x N tot ) to be introduced as a variationalapproximation to the full N tot -body problem. We take the same Slater ansatz used in the previous section:Ψ ( x , . . . , x N tot ) = ˆ A N tot Y k =1 u k ( x ) (61)with u k ( x + mL tot ) = u k ( x ). Because of the periodicity of the one body density, each sub-cell of linear dimension L can be exactly associated to a finite number N = N tot /N r particles, and the wave function can be written introducinga sub-cell index m as: Ψ ( x , . . . , x N tot ) = ˆ A N Y k =1 N r Y m =1 u k,m ( x ) (62)Moreover, the periodicity of ρ can be written as ρ ( x + L ) = N r X m,m ′ =1 N X k,j =1 u ∗ k,m ( x + L ) u j,m ′ ( x + L ) B − jk = N r X m,m ′ =1 N X k,j =1 u ∗ k,m ( x ) u j,m ′ ( x ) B − jk = ρ ( x ) (63)This condition can be satisfied introducing only N independent one body wave functions u k,m ( x ) = u k ( x − mL ) wherein principle u k are arbitrary (non periodic) functions. If these functions would form an orthonormal basis < j, m | k, m ′ > = δ jk δ mm ′ (64)the problem would be reduced to a N body problem as with the ansatz eq.(21). Indeed in the calculation (57) of anarbitrary one-body observable only the u k , k = 1 , . . . , N functions would be needed < ˆ O > L = 1 N r N tot X k =1 o kk = N X k =1 o kk . (65) B. Application to non-orthogonal single-particle basis
For practical applications we will not work with an orthonormal basis, but we want to take the FMD/AMD gaussianwave packet variational ansatz eq.(23), namely u k ( x − ml ) = g Z k + ml ( x ) (66)2Since the g Z k + ml do not form an orthonormal basis, we are left with a problem of computational size N tot , that isimpossible to solve (recall N tot → ∞ ).A possibility would be again to impose the simplified expression:Ψ ( x , . . . , x N tot ) = N r Y m =1 ˆ A N Y k =1 g Z k + ml ( x ) = N r Y m =1 Ψ ( x , . . . , x N ) (67)which allows to work with only N functions.The conceptual problem with this simple ansatz eq.(67) is that it neglects the antisymmetrization correlationsamong different cells. We know that such correlations are important as they are at the origin of the band structure.For this reason we need to keep the antisymmetrization correlation among all the N tot = N r · N particles.By introducing a linear transformation ~u k = ˆ U ~g k (68)defined by u k,n ( x ) = 1 √ N r N r / X m = − N r / exp (cid:18) im πnN r (cid:19) g Z k + mL ( x ) , (69)we can reduce the computation size to N r N keeping the antisymmetrization correlations among all the N tot = N N r particles. It is interesting to remark that, in the limit of N r → ∞ , this transformation coincides with the definitionof Wannier functions[1]. This transformation preserves the determinantˆ A N r Y m =1 g Z k + mL ( x ) = ˆ A N r Y n =1 u k,n ( x ) (70)so that we can write Ψ ( x , . . . , x N tot ) = ˆ A N r Y n =1 N Y k =1 u k,n ( x ) (71)The advantage of the transformation (69) is that matrices are block-diagonal, that is functions corresponding todifferent phases γ = (2 n − N r ) π/N r are orthogonal < u j,n | u k,n ′ > = δ n,n ′ B jk ( n ) . (72)We have thus shown that the decoherence hypothesis eq.(19) of the TABC method is exactly verified in the case ofvariational approaches based on Slater determinants.Having recovered the orthogonality on the level of the replica quantum number we can write < ˆ O > L Ψ = 1 N r N r X n =1 N X j,k =1 < u j,n | ˆ o | u k,n > B − kj ( n ) (73)The matrix element of any operator is given by o jk ( n ) = < u j,n | ˆ o | u k,n > = Z L/ − L/ dx ∞ X m,m ′ = −∞ g ∗ Z j + mL ( x ′ ) o ( x, x ′ ) g Z k + m ′ L ( x ) exp (cid:18) − i ( m − m ′ ) 2 πnLN r (cid:19) (74)and explicitly depends on the n quantum number associated to the antisymmetrization among the different subcells.In particular the overlap matrix reads B jk ( n ) = < u j,n | u k,n > = Z L/ − L/ dx ∞ X m,m ′ = −∞ g ∗ Z j + mL ( x ) g Z k + m ′ L ( x ) exp (cid:18) − i ( m − m ′ ) 2 πnLN r (cid:19) (75)Using the same argument that we had for the expectation (36), the series can be reduced to a finite sum B jk ( n ) = M X m,m ′ = − M Z L/ − L/ dxg ∗ Z j + mL ( x ) g Z k + m ′ L ( x ) exp (cid:18) − i ( m − m ′ ) 2 πnLN r (cid:19) (76)3In practical application, a large value of M ≈
10 is needed to describe non-interacting particles, but M = 1 turnsout to be sufficient in the presence of the nuclear and Coulomb interaction for the particle densities relevant for theneutron star crust. From the computational viewpoint, M = 1 is equivalent to calculate a single gaussian overlapbetween each pair if periodic boundary conditions are applied to the gaussian, meaning that the simulation of theinfinite system is not more expensive numerically than the one for the finite system. The analogous expression forthe matrix elements reads < u j,n | ˆ a | u k,n > = M X m,m ′ = − M Z L/ − L/ dxg ∗ Z j + mL ( x ) a ( x ) g Z k + m ′ L ( x ) exp (cid:18) − i ( m − m ′ ) 2 πnLN r (cid:19) (77)Going to the N r → ∞ limit we can write lim N r →∞ πnLN r = γ (78)with − π < γ ≤ π . The observables read < ˆ O > L Ψ = 12 π Z π − π dγ N X j,k =1 < u j,γ | ˆ o | u k,γ > B − kj ( γ ) (79)which is expected to give equivalent results to the implementation where TABC are applied. E N [ M e V ] N r FIG. 2: Energy associated to a system of two free FMD particles in one dimension at a density ρ = . π as a function of thenumber of replicas. Dash-dotted line: free Fermi gas in one dimension. The FMD ground state energy for a system of two free particles in one dimension is shown in Figure 2 as a functionof the number of replicas, with the Wannier choice for the wave packets. For each calculation, the number M waschosen high enough to have convergent results. We have already observed that the FMD variational space is sufficientlyflexible to describe a free-particles system. The FMD solution for each value of N r thus represents the exact energyof a finite system with periodicity L = 2 /ρ and size L tot = L · N r . The thermodynamic limit, corresponding to theanlytical Fermi gas result, is obtained for about ten replicas, corresponding to ten values for the phase γ of eq.(78). IV. COMPARISON OF THE TWO METHODS
In the previous sections, we have introduced the Bloch and the Wannier representation as two alternative ways toaddress the thermodynamical limit in variational theories. From a principle point of view the two approaches are very4similar, however they lead to very different expressions in the calculation of observables. In particular, the fact ofconsidering only diagonal terms in the phase quantum number naturally emerges as a consequence of the orthogonalityof the Wannier basis in the replica method, while it appears as a simple working hypothesis in TABC. Moreover,the equivalence is expected in the ideal case of an infinite number of twists and under the condition that the wavefunction is an exact eigenvalue of the Hamiltonian. This first condition is numerically expensive, while the second isnot assured in variational methods. For all these reasons, the equivalence of the two formalisms is not a-priori clearand will be discussed in the present section.
A. Finite number of replicas and finite number of twists
We now show that, in the hypothesis that the variational theory produces the exact eigenvectors of the many-bodyHamiltonian, the two methods to treat the antisymmetrisation correlations can be mapped on each other for anysystem size.In the Bloch formalism, if we impose that the condition of periodic boundaries eq.(14) is verified for a finitetranslation of size N r L , the original problem of an infinite system with periodicity L is transformed into the problemof a finite system of size N r L , with periodic boundary conditions.The allowed values for θ are then discrete θ = n − N r N r π , where n is an integer. This amounts to discretize theintegral (20) with ∆ θ = πN r . We expect then the two methods to give exactly the same results. This expectation isconfirmed by Figure 3, which compares for the case of the non-interacting one-dimensional system the energy obtainedconsidering a periodic system of N tot = ρN r L particles, with the calculation of a reduced N = ρL system averagedover N r phases with the two (Bloch and replica) proposed methods. We can see that the agreement is perfect forall sizes, and all calculations converge to the free Fermi gas result with increasing number of replicas (respectively:twists). E N [ M e V ] N r L tot [ fm ]FIG. 3: Energy of a free periodic one-dimensional system at a density ρ = . π as a function of the number of replicas. Theupper abscissa represents the total linear size in the case of simple periodic boundary conditions (full line), the lower one is anumber of replicas for the replica method (crosses), and a number of twists for the Bloch method (circles). Dash-dotted line:free Fermi gas in one dimension. This shows the equivalence of the different methods when the wave vector is an eigenvalue of the Hamiltonian.5
B. Limitations of the Bloch method in variational applications
In realistic applications of interacting fermion system, the variational solution is always an approximation of theexact eigenvalue. This may create a bias in the different method that we turn to discuss. For eq.(4) to correctlyaccount for the correlations due to the particles outside the cell L , it is necessary that the wave function Ψ θ is aneigenfunction of the hamiltonian. Conversely the Wannier representation can be obtained from the solution of theinfinite periodic system via a simple linear transformation eq.(68). Of course if the variational ansatz is not adequateto describe the exact state both representations will be false. However the error due to the inadequacy of the solutionin the cell L will be propagated in an uncontrolled way in the Bloch method for the description of the global system,and this will not be the case for the replica method that directly addresses the replicated system of size N r L .To illustrate this difference, we can once again take the simple example of the non-interacting Fermi gas. There isa popular microscopic theory in nuclear physics, the so-called AMD model[27], where the single particle wave packetsare gaussians of fixed width. With this ansatz the variational space does not contain the exact solution of this problemat the bulk limit, that is an antisymmetrized combination of plane waves. In particular, if we choose a small value forthe width, as it is variationally obtained to describe finite nuclei[27], it is clear that the states will be very far fromthe exact eigenvectors of the kinetic hamiltonian. The comparison of the two methods in this model case is shown inFigure 4. E / N [ M e V ] L [ F m ]13.213.614 5 10 15 20FIG. 4: Comparison of the different boundary conditions for the AMD model in the case of a non interacting system as afunction of the system size. Dash-dotted line: free Fermi gas. Dashed line: AMD ground state energy for a system size largeenough for finite size effects be negligible ( L = 150).Circles: finite system with Bloch boundary conditions. Cross: replicamethod. Stars: periodic boundary conditions. We can see in Figure4 that the AMD ansatz is not adapted for this problem, from the fact that asymptoticallythe energy of the free Fermi gas (full line) is not recovered by the AMD calculation (dashed line). This asymptoticvalue is perfectly reproduced by the replica method (crosses) for any system size, while the Bloch method (opencircles) produces artificial fluctuations and attains convergence only when the Bloch phase is negligibly small, and thecalculation is perfectly equivalent to simple periodic boundary conditions (stars).The non-interacting system is certainly not the ideal application ground of molecular dynamics model. As a secondmodel case which is closer to the physical case of the neutron star crust we take a one-body periodic potential, whichis obtained as a periodic self-consistent mean field in realistic applications of mean-field variational models [7–10, 12] ifthe two-body nuclear and Coulomb interaction is added, including the interaction with a uniform electron background.For this model one-dimensional application we only wish to discuss the effect of the boundary conditions, thereforewe restrict to the simpler case of an external periodic potential. In particular, the choice of a gaussian potential has6the advantage of making all calculations analytical. Let us take a potential form: V gaus (ˆ x ) = V N X n = − N ∞ X m = −∞ exp (cid:20) − a (ˆ x − b n − mL ) (cid:21) (80)where V , a, b n are parameters. In particular, the choice a = ρ , b n = nρ for a system with average density ρ , correspondsto having a single particle per potential well. The average energy is readily calculated as D ˆ V gaus E Ψ = V lim N r →∞ N r N X i,j =1 B − ij s π a i a ∗ j aa ∗ j a + a i a + a i a ∗ j ∞ X m ,m ,n = −∞ exp " − a ∗ j ( Z i − b n + m L ) + a i ( Z ∗ j − b n + m L ) + a ( Z i − Z ∗ j − ( m − m ) L ) a ∗ j a + a i a + a i a ∗ j (81)Results from the Bloch method for a case of 3 particles with an average density ρ = . π are presented in Figure 5.At variance with the free particle case, the periodic part of the wave function Φ now explicitly depend on the θ quantum number, leading to a θ dependence of the one body density.The comparison between the different methods of treating the boundary conditions is shown in figure 6.This model is less trivial than the free Fermi gas, although still very schematic. The FMD model is expected tobe a good approximation of this system, but there is no guarantee the FMD solution should be exact. Similar to theprevious application, the replica method gives by construction the asymptotic result for any finite size, provided asufficient number of phases is considered. Conversely for the Bloch method, which is calculated with the same numberof phases as the replica method for the application of Fig.6, a convergence is reached only when the effect of the phasecan be neglected.This model example shows that in realistic applications, where any variational ansatz might be simply an approxi-mation of the exact energetics of the system, the different boundary conditions are not equivalent, and the employ ofBloch corrections to the single particle basis may increase the deviation respect to the exact solution. V. CONCLUSIONS
In this paper we have introduced and critically discussed two different ways of addressing the boundary conditionsin finite fermionic system in order to account for the infinite range antisymmetrisation correlations present at thethermodynamical limit. In the case of state vectors which are eigenstates of the many-body Hamiltonian, we haveshown that the use of Bloch wave functions with twist averaged boundary conditions physically corresponds exactly toa bigger system constituted of a number of replicas of the system under study equal to the number of phases. In turn,this last formalism can be interpreted as the use of localized Wannier states. While these different representations areequivalent when dealing with eigenvectors, the same is not true in the case of approximate solutions of the many-bodyproblem, as it is the case in the variational approach. In this case, the approach of the thermodynamic limit is notproperly treated by the Bloch method, while the replica method provides a very accurate evaluation of the Fockenergy. Even in the case of a non-orthogonal single particle basis this method can be numerically implemented witha moderate computational cost, opening the possibility of a completely quantum-mechanical treatment of nuclearmatter with molecular dynamics approaches.
Acknowledgments
This paper has been partly supported by ANR under the project NEXEN.Discussions with Klaas Vanternhout and Karim Hasnaoui are gratefully acknowledged. [1] N. W. Ashcroft and N. D. Mermin, ’Solid State Physics , Holt-Saunders (1976).[2] C. Lin, F.-H. Zong and D. M. Ceperley, Phys.Rev.
E64 , 026502 (2011).[4] G. Rajagopal, R.J.Needs, A.James, S.D.Kenny, W.M.C.Foulkes, Phys.Rev. B 51 -0.8-0.6-0.4-0.200.20.40.60.8 -4 -3 -2 -1 0 1 2 3 4 -60-40-200204060 ρ [ f m − ] V ( x ) [ M e V ] x [ F m ] -2 -1 0 1 2 -60-40-200204060 < k | ˆ ρ | k > k F k Nk F -1-0.500.51 -1 0 1 h ˆ K i k F k Nk F -1 0 1 -10.8-10.6-10.4-10.2-10-9.8-9.6 h ˆ H i A k Nk F FIG. 5: FMD calculation of a system of three particles in a one dimensional box of size L subject to a periodic externalpotential with Bloch boundary conditions. Upper part: behavior of the total linear momentum (left) and total energy (right)as a function of the Bloch wave number k = θ/L Full lines: total expectation values h ˆ K i Ψ , h ˆ H i Ψ ; dash lines: averagesover the variational part of the wave function h ˆ K i Φ , h ˆ H i Φ without the contribution of the Bloch phase. Dashed-dotted lines:expectation values including the phase average. Lower part: dash-dotted lines: spatial (left) and momentum (right) density;full line: gaussian potential.[5] J. M. Lattimer and M. Prakash, Science Vol.304 no. 5670, 536 (2004).[6] P.Haensel,A.Y.Potekhin,D.G.Yakovlev, ’Neutron stars: equation of state and structure’ , Springer, Berlin (2007).[7] N. Sandulescu, N.V. Giai, and R.J. Liotta, Phys, Rev.
C 69 , 045802 (2004).[8] E. Khan, M. Grasso, J. Margueron, N.V. Giai, Nucl. Phys.
A 800 , 37 (2008).[9] M. Grasso, E. Khan, J. Margueron, N. Van Giai, Nucl. Phys.
A 807 , 1 (2008).[10] C. Monrozeau, J. Margueron, and N. Sandulescu, Phys. Rev.
C 75 , 065807 (2007).[11] S. S. Avancini, S. Chiacchiera, D. P. Menezes, C. Providencia, Phys.Rev.
C 82 , 055807 (2010).[12] F.Douchin, P.Haensel, Phys.Lett.
B 485 , 107 (2000).[13] J.W. Negele, D. Vautherin, Nucl. Phys.
A 207 , 298 (1973).[14] C.J. Pethick and D.G. Ravenhall, Annu. Rev. Nucl. Part. Sci. (1995).[15] M. Baldo, E.E. Saperstein, and S.V. Tolokonnikov, Nucl. Phys. A 775 , 235 (2006).[16] B. Carter, N. Chamel, P. Haensel, Nucl. Phys.
A 748 , 675 (2005).[17] N. Chamel,S. Naimi, E. Khan, and J. Margueron, Phys. Rev.
C 75 , 055806 (2007).[18] N. Chamel, J. Margueron, and E. Khan, Phys. Rev.
C 79 , 012801(R) (2009).[19] K. Hasnaoui, PhD Thesis, Caen University, 2008,http://tel.archives- ouvertes.fr/tel-00337606/fr. -10.8-10.7-10.6-10.5-10.4-10.3-10.2-10.1 5 10 15 20 25 02468101214 E N [ M e V ] E F G N L [ F m ]FIG. 6: Comparison of the different boundary conditions for the FMD model in the case of a system subject to an externalgaussian potential as a function of the system size. Dash-dotted line: free Fermi gas. Dashed line: FMD ground state energy fora system size large enough for finite size effects be negligible ( L = 30). Circles: finite system with Bloch boundary conditions.Cross: replica method. Stars: periodic boundary conditions.[20] C.J. Horowitz, M.A. Perez-Garcia, D. K. Berry, J. Piekarewicz, Phys.Rev. C 72 , 035801 (2005).[21] G. Watanabe, H. Sonoda, T.Maruyama, K. Sato, K. Yasuoka and T. Ebisuzaki, Phys. Rev. Lett. , 121101 (2009).[22] W.G.Newton and J.R.Stone, Phys. Rev.
C 79 , 055801 (2009).[23] F. Sebille, S. Figerou, V. de la Mota, Nuclear Physics
A 822 , 5173 (2009).[24] K.Vantournhout, ’Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermionsystems’ , Ph.D. dissertation, Ghent University (2009); K.Vantournhout, T.Neff, H.Feldmeier, et al., Prog. Part. Nucl.Phys., to be published (2011).[25] R.M.Lee, N.D.Drummond, Phys. Rev.
B 83 , 245114 (2011).[26] Feldmeier H, Schnack J, Rev. Mod. Phys. , 655688 (2000).[27] Ono A, Horiuchi H, Prog. Part. Nucl. Phys.53