Energy minimization of paired composite fermion wave functions in the spherical geometry

We perform the energy minimization of the paired composite fermion (CF) wave functions, proposed by M\"oller and Simon (MS) [PRB 77, 075319 (2008)] and extended by Yutushui and Mross (YM) [PRB 102, 195153 (2020)], where the energy is minimized by varying the CF pairing function, in the case of an approximate model of the Coulomb interaction in the second Landau level for pairing channels $\ell = -1, 3, 1$ which are expected to be in the Pfaffian, anti-Pfaffian and particle-hole symmetric (PH) Pfaffian phases respectively. It is found that the energy of the $\ell = -1$ MS wave function can be reduced substantially below that of the Moore-Read wave function at small system sizes, however, in the $\ell = 3$ case the energy cannot be reduced much below that of the YM trial wavefunction. Nonetheless, both our optimized and unoptimized wavefunctions with $\ell=-1,3$ extrapolate to roughly the same energy per particle in the thermodynamic limit. For the $\ell = 1$ case, the optimization makes no qualitative difference and these PH-Pfaffian wave functions are still energetically unfavourable. The effective CF pairing is analyzed in the resulting wave functions, where the effective pairing for the $\ell = -1, 3$ channels is found to be well approximated by a weak-pairing BCS ansatz and the $\ell = 1$ wave functions show no sign of emergent CF pairing.


I. INTRODUCTION
Well over 30 years since its discovery [1], the precise nature of the fractional quantum Hall effect in the halffilled second Landau level in GaAs heterostructures still remains elusive.The three leading candidate phases of matter, all of which could potentially host non-abelian anyons, are the Pfaffian [2], the anti-Pfaffian [3][4][5] and the particle-hole symmetric Pfaffian (PH-Pfaffian) [6].All three of these phases can be understood as paired composite fermions (CFs) at pairing channels ℓ = −1 [7][8][9], ℓ = 3 [10] and ℓ = 1 [10,11] respectively.Numerical studies have pointed toward either the Pfaffian or anti-Pfaffian as being the most likely to occur in this setting [12][13][14][15][16][17][18][19].However, experimental measurements of the heat conductivity along the edges of these systems [20] and measurements of noise along interfaces of these systems with other quantum Hall states [21] are more consistent with the PH-Pfaffian state.We are thus left with an apparent contradiction.
Several proposals have been made in an attempt to resolve this contradiction with one possibility being that the PH-Pfaffian phase could be induced by disorder [22][23][24] and another being the incomplete thermal equilibration of the edge modes affecting the measurements of heat conductivity along the edge [25][26][27][28].While these lines of inquiry might be physically relevant, neither approach has yet provided a fully acceptable explanation of all the experimental facts.
At the same time, it is worth noting that the numerical work that has been performed on the 5/2 state is necessarily limited to finite-sized systems.Exact diagonalization, and even DMRG, are limited to fairly modest system sizes.While it seems a bit unlikely that the results obtained at small sizes (particularly that the PH-Pfaffian is unfavored) will change at larger system sizes, given the enduring conflict with experimental observation, it is worth exploring this possibility.To fully resolve this issue it would then be useful to be able to perform numerical simulations at larger system sizes, than previously studied, with accurate numerically tractable trial wave functions, which is what we shall pursue in the current work.
Whilst the Moore-Read (MR) wave function, which is a representative of the Pfaffian phase, is numerically tractable to large system sizes, it does have a relatively low overlap with the exact Coulomb ground state, in the second Landau level (LL), compared with other trial wave functions for other states in the lowest Landau level (LLL) [12,13].Furthermore, both the representative trial wave functions for the anti-Pfaffian, taken to be the particle-hole conjugate of the MR wave functions denoted MR, and the PH-Pfaffian, proposed in Ref. [11], are only numerically tractable at small system sizes.There is also mounting evidence that the current representative trial wave function for the PH-Pfaffian does not in fact represent a gapped phase of matter [29][30][31].
Motivated in part by this poor overlap of the MR wave function with the exact ground state, Möller and Simon (MS) proposed a numerically tractable trial wave function for the Pfaffian phase, denoted MS −1 , which can be interpreted as a wave function of paired CFs at pairing channel ℓ = −1, where the pairing function can be varied [32].By varying the pairing function so as to minimize the energy of the wave function MS where then able to obtain a more accurate approximation to the exact ground state.This MS construction was later extended by Yutushui and Mross (YM) to CF paring channels ℓ = 3 (denoted MS 3 ) and ℓ = 1 (denoted MS 1 ) which are expected to be in the anti-Pfaffian and PH-Pfaffian phases respectively.YM studied these wave functions with fixed pairing functions (i.e. with the variational pa-rameters fixed), where they found the ℓ = 3 wave function to be an accurate representative of the MR wave function and found the ℓ = 1 wave function to show no sign of emergent paring of CFs in the density-density correlation function.In their work, YM actually proposed two versions of these wave functions: the single-particle projected and so-called "pair-projected" wave functions.
Here we shall only be interested in the single-particle projected wave functions as these are more numerically tractable.We denote the single-particle projected fixed parameter wave functions proposed by YM at pairing channel ℓ by YM ℓ .
A comparison of the energetics of the MS ℓ versus the YM ℓ wave functions have not yet been performed, and it is not known if the MS ℓ wave functions with their energy minimized can offer significantly better approximations to the corresponding ground states at ℓ = 1, 3. In particular, minimising the energy of the MS 1 wave functions could, in principle, produce PH-Pfaffian trial wave functions with energies significantly closer to the ℓ = −1, 3 trial wave functions.Performing such optimizations comes with a practical challenge in that they must be done in an efficient way in order to access larger system sizes.
There has also been some recent interest in understanding if the effective CF pairing in these phases can be approximated by some weak pairing BCS-type description, where the BCS gap parameter can be estimated [33].As well as offering some physical insight into these phases of matter, it also has some practical use, by allowing, for example, to observe when the system is transitioning from CF pairing to the CF Fermi liquid [32].The precise BCS weak pairing description has only so far been studied for the ℓ = −1 pairing [32,33] and it is not known if such a description is accurate for the ℓ = 3, 1 cases.
In this work, we will show how the energy of the MS ℓ wave functions can be minimized at larger system sizes, in the spherical geometry, where we will present the results of this optimization for the case of an approximate model of the Coulomb interaction in the second LL in the absence of LL mixing.At small system sizes, we find, in agreement with Ref. [32], that the energies of the MS −1 wave functions can be reduced substantially below the corresponding energy of the MR wave function, with the optimized MS −1 wave functions showing a muchimproved overlap with the corresponding exact ground state.It is further observed that the optimized MS 3 wave functions offer considerably less energy reduction over the corresponding zero-parameter trial wave functions in comparison to the MS −1 wave functions, although the optimized MS 3 still show a noticeable improvement in the overlap with the corresponding exact ground state compared with the YM 3 wave function.Furthermore, the amount by which the energy of the MS 1 wave functions can be minimized compared with the YM 1 wave functions is negligible where we find both wave functions to be energetically unfavourable.We demonstrate that the effective CF pairing in the optimized MS ℓ and YM ℓ wave functions, for ℓ = −1, 3, can be well approximated by a weak pairing BCS type description, where, for finite size systems, the MS ℓ wave functions show stronger pairing than the corresponding YM ℓ wave functions.However, whether we consider YM ℓ or MS ℓ with ℓ = −1, 3 or MR wave functions, within our numerical error these all extrapolate to roughly the same energy per particle in the thermodynamic limit.Finally, further pathologies of the MS 1 wave functions are found where no evidence of emergent pairing between the CF (at this pairing channel) is found even after some optimization.
The MS ℓ and YM ℓ wave functions are introduced in Sec.II.Then in Sec.III we present the approximate LLL model of the Coulomb interaction in the first excited LL.Finally, in Sec.IV the results of these energy minimizations are presented in the case of the approximate model of the Coulomb interaction in the second LL, where it is also shown how the effective weak-pairing BCS description can be extracted from these wave functions.The optimization algorithm used in this work is detailed in Appendix B.
Throughout this work, we will assume the ν = 5/2 system to be such that the LLL is completely full for the spin-up and spin-down electron orbitals and that the second LL is at half-filling with the electrons being spin polarized.This second LL system is then mapped to the LLL in the usual way (see Sec. III).Note that we will not include any Landau level mixing in the calculations presented in this work.

II. PAIRED CF WAVE FUNCTIONS ON THE SPHERE
In the BCS theory of superconductivity [34] one typically starts with the mean-field Hamiltonian for a system of fermions, which we will take to be in two spatial dimensions, which takes the form and is assumed to be a reasonable approximation to the actual Hamiltonian of the system at low energies, where ∆ k is known as the gap function and ε k is the kinetic energy relative to the Fermi level of a single particle state labelled by k (i.e.ε k = E k −µ with µ the chemical potential).The unnormalised ground state of this Hamiltonian is given by, where Note that for this case of spinless fermions g k must be an odd function of k, g −k = −g k .In real-space this can be expressed as (3) |Ψ BCS ⟩ is physically interpreted as a paired state, where particles near the Fermi level become bound into pairs with g(r) being referred to as the pairing function.
The average occupation of the orbital labelled by k, n k , can be expressed as ).This can be expressed in terms of ε k and ∆ k as For a rotationally symmetric microscopic Hamiltonian, the gap function is some eigenstate of the rotation operator.In two dimensions this would mean that under a rotation by angle ϑ the gap function transforms as ∆ k → e iℓϑ ∆ k , where ℓ is known as the pairing channel.For spinless fermion cases we are considering here ∆ k is an odd function of k, as g k is, and a typical ansatz for ∆ k is ∆ k = ∆|k|e iℓθ , where θ is the angle from the x-axis in k-space.For |r| ≫ ∆/µ the pairing function corresponding to this gap function is g(r) ∝ e iℓθ |r| (see Appendix A of Ref. [35]).When this long-distance form of g occurs the fermions are said to be in a weak-pairing phase [9].
This BCS mean-field wave function can be used to create a trial wave function for a fixed number of fermions N , by projecting it to the space of states with N particles, with the projector written as P N .We then define |Ψ N ⟩ ≡ P N |Ψ BCS ⟩.The wave function of |Ψ N ⟩ is given by Ψ Hence, the average orbital occupations are given by n ⟩. Now we move to a system of N spinless fermions moving on a sphere with N ϕ = 2Q = 2(N − 1 + q) flux quanta passing through its surface, where q is on the order of one (i.e.does not scale with N ) and it is presumed that the magnetic field is strong enough that the fermions are confined to the LLL.Thus, in the thermodynamic limit, this system is at filling fraction Let us then assume they form CFs with each fermion being bound to 2 wave function vortices.The effective flux that the CFs experience is then N * ϕ = 2q.As the effective magnetic field is negligible in the thermodynamic limit, if the effective interaction between the CFs is weakly attractive we then expect them to form some weakly paired BCS state.Let u i , v i be the spinor coordinates for i th fermion.On the sphere, the flux attaching Jastrow factor is From standard CF theory we then write the "ideal" trial wave function as 2 , where, by imposing rotational invariance, the pairing function takes the general form g(r i − r j ) = lm (−1) m+|q| g l Y qlm (Ω i )Y ql(−m) (Ω j ) for some unspecified g l ∈ C with Y qlm (Ω) being the monopole harmonics [36].In principle, one can then extract the pairing physics of the CFs by finding the g l that minimizes the energy of this wave function.This is, however, numerically intractable for N ⪆ 10 due to the projection to the LLL.To create a numerically tractable trial wave function MS proposed using the Jain-Kamila [37] procedure where we can produce a LLL wave function simply by replacing the single-particle orbitals Y qlm (Ω i ) by the corresponding CF "orbitals" defined by, where ) and P LLL projects particle i to the LLL with magnetic flux 2q +N −1.The resulting family of pairing wave functions then defines MS paired CF wave functions [32], In particular, we denote the MS family of wave functions at effective flux q by MS 2q .One can then vary the g l to minimize the energy.The YM wave function at effective flux 2q, denoted YM 2q , is defined to be the MS 2q wave function with g l = 1 2l+1 .Note that as the interaction potential V (r) is real-valued, we must have that ⟨Ψ MS | V |Ψ MS ⟩ is invariant under time reversal where we simply replace the wave function by its complex conjugate.Furthermore, by expressing ⟨Ψ MS | V |Ψ MS ⟩ as an integral one can perform a change of variable where we change the azimuthal angle ϕ → −ϕ which is the same as the transformation u → u * and v → v * .Combining these two transformations we must have that Assuming that the minimum energy solution is unique up to multiplying all g l by the same complex number, it then follows that the minimum energy wave function can be expressed with g l all being real numbers.We will take g l ∈ R from now on.
From the usual CF theory, we expect that this can be physically interpreted as being analogous to a BCS state, of the CFs, of the form where c † qlm is the creation operator for the orbital Y qlm (Ω).The pairing wave function of this BCS state is g( ).This expansion of the pair wave functions was shown to describe pairing in the ℓ = 2q channel [38][39][40], as follows from the properties of the Dirac monopole harmonics on the sphere [41] by considering relations for the complex conjugation of these functions (see Eq. A5 in [40]).Alternatively, using Eq.B4 of Ref. [35] the pair wave function can be expressed as g l−q (cos θ), where θ is the angle between particle i and particle j and P (2q,0) n are the Jacobi polynomials, where it should be noted that if q < 0 then one should replace (u, v) → (u * , v * ).Here, the Jastrow factor (u i v j − v i u j ) 2q allows one to read off the pairing channel explicitly, as the remainder of the expression is real, confirming the result ℓ = 2q [42].Also from Eq. B6b of Ref. [35] we have , where again negative q requires complex conjugation.This pairing function scales with the distance between particle i and particle j as g ∼ 1/r.Thus the analogous weak pairing on the sphere is given by g l ∼ 1 2l+1 at small l, the form assumed by YM [35].Now let )) where R is the radius of the sphere, m CF the the CF effective mass and l F is the Fermi "shell".We can express g l = (ε l − ε 2 l + ∆ 2 l )/∆ * l where ∆ l is the gap function on the sphere.We can then make the analogous weak pairing ansatz for the gap function ∆ l = ∆ l+1/2 R .By matching the kinetic energy for l ≫ 1 one can see that the correspondence between l and the wave vector k is k = l R for large l.This then reproduces the previous ansatz for the gap function in terms of k for large l, ∆ l ≈ ∆k.By symmetry, the occupation of all orbitals in the l shell must be the same n l , which can then be expressed as For the weak pairing ansatz, this gives, Whilst this is perhaps an appealing way to physically interpret the MS wave function, it is not obvious if the CF orbitals are "normalised" in such a way that the g l of the MS wave function should be the same as the g l of the effective BCS state.As discussed by MS in Ref. [32], one can get around this issue by defining the effective CF occupations as This can also be expressed as, where the expectation value is taken with respect to |Ψ MS ⟩.From their definition, these clearly do not depend on the normalisation of g l .These effective occupation probabilities for the minimum energy g l were empirically found by MS to behave as the occupation probabilities of an actual fermion system.We will then use these n CF l in this work to relate optimized MS wave functions to their effective CF-BCS description.
As the pairing channel is given by ℓ = 2q, we expect that the MS wave functions MS −1 , MS 3 and MS 1 to be in the Pfaffian, anti-Pfaffian and PH-Pfaffian phases respectively.In what follows we will also use yet another variational wave function for the Pfaffian phase which we will denote MR* and can be expressed as, Roughly speaking, the MR* wave function allows us to "perturb" around the MR wave function, where if we take the g l to zero we recover exactly the MR wave function.By varying the g l parameters, one may expect that an MR* wave function can always be found with lower energy than the corresponding MR wave function (at the same system size).We will thus use the MR* wave function as a benchmark for the optimization of the standard MS −1 wave function.

III. MODEL INTERACTION
As is standard in the FQHE literature we model a system of electrons confined to the first excited Landau level (LL), and without LL mixing, by an effective description of electrons in the LLL.In the spherical geometry, where the actual system of interest has magnetic flux 2Q passing through the sphere's surface, the effective interaction V eff for the LLL system is defined by, where |Q, l, m 1 m 2 ⟩ is a two-particle state (with the particles not identical) with particle 1 and 2 in the Y Qlm1 and Y Qlm2 orbitals respectively and V is the interaction of the original problem.This definition of V eff is equivalent to requiring all the Haldane pseudo-potential coefficients V eff in the LLL to match the corresponding pseudo-potential coefficients of the actual interaction V in the second LL.Note that for the second LL system, we will take the radius of the sphere to be R = l B √ Q and for the LLL system we will take R * = l B √ Q + 1, with l B being the magnetic length.
In this work, we are interested in modelling the Coulomb interaction in the second LL.In particular, as the trial wave functions we are using (for the effective system) are all expressed in real-space (position) representation we then require a real-space form of the effective interaction.An ideal real-space form of the effective interaction would be such that all its pseudo-potential coefficients matched those of V eff (where it should be noted that there exist many ideal real-space potentials for a given magnetic flux).For a fixed magnetic flux, one can in principle find an effective real-space potential by expressing the effective interaction as a sum of 2Q + 4 real-space potentials whose vectors of pseudo-potential coefficients are linearly independent.In practice, this can make evaluating the real-space potential, which must be performed many times when using Monte Carlo methods, computationally expensive.Instead in this work, we use an approximate effective interaction, whose realspace form is simple to evaluate.
The model interaction we use to approximate the ideal effective interaction of the second LL Coulomb interaction is, where r is the chord length between the two particles on the sphere, and a i and α i are parameters that must be fit.This has been shown in previous works [8,33] to provide a good approximation to the desired effective interaction.In this work, we allow for the parameters a i , α i to vary with the system size, where the parameters are determined by minimising the sum of squared differences between the pseudo-potential of this interaction in the LLL and those of the Coulomb interaction in the second LL for L = 0, 1, 2, . . ., 2Q + 2 (with 2Q being the magnetic flux for the system in the second LL).
using Eq.2.37 of Ref. [43], where the pseudo-potential coefficients are expressed in terms of Wigner 3-j and 6-j symbols, and V k which are the coefficients in the expansion of V (r) in terms of Legendre polynomials, V (r) = k V k P k (cos θ) (with θ being the angle between the two particles on the sphere).The results of this fitting procedure for the 2Q = 35 case can be seen in Fig. 1, where it can be seen V eff can accurately reproduce the Haldane pseudo-potential coefficients for the Coulomb interaction in the 2nd LL with only slight deviations at intermediate m.The parameters used for V eff for the Coulomb interaction in the second LL for the system sizes used in this study are given in Appendix A.

A. Energetics and overlaps
Throughout this section, we will use the shorthand notation ED ℓ to denote the exact ground state of the approximate model of the Coulomb interaction in the second Landau level (see Sec. III) at L z = 0 and for total magnetic flux N ϕ = 2Q = 2N − 2 + ℓ.Note that at fixed N the ED ℓ states for ℓ = −1, 1, 3 occur at the same total magnetic fluxes as the MS ℓ for ℓ = −1, 1, 3 respectively.
Shows the Ẽ/N (see Eq. 12) for the optimized wave functions as well for the MR, the particle-hole conjugated wave functions of these, denoted with a bar, and the exact ground states ED ℓ , which are at Lz = 0 with total magnetic flux N ϕ = 2N − 2 + ℓ, for the approximate model (see Sec. III) of the Coulomb interaction in the 2nd LL.Lines show the resulting polynomial fits in 1/N (see main text for full details).Whilst error bars are included for the Monte Carlo estimated energies their sizes are similar to those of the markers.
interaction in the 2nd LL (see Sec. III), where a bar denotes the particle-hole conjugate of a wave function.The energies of the MS ℓ and MR* wave functions have been minimized using the optimization algorithm of Appendix B. The adjusted energy is defined by, In the case of the Coulomb interaction, the multiplicative factor adjusts the energies so that the particle density is kept constant, by rescaling the radius of the sphere N , and the Q term is the electrostatic energy of a uniformly charged sphere of radius l B √ Q with total charge N e.As discussed in Ref. [44] this is used to improve the estimates of the thermodynamic energy per particle, which we achieve by fitting a quadratic polynomial in 1/N to the Ẽ/N of the MR, YM ℓ , MS ℓ and MS ℓ , with ℓ = −1, 3, wave functions and a linear function of 1/N for the Ẽ/N of the YM 1 and MS 1 wave functions.The resulting estimated thermodynamic energies per particle can be found in Tab.I, where we have included the thermodynamic energy estimates for two different fitting procedures: Fit1 includes all data points to fit the given polynomial and Fit2 which includes all data points except the smallest system size to fit the given polynomial.This is to demonstrate that the thermodynamic energy estimates are somewhat robust.However, it should be noted that the quoted errors in the thermodynamic estimates, given from the output of the curve fitting program, are evidently larger which can be seen from the differences in Fit1 and Fit2.The curves in Fig. 2 have been fit using the Fit1 procedure.
For those wave functions whose real-space form is known exactly, the energies at each system size have been estimated using ∼ 5 × 10 9 Monte Carlo samples.The energies of the particle-hole conjugated wave functions have been calculated using the result Ref. [45], where for a rotationally symmetric interaction under a particle-hole transform the energy transforms as where E filled is the energy of the corresponding filled Landau level.
As can be immediately seen from Fig. 2 despite allowing for some optimization, the MS 1 wave functions, which are expected to be in the PH-Pfaffian phase, are energetically unfavourable in comparison with the trial wave functions at the other pairing channels tested in this work in the case of the (approximate) Coulomb interaction in the 2nd LL.In fact, the amount by which  II. Showing overlaps of the YM1 and MS1 wave functions, expected to be in the PH-Pfaffian phase, with the exact ground state, ED1, at the corresponding shift and at Lz = 0, in the case of the approximate model interaction for the Coulomb interaction in the first excited Landau level (see Sec. III).
one can reduce the energy of these ℓ = 1 wave functions is negligible in comparison with the energy scales in Fig. 2.
One can also see that the energies of the YM 1 and MS 1 wave functions are far higher than the energies of corresponding ED 1 's.One can also see from Tab. II that the overlaps of the YM 1 and MS 1 wave functions with the exact ground state for the chosen model interaction, ED 1 , at the corresponding shift, are nearly zero.Thus, the YM 1 and energy minimized MS 1 wave functions do not offer a good approximation to the exact ground state.III.Showing overlaps of the MR, YM−1 and the energy minimized MS−1 wave functions, expected to be in the Pfaffian phase, with the exact ground state, ED−1, at the corresponding shift and at Lz = 0, in the case of the approximate model interaction for the Coulomb interaction in the first excited Landau level (see Sec. III).
Another observation that one can make from Fig. 2 is that the energy minimized MS −1 wave functions offer substantial energy reduction over the MR and YM −1 wave functions at small system sizes, with energies close to those of the corresponding exact ground states ED −1 , however still with a notable error.This is further detailed in Tab.III, where it can be seen that the energy minimized MS −1 wave functions have significantly higher overlaps with the corresponding ED −1 states compared with those of the YM −1 and MR wave functions, at small system sizes.
On the other hand, in the ℓ = 3 sector it can also be seen that the energy minimized MS 3 offer noticeably less  IV.Showing overlaps of the YM3 and the energy minimized MS3 wave functions, expected to be in the anti-Pfaffian phase, with the exact ground state, ED3, at the corresponding shift and at Lz = 0, in the case of the approximate model interaction for the Coulomb interaction in the first excited Landau level (see Sec. III).
energy reduction relative to the other trial states relative to the energy gain of the MS −1 wave functions in the l = 1 sector.The particle-hole conjugated MS −1 wave functions, MS −1 , in fact, offer a better approximation to the energies of the corresponding exact ground states ED 3 .As can be seen from Tab. IV the MS 3 wave functions have noticeably better overlaps with ED 3 compared with YM 3 , although these improvements in the exact ground state overlaps are still not as significant as those of the MS −1 wave functions.
Furthermore, as can be seen from Fig. 2 and Tab.I the energies of the MR, MR, YM ℓ , with ℓ = −1, 3, and the optimized MS ℓ , along with their particle-hole conjugates MS ℓ , wave functions all converge in the thermodynamic limit.This is perhaps surprising as one would typically expect that allowing for some energy minimization away from the zero-parameter trial wave functions would allow for an improved estimate of the energy of the actual ground state in the thermodynamic limit.Even the MR* wave functions do not offer a better thermodynamic energy estimate, even though we expect these to always have lower energy than the corresponding MR wave function, which can be seen in Fig. 2 where they have energies that are not discernable from those of the optimized MS −1 wave functions.
There are several possibilities at this point.Firstly, it may, in fact, be the case that the amount by which the energy can be reduced by optimizing the MS −1 and MS 3 wave functions compared with the MR and YM 3 wave functions falls to zero in the thermodynamic limit.In short, these wave functions may just not offer enough variational freedom at larger system sizes.One may expect, however, that the MS −1 wave functions should offer better thermodynamic estimates given that at smaller system sizes their energies are comparable with more exact methods.On the contrary, it should be noted that the polynomial used to extrapolate MS −1 has a noticeable curvature, which implies even slight differences between the exact energies and the MS −1 can result in a large difference between their thermodynamic extrapolations [46].Of course, one could object to the extrapolation method used here.Whilst it can never be definitively known if the extrapolation method is correct, we can at least verify it is a "reasonable" method by the fact that the thermodynamic energy extrapolations of each wave function converge precisely with its corresponding particle-hole conjugate, when the conjugate wave function has been included.
Another possibility is that the optimization algorithm is getting stuck at a local minimum.This could be alleviated by starting the optimization at randomized g l , however, this would come at an increased computational cost as the algorithm must explore a larger area of parameter space to find a minimum.Indeed, the computational cost of optimizing these wave functions using the procedure outlined in Appendix B should be emphasised.For example, at N = 26 particles the optimization of a wave function with 8 g l using 300 computing nodes can take around a week, at current computer standards, from the beginning of the optimization to obtaining an accurate estimate of the energy of the optimized wave function.In the worst case, if several fine-tuning phases are required, obtaining the optimum energy at the desired level of accuracy can take on the order of a month.As can be seen from Fig. 2 this is partly due to the fact that many Monte Carlo samples are required at each iteration of the optimization in order to resolve the rather small differences in the energies of the various wave functions.Thus, although the algorithm could be started from randomized g l this would come at an increased cost which would render this procedure impractical for most interesting use cases.Whilst, we have checked in Fig. 5 if using extra g l in the optimizations performed makes very little difference in the energies, it is possible that adding these extra g l could make a difference to the thermodynamic extrapolation, particularly for the MS −1 where the extrapolating polynomial has noticeable curvature.However, with each new g l the computational cost increases as many more CF orbitals need to be computed.In short, this is perhaps not a practical method for obtaining better thermodynamic estimates in comparison with other numerical methods.
Finally, we would also like to emphasise that at small system sizes the YM −1 and YM 3 are very good approximations to the MR and MR wave functions respectively.From Tab.V one can see that the YM 3 wave function has an overlap with the MR wave function that is at least 98% for system sizes of 10-14 particles and from Tab. VI that the YM −1 wave function has an overlap with the MR wave function of around 98% for system sizes of 12-16 particles.Such overlaps had already been reported by YM between the YM 3 and MR wave function, at smaller system sizes, and between the YM −1 and MR wave functions at larger system sizes [35].The fact that these YM wave functions are good approximations to the MR wave function and it's particle-hole conjugate can also be observed from Fig. 2 where the energies of the YM 3 and YM −1 wave functions are very close to those of the MR and MR wave functions respectively.Whilst the YM −1 wave function is clearly not as numerically tractable as the MR wave function, the YM 3 does provide a more numerically tractable approximation to the MR wave func-tion, at least at small and intermediate system sizes.

B. Effective CF pairing
Whilst the optimized MS ℓ wave functions do not appear to offer any better thermodynamic energy estimates compared with the YM ℓ , they do offer more insight at finite-size systems through the effective paired CF description.We have estimated, using Monte Carlo, the effective CF occupation probabilities n CF l (Eq.8) for the optimized MS ℓ wave functions and for the corresponding YM ℓ .We have then fitted the BCS type weak pairing ansatz (WPA) of Eq. 7 for the various estimated n CF l , where we allow for the Fermi-level, l F , to be a continuous parameter that can be varied in the fit and we take the radius of the sphere used in Eq. 7 to be R = l B √ N so as to keep the particle density constant.This then allows us to give an estimate for the dimensionless parameter , where m CF is the effective CF mass which we will assume to be constant so that we can take mCFl B ∆ ℏ 2 as a measure of the effective CF pairing strength.
Fig. 3 shows the estimated n CF l and corresponding fitted WPA for the optimized MS ℓ and fixed parameter YM ℓ wave functions at ℓ = −1, 3, 1 for N = 20, 18, 20 respectively.The n CF l of the MS ℓ and YM ℓ for ℓ = −1, 3 can be fit reasonably well to the WPA.
In Fig. 4 shows the fitted m CF l B ∆/ℏ 2 parameter of the optimized MS ℓ and YM ℓ wave functions for ℓ = −1, 3 as a function of 1/N , where it can generally be seen that the effective CF pairing in the optimized MS ℓ wave functions is higher than that of the YM ℓ .Interestingly we find the effective pairing strength for the optimized MS ℓ for ℓ = −1, 3 wave functions to be roughly the same at the same number of particles N , where it should be emphasised that at the same N the MS −1 and MS 3 wave functions occur at different total magnetic flux through the sphere (i.e. have different shifts S).
The data from Fig. 4 can also offer more insight into some of the observations one can make from Fig. 2. One can see that from Fig. 2 that the amount by which the energy of the MS −1 wave functions can be lowered compared with the corresponding MR wave functions is generally larger than the amount by which the energy of the MS 3 wave function can be lowered in comparison with the YM 3 wave functions.Taking the YM −1 wave function as an approximation of the MR wave function, this has a simple interpretation in the effective pairing description in that the effective pairing strength of the YM −1 wave function is generally lower than that of the YM 3 wave functions.Thus the YM 3 wave functions have an effective pairing closer to the optimum compared with the YM −1 wave functions which gives some explanation as to why the energy of the YM 3 wave functions are already close to optimal MS 3 wave functions, whereas at intermediate system sizes the optimum MS −1 wave functions have energy noticeably lower than the corresponding MR wave function.(see Eq. 8) for various YM ℓ and optimized MS ℓ wave functions (dashed lines) along with a corresponding fit to the weak-pairing ansatz of Eq. 7 (solid lines), where lF is taken to be an adjustable continuous parameter.The fitted parameters are indicated in the bottom left corner of the corresponding plot, where mCF is the effective CF mass.Note that mCFlB∆/ℏ 2 parameter has been estimated where we take the radius of the sphere to be R = lB √ N , which keeps the particle density constant.
Finally, as can also be seen in Fig. 3 the n CF l of the optimized MS 1 and YM 1 wave functions are noticeably larger than one in some cases and are negative in some other cases, which is clearly inconsistent with interpreting these n CF l as occupation probabilites.This was found to occur at all other tested system sizes.This interpretation of the n CF l is based on the assumption that the Jastrow factor of the wave function can be approximated in some "mean-field" way, which is usually assumed when interpreting generic CF wave functions.Combining this with the observations of YM, that the YM 1 wave function shows no sign of emergent pairing in the pair correlation function, and with the evidence that the "ideal" PH-Pfaffian wave function may, in fact, represent a gapless phase of matter [31], indicates the possibility that this usual "mean-field" interpretation of the Jastrow factor might break down for the current candidate PH-Pfaffian wave functions and so may in fact not be a representative for the phase of matter predicted by Son [6].It has also recently been argued by Haldane [47], based on a conjecture that FQH states must have a non-zero so-called "guiding centre quadrupole moment", that particle-hole symmetric states can never be FQH states, which may give an explanation for these observed pathologies.Despite this, we can still see that the n CF l of the optimized MS 1 wave functions can still be roughly fit to the WPA with ∆ = 0, which would correspond to the gapless CF Fermi-liquid.In summary, even after some optimization, we see no evidence of an effective paired CF description for the PH-Pfaffian wave functions.

V. CONCLUSION
In this work, we have shown how the energy of the paired CF wave functions, MS ℓ , proposed by MS [32] and extended by YM [35], can be minimized in a practical manner by varying the pairing function for the pairing for the optimized MS ℓ and the fixed parameter YM ℓ wave functions for ℓ = −1, 3 over a variety of system sizes, where mCF is the unknown effective CF mass.Note that mCFlB∆/ℏ 2 parameter has been estimated where we take the radius of the sphere to be R = lB √ N , which keeps the particle density constant.
channels ℓ = −1, 3, 1 which are believed to represent the Pfaffian, anti-Pfaffian and PH-Pfaffian topological orders respectively.We have presented the result of such optimizations in the case of an approximate model of the Coulomb interaction in the second Landau level (LL).
For pairing channel ℓ = −1 we found that the energy can be reduced substantially below that of the MR wave function at intermediate system sizes, with a noticeable improvement in the overlaps with the corresponding exact ground state.For the ℓ = 3 pairing channel, however, we find the energy cannot be reduced as much, although the resulting improvement in the overlaps with the corresponding exact ground state are still notable.We find that optimizing the pairing channel ℓ = 1 wave functions makes no qualitative difference and these PH-Pfaffian wave functions are still very energetically unfavourable compared with the ℓ = −1, 3 pairing channels.We have further emphasised that the fixed parameter versions of these wave functions YM ℓ , proposed by YM [35], at pairing channels ℓ = −1 and ℓ = 3 have very high overlaps (98%-99%) with the MR and MR wave functions respectively, at system sizes of around 10-16 particles.
The effective CF pairing of these MS ℓ and YM ℓ wave functions was then investigated by studying the effective CF orbital occupation probabilities.For the pairing channels ℓ = −1, 3 we found both the energy minimized MS ℓ and YM ℓ wave functions show effective CF pairing which can be well approximated by some weak paring BCS type ansatz where the pairing parameters could be extracted, with the MS ℓ wave functions showing a stronger pairing strength than the YM ℓ wave functions.For the YM 1 and energy minimized MS 1 wave functions showed no sign of effective CF pairing, which adds to the growing list of pathologies found for current trial wave functions for the PH-Pfaffian phase.
It is not clear from this work if the observed convergence of the thermodynamic energies of both the energy minimized and fixed parameter wave functions in the Pfaffian and anti-Pfaffian phases is unique to the Coulomb interaction in the second LL.It would be useful to perform these optimizations at slightly different interactions to see if this convergence is generic.Although one should be cautious of the observed charge density wave phase in the vicinity of the Coulomb interaction in the second LL [13].It would also be interesting to see if minimizing the energy of the MS 1 wave function for some interaction near the Coulomb interaction will result in an ℓ = 1 wave function which shows a definite signature of emergent CF pairing.given set of parameters is computed exactly.This is, however, computationally expensive and would restrict one to only working with smaller system sizes.For the energy minimization of wave functions of large systems, there exists a large class of algorithms under the name variational quantum Monte Carlo [48], where Monte Carlo methods are used to estimate the energy and energy gradients.Whilst there are many specific algorithms to choose from, for the problem at hand we have used a combination of the Stochastic Reconfiguration (SR) algorithm [49] and the Adam optimizer [50], which we will now describe.
First, let us describe the Adam optimizer.Let E(g) be the expectation value of the Energy as a function of the wave function parameters g.We will then impose a yet unspecified geometry on this parameter space defined by some metric tensor S which can vary over the parameter space.At each iteration t of the optimization, we have three vectors, f t ≡ ∂E ∂gt (gradient vector), m t (momentum) and mt (bias-corrected momentum).We also have the real numbers v t and vt the bias-corrected version.The algorithm has four hyper parameters: γ (learning rate), β 1 , β 2 and ϵ.
Using the momentum m t allows the algorithm to smooth out the noise in estimating gradients and v t ensures we move roughly the same distance in parameter space at every iteration.One can view this algorithm as a particle moving in parameter space with friction and noise in the potential given by E. β 1 can then be thought of as setting our inertia.β 2 sets the time scale over which we average the sizes of S −1 g t .The ϵ parameter is simply a distance cut-off.
At each iteration, everything is updated by, Let O l be an operator which is diagonal in the position space representation, which is defined by O l ≡ ∂g l ΨMS ΨMS .The energy gradients can then be expressed as, where ℜ denotes the real part.We compute these expectation values using the Hastings-Metropolis Monte Carlo algorithm, which introduces stochastic noise into the optimization.
The SR algorithm can be interpreted as a standard gradient descent algorithm where one uses the Fubini-Study metric to define the geometry of the parameter space of the wave function.To combine this with the Adam optimizer we simply take S to be the Fubini-Study metric.This metric is given by defining the distance between two nearby points in parameter space to be the Hilbert space norm of the difference between the normalised states at the two coordinate points.This can be straightforwardly shown to be given by, Using this metric has the advantage that it removes the ambiguity over how the parameters are normalised.At each iteration of the algorithm, we again use the Hastings-Metropolis Monte Carlo algorithm to estimate S.
Both the gradients f and S are estimated using ∼ 10 6 Monte Carlo samples, with the actual number of samples used increasing with the system size.To lessen the computation time sampling at each iteration was run in parallel using around 10 2 computing nodes.
In practice, we regularise this metric as the Ψ MS is invariant under multiplying all the g l by a constant which will mean S will always have determinant zero.We regularise by, S l1l2 → (1 + εδ l1l2 )S l1l2 , for some small ε which is kept constant throughout the optimization.
Throughout the optimizations performed in this study we use the recommended hyperparameters from Ref. [50] of β 1 = 0.9, β 2 = 0.999 and ϵ = 10 −5 .We found that a suitable learning rate for this optimization problem is around γ ∼ 0.005 and we use the same ε to regulate S as given in Ref. [51] with ε = 10 −3 .
The issue of how to pick the number of g l to use is addressed as follows.For the MS −1 case, we ran the optimization at N = 12 with the first 7 g l , where it was found that the broadening of the Fermi surface, as seen through the n CF l , was around 2 in l space (i.e.n CF l were found to be sufficiently close to zero for l > l f + 2).We expect this broadening in l space to scale as √ N as l ∼ kR.Thus, for system sizes from 12-20 particles, we used the first 7 g l to optimize the MS −1 wave functions and at 26 particles we used the first 8.To demonstrate this number of parameters is sufficient we estimated the energy as a function of the 9 th g l around the optimum set when using the first 8 g l for the N = 26 MS −1 optimization with the approximate model of the Coulomb interaction in the second LL.The result of this parameter  scan can be seen in the plot of Fig. 5, where it can be seen that the possible energy reduction per particle that can be achieved by varying this parameter is small in comparison to the difference between the previous optimum energy per particle and the thermodynamic energy per particle of the Moore-Read wave function (see Fig. 2 and Tab.I for details of thermodynamic estimates where we use the Fit1 estimate for the thermodynamic energy per particle of the MR wave function).Similarly, for both the MR* and the MS 1 wave functions we used the first 7 g l for 12-20 particles and used the first 8 for 26 particles.By the same method, we found that using the first 6 and 7 g l was sufficient to optimize the MS 3 wave functions for 10-18 and 24 particles respectively.
When optimizing a given MS ℓ wave function we start the algorithm with g l = 1 2l+1 for the parameters that are actually varied and all other g l are set to zero, g l = 0.This is the approximate YM ℓ wave function which we expect to be close to the minimum energy.
The plot of Fig. 6 shows the n CF occupation probability at each iteration t for the optimization of the N = 20 MS −1 wave function with the approximate model of the 2nd LL Coulomb interaction (note that l = 9  2 is just at the Fermi level).Initially, the algorithm moves sharply towards a minimum where it plateaus near the minimum energy at around t = 4000.We found that this algorithm converged after around t ∼ 4000 iterations for all cases considered using the hyperparameters given above.One can also see that the apparent noise in the n CF (t) path increases as we approach the minimum.At larger system sizes this sometimes required a fine-tuning stage with around 100 iterations where the number of samples and the learning rate are increased to obtain the minimum energy solution with the desired accuracy.
FIG.1.Showing the Haldane pseudo-potential coefficients, Vm, of both the Coulomb interaction in the second LL along with the Vm of the fitted LLL V eff , where the total number of magnetic flux quanta is 2Q = 35 for the system in the second LL.Note that VL = Vm=2Q+2−L.

FIG. 3 .
FIG. 3. Showing the CF orbital occupation probabilities n CF

FIG. 5 .
FIG.5.Adjusted energy Ẽ (see Eq. 12) as a function of the 9th g l , around the optimum solution for the first 8 g l for the N = 26 MS−1 wave function with the approximate model of the 2nd LL Coulomb interaction, where ĒMR is the thermodynamic Ẽ/N of the MR wave function (see Tab.I where we use the Fit1 value).

2 N = 20 MS − 1 FIG. 6 . 2 (
FIG. 6. Showing n CF 9 2 (see Eq. 8) at each iteration t of the optimization of the N = 20 MS−1 wave function with the approximate model of the 2nd LL Coulomb interaction.
4IG.4.Showing the fitted mCFlB∆/ℏ 2 of the weak pairing ansatz, Eq. 7, for the effective CF occupations n CF l

TABLE VII .
Fitted V eff parameters for the 2nd LL Coulomb interaction at all 2Q (of the second LL system) considered in this study.