Crystalline Solutions of Kohn-Sham Equations in the Fractional Quantum Hall Regime
CCrystalline Solutions of Kohn-Sham Equations in the Fractional Quantum HallRegime
Yayun Hu, ∗ Yang Ge, ∗ Jian-Xiao Zhang, and J. K. Jain
Physics Department, 104 Davey Laboratory, Pennsylvania State University, University Park, PA 16802
A Kohn-Sham density functional approach has recently been developed for the fractional quan-tum Hall effect, which maps the strongly interacting electrons into a system of weakly interactingcomposite fermions subject to an exchange correlation potential as well as a density dependentgauge field that mimics the “flux quanta” bound to composite fermions. To get a feel for the roleof various terms, we study the behavior of the self-consistent solution as a function of the strengthof the exchange correlation potential, which is varied through an ad hoc multiplicative factor. Wefind that a crystal phase is stabilized when the exchange correlation interaction is sufficiently strongrelative to the composite-fermion cyclotron energy. Various properties of this crystal are examined.
I. INTRODUCTION
The Kohn-Sham (KS) density functional theory (DFT)treats an interacting electronic system by mapping it intoa system of non-interacting electrons moving in an ef-fective single-particle KS potential [1]. The validity ofthe KS scheme relies on the assumption of the formof the universal Hohenberg-Kohn energy functional, orthe exchange-correlation (XC) functional, the search forwhich has motivated extensive studies. Despite the suc-cess of DFT, for example at the level of local densityapproximation (LDA), its application to strongly inter-acting systems has been more challenging and requiresmore sophisticated treatments [2–5].The system of interest to us is the fractional quan-tum Hall effect (FQHE) [6], which occurs when electronsin two dimensions are subjected to a strong perpendic-ular magnetic field which quenches their kinetic energyand as a result enhances the effects of Coulomb inter-action. Very few papers [7–12] have been written ap-plying DFT to the FQHE since its discovery over thepast four decades. The difficulty of applying DFT toFQHE traces back to the construction of KS DFT bymapping into non-interacting electrons, because the non-interacting system possesses a large degeneracy. Thatcan be seen by considering the canonical Hamiltonianthat describes the bulk of the FQHE sample: H LLL = V ee = N (cid:88) i FIG. 1. Schematic illustration of our KS approach for theFQHE. (a) The real system under consideration, taken to beat 1 / V ext . In thebulk, each orbital is occupied, on average, by 1/3 of an elec-tron. The fractional occupation is illustrated by the partialcoloring of the otherwise empty circles that represent unoccu-pied orbitals. (b) The KS spectrum of an auxiliary problemof non-interacting electrons in a KS potential V KS , with therenormalized LLs illustrated by purple lines. The orbitals areeither fully occupied or empty. (c) The occupation configura-tion of the auxiliary system of emergent composite fermions(shown as electrons plus two flux quanta) in the composite-fermion Landau level (called Λ level) spectrum. The integerquantum Hall effect of CFs in a KS ∗ potential V ∗ KS can repro-duce the uniform density of the real system in (a). The DFTtreatment in this paper maps the real system of FQHE to anauxiliary system in (c). fermions. It is therefore unclear how our results relate toprevious studies or to experiments.The plan of the manuscript is as follows. In Sec. II wereview the Kohn-Sham equations. In Sec. III we explainin detail how to numerically solve the Kohn-Sham equa-tions using a finite difference method on a square lattice.As an application of our method, in Sec. IV, we studythe ground state density as a function of the strength ofthe exchange-correlation energy. The effects of tempera-ture, weak disorder, the form of XC energy, and systemsize are also considered. We summarize our findings inSec. V. II. KS ∗ EQUATIONS FOR COMPOSITEFERMIONS We consider the following Hamiltonian of a FQHE sys-tem in a 2D x - y plane with an external potential V ext : H = T ee + V ee + V ext , (2)where T ee = (cid:80) Ni =1 12 m b ( p + ec A ) is the kinetic energyoperator, m b is the band mass of electrons and A is thevector potential due to a uniform external magnetic field B = B e z = ∇ × A ( r ) along the z direction.Following the magnetic-field DFT (BDFT) [20–22], thetotal energy functional E tot of the FQHE system can beexpressed as a functional of the ground state density ρ ( r )as: E tot [ ρ ] = E K [ ρ ] + E xc [ ρ ] + E H [ ρ ] + (cid:90) d r V ext ( r ) ρ ( r ) . (3)Here the Hartree energy E H [ ρ ] takes the standard form: E H ( ρ ) = 12 (cid:90) d r d r (cid:48) ρ ( r ) ρ ( r (cid:48) ) (cid:15) | r − r (cid:48) | . (4)In order to study the FQHE in the LLL, throughoutthis paper we have defined the non-interacting kinetic en-ergy functional to be E K [ ρ ] ≡ (cid:126) ω B (cid:82) d r ρ ( r ) = N (cid:126) ω B ,where ω B = eBm b c is the cyclotron frequency. The elec-tron XC energy functional E xc is defined by Eq. (3). Itcan be equivalently defined through a constrained searchformalism [23, 24] E xc [ ρ ] ≡ min Ψ → ρ ( r ) (cid:104) Ψ | T ee + V ee | Ψ (cid:105) − E K [ ρ ] − E H [ ρ ] , (5)which is further simplified by using (cid:104) Ψ LLL | T ee | Ψ LLL (cid:105) = N (cid:126) ω B in the LLL: E xc [ ρ ] = min Ψ LLL → ρ ( r ) (cid:104) Ψ LLL | V ee | Ψ LLL (cid:105) − E H [ ρ ] , (6)where the many-body wave function Ψ LLL searches foran energy minimum of V ee within the LLL Hilbert space.We adopt Eq. (6) in the following.The above formulation is the standard version of theHohenberg-Kohm (HK) theorem. Here E xc depends onthe external magnetic field B in the BDFT formalism,but is otherwise a universal functional of density thatdoes not depend on the external potential.We next construct an auxiliary KS* system of CFs(rather than electrons), which is not typical of the stan-dard KS scheme, but is within the formulation of the gen-eralized KS scheme [25] and the concepts in the standardKS scheme apply as usual. We imagine that there existsa reference system that consists of non-interacting CFs,whose ground state density is the same as the groundstate density of the FQHE system and is expressed as the sum of the contribution from the occupied KS or-bitals as ρ ( r ) = (cid:88) α c α | ψ α ( r ) | , (7)where c α is the occupation number of the KS orbital la-beled by α . We then define the “non-interacting” kineticenergy T ∗ s [ ρ ] of CFs as T ∗ s [ ρ ] = (cid:88) α (cid:104) ψ α | T ∗ | ψ α (cid:105) , (8)where the kinetic energy operator T ∗ of CFs is T ∗ = 12 m ∗ (cid:16) p + ec A ∗ ( r ; [ ρ ]) (cid:17) . (9)The important physics of CFs is incorporated throughthe density-dependent effective vector potential A ∗ ( r ),or effective magnetic field B ∗ , through ∇ × A ∗ ( r ) = B ∗ ( r ) e z = [ B − ρ ( r ) φ ] e z . (10)As in the standard KS scheme, a further connection be-tween the FQHE system and the KS system of CFs ismade by rewriting E xc as E xc [ ρ ] = T ∗ s [ ρ ] + E ∗ xc [ ρ ] , (11)which also defines E ∗ xc [ ρ ] as the XC energy of CFs. Toclarify, despite the absence of the kinetic energy of elec-trons, the kinetic energy of CFs arises from Coulombinteraction and is included as part of E xc [ ρ ]. In Ref. 12,the CF XC energy was approximated in LDA as: E ∗ xc ( ρ ) = ς (cid:90) d r (cid:104) aν / + ( b − f / ν + g (cid:105) ρ ( r ) , (12)with parameters a = − . b = 0 . f = 0 . g = − . 050 in units of e (cid:15)l B , and ν ( r ) = 2 πl B ρ ( r ) is the localfilling factor, where l B = (cid:113) (cid:126) ceB is the magnetic length.The parameter ς will be used to control the strength ofthe XC potential; ς = 1 corresponds to the choice inRef. 12. The first term aν / in E ∗ xc is chosen to matchwith the known classical value of energy of the Wignercrystal in the limit ν → ν are chosen to fit the electronic XCenergies that are obtained using trial wave functions atthe Jain fillings ν = n/ (2 n +1). This XC form is suited forthe filling factor range 1 / < ν < / 2, but we will use ituncritically for arbitrary filling factors below. The valueof g gives a constant energy offset and does not affectthe KS orbitals or the ground state densities. (We notethe a similar multiplicative factor to tune the strengthof the XC potential has been used in other contexts, forexample for a model hydrogen molecule [27].)Minimization of E tot is achieved by variations with re-spect to the KS orbitals, which leads to the KS* equation: H ∗ ψ α ( r ) = [ T ∗ + V H ( r ) + V ext ( r ) + V ∗ xc ( r ) + V ∗ T ( r )] ψ α ( r ) = (cid:15) α ψ α ( r ) , (13)where the Hartree potential takes the standard form V H ( r ) = e (cid:15) (cid:90) d r (cid:48) ρ ( r (cid:48) ) | r − r (cid:48) | . (14)The CF XC potential is obtained through V ∗ xc ( r ) ≡ δE ∗ xc /δρ ( r ) as: V ∗ xc ( r ) = ς (cid:20) aν / ( r ) + (2 b − f ) ν ( r ) + g (cid:21) . (15)In the KS potential experienced by the CFs, there is anon-standard term V ∗ T that is defined as: V ∗ T ( r ) = (cid:88) α c α (cid:104) ψ α | δT ∗ δρ ( r ) | ψ α (cid:105) , (16)which comes from the density-dependence of the vectorpotential A ∗ inside the kinetic energy operator T ∗ ofCFs. This term is typically much smaller than the otherterms in the KS potential. In particular, the effect of V ∗ T is irrelevant for the topological properties [12, 17]. Itis worth emphasizing that while the kinetic energy op-erator of electrons, T ee , is absent in the KS* equation,the kinetic energy of CFs enters the KS* equation, andincorporates the non-perturbative effect of the Coulombinteraction.At finite temperatures, { c α } and the chemical potential µ are determined by c α = 11 + exp[( (cid:15) α − µ ) /k B τ ] , (17) N = (cid:88) α c α . (18)The occupation number c α reduces to either 0 or 1 in thelimit of zero temperature. III. NUMERICAL PROCEDURE FOR SOLVINGTHE KOHN-SHAM EQUATIONS In this section, we outline the numerical procedureadopted to find the KS solutions. We show how the finite-difference method is implemented on a discretized lattice.We also discuss the algorithm applied in successive iter-ations to achieve convergence. A. Choice for the magnetic vector potential A ∗ We consider a rectangular 2D system of sides L x and L y . We discretize the system into a lattice and labeleach lattice point as r = ( x, y ), where x = ia x , y = ja y ,with i = 1 , , · · · , N x and j = 1 , , · · · , N y and the lattice constants are a x = L x /N x and a y = L y /N y , respectively.This allows a discretization of the physical quantities.The effective magnetic field for CFs is given by B ∗ ( r ) = B ∗ ( r ) e z = [1 − ν ( r )] B e z , (19)which is equivalent to Eq. (10). In order to write downthe KS Hamiltonian explicitly, we pick the symmetricgauge for the bound flux of CFs. The vector potential A ∗ ( r ) reads A ∗ ( r ) = (cid:90) d r (cid:48) B ∗ ( r (cid:48) )2 π | r − r (cid:48) | e z × ( r − r (cid:48) )= (cid:88) r (cid:48) (cid:54) = r B ∗ ( r (cid:48) ) a x a y π | r − r (cid:48) | ( y (cid:48) − y, x − x (cid:48) ) . (20)The sum over r (cid:48) extends over all space. This choice sat-isfies the Coulomb gauge condition ∇ · A ∗ = 0, whichcan be checked explicitly and implies the commutationrelation [ p , A ∗ ] = 0. B. The discretized Hamiltonian The discretized form of the KS Hamiltonian is straight-forward. We show here the form explicitly for the V T term, which is non-standard and also the most complex.We proceed as follows V T ( r ) = (cid:88) α c α (cid:104) ψ α | m ∗ δ (cid:0) p (cid:48) + ec A ∗ ( r (cid:48) ) (cid:1) δρ ( r ) | ψ α (cid:105) (21)= (cid:88) α c α (cid:104) ψ α | em ∗ c δ A ∗ ( r (cid:48) ) δρ ( r ) · (cid:16) p (cid:48) + ec A ∗ ( r (cid:48) ) (cid:17) | ψ α (cid:105) (22)= (cid:126) eB m ∗ c (cid:88) α c α (cid:90) d ¯ r (cid:48) ψ ∗ α (¯ r (cid:48) ) (cid:20) ¯ y − ¯ y (cid:48) | ¯ r (cid:48) − ¯ r | (cid:18) − i ∂∂ ¯ x (cid:48) + ¯ A ∗ x (¯ r (cid:48) ) (cid:19) + ¯ x (cid:48) − ¯ x | ¯ r (cid:48) − ¯ r | (cid:18) − i ∂∂ ¯ y (cid:48) + ¯ A ∗ y (¯ r (cid:48) ) (cid:19)(cid:21) ψ α (¯ r (cid:48) ) , (23)where (cid:104) ψ α | O ( r (cid:48) ) | ψ α (cid:105) ≡ (cid:82) d r (cid:48) ψ ∗ α ( r (cid:48) ) O ( r (cid:48) ) ψ α ( r (cid:48) ), ¯ r = r /l B , and ¯ A ∗ = e A ∗ l B /c . In Eq. (22), (cid:16) p (cid:48) · δ A ∗ ( r (cid:48) ) δρ ( r ) (cid:17) = 0is used, which can be checked explicitly by noticing that ∂∂x (cid:48) δA ∗ x ( r (cid:48) ) δρ ( r ) = − ∂∂y (cid:48) δA ∗ y ( r (cid:48) ) δρ ( r ) = 2 φ π ( y − y (cid:48) )( x (cid:48) − x ) | r (cid:48) − r | , (24)where we have used δ [ B ∗ ( r (cid:48) )] /δρ ( r ) = − φ δ ( r − r (cid:48) ) , (25)and δ A ∗ ( r (cid:48) ) δρ ( r ) = φ π | r (cid:48) − r | ( y (cid:48) − y, x − x (cid:48) ) . (26)The conversion (cid:126) eB m ∗ c = 0 . e α m ∗ (cid:15)l B is used in our nu-merical calculation, where α m ∗ relates the CF mass toelectron mass m e by m ∗ = α m ∗ (cid:112) B [ T ] m e . We take α m ∗ = 0 . 08, which is a good approximation for theo-retical transport gaps [14].The Hartree potential is calculated using a discretizedform of Eq. (14). To avoid the singularity in the self en-ergy at r = r (cid:48) , we replace point charge by a uniformlydistributed charge on a square region of size a x × a y cen-tered around r .In this paper, we consider an external potential V ext generated by a uniform positive background charge in-side a circular region around the origin. For a systemof N electrons, the background charge density is chosenas ρ b = ν b πl B with a radius R b = (cid:113) Nν b l B , where ν b isthe average ion filling factor ν b = 2 πl B ρ b . We make surethat the rectangle L x × L y is chosen to be large enoughso as to comfortably enclose the electron system. C. Numerical procedure for iterations We obtain the self-consistent solution of Eq. (13) us-ing the following iterative procedure. (i) We start withan input density ρ in . (ii) We obtain T ∗ and V ∗ KS ( r ) onthe left-hand side of Eq. (13), diagonalize the Hamil-tonian to obtain the KS ∗ orbitals, and determine theoutput density ρ out = (cid:80) α c α | ψ α | according to Eq. (7),(17) and (18). Note that we work with a fixed parti-cle number, and therefore need to adjust the chemicalpotential suitably in each iteration. (iii) The relativedifference ∆ NN between the input and output ρ , where∆ N = (cid:82) | ρ in − ρ out | d r , is called the absolute difference.We accept ρ out as converged if the relative differences be-tween any of the two output densities for 2000 successiveiterations satisfy ∆ NN < . ρ has not converged, we preparenew input density ρ in by mixing some output density intothe previous input: ρ in → ηρ in + (1 − η ) ρ out , where themixing coefficient is η ≥ . 9. The choice of η close toone helps avoid the so-called occupation sloshing , whichcan occur due to the large degeneracy in our system [28].We iterate the process until convergence is reached. Itis worth mentioning that the calculation of V T requiresthe information of the KS orbitals. V T is set to zero inthe initial input, but in later iterations, it is necessary toalso mix the input and output V T in the same way as themixing of density in each iteration in order to ensure con-vergence. We always start at a sufficiently high tempera-ture, where convergence is straightforward, and slowly goto lower temperatures, while using the converged densityof the previous temperature as the input.The KS Hamiltonian needs to be updated in each it-eration. We notice that a direct calculation of the A ∗ , V T and V H terms on each lattice site using for-loops can be time consuming. To increase efficiency, we have uti-lized the convolution algorithm that is available in theIntel ® Math Kernel Library (MKL) to calculate theseterms. For the diagonalization of the Hamiltonian, weuse the Feast algorithm [29], which is also available inthe MKL and can take a guess of eigenstates as input toincrease efficiency. This is suitable for our purpose be-cause the KS orbitals from the previous iteration serveas a good guess of eigenvectors for the new diagonal-ization. In this paper, the typical system we considerhas N = 40, L x = L y = 35, N x = N y = 210. Thecorresponding Hamiltonian is a sparse matrix with a di-mension of 44100. With η = 0 . 95, the convergence takesseveral thousand to tens of thousands of iterations de-pending on whether the converged density is liquid-likeor crystal-like. Liquid-like solutions are largely uniformin the bulk and converge quickly. In contrast, crystal-likesolutions require significantly larger number of iterationsfor convergence, in order to adjust the position and shapeof the crystalline sites. The corresponding computationtime can range from half a day for liquid to one week forcrystal solutions respectively for the above typical systemsize in a single cluster node with 10 cores. IV. RESULTS The validity of our DFT results depends on the accu-racy of the choice of the XC energy for CFs. We restrictthe approximation of E ∗ xc to the level of LDA and theform of XC energy in Eq. (12) is obtained by fitting to theground state energies of the uniform systems in the fill-ing factor range 1 / < ν < / 2, which is the range wherecomposite fermions carrying two flux quanta are relevant.Since the CFs are weakly interacting, it is reasonable touse a smooth fitting curve for the XC energy, althoughthere has been no study of the exact constraints [30, 31]on the proper choice of the CF XC energy. Eq. (12)is not unique and a different form has been applied inRef. 11. These slightly different fitting forms do not in-fluence the ground state energy as well as the topologicalproperties when the system is in the filling factor range1 / < ν < / 2. However, we will use this form for theXC energy uncritically for all filling factors, and all ofour results are subject to this approximation.An important point for our discussion below is that theform of the exact E ∗ xc also depends on various physicalfactors [16]. For example, one possible factor is LL mix-ing, which is absent in the theoretical limit of very strongmagnetic fields but is relevant for typical magnetic fieldsand can be quite significant. With LL mixing, one canreformulate the problem in terms of electrons still resid-ing in the LLL but with an effective interaction, whichis less repulsive than the Coulomb interaction at shortdistances. A similar correction arises due to finite widthof the quantum well. In principle, one then needs to eval-uate the CF cyclotron energy and CF XC energy for theeffective interaction, which is likely to change the rela- FIG. 2. The ground state density ( ρ ) and the density of states (DOS) of KS* orbitals for a system of N = 40 electrons withaverage filling ν = 1 / k B τ = 0 . e (cid:15)l B . The system size is L x = L y = 35 l B , and N x = N y = 210. Panels(a-d) depict how the density varies as a function of the exchange-correlation (XC) potential, whose strength is tuned by theprefactor ς ; these panels correspond to ς = 0 , , . , 2. The density is quoted in units of (2 πl B ) − . The FQHE liquid evolvesinto a crystal with increasing XC energy. Panels (e-f) show the corresponding density of states, ρ E . The cumulative state count N ( E ) gives the number of Kohn-Sham orbitals below energy E . Only the lowest 85 orbitals are shown. The dashed line marksthe location of the Fermi energy. All the energies are in units of e (cid:15)l B . tive importance of the two terms. We have not made arealistic determination of these effects.We will tune the strength of the XC energy E ∗ xc inEq. (12) by varying ς . Notice that E ∗ xc is negative, so theXC potential increases in magnitude when ς > 1. Werefer to ς → ∞ and ς → A. Appearance of a crystal phase We consider a problem with rotational symmetry bychoosing the external potential to be generated by a uni-form positive charge in a circular region around the ori-gin, as explained in Sec. III C. The ground state densitiesfor a system of N = 40 are shown in Fig. 2(a-d) for cer-tain choices of ς . The corresponding density of states(DOS) are shown in Fig. 2(e-h). Here we assume a smalltemperature of k B τ = 0 . e (cid:15)l B , which is less than 10%of the cyclotron gap of ΛLs and is useful for finding con-verged solutions. (Temperature dependence of KS solu-tions is discussed later.) In the weak XC energy limit( ς = 0), the ground state density of electrons almost per-fectly screens the background density, due to the dom-inance of the Hartree term and the external potential.The spectrum of the KS solutions shows the formationof ΛLs, where the positions of the lowest two ΛLs can beseen from the two peaks in the DOS plot in Fig. 2(e). In order to obtain a smooth curve for DOS, we havereplaced the δ ( E − (cid:15) i ) in the standard definition of DOS, ρ E = (cid:80) i δ ( E − (cid:15) i ) (where (cid:15) i is the eigenvalue of KSorbitals sorted by (cid:15) < (cid:15) < . . . , increasing with i ), by aGaussian to define: ρ E = (cid:88) i exp[ − ( E − (cid:15) i ) /σ ] / √ πσ , (27)where σ = 0 . e (cid:15)l B throughout this paper. The discretepoints of { (cid:15) i } can be seen from the plot of the cumulativestate count N ( E ) = (cid:90) E −∞ (cid:88) i δ ( E (cid:48) − (cid:15) i ) dE (cid:48) . (28)When ς = 1 in Eq. (12), the total density showsstronger oscillations near the edge but still respects arotational symmetry. In particular, the density profilenear the edge of the system first shoots up before itcomes down to zero, which is also seen in results fromexact diagonalization (ED) in a rotationally symmetricsystem [32]. For the stronger XC potential of ς = 1 . ς = 2. Theformation of crystalline structures breaks the rotationalsymmetry of the system, which is allowed in our numeri-cal method where no symmetry is assumed. (In contrast,the calculations in Ref. 12 choose the angular momentumas a good quantum number and reduce the 2D system toeffectively a 1D system along the radial direction; theresults therein are rotationally symmetric by construc-tion regardless of the strength of the XC potential or thechoice of the CF mass. Rotation symmetry is also im-posed in the ED calculation in the disk geometry [32].)The results are stable and driven by the XC energy. Wediscuss in Appendix A that our results, in particular theappearance of a crystal phase, are not a numerical arti-fact of discretization and lattice configuration. B. Nature of the crystal phase The competition between the correlated Wigner crys-tal state and the liquid state has been studied theo-retically in many articles [19, 33–39]. In particular,Ref. 19 studies the role of LL mixing and finds that theFQHE liquid yields to a crystal when the LL mixing islarge. However, there is an important difference betweenthe crystal found in that study and that in our study.In Ref. 19, the n/ (2 n + 1) FQHE liquid of compositefermions carrying two vortices freezes into an electron crystal with increasing LL mixing. In our study, on theother hand, we obtain a CF crystal of composite fermionscarrying two vortices [38, 39].The Wigner crystallization of a system of 2D electronsin a circularly symmetric confining potential has beenstudied in Ref. 37. They find that crystallization oc-curs in two stages: first in the radial ordering, and thenin the angular ordering. The situation is similar to ourfindings, though we have CF crystals rather than electroncrystals. We calculate the evolution of variances in bulkdensity along the radial and the azimuthal directions re-spectively as we increase ς (results not shown). For asmall ς ( ς < ς is larger than a threshold value,the radial variance first increases significantly; beyond agreater threshold, the azimuthal variance also increasesabruptly. This indicates a two-stage crystallization inour results.It is interesting to ask if our crystal is an example ofthe so-called Hall crystal [40–42], which is the quantumHall effect counterpart of the putative supersolid phaseof He atoms. Because our crystal is a correlated crystalof composite fermions, it may appear to be a promisingcandidate for the Hall crystal phase. One feature thatmay distinguish the Hall crystal from the Wigner crystalis that, in the former, the number of particles per unitcell is not necessarily an integer [43]. We find, for allcases we have studied, that the number of crystal sitesin our KS solution at the smallest temperature is equalto the number of composite fermions (which also justi-fies the term crystal rather than a charge density wave).Another character of the Hall crystal is that, similarlyto the Hall liquid state, it hosts chiral edge states. Wefind that for our crystal, there are no gapless edge states; this is indicated by the presence of a gap at the chemi-cal potential. We thus conclude that our crystal phase isgenerically not a Hall crystal, but we do not rule out thepossibility that the Hall crystal state could be stabilizedfor some forms of the XC energy. C. Effects of temperature, disorder, form of XCinteraction and system size We ask how temperature influences the density in ourcalculation. Fig. 3 depicts the evolution of a system with ς = 1 . k B τ = 0 . e (cid:15)l B , the system is crystalline(panel d). As the temperature is raised, the system meltsfrom the edge into the bulk and becomes a liquid-likestate when k B τ = 0 . e (cid:15)l B , which is approximately thevalue of the ΛL gap ( ≈ . e (cid:15)l B ), as can be seen fromthe DOS plot in Fig. 3(e). It has been proposed thatthe re-entrant of a solid state can occur in certain pa-rameter regimes when the temperature of a liquid stateis raised [36]; we have not explored that physics in ourcalculations.Next, we test the stability of the results against a weakdisorder. We consider onsite disorder by adding to theexternal potential a term (cid:80) i δV ext ( r i ), where V ext ( r i ) israndomly chosen according to a uniform distribution inthe range [ − W, W ], where W is the strength of disorder.In the ν = 1 / W = 0 . e (cid:15)l B , although we expect thatthe phase boundary will be slightly modified by disor-der [35].One may ask how the detailed form of the XCenergy/potential influences the results. For that weconsider another form of the XC energy E ∗(cid:48) xc = ς (cid:82) d r [ − . ν . − . ν ] ρ ( r ), which gives an XC po-tential V ∗(cid:48) xc = ς ( − . ν . − . ν ). (These forms are dif-ferent approximations for the exact energies in the range1 / < ν < / 2, but have significant differences outsidethis range.) The qualitative behavior, namely a liquid forsmall ς and a crystal at large ς is also seen for the newXC potential. However, the phase boundaries are differ-ent; for example, the low-temperature KS solution is acrystal for V ∗(cid:48) xc with ς = 1. (We note that both forms ofXC energy produce a uniform liquid state in the bulk for ς = 1 when the system is constrained to be rotationallysymmetric [12].)We have also studied the effects of the CF mass m ∗ ,and the density rings and crystal sites emerge for a large m ∗ that decreases the ΛL gap. This suggests that theformation of ΛLs is also important for the stability ofthe liquid phase.We have also investigated the behavior as a functionof the system size. Away from the transition region, wefind that the nature of the ground state is not sensitive tosystem size. This is illustrated in Fig. 4. Here, the quali- FIG. 3. The temperature dependence of the ground state density ( ρ ) and the density of states (DOS) of KS* orbitals for asystem with N = 40 particles at average filling ν = 1 / 3. We choose the XC potential with ς = 1 . 5. The density is quoted inunits of (2 πl B ) − . The panels (a-d) correspond to temperatures k B τ = 0 . , . , . , . 001 in units of e (cid:15)l B . The panels (e-f)show the corresponding density of states and the cumulative state count for Kohn-Sham solutions in (a-d). Only the lowest100 orbitals are shown. Other parameters are the same as those in Fig. 2. tative features of the solution in Fig. 2(c) are retained insmaller systems. For the systems with the same numberof particles in Fig. 4, the solutions remain liquid-like for ς = 1 and crystal-like for ς = 2 (results not shown). V. CONCLUSIONS We have studied how the strength of the XC potentialbetween composite fermions dictates their state. For thispurpose we develop a numerical procedure to solve theKS equations of CFs in a fashion that allows for crys-talline solutions. Our primary finding is that the stateevolves from a liquid-like state to a crystal-like state asthe strength of XC energy increases.We mention again that our study is not to be taken asa quantitative treatment of the physics of crystals in theFQHE regime. A notable limitation is that we only con-sider states of composite fermions carrying two vortices,and do not consider the possibility of a crystal or a liquidof electrons, or of composite fermions with greater num-ber of attached vortices (as might be relevant in regionsof small densities).An obvious direction for future study will be to buildbetter XC potentials that apply to a larger range of fill-ing factors and also include the effects of finite thicknessand Landau level mixing. It is possible that the strengthof the XC energy of composite fermions relative to theircyclotron energy may also depend on the filling factor,which may be relevant to the formation of a crystal atlow fillings. It would be interesting to apply the DFT FIG. 4. The ground state density ρ for different particle num-bers N = 5 , , , , , 30 with ς = 1 . ν b = 1 / 3. The qualitative featuresof a crystalline structure in the bulk and a liquid-like ringalong the edge are consistent with those in the large systemof N = 40 shown in Fig. 2(c). We have used L x = L y = 30 l B , N x = N y = 180, k B τ = 0 . e (cid:15)l B . Other parameters are thesame as those in Fig. 2. method to study the edge structure, the effect of dis-order and/or anisotropy, spin physics, screening, and offractional quantum Hall effect in mesoscopic devices. FIG. 5. This figure shows the wave functions of single-particle Kohn-Sham (KS) orbitals for ς = 1 , . , 2. For each value of ς ,we show eight orbitals of increasing energy, with state indices i = 1 , . . . , 57. The numbers in the parenthesis represent the pair( ς, i ). The height and color in the plot represent the magnitude and phase the KS orbitals, respectively. The broken rotationalsymmetry in the magnitude of a wave function indicates that angular momentum L z is no longer a good quantum number. Inthe liquid-like phase in (a-h), the expectation value of L z (which can be estimated from the number of phase windings overthe azimuthal angle) is positively correlated with the average radius of an orbital. This correlation is absent in the crystal-likephases, where KS orbitals are delocalized over a few crystalline peaks. We choose the background charge at ν b = 1 / N = 40and k B τ = 0 . e (cid:15)l B . Other parameters are the same as those in Fig. 2. ACKNOWLEDGEMENTS Y. H. thanks Junyi Zhang and Jiabin Yu for helpfuldiscussions. J.K.J. thanks Steve Kivelson for an insight-ful discussion. The work at Penn State was made possibleby financial support from the US Department of Energyunder Grant No. DE-SC0005042. Y. H. acknowledgespartial financial support from China Scholarship Coun-cil. Computations for this research were performed onthe Pennsylvania State University’s Institute for Com-putational and Data Sciences’ Roar supercomputer. Appendix A: Discussion of numerical stability of theresults In exact diagonalization studies, the ground state ofFQHE on a disk geometry is assumed to be an eigenstateof angular momentum, which therefore is rotationallysymmetric. The liquid or crystal nature of the groundstate can be seen by studying the structure of the paircorrelation function [39]. However, once the rotationalsymmetry is spontaneously broken, the system can pickone of the infinitely many possible ground state configu-rations that no longer preserves angular momentum as agood quantum number. These configurations are relatedto each other by a rotation around the origin. Ideallythis is expected to be the case for our KS solutions. (We refer the readers to Ref. 44 for discussions of symmetry-breaking solutions in DFT calculations.) However, ourchoice of the square open boundary and our lattice dis-cretization effectively break the rotational symmetry. Wediscuss in this appendix various numerical tests to showthat these do not affect the nature of the state in anessential manner.We have examined the effect of the square open bound-ary. We have found that adding a circular potential wallof infinite height that is internally tangent to the squareboundary, or expanding the size of the square boundaryby a factor of two, leaves the essential features of thesolutions unchanged.Discretization reduces the rotational symmetry ofspace into a four-fold rotation C and the mirror reflec-tion, which are exact symmetries regardless of the latticeconstant. Interestingly, these symmetries can be brokenin the KS solutions, as is evident in the density plots ofcrystalline solutions [for example, see Fig. 3(c)]. Even inthe liquid solutions, the C symmetry is broken by the KSorbitals, as shown in Fig. 5. In fact, the KS orbitals canform rather complex structures. In the liquid-like phase,each KS orbital retains a ring-like shape around the ori-gin but no longer preserves the angular momentum as agood quantum number. In the crystal-like phase, the KSorbitals typically do not respect any of the symmetries.We note that the KS orbitals are delocalized over severalcrystalline sites, hinting at a highly correlated nature of0the crystal.We have tested that for any solution of the KS equa-tion, the densities related by the C or the mirror sym-metries are also valid solutions. Which of the degeneratesolutions is obtained depends on the initial input. Fur-thermore, for the crystal phase, depending on the initialconditions, we can also obtain solutions that are not re-lated to one another by rotation, suggesting the presenceof many nearly degenerate solutions in the continuumlimit. Nonetheless, all of the solutions for a given set of parameters are in the same phase.We need to make sure that our lattice is fine enough tocapture the physics in the continuum limit. To reduce thenumerical expense, the results above are obtained with alattice resolution of a x = a y = l B / 6. We have tested thatgoing to a resolution of l B / 10 or l B / 15 does not alter theresults appreciably. In the crystalline phase, the numberof nearly degenerate solutions increases as we go to finerlattices, but for all cases that we have studied, going to afiner lattice does not change the liquid or crystal natureof the solution. [1] G. Giuliani and G. Vignale, Quantum Theory of the Elec-tron Liquid (Cambridge University Press, The EdinburghBuilding, Cambridge CB2 2RU, UK, 2008).[2] M. Seidl, J. P. Perdew, and M. Levy, Phys. Rev. A ,51 (1999), URL https://link.aps.org/doi/10.1103/PhysRevA.59.51 .[3] P. Gori-Giorgi, M. Seidl, and G. Vignale, Phys. Rev. Lett. , 166402 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.166402 .[4] F. Malet and P. Gori-Giorgi, Phys. Rev. Lett. ,246402 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.246402 .[5] H. S. Yu, S. L. Li, and D. G. Truhlar, The Journal ofchemical physics , 130901 (2016).[6] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Phys. Rev.Lett. , 1559 (1982), URL http://link.aps.org/doi/10.1103/PhysRevLett.48.1559 .[7] M. Ferconi, M. R. Geller, and G. Vignale, Phys. Rev. B , 16357 (1995), URL http://link.aps.org/doi/10.1103/PhysRevB.52.16357 .[8] O. Heinonen, M. I. Lubin, and M. D. Johnson, Phys. Rev.Lett. , 4110 (1995), URL http://link.aps.org/doi/10.1103/PhysRevLett.75.4110 .[9] Y.-H. Zhang and J.-R. Shi, Phys. Rev. Lett. ,016801 (2014), URL http://link.aps.org/doi/10.1103/PhysRevLett.113.016801 .[10] Y.-H. Zhang and J.-R. Shi, Chinese Physics Let-ters , 037101 (2015), URL http://stacks.iop.org/0256-307X/32/i=3/a=037101 .[11] J. Zhao, M. Thakurathi, M. Jain, D. Sen, and J. K. Jain,Phys. Rev. Lett. , 196802 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.196802 .[12] Y. Hu and J. K. Jain, Phys. Rev. Lett. ,176802 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.176802 .[13] J. K. Jain, Phys. Rev. Lett. , 199 (1989), URL http://link.aps.org/doi/10.1103/PhysRevLett.63.199 .[14] J. K. Jain, Composite Fermions (Cambridge UniversityPress, New York, US, 2007).[15] B. I. Halperin and J. K. Jain, eds., Fractional QuantumHall Effects New Developments (World Scientific, 2020),https://worldscientific.com/doi/pdf/10.1142/11751,URL https://worldscientific.com/doi/abs/10.1142/11751 .[16] R. Price and S. Das Sarma, Phys. Rev. B ,8033 (1996), URL https://link.aps.org/doi/10.1103/PhysRevB.54.8033 . [17] Y. Hu, G. Murthy, S. Rao, and J. K. Jain, Phys. Rev. B , 035124 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.035124 .[18] A. C. Archer, K. Park, and J. K. Jain, Phys. Rev. Lett. , 146804 (2013).[19] J. Zhao, Y. Zhang, and J. K. Jain, Phys. Rev. Lett. , 116802 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.116802 .[20] C. J. Grayce and R. A. Harris, Physical Review A ,3089 (1994).[21] W. Kohn, A. Savin, and C. A. Ullrich, International jour-nal of quantum chemistry , 20 (2004).[22] E. I. Tellgren, S. Kvaal, E. Sagvolden, U. Ekstr¨om,A. M. Teale, and T. Helgaker, Phys. Rev. A ,062506 (2012), URL https://link.aps.org/doi/10.1103/PhysRevA.86.062506 .[23] M. Levy, Proceedings of the National Academy of Sci-ences , 6062 (1979).[24] E. H. Lieb, Int. J. Quantum Chem. , 243 (1983).[25] A. Seidl, A. G¨orling, P. Vogl, J. A. Majewski, andM. Levy, Phys. Rev. B , 3764 (1996), URL https://link.aps.org/doi/10.1103/PhysRevB.53.3764 .[26] L. Bonsall and A. A. Maradudin, Phys. Rev. B , 1959(1977).[27] M. Holst, H. Hu, J. Lu, J. L. Marzuola, D. Song, andJ. Weare, arXiv preprint arXiv:1902.03497 (2019).[28] N. Woods, M. Payne, and P. Hasnip, Journal of Physics:Condensed Matter , 453001 (2019).[29] E. Polizzi, Phys. Rev. B , 115112 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.79.115112 .[30] S. Pittalis, C. R. Proetto, A. Floris, A. Sanna, C. Bersier,K. Burke, and E. K. U. Gross, Phys. Rev. Lett. ,163001 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.107.163001 .[31] J. W. Dufty and S. B. Trickey, Phys. Rev. B ,125118 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.125118 .[32] E. V. Tsiper and V. J. Goldman, Phys. Rev. B ,165311 (2001), URL https://link.aps.org/doi/10.1103/PhysRevB.64.165311 .[33] P. K. Lam and S. M. Girvin, Phys. Rev. B , 473 (1984).[34] D. Levesque, J. J. Weis, and A. H. MacDonald, Phys.Rev. B , 1056 (1984).[35] R. Price, P. M. Platzman, and S. He, Phys. Rev. Lett. , 339 (1993).[36] P. M. Platzman and R. Price, Phys. Rev. Lett. , 3487(1993). [37] A. V. Filinov, M. Bonitz, and Y. E. Lozovik, Phys.Rev. Lett. , 3851 (2001), URL https://link.aps.org/doi/10.1103/PhysRevLett.86.3851 .[38] G. S. Jeon, C.-C. Chang, and J. K. Jain, Journal ofPhysics: Condensed Matter , L271 (2004).[39] G. S. Jeon, C.-C. Chang, and J. K. Jain, Phys. Rev.B , 241304 (2004), URL https://link.aps.org/doi/10.1103/PhysRevB.69.241304 .[40] S. Kivelson, C. Kallin, D. P. Arovas, and J. R. Schrieffer,Phys. Rev. Lett. , 873 (1986), URL https://link.aps.org/doi/10.1103/PhysRevLett.56.873 .[41] B. I. Halperin, Z. Teˇsanovi´c, and F. Axel, Phys. Rev. Lett. , 922 (1986), URL https://link.aps.org/doi/10.1103/PhysRevLett.57.922 .[42] E. Fradkin and S. A. Kivelson, Phys. Rev. B ,8065 (1999), URL http://link.aps.org/doi/10.1103/PhysRevB.59.8065 .[43] Z. Teˇsanovi´c, F. Axel, and B. Halperin, Physical ReviewB , 8525 (1989).[44] J. P. Perdew, A. Ruzsinszky, J. Sun, N. K. Nepal, andA. D. Kaplan, Proceedings of the National Academy ofSciences118