Finitesize scaling for quantum criticality using the finiteelement method
Abstract
Finite size scaling for the Schrödinger equation is a systematic approach
to calculate the quantum critical parameters for a given Hamiltonian. This
approach has been shown to give very accurate results for critical parameters
by using a systematic expansion with global basistype functions. Recently, the finite element method
was shown to be a powerful numerical method for ab initio electronic structure calculations
with a variable realspace resolution. In this work, we demonstrate how to obtain quantum critical
parameters by combining the finite element method (FEM) with finite size scaling (FSS) using
different ab initio approximations and exact formulations.
The critical parameters could be atomic nuclear charges, internuclear distances, electron density, disorder,
lattice structure, and external fields for stability of atomic, molecular systems and quantum
phase transitions of extended systems.
To illustrate the effectiveness of this approach we provide detailed
calculations of applying FEM to approximate solutions for the twoelectron
atom with varying nuclear charge; these include HartreeFock, density functional theory under
the local density approximation,
and an “exact” formulation using FEM. We then use the FSS approach to determine its critical
nuclear charge for stability; here, the size of the system is related to the number of
elements used in the calculations. Results prove to be in good agreement with previous Slaterbasis set
calculations and demonstrate that it is possible to combine finite size scaling with the finiteelement method
by using ab initio calculations to obtain quantum critical parameters. The combined approach provides a promising
firstprinciples approach to describe quantum phase transitions for materials and extended systems.
PACS numbers: 02.70.Dh 64.60.an 05.70.Jk 82.60.s
I Introduction
The theory of finite size scaling (FSS) in statistical mechanics can provide us with numerical methods fisher (); widom (); barber (); privman (); cardy (); nightingale1 (); Peter1 (); Peter2 () capable of obtaining accurate results for infinite systems by studying the corresponding small systems. In the past, basis set approximation to manybody systems have been combined with FSS to obtain critical parameters, where calculations for these systems were done by expanding the wave function in Slatertype basis sets or Gaussiantype basis functions. However, a generalization of larger atomic and molecular systems proved itself difficult using these type of functions MSK (); PSK ().
More recently, we have combined FSS with the finite element method (FEM) in order to obtain critical parameters in a potential that admits a finite number of bound states and we obtained excellent results by using both finiteelement and Slatertype basis functions AMWK (). For this potential we also established that the finitedifference (FD) method can be combined with FSS. The finite difference method, however, is not a variational method since it approximates the operators (i.e. the derivative) and therefore can result in errors of either sign TT2 (). As a result, we observed non monotonic behavior of the calculated critical parameters with respect to the number (N) of grid points used Daniel (). Nevertheless, for a large enough number of grid points N one can recover the correct asymptotic behavior in the large N limit. For this particular example we used FD to fourth order and in FEM we used continuous basis.
In this paper, we present a method to combine FEM with FSS for electronicstructure calculations near regions of criticality. We note that any method where the wavefunction can be expanded in a complete basis (or set of grid points) should also work, provided that in the infinite basis limit () it becomes arbitrarily close to the the exact wavefunction. The main advantage of finite element methods, over other global basis methods, resides in the fact that matrices obtained by these numerical methods are often banded and tend to be sparse since variables are not coupled over arbitrarily large distances. This is ideal for describing extended systems since it facilitates parallelization pask (). Moreover, due to the polynomial nature of the basis set most integrals involved can be evaluated easily and accurately with a better convergence with respect to the basis set size.
We examine the two electron system with different levels of approximations to the ground state energy, by solving the system using FEM both exactly and using meanfield equations and demonstrate that this formalism is general and that it has potential for more complex systems by solving HartreeFock equations or KohnSham equations in densityfunctional theory. The paper is organized as follows: In Sec. (II) we discuss FEM formalism, in Sec. (III) we discuss the FSS approach to quantum systems, in Sec. (III) and (IV) we combine FSS with FEM to compute calculations of critical parameters for two electron atoms, ,and, lastly in Sec. (V) we discuss the results and applications to extended systems.
Ii Finite Element Method in Electronic Structure Calculations
Computational methods in quantum chemistry typically make use of Slater or Gaussian basis sets to approximate solutions in atomic or molecular systems. The advantages of these methods are that some of the integrals involved in the calculations are analytical and very accurate results can be obtained for these integrals. However, there are problems in extrapolating to the complete basis limit. Plane wave methods, having a complete basis, present some difficulties if nonperiodic boundary conditions are to be used such as in molecule calculations. A common trait of all these methods is that they employ global basis, which can result in some difficulties for the formulation of large problems. First, the matrix formulation becomes dense using global basis and therefore illconditioned for parallelism pask (). While for some complex objects, convergence issues may arise if there are large variations of potential terms in small regions which cannot be captured by the traditional Gaussian basis set that are centered on the nuclei.
The FEM presents an alternative solution to these problems. This method uses a variational method to solve differential equations by discretizing a continuous solution to a set of values in subdomains, called elements, and employs a local basis expansion in terms of polynomials with variable realspace resolution, allowing for convergence to be controlled systematically. The locality of the bases in real space results in sparse and banded matrices, where the number of operations for matrixvector multiplication can be reduced from for dense matrices to for sparse matrices, where N is the dimension of the matrix pask (). Some of the first applications of FEM in engineering date back to 1950, in problems such as elasticity and structural analysis in civil engineering and aeronautical engineering, respectively. In the 1970s, some applications to quantum mechanical problems appeared, but were limited to small systems due to the storage limitation FEMBook1 (); FEMBook2 () Recently, advances in highperformance computing make FEM a viable alternative to the traditional approaches of electronicstructure calculations FEMReview1 (); FEMReview2 (). The sparsity of the global matrices resulting from the FEM makes the parallel formulation more affordable. Alizadean et al, have recently applied a divideandconquer method to the FEMHartreeFock approach for electronic calculations showing a facilitation on parallelization and reduced scaling for larger systems AHM ().
We briefly describe the FEM formulation for a given Hamiltonian to be used in later sections. Let us consider the following oneparticle Hamiltonian:
(1) 
where and are the energy eigenvalue and eigenvector respectively, and V is an arbitrary potential that will represent a sum of terms such as the electronnuclear potential (), Hartree potential () and exchangecorrelation potentials (). We can reformulate this problem into a weak form by taking an inner product of an arbitrary “test function” with Eq (1):
(2) 
This relaxes the problem into a variational form where instead of finding an exact solution in the domain , we find a solution that satisfies Eq.(1) in an average sense everywhere. We then arrive at an equivalent equation by integrating by parts (in real space):
(3) 
where is a boundary term and is the outward unit normal at r. In the absence of external fluxes we can employ Neumann boundary conditions []. Furthermore, in radial coordinates () we can write the surface term in the following form:
(4) 
This surface term vanishes trivially at , while for , it can be removed if we choose a test function that vanishes at . Since Eq (3) is satisfied for any , we can choose and in this manner it reduces the problem to an eigenvalue problem. The wavefunction can be then discretized using the Garlekin approximation FEMBook1 (); FEMBook2 ():
(5) 
where are piecewiselinear continuous functions (i.e. continuous but not necessarily smooth) and are the nodal degrees of freedom. Smoother bases can be constructed by using higher order polynomial ( continuous, ) and by extending the nodal space to include the derivative of the wavefunction at the node, viz, ; this ensues continuity of the wavefunction between the nodes. In the nodal basis, the generalized eigenvalue problem reads
(6) 
where the above terms can be written in terms of local basis representation, viz
The finite element method takes advantage of the strict locality of its basis to yield matrices that are sparse and structured and due to the polynomial nature of the basis most integrals can be easily evaluated, as in global basis approaches. Moreover, since the method is variational, all errors are of the same sign, leading to monotonic convergence often from below pask2 ().
Iii Finite Size Scaling
We now discuss the role of finite size scaling in systems exhibiting critical behavior. Phase transitions are associated with singularities of the free energy that occur only in the thermodynamic limit yanglee1 (); yanglee2 (). Near criticality, fluctuations in the system become correlated over distances of the order of the correlation length , while for second order phase transitions, this correlation length diverges.
The theory of finite size scaling can be used to extract universal parameters describing the transition by studying systems of finite size widom (); barber (); privman (); cardy (); nightingale1 (); Peter1 (); Peter2 (). Any physical quantity () that would diverge in the thermodynamic limit, becomes now “rounded off” in the finite systems (of size L) by a scaling function cardy (); nightingale1 ()
(7) 
where is a reduced critical parameter describing the phase transition and is the correlation length of the system in the thermodynamic limit.
While classical phase transitions are governed by the fluctuation of temperature that tend to minimize the free energy, at zero temperature, these fluctuations are driven by Heisenberg’s uncertainty principle sachdev (). A signature of a phase transition is nonanalyticity of the ground state energy with respect to a parameter in the Hamiltonian. The nonanalyticity can be the result of avoided crossing, or actual level crossing. In the first case, the configuration of the system that minimizes the energy maintains the energy gap () between the first excited energy state and the ground state as the parameter is varied. In the latter case, the energy levels cross so that the energy is now minimized by a configuration that previously corresponded to the first excited state of the system.
To pin down the transition point involves computing the energy accurately near the critical point. However, this can be quite a challenge since most real systems lack an analytical solution to the ground state energy and approximation methods, such as perturbation theory or meanfield methods (HartreeFock , DFT) which are viable options in stable systems, in general are not well defined or become illconditioned near critical points.
This is where the theory of finite size scaling can be recast in the context of the Hilbert space. Any given approximation to the ground state energy will result in a “rounding” effect from the “true” critical behavior, since any approximation to a wavefunction can be expanded in some basis and truncated to some order N; that is,
(8) 
FSS has been used in quantum systems neirotti0 (); serra2 (); kais (); snk1 (); snk2 (); nsk (); qicun (); kais1 (); review (); adv (). In this approach, the finite size corresponds not to the spatial dimension but to the number of elements (N) in a complete basis set used to expand the exact eigenfunction of a given Hamiltonian dipole (); quadrupole (); qicun1 (); qicun2 (); FerronLast (). In the variational approach singularities in different mean values will occur only in the limit of infinite number basis function, viz.,
(9) 
The FSS ansatz in the Hilbert space takes the following form (for expectation values of operators acting on the Hilbert space). There is a scaling function such that
(10) 
where is a critical parameter in the Hamiltonian, is a quantum mechanical operator near criticality, and is a scaling function for each different operator but with a unique scaling exponent. At the critical point, the expectation value is related to as a power law, . In this manner the scaling function will remove the divergence as it is computed in the finite basis.
The asymptotic behavior of operators near the critical point, in this case the ground state energy, can be expressed using critical exponents similarly to classical phase transitions, viz
(11) 
where is a universal critical exponent. In addition to the vanishing energy scale, a diverging length scale also occurs which is associated with the exponential decay of equaltime correlations in the ground state of the system sachdev (). For our case, we focus on the behavior of the wavefunction at threshold, where the correlation length is defined by:
(12) 
In order to apply finite size scaling method to quantum system, we assume a Hamiltonian of the form
(13) 
where is the coupling parameter showing a phase transition at . For potentials linear in ( ), Simon simon () showed that the critical exponent is equal to one if and only if the Hamiltonian has a normalizable eigenfunction with eigenvalue equal to zero. The existence or absence of a bound state at the critical point is related to the type of singularity in the energy. Using statistical mechanics terminology, we can associate a“first order phase transition” with the existence of a normalizable eigenfunction at the critical point. The absence of such a function could be related to a “continuous phase transition.” To obtain the critical parameters, we define the following function that will remove the Ndependence at a critical point:
(14) 
For the energy operator , and using the customary Greek letter for the corresponding exponent we have
(15) 
To arrive at an expression that isolates the critical exponent from numerical calculations, we make use of the HellmannFeynman theorem landau ()
(16) 
At this point it becomes convenient to define another function that will be used in the following sections.
(17) 
At the critical point this function is independent of and and takes the value of . Namely, for and any values of and we have
(18) 
Iv Implementation
In the following, we illustrate the combined approach (FEM+FSS) for the twoelectron atom, but emphasize that the method is general and can be applied to more complex atomic and molecular systems. We seek the minimum charge (critical charge) needed to bind both electrons to the nucleus with charge Z. It is well known that ignoring the electronelectron energy correlation will not predict the stability of chandr (), since the critical charge obtained by such assumption is . Stillinger (1974) Stillinger () took into account the ee correlation by using a nonlinear wavefunction in order to find the minimum energy as the nuclear charge Z was varied. He found a critical charge around . Baker et al (1990) Baker () performed a thorough 401order perturbation calculation to resolve the controversy for the critical charge; they obtained a critical charge of . This is the value we use as reference. In the following sections we illustrate the meanfield equations and exact formulations for the ground state of the twoelectron system, while in the next section we apply FSS to these approximations. Let us start with the Hamiltonian for the twoelectron atoms; it is given by
(19) 
where and are the coordinates of the two electrons referenced from the nuclei. We use Hartree units () throughout in our calculations.
Hartree Fock (HF) Approximation:
In the Hartree approximation the total wavefunction is approximated as a product of single particle wavefunctions is . We can arrive at a selfconsistent form by multiplying the Hamiltonian by the complex conjugate of the aforementioned wavefunction and integrating one coordinate out, we arrive at
(20) 
Using the FEM (see Sec. II) this expression can be reduced to a generalized eigenvalue problem . The Hamiltonian is written in terms of the single electron energy term (i.e. ) and a direct coulomb term defined below, viz
(21) 
Since electrons are fermions, the product of the single particle states needs to be antisymmetric with respect to exchange of the coordinates on the electrons. This can be achieved by writing the wave function as a Slater determinant thijssen ():
(22) 
Searching for the best Slater determinant selfconsistently in the Hamiltonian [Eq.19] leads to the HatreeFock equations thijssen (); Schw (); GG ().
(23) 
where the additional term on the right is known as the exchange Coulomb term, denoted by , and tends to repel states of the same (parallel) spin. We can write this term in the local basis form by looking at the direct Coulomb term and noting that the indices and similarly indices . The operator is obtained from by switching indices (), which effectively switches the coordinates and on two of the wavefunctions and makes this operator nonlocal.
The ground state energy (for electrons) is determined by solving for the ground state eigenfunction for the eigenvalue problem (Eq.23) and computing the expectation value below (to avoid doublecounting terms),
(24) 
Note that for our case, the exchange term does not arise since the ground state of the helium atom does not have parallel spins. In order to improve upon the HartreeFock approximation we can take into account the correlation energy from an approximation to the free electron gas due to Wigner. Here the “correlation” potential can be obtained from Refs. wigner1 (); wigner2 (). Note that the references cited use Rydbergs units, and reference wigner1 () additionally uses . In Hartrees the correlation potential is given by:
(25) 
where is the local density parameter for the atom, that is
(26)  
(27) 
The last expression follows for helium since in the ground state the two electrons have opposite spin, but share the same spatial orbital. In this manner the correlation potential is included selfconsistently and the density is obtained from the sum of the densities for each electron i.e., . Finally, the correlation energy is obtained as
and the value for the total energy is given by:
(28) 
Density Functional Theory
In the HartreeFock approximation, the manybody wavefunction is written in the form of a Slater determinant while the correlation energy is treated separately. This can be done for small systems, but for larger systems this is not an easy task. The theory of Kohn and Sham developed in 1965 simplifies both the exchange and correlation effect of electrons by reducing their effect to a effective potential that is in principle exact. In their seminal paper KS (), they showed that to determine the ground state energy it is sufficient to determine the ground state electron density . By constructing a simple system of noninteracting electrons which yields the same density as the “real” fully interacting system, the ground state energy can be obtained.
The KohnSham equations for two electrons are:
(29) 
These equations are selfconsistent in the density and differ from the HartreeFock equations only in the exchangecorrelation potential which is a function of the total local density:
(30) 
While the exact form of the exchangecorrelation functional is not known, the simplest approximation is the local density approximation (LDA)
(31) 
(32) 
(33) 
It follows from Eq.(30) that the potentials are given by:
(34) 
(35) 
and the ground state energy of the system is obtained from thijssen ():
(36) 
LDA  Total Energy  

FEM  FEM  PerdueWang  FEM  FEM  Exact  
(N=100)  (N=200)  UG (); PW ()  (N=100)  (N=200)  UG ()  
2.813 661  2.821 852  2.834 455  2.886 297  2.904 925  2.903 724  
2.703 352  2.726 064  2.767 389  2.806 848  2.855 856  2.867 082  
6.538 681  6.571 185  6.624 884  6.664 319  6.737 783  6.753 267  
1.973 543  1.978 456  1.995 861  1.020 329  1.026 754  1.024 568  
0.850 814  0.854 086  0.861 740  0.000 000  0.000 000  0.000 000  
0.101 061  0.101 092  0.111 080  0.049 135  0.049 814  0.042 107  
0.563 621  0.565 805  0.570 256  0.957 555  0.963 991  0.903 724 
The KohnSham equations are in general simpler than HF equations, while including correlation effects beyond the HF approximation. Furthermore their computing time grows as a power of the number of the Nelectrons, whereas the exact solution
for the Nelectron Schrodinger equations is exponential in N GG ().
Let us comment on the difference between the two approaches, even though both HartreeFock and KohnSham equations have similar potential terms: [see Eqs.(23) and (29)]. In HF the energy is calculated using the wavefunction explicitly, whereas in DFT the energy is a functional of the density. For example, for the correlation energy we have in the HF approximation:
while for the DFT approximation we have:
Care should be exercised when computing the LDA exchange potential since this term is nonlinear in the density and as a result the normalization does not trivially cancel out when solving the eigenvalue problem. For example, in this work we use the normalization so that a factor of is absorbed into the wavefunction . This in turn changes the density to , and therefore the exchangepotential in LDA becomes
In order to satisfy the boundary conditions as , we use a uniform mesh with a cutoff radius at . The radial equations where solved using continuous polynomials, that is, and ; where the function interpolates the wavefunction from the left node to the right node and similarly for . We solved all the integrals numerically using 10point Gaussian quadratures and in order to solve the generalized eigenvalue problem we have used a fortran subroutine EWEVGE EWEVGE ().
In Table 1 we show the energy components of the KohnSham equations under LDA approximation
and the HartreeFock equations using the Wigner correlation potential for the helium atom.
We note that the reference values for the Total Energy (HF+Wigner correlation) are exact, therefore even
in the infinite basis limit the expression will not converge to the value shown.
We also point out that the PerdueWang values on the correlation energy are obtained by using the
randomphase approximation (RPA) for the free electron gas PW ().
Our convergence value is be limited by how well we can approximate this value.
We have checked our level of convergence by increasing the number of elements to N=1000 in a linear mesh of size
with continuous basis (linear polynomials).
We obtain a total energy of and correlation energy of , and when
compared to the values by PerdueWang (see Table 1 ), the errors are and respectively.
We see that the error in the total energy is dominated by the correlation term arising from using Wigner correlation rather than RPA description. However, taking this difference into account we can get an estimate for the error
of all other terms, viz . In FEM, the error on integrals such as the energy is ,
where h is the element width and p is the order of the polynomial used. Thus, for our example using N=1000
in a uniform grid size with a.u. we have , and since we employ a linear basis, the error we
expect is . In this manner, if we take into account the intrinsic error due to the different
energycorrelation terms we can recover the level of convergence expected by using FEM.
Exact Formulation
We can compare the oneparticle formalism above with an exact formulation of the system. In 1930 Breit showed that the Schrödinger equation for the twoelectron atom becomes separable if one of the coordinate axes is aligned with one of the radius vector Breit (). For the s state, the exact wavefunction depending on six total variables can be reduced to only three variables: the positions of the two electrons and , and the angle between them , viz
(37) 
where
We approximate the solution to the above Hamiltonian using FEM by discretizing each variable above independently as in Levin (). The boundary conditions implemented make the wavefunction vanish at a cutoff radius of a.u, and the angular range is between . Furthermore an adaptive grid was refined closer to the nuclei as shown in Fig.(1) to speed up calculations. In our implementation we use continuous basis, that is Hermite polynomials and to interpolate the left and right nodal values of the wavefunction and additionally and to interpolate the derivatives of the wavefunction, where is the element width and . In this manner, the continuity of the wavefunction is enforced by using this basis. All the integrals involved in constructing the matrices were solved numerically using 10point Gaussian quadratures, and in order to solve the generalized eigenvalue we used a sparse matrix solver in ARPACK++ ARPACK ().
V Results and Discussion
In order to find the critical parameter we approximate the energy gap in the ground state numerically as
where is the binding energy obtained for two electrons and is binding energy for a single electron. The theory of finite size scaling is applied to a heliumlike Hamiltonian of the form,
(38) 
where and are the coordinates of the two electrons referenced from the nuclei, and Z is the charge of the nucleus, this will play the role of the critical parameter. For reference, we note that the transformation takes Eq. (38) into:
(39) 
this form is often used in the literature to obtain the critical parameter (see Stillinger (); Baker ()). Below we outline a general procedure that can be used to combine FSS with numerical methods that expand a given wavefunction in the Hilbert space in any complete basis at a truncated level N.
General Procedure:

Compute the expectation value of and for various values of N (number of elements). In this case the HellmannFeynman theorem (see Eq. 38) becomes useful since in this case is linear in and the derivative is quite simple .

Construct the gamma function (see Eq. 17) for several values of vs . At criticality, the gamma function becomes independent of N and N’ and the curves will cross. In Fig.2 we show the gamma function we obtain by using different number of elements and how it begins to cross for various approximations to the energy.

In order to go to the complete basis set limit , we have to let . The best choice is to keep minimal nightingale1 (); privman (). Then, we proceed as follows: We solve for (N1, N , N+1). This provides us with two gamma curves: and , with minimum . The Crossing of these two curves gives us the pseudocritical parameters and with a truncated basis set number N. Fig. 3 illustrates the N dependence of the crossing.
(a)  (b) 
(c)  (d) 
While finding values and and extrapolating to the limit some numerical fluctuations can arise depending on the approximation used. We discuss some ways of dealing with these fluctuations. First, one can improve the basis used from continuous polynomials to continuous polynomials in order to obtain greater accuracy in the calculations, however, matrices double in size increasing the memory cost as well as the computing time. Second, one can construct a gamma function from a potential in the Hamiltonian that explicitly takes into account all variables. For example in the three dimensional (3D) case we found that it is better to use the transformed Hamiltonian (Eq.39) where the perturbation is given by rather than . This does not change the physics of the system since both expressions can be used to find the critical value: either or , which corresponds to one of the electrons going from a bound state into the continuum. The first perturbation takes into account all three variables when the expectation value is computed whereas the latter introduces numerical errors in the expectation value.
HF  LDA  Total Energy  Exact  Reference Baker (); PSK ()  

critical charge  
1.03114  0.92808  0.90946  0.91857  0.91103  
critical exponent  
0.99971  1.00011  1.00009  1.00082  1.00000 
(a)  (b)  (c) 
Last, another method is to increase the minimum difference used in finding a crossing between and until a systematic crossing is obtained for all values of N; in the LDA approximations we used , since the Gamma curves change very little as the number of elements increases (see Fig.2c), whereas in the 3D formulation we used since in this case the number of elements increases as , and was kept fixed at three. In solving the selfconsistent equations [ Eqs.(23) and (29)], we found that the region below became problematic when updating the selfconsistent equations. This is because as the nuclear charge decreases, the wavefunction decay becomes slower with respect to r and the cutoff radius might lead to unphysical solutions (i.e. solutions which depend on ). To get around this problem, we found it is useful to perturb the wavefunction with the arithmetic average of the current and predicted wavefunctions, rather than just the predicted wavefunction. In this manner, we can avoid converging to a local minima that corresponds to an unphysical solution.
In phase transition theory it is well known that using mean field theory (for dimensionality ) does not give correct values near a critical point. This is indeed the case by using only the HartreeFock approximation (see Table 2 ), however we see that incorporating the exchangecorrelation energy within the LDA approximation improves substantially the critical charge within 2% of the accepted value and within 1% using Total Energy approximation and exact formulation. It is interesting to note that all approaches give a consistent estimate for the critical exponent, regardless of the approximation used. This study demonstrates that Finiteelement method is a viable tool that can be applied in combination with finite size scaling in order to obtain good estimates of critical parameters starting with ab initio calculations for electronic calculations such as DFT and HF and compare well with other methods that are deemed more precise.
We now check the consistency of our approach. In the theory of continuous phase transitions, a quantity diverging in the thermodynamic limit (), becomes a regular function for finite N, viz . We can check that this assumption is indeed satisfied for any value N. In the infinite basis limit we have the following form near criticality: (Eq.11). In this manner, can obtain the derivative of the energy near criticality given that we have a limiting expression for the energy, viz
(40) 
In a finite basis the above expression becomes a power law in the number of basis with the following characteristic exponents.
(41) 
Excluding the trivial case N=0, one can therefore use Taylor expansion around .
(42)  
(43)  
(44)  
(45) 
In this manner we obtain a universal scaling function . Plotting vs the argument of the function would confirm whether the scaling function depends on N or not. Fig. 6 shows that all quantities obtained from different N values do in fact collapse into one curve, and this confirms the scaling hypothesis for FSS with FEM (Eq.10). In order to make the data collapse, it should be emphasized that the values and used to produce Fig.6 correspond to those obtained in Table 2. If other values are used, the curves will not collapse and they will be distributed in a family of curves (N dependent) all deviating away from the unique curve shown in Fig.6. The correlation length exponent was varied until the best collapse curve was found. This occurred around for all implementations.
(a)  (b) 
(c)  (d) 
The solution to a set of HartreeFock or KohnSham equations for large number of atoms becomes the limiting step, and it is here that the FEM can extend the range of sizes that can be investigated by other current means. The efficiency of FEM in the context of largescale calculations derives from the strict locality of the FEM basis in real space, thus minimizing the need for extensive interprocessor communications on massively parallel computational architectures pask (); pask2 (). The work of Pask et al pask (); pask2 () has formulated a realspace DFT approach for condensed matter applications without the need of using expensive Fourier transforms. Impressive results are obtained for GaAs calculations using this approach on DFT+FEM that systematically converges, as the number of elements is increased, to the values obtained using the conventional plane wave calculations, employing Fourier transforms and Ewald summations.
The scope of the dimension attainable with FEM has also been demonstrated by Sterne et al Sterne (). They have employed FEM to compute positron distributions and lifetimes for various defect structures. Because the distribution is nonelectronlike (concentrating away from the nuclei rather than around them), the generality of the FEM basis makes it particularly well suited for such problems. In order to study defects they have modeled a 5488 atom cell; as far as we know, this is one of the largest such calculations reported by any method to date.
Another applications of FEM to largescale ab initio calculations is the work of Tsuchida et al on molecular dynamics simulations of liquids TT (). They apply FEM study of liquid formamide (HCONH2) using a generalized gradient approximation to the exchange and correlation terms to improve the modeling of hydrogen bonds. The use of adaptive coordinate transformations helped to double the resolution of the basis at the C, N and O centers. The liquid was modeled by 100 molecules in a cubic supercell of side 35.5 a.u., for a total of 600 atoms. The scope of this simulation is of the order of the largest one performed by standard PW methods to date and demonstrates the capacity of FEM for such largescale computations already.
A new promising approach developed recently by Mazziotti Mazziotti1 () is the largescale semidefinite programming for manyelectron quantum mechanics. In this approach, the energy of a manyelectron quantum system can be approximated by a constrained optimization of the twoelectron reduced density matrix (2RDM) that is solvable in polynomial time by semidefinite programming (SDP). The developed SDP method for computing strongly correlated 2RDMs is 1020 times faster than previous methods Mazziotti2 (). This approach was illustrated for metaltoinsulator transition of H, with about variables.
These examples along with our current work, suggest that we can combine the Finiteelement method with finite size scaling to obtain critical parameters in extended systems. One possibility is to calculate the phase transition at zero temperature in the 2dimensional electron gas. The ground state energy of the system favors a crystallize state as the density of the electrons is lowered. Tanatar and Ceperey TC () performed a variational Monte Carlo simulations to the total energy using HartreeFock approach
(46) 
They find the gastocrystal transition at . This system would be ideal to test FEM with finite size scaling for quantum criticality in extended systems.
Vi Conclusion
The theory of finite size scaling (FSS) has been demonstrated to be of great utility in analyzing numerical data of finite systems fisher (); widom (); barber (); privman (); cardy (); nightingale1 (); Peter1 (); Peter2 (). Moreover, it has been demonstrated that the theory also works for quantum systems, where the formalism applies to a complete basis expansion in the Hilbert space neirotti0 (); serra2 (); kais (); snk1 (); snk2 (); nsk (); qicun (); kais1 (); review (); adv (). In this study we implemented the finiteelement method which uses a local basis expansion consisting of arbitrary polynomials in real space and with variable resolution. It was demonstrated that this method can be used to compute matrix elements (used in the finite size scaling formalism) accurately and with monotonic convergence with respect to the system size N. We have implemented the combined approach on different approximations to electronic calculations such as HartreeFock approximation (+ electronelectron correlation) and densityfunctional theory (DFT) under the local density approximation and found that even at the simplest level in DFT, the results we have obtained are in good agreement with other more accurate estimates such as large order perturbation theory Baker ().
The finiteelement method has already shown great promise for extended systems since it uses local basis expansion to yield sparse matrices that lend the methodology easily to parallel architectures pask (); AHM (). Furtermore, the KhonSham equations in DFT provides a general method for extended system. Therefore we expect this combined approach can be applied to electron structure in atoms, molecules, and extended systems without loss of generality. Nevertheless, more studies need to be done to confirm the applicability of FSS with FEM for larger systems, where this combined approach can be used to shed some light into the nature of phase transitions.
Acknowledgments
We would like to acknowledge the financial support of ARO and B. W.K. thankfully acknowledges support form the NSF under the grant PHY0969689. We would also like to thank the referee for his helpful feedback and for pointing out the form for the convergence error in FEM.
References
 (1) M.P. Nightingale, Physica 83A, 561 (1976).
 (2) J.L. Cardy, FiniteSize Scaling, (Elsevier Science Publishers B.V. New York 1988).
 (3) M.E. Fisher, in Critical Phenomena, Proceedings of the 51st enrico Fermi Summer School, Verena’s, Italy, edited by M.S. Green (Academic Press, New York, 1971); M.E. Fisher and M.N. Barber, Phys. Rev. Lett. 28, 1516 (1972)
 (4) B. Widom, Critical Phenomena in Fundamental Problems in Statistical Mechanics Ed: E.G.D. Cohen, (Elsevier Publishing Company, NY 1975).
 (5) M.N. Barber, Finitesize Scaling, in Phase Transitions and Critical Phenomena Vol. 8. C. Domb and J.L. Labors eds. (Academic Press, London 1983).
 (6) V. Privman, Finite Size Scaling and Numerical Simulations of Statistical Systems, ( World Scientific,Singapore 1990).
 (7) P. J. Reynolds, H. E. Stanley, and W. Klein, J. Phys. A 11, L199 (1978).
 (8) P. J. Reynolds, H. E. Stanley, and W. Klein, Phys. Rev. B 21, 1980 (1223).
 (9) W. Moy, S. Kais, and P. Serra, Molecular Physics 106, 203 (2008).
 (10) J. P. Neirotti, P. Serra, and S. Kais, Phys. Rev. Lett. 79, 3142 (1997).
 (11) E. Antillon, W. Moy, Q. Wei and S. Kais, J. Chem. Phys. 131, 104105 (2009).
 (12) E. Tsuchida and M. Tsukada, Phys. Rev. B 54, 7602 (1996).
 (13) D. Stiles (private communication).
 (14) J.E. Pask and P.A. Sterne, Modelling Simul. Mater. Sci. Eng. 13, R71 (2005).
 (15) L. R. RamMohan, Finite Element and Boundary Element Applications in Quantum mechanics, (Oxford Uinv. Press, 2002).
 (16) D.W. Pepper and J.C. Heinrich, the Finite element method, (Taylor and francis group, New York 2006).
 (17) O. Taisuke; T. Masayuki, Comp. Phys. Commun., 182, 1245 (2011).
 (18) J.E. Bylaska; M. Holst; J.H. Weare. J. Chem. Theor. Comp. 5, 937 ( 2009).
 (19) R. Alizadegan, K. J. Hit, and T. J. Martinez J. Chem. Phys. 132, 034101 (2010).
 (20) J.E Pask, B.M. Klein, C.Y. Fong, P.A Sterne, Phys Rev B 59 12352 (1999).
 (21) C.N. Yang and T.D. Lee Phys. Rev. 87 404 (1952).
 (22) T.D. Lee and C.N. Yang Phys. Rev. 87 410 (1952).
 (23) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 1999).
 (24) J.P. Neirotto, P. Serra, and S. Kais,Phys. Rev. Lett., 79, 3142 (1997).
 (25) P. Serra, J.P. Neirotti, and S. Kais, Phys. Rev. Lett. 80, 5293 (1998).
 (26) S. Kais, J.P. Neirotti, and P. Serra, Int. J. Mass Spectrometry 182/183, 23 (1999).
 (27) P. Serra, J.P. Neirotti, and S. Kais, Phys. Rev. A 57, R1481 (1998).
 (28) P. Serra, J.P. Neirotti, and S. Kais J. Phys. Chem. A, 102 47, 9518 (1998)
 (29) J.P. Neirotti, P. Serra, and S. Kais, J. Chem. Phys. 108, 2765 (1998).
 (30) Q. Shi and S. Kais, Mol. Phys. 98, 1485 (2000)
 (31) S. Kais and Q Shi, Phys. Rev. A62,60502 (2000).
 (32) S. Kais, and P. Serra, Int. Rev. Phys. Chem. 19, 97 (2000).
 (33) S. Kais, and P. Serra, Adv. Chem. Phys., 125, 1 (2003).
 (34) P. Serra, and S. Kais, Chem. Phys. Lett. 372, 205209 (2003).
 (35) A. Ferron, P. Serra and S. Kais, J. Chem. Phys. 120, 84128419 (2004).
 (36) Q. Shi, and S. Kais, Mol. Phys. 98, 14851495 (2000)
 (37) Q. Shi, and S. Kais, Int. J. Quantum Chem. 85, 307314 (2001).
 (38) A. Ferron, P. Serra, and S. Kais, J. Chem. Phys. 128, 044307 (2008).
 (39) B. Simon J. Functional Analysis 25, 338 (1977).
 (40) L.D. Landau and E.M. Lifshitz, Quantum Mechanics: NonRelativistic Theory. Vol. 3 (3rd ed.) (ButterworthHeinemann, London, 1981).
 (41) S. Chandrasekhar, Astrophysical Journal, 100, 176 (1944).
 (42) F. H. Stillinger and D. K. Stillinger, Phys. Rev. A 10, 1109 (1974).
 (43) J. D. Baker, D. E. Freund, R. D. Hill, and J. D. Morgan, III. Phys. Rev. A 41, 1247 (1990)
 (44) J.M. Thijssen, Computational Physics. (Cambride University Press 1999).
 (45) W. Schweizer. Numerical Quantum Dynamics (Kluwer Academic Publishers, Dordrecht, 2001).
 (46) G. Giuliani, G. Vignale, Quantum Theory of the Electron Liquid. (Cambride University Press 2005)
 (47) H. Mitler, Phys. Rev. 99, 1835 (1955).
 (48) C. W. Ufford,J.G. Thomas Phys. Rev. 133, A121 (1964).
 (49) W. Kohn, L. Sham, Phys. Rev. 140, A1133 (1965).
 (50) C. J. Umrigar and X. Gonze, Phys. Rev. A 50, 3827 (1994)
 (51) J.P. Perdew, Y. Wang, Phys. Rev. B 45, 13244 (1992)
 (52) D. Porezag, http://beam.acclab.helsinki.fi/akrashen/escalc/mathsub.f
 (53) G. Breit, Phys. Rev. 35, 569 (1930)
 (54) F. S. Levin and J. Shertzer, Phys. Rev. A 32, 3285 (1985).
 (55) F. M. Gomes http://www.ime.unicamp.br/chico/arpack++
 (56) C. Geuzaine and J. Remacle. http://geuz.org/gmsh/
 (57) R. Bulirsch, J. Stoer, Num. Math, 6, 413 (1964)
 (58) P.A. Sterne, J.E. Pask and B.M. Klein Appl. Surf. Sci. 149, 238 (1999).
 (59) E. Tsuchida , J. Chem. Phys. 121, 4740 (2004).
 (60) D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011).
 (61) D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004).
 (62) B. Tanatar and D. M. Ceperley, Phys. Rev. B 39, 5005 (1989).