Study of nuclear pairing with Configuration-Space Monte-Carlo approach
aa r X i v : . [ nu c l - t h ] M a r Study of nuclear pairing with Configuration-Space Monte-Carlo approach
Mark Lingle and Alexander Volya
Department of Physics, Florida State University, Tallahassee, FL 32306, USA (Dated: July 31, 2018)Pairing correlations in nuclei play a decisive role in determining nuclear drip-lines, binding ener-gies, and many collective properties. In this work a new Configuration-Space Monte-Carlo (CSMC)method for treating nuclear pairing correlations is developed, implemented, and demonstrated. InCSMC the Hamiltonian matrix is stochastically generated in Krylov subspace, resulting in theMonte-Carlo version of Lanczos-like diagonalization. The advantages of this approach over othertechniques are discussed; the absence of the fermionic sign problem, probabilistic interpretation ofquantum-mechanical amplitudes, and ability to handle truly large-scale problems with defined pre-cision and error control, are noteworthy merits of CSMC. The features of our CSMC approach areshown using models and realistic examples. Special attention is given to difficult limits: situationswith non-constant pairing strengths, cases with nearly degenerate excited states, limits when pairingcorrelations in finite systems are weak, and problems when the relevant configuration space is large.
PACS numbers: 21.60.Ka, 02.70.Ss, 21.60.Cs
I. INTRODUCTION
Pairing correlations are a salient component of the nu-clear many-body dynamics which has a profound impacton most of the nuclear properties, and on the nuclearlandscape in general. The recently published volume
Fifty Years of Nuclear BCS [1] offers a unique overview ofmore than fifty years of research in this area. Advancesin experimental techniques and emergence of many newfacilities made it possible to explore nuclear systems atthe edge of stability. Interest in pairing is reinvigoratedby extraordinary effects observed recently in near-dripline nuclei. This includes the so-called nuclear halo effectin neutron rich nuclei, such as the case of borromean nu-cleus Li which is bound only due to the specifics of thepairing dynamics of the two neutrons above the Li core.Recent observations of di-neutron decay [2] and other re-markable manifestations of pairing, seen in structure andreactions with exotic nuclei [3, 4], all encourage theoret-ical effort to be continued.The pairing interaction involves pairs of time-conjugate single-particle states. We use p † k , p k , and ˆ n k to denote the corresponding pair creation, pair annihila-tion, and number-of-pairs operators; index k is used toidentify various distinct pair-states in the system. Thepairing Hamiltonian of interest is defined as: H = 2 X k ǫ k ˆ n k − X k,k ′ G kk ′ p † k p k ′ , (1)Here ǫ k are the single-particle energies and G kk ′ are thematrix elements of pairing interaction. In this work welimit our discussion to fully paired systems of n pairs in ω pair-states, which corresponds to 2 n fermions within thetotal particle capacity 2 ω of the valence space. Workingunder the assumption of a fully paired state is completelygeneral: any unpaired nucleons remain untouched by theHamiltonian (1); these nucleons effectively block somepart of the valence space so that the problem is thenreduced to a fully paired state in a reduced space. Thus, the configuration space of interest spans over all ω choose n basis states | n i = | n , n , . . . n ω i . (2)Here we use occupation representation where for eachpair-state n k = h n | ˆ n k | n i = 1 or 0 depending on whetherthe pair-state is occupied or not. Clearly, the totalnumber of pairs n = P k n k . Any state can be rep-resented as a linear combination of the basis states, | Φ i = P n h n | Φ i | n i . Pairing correlations have been traditionally exploredwith the help of the BCS theory of superconductivity[5]. This variational technique, which is formally exact inthermodynamic limit, is very well integrated into moregeneral mean-field approaches and into techniques be-yond mean-field. Starting from pioneering works [6, 7],the BCS theory has been applied in nuclear physics withgreat success. However, non-conservation of the parti-cle number and difficulty in handling limits where pair-ing is weak as compared to the characteristic mean-fieldsingle-particle level spacing, have proven to be significantdrawbacks in applications of BCS to finite nuclear sys-tems [1, 8–11]. Over the years a number of remedies havebeen proposed to overcome these drawbacks. For exam-ple, the issue of the particle number non-conservation hasbeen addressed with a variety of techniques proposed inRefs. [12–19].With theoretical and computational advances a grow-ing number of pairing problems in mesoscopic systems,such as atomic nuclei, can be treated exactly; thus avoid-ing the BCS and its drawbacks. There are several ma-jor groups of exact methods. Symmetry-based algebraicmethods were introduced by Racah [20–22] even beforethe BCS theory. These methods found wide applicabilityboth independently [23–27] and as components of othertechniques [28, 29].Presented more than 40 years ago by Richardson [30–34], an exact solution that reduces the pairing eivenvalueproblem to a set of non-linear equations have been suc-cessfully generalized and applied [35–38] in multiple sit-uations. Some generalizations and interpretations, suchas those related to electrostatic analogies [39], are of par-ticular theoretical interest [36].Computational advances and iterative sparse matrixdiagonalization algorithms allowed for direct diagonal-ization methods to emerge as extremely simple, stable,and robust alternatives [40–42]. Nevertheless, the dimen-sion of the Hamiltonian matrix in the relevant basis spacegrows exponentially with the number of pairs, eventuallyrendering these methods computationally impractical es-pecially for model spaces required for problems with con-tinuum of scattered states. Mote Carlo approach, whichis the main subject of this work, so far appears to pro-vide the only reasonable technique that overcomes, in acontrolled way, the exponentially growing computationaldifficulty. There exist numerous variations of the MonteCarlo approach, varying in philosophy and implementa-tion. Many of these methods can be found in the text-book [43]. The Shell Model Monte-Carlo [44, 45] andQuantum Monte-Carlo involving variational, auxiliary-field and Green’s function versions (see review [46]) areamong well-known successful examples used in low en-ergy nuclear physics. Random Monte-Carlo sampling,either for variational purposes or in order to evaluatemulti-dimensional integrals, such as those emerging inHubbard-Stratanovich transformation is at the center ofthese techniques.The pairing Hamiltonian (1) is equivalent to a Hamil-tonian describing ω spin-1 / n k = 1 and 0 for up and down spin orientations,respectively, [10, 41, 42]. This offers opportunities for abroad class of Monte-Carlo methods known in spin sys-tems [43] to be applied. The idea of using the connectionbetween spin physics in condensed matter and quasispinin pairing problems was originally explored by Cerf andMartin [47, 48]. In their approach the ground state | Ψ i is found as an asymptotic state that emerges from an ar-bitrary initial state Φ as a result of evolution along theimaginary time: | Φ i ≃ e − τE h Ψ | Φ i| Ψ i , as τ → ∞ . (3)In order to propagate the initial state along the imaginarytime, Cerf and Martin proposed breaking the Hamilto-nian into the two non-commuting parts H and H , cor-responding to one-body and two-body terms. Then thepropagation can be done in small steps ∆ τ using Trotter-Suzuki operator decomposition [49, 50], e − ∆ τ ( H + H ) = e − ∆ τ H e − ∆ τH e − ∆ τ H + O (cid:0) ∆ τ (cid:1) . (4)The principal advantage of the technique is that H isdiagonal in the basis states | n i and the correspondingexponent can be easily evaluated. Assuming a constantpairing, where G kk ′ ≡ G, Cerf and Martin proposed toevaluate e − ∆ τH stochastically by breaking the exponentinto a Taylor series and taking advantage of the fact thatfor constant pairing strength the probability of a walk in configuration space to have a given number of steps isexactly Poissonian. Given that for constant pairing H can be diagonalized analytically using quasispin algebra,it may be possible to use the Trotter-Suzuki propagationwithout involving Monte-Carlo.With some degree of success, the method of Cerf andMartin was picked up recently by other research groups[51, 52]. The algorithm has so far been applied mainlyto problems with constant pairing matrix elements. Thisapparent limitation is not well addressed in the literature,but appears to be related to unknown quality and relia-bility of the stochastic evaluation of exponential operatorof the two-body interaction with non-constant matrix el-ements.In the Configuration-Space Monte-Carlo (CSMC) al-gorithm presented in this work we completely avoid theimaginary time evolution and Trotter-Suzuki decompo-sition; instead, using a stochastic process of nucleon-pair diffusion through the configuration space, we built aKrylov subspace which contains the set of lowest eigen-values. The resulting algorithm can be seen as a Monte-Carlo version of the well-known Lanczos algorithm. Thisclass of algorithms are often referred to as projector al-gorithms [43] since repeated application of the Hamilto-nian operator to a random state eventually amounts tothe ground state being projected out. Excited states canbe obtained as well by storing the wave functions and byenforcing orthogonality. II. CONFIGURATION SPACE MONTE-CARLO
In this section we present the Configuration SpaceMonte-Carlo method. It should be mentioned that themethod is generic, and is not limited to pairing, but thespecifics of the pairing Hamiltonian offer some big advan-tages, which is discussed in Sec. II C and demonstratedin Sec. III.
A. CSMC formalism
Let us consider a sequence of states | Φ L i ≡ V L | Φ i (5)which is generated by a repeated application of theHamiltonian H = − V onto a random initial vector | Φ i .These states span over the Krylov subspace. Eigenval-ues of the Hamiltonian matrix in this subspace converge,after enough iterations, to the eigenvalues (greatest inabsolute value) of the Hamiltonian in the entire space.Since we are interested in the lowest, most negative,states it is convenient to carry out this discussion us-ing V = − H. The repeated application of the opera-tor V can be written as a summation over all possible L + 1 intermediate states which are given by the sets { n } L ≡ { n , n , . . . n L } , | Φ L i = X { n } L | n L i A ( { n } L ) , (6)where the amplitude is A ( { n } L ) ≡ h n L | V | n L − ih n L − | V | n L − i . . . h n | V | n ih n | Φ i . (7)One advantage of evaluating powers of the Hamiltonianoperator is that the summation in Eq. (6) is restrictedto all possible paths n → n · · · → n L where each con-secutive configuration is connected to the previous oneby the matrix element of the interaction V. Therefore, inwhat follows { n } L denotes a connected L -step long path n → n · · · → n L . The path summation can be performed using Monte-Carlo sampling. To be more specific, if one generates N paths { n } ( s ) L ≡ n ( s )0 → n ( s )1 · · · → n ( s ) L , labeled here withsuperscript s = 1 . . . N, then | Φ L i ≈ N N X s =1 | n ( s ) L i B ( { n } ( s ) L ) (8)where B ( { n } ( s ) L ) ≡ A ( { n } ( s ) L ) P ( { n } ( s ) L ) (9)is the amplitude for the s -th random path weightedby the inverse of the probability to generate this path P ( { n } ( s ) L ). In applications where sampling is done withuniform probability P ( { n } ( s ) L ) ≡ P , the common term1 / P is just the total number of all possible paths { n } L which equals to the number of terms in the sum inEq. (6).Each sampling path can be generated as a randomwalk. The amplitude in Eq. (7) is subject to the recursionrelation A ( { n } L +1 ) = h n L +1 | V | n L i A ( { n } L ) , (10)where A ( { n } ) = h n | Φ i . Similarly, the probabilityfor a path is the product of probabilities for each step.Therefore, starting from the probability to pick the firstconfiguration P ( n ) ≡ P ( { n } ) , the probability for theentire path is generated recursively as P ( { n } L +1 ) = P ( n L → n L +1 ) P ( { n } L ) ; (11)here P ( n → n L +1 ) is the conditional probability tomove to configuration n L +1 given the current positionat n L . Thus, while going along a random path the coef-ficients B are generated recursively, B ( { n } ) = h n | Φ iP ( n ) and (12) B ( { n } L ) = h n L | V | n L − iP ( n L − → n L ) B ( { n } L − ) . (13)The probability distribution for selecting an initial posi-tion in configuration space P ( n ) and the distribution ofconditional probabilities P ( n L − → n L ) describing thedirection in which each next random step is to be taken,are both arbitrary user-supplied functions. Strategies forselecting these functions are discussed in what follows.It is important that the probability of taking a certainstep depends only on the current position and not onthe preceding history, therefore the process represents aMarkov chain [43]. The computational implementation ofthe Markov Chain Monte-Carlo methods is a well studiedsubject; see Ref. [43] and references therein.The Configuration Space Monte-Carlo approach, de-fined by Eq. (8), is implemented using an ensemble of N “walkers” starting from configurations n ; the initial con-figurations are generated with the probability distribu-tion P ( n ) . Then each walker independently takes L ran-dom steps; the probability distribution P ( n L → n L +1 ) isused to generate steps. We envision that each walker car-ries a “bag” B that is initialized and modified along thepath following Eqs. (12) and (13). Contributions fromthe bags of all walkers arriving to a given configuration n L comprise the component h n L | Φ L i as shown by Eq.(8).As the most straightforward application of the method,one could assume the probabilities for steps in all “di-rections” to be equal, then the conditional probabil-ity P ( n → n ′ ) depends only on the initial configuration n and the inverse of it equals to the number of con-figurations connected to n . In most cases the numberof connected configurations is the same for all states,which makes the conditional probability for each stepbeing an absolute constant, i.e., independent of initialand final positions. For example, for any paired con-figuration with n pairs and ω pair-spaces the pairingHamiltonian can generally move one of the n pairs ontoone of the ω − n + 1 unoccupied pair-states (this in-cludes diagonal move back to the same pair-state). Thus,in the pairing case, for equiprobable steps the condi-tional probability becomes a configuration-independentconstant P ( n → n ′ ) = ( n ( ω − n + 1)) − and the result-ing random paths are all generated with equal probabil-ity. This amounts to uniform Monte-Carlo sampling ofterms in sum (6). B. Importance sampling
Uniform sampling is convenient and effective when con-tributions from most paths are nearly equal; constant-strength pairing Hamiltonian discussed by Cerf and Mar-tin in Refs. [47, 48] is a good example of this situa-tion. However, sampling uniformly can be extremely in-effective if certain amplitudes A ( { n } L ) are very small orequal to zero; importance sampling can be introduced asa remedy. In the CSMC the contributions from differentsampling paths can be made comparable in magnitudeif steps are generated with probabilities proportional tothe magnitude of the corresponding matrix elements, P ( n L − → n L ) ∝ |h n L | V | n L − i| . (14)This way the scaling factor in Eq. (13) would not dependon the direction of the step. It should be emphasized,that satisfying the proportionality (14) exactly, whichmay be computationally expensive, is not necessary. Anyprobability distribution that in some general way followsthe distribution of the matrix elements is sufficient.The approach described here, referred to as Configura-tion Space Monte Carlo, allows one to build stochasticallythe Krylov subspace and find the eigenstates and eigen-values of the Hamiltonian using steps similar to those inLanczos approach. Clearly, the method is applicable toany Hamiltonian; however, different signs of matrix ele-ments h n L +1 | V | n L i and thus different signs of the am-plitudes can lead to poorly convergent sums. This issue,commonly known as the Monte-Carlo sign problem, is notpresent in applications of the CSMC to pairing problemsthat are discussed next. C. Features of the pairing Hamiltonian
Let us summarize some of the important features ofthe pairing problem that boost the effectiveness of theCSMC method. (i)
For fully paired systems the diagonal pairing matrixelements G kk are equivalent to the single-particle ener-gies. Thus, by redefining the diagonal pairing matrixelements as G kk → G kk − ǫ k / , the pairing Hamiltoniancan be written in the following form V ≡ − H, where V = X k,k ′ G kk ′ p † k p k ′ . (15) (ii) The pairing interaction is attractive. Therefore, withthe proper choices of phases all off-diagonal matrix ele-ments G kk ′ can be made non-negative. Without any lossof generality the single-particle energies can be measuredrelative to some chemical potential µ ; and the constant µ can be selected so that all diagonal many-body matrixelements h n | V | n i are also positive. (iii) The nucleon pairs are the only degrees of freedom,and the entire dynamics is represented by the “hopping”of the pairs between available states. Each pair hoppingleads to a step in configuration space where h n ′ | V | n i = h ...n k , ...n k ′ , ... | V | . . . n k , ...n k ′ , ... i = G kk ′ ≥ k = k ′ , and h n | V | n i = P k G kk n k > n ( ω − n + 1) different finalconfigurations that can be reached in one step. (iv) Given that the matrix elements of V are all positive,the Monte-Carlo sign problem does not appear. (v) Positive matrix elements of V imply that if in the ini-tial wave function Φ all components h n | Φ i ≥ , whichcan always be accomplished by defining phases of thebasis states | n i , then all components of any Φ L are non-negative, that is, h n | Φ L i ≥ L and any n . Thisalso allows one to introduce a linear L norm || Φ L || ≡ X n h n | Φ L i . (16) (vi) Asymptotically, as L → ∞ , | Φ L i ≃ ( − E ) L h Ψ | Φ i| Ψ i , (17)where | Ψ i is the ground state wave function and E isthe ground state energy. Since E < , in fact for thenegative-definite Hamiltonian, all eigenvalues are nega-tive, and the phase of the ground state wave functioncan be selected so that h Ψ | Φ i > , and all componentsof the ground state wave function are also non-negative: h n | Ψ i ≥ . (vii) Given that all many-body states that span theKrylov subspace have positive-definite amplitudes rel-ative to the basis states | n i , these amplitudes can betreated as probabilities. Therefore in the ideal limit ofthe importance sampling Monte-Carlo, when Eq. (14) issatisfied exactly, all walkers’ bags are equal. In this limitthe number of walkers arriving to a certain many-bodyconfiguration n on step L is proportional to h n | Φ L i . Lin-ear norm can be used to normalize the wave function.
III. NUCLEAR PAIRING WITHCONFIGURATION SPACE MONTE CARLO
In what follows we demonstrate the CSMC applica-tions to pairing problems. We organize our presentationby progressing from simple to more elaborate applica-tions, we use model and realistic examples to highlightthe CSMC and to address technical details.In the following subsections A-D we consider a systemconsisting of ω double-degenerate equally-spaced singleparticle orbitals, ǫ k = ǫ k with k = 0 , . . . ω − , wherethe single-particle level spacing ǫ defines the unit of en-ergy. This system, often referred to as a picket-fenceor ladder model, is commonly used for testing variousapproaches to pairing [42]. The model has a minimalsymmetry, time-reversal only, making it the most com-putationally challenging one. For the set of studies usingthe ladder model we assume constant pairing interaction,where G kk ′ = G in Eq. (1). This choice is not related toany limitation, this merely minimizes the number of pa-rameters and is convenient for comparison with numerousprevious studies where the same model was used. For thehalf-occupied ladder model the critical BCS strength isapproximately G cr = ǫ/ ln(2 ω ) , see Ref. [42].Realistic examples with non-constant pairing strengthand large scale application are discussed in subsectionsE and F. A. Linear norm
We start with a very simple and quick technique thatcan be used to determine the ground state energy andrequires no storage for wave functions. Despite certainlimitations, the method is elegant, simple in implementa-tion, very computationally efficient, and is an importantcomponent in the general approach.In the implementation of the CSMC through multi-ple walkers in configuration space the construction of thewave function | Φ L i is the most challenging task because,according to Eq. (8), it requires organizing walkers basedon their arrival locations. Even if performed in parallel,this is still a daunting task when the number of contribut-ing configurations becomes large. As we show next, forcertain observables, and for the ground state energy inparticular, this task can be avoided; thanks to the prop-erties of pairing interaction outlined in Sec. II C.According to Eq. (8) the average of all bags gives thelinear norm (referred to as L ) of the wave function1 N N X s =1 B ( { n } ( s ) L ) ≈ X n h n | Φ L i ≡ || Φ L || ; (18)computing the bag average is a simple and fast operation.Following Eq. (17), the bag averages for two consecu-tive values of L as L → ∞ give an estimate for the groundstate energy as E ≃ E ( L ) ≡ − P n h n | Φ L +1 i P n h n | Φ L i = − || Φ L +1 |||| Φ L || . (19)Clearly, any procedure on the Hamiltonian leading to theground state can be subjected to the linear norm. Forexample, using projection (3) one could evaluate energyin the limit τ → ∞ as E ≃ E ( τ ) ≡ || He − τH Φ |||| e − τH Φ || , (20)where exponents are evaluated in CSMC approach usingTaylor series || e − τH Φ || = ∞ X L =0 τ L L ! || Φ L || . (21)Obviously, it is possible to compute the linear norm forany operator || O Φ || ≡ P n h n | O | Φ i , but unfortunately inmost situations this linear norm does not have a trans-parent physical meaning.This quick linear-norm-based technique for evaluatingthe ground state energy within CSMC algorithm is il-lustrated in Fig. 1, and is compared to the general ap-proach discussed in the following subsection. A ladder model with ω = 18 levels and n = 9 nucleon pairs, where G = 1 , is used in this example. Two different startingwave functions | Φ i are considered. In the Fermi state alllowest single-particle levels are occupied up to the Fermisurface, | Φ (Fermi)0 i = n Y k =1 p † k | i = | , . . . | {z } n spaces , , . . . i . (22)The Fermi state is an exact ground state of non-interacting fermions ( G = 0 limit).The BCS solution offers a second convenient startingwave function Q ωk =1 ( u k + v k p † k ) | i . In our applicationsit is projected onto an appropriate number of particles;therefore | Φ (BCS)0 i is defined via amplitudes as h n | Φ (BCS)0 i = ω Y k =1 ( u k δ n k , + v k δ n k , ) . (23)The coefficients u k and v k are determined by solving theusual BCS equations. Our procedure does not requirethe starting wave function to be normalized. Given aproduct form of Eq. (23), it is efficient to generate theBCS based initial state stochastically by selecting | n i ina process where each randomly selected state k is chosento be occupied or empty with a probability proportionalto the corresponding v k and u k ; the process is stoppedonce a desired number of occupied states given by thetotal number of particles is reached. It is important thatvariations in implementation of the projected BCS or notfollowing Eq. (23) exactly, are not essential since thestarting wave function can be arbitrary.In Fig. 1(a) the convergence of energy as a function of L, following Eq. (19), is shown for the two initial wavefunctions. The method just outlined is based on eval-uation of the linear norm, the sum of all walkers’ bags,and since the actual wave functions are never constructedthese are labeled as “no-wf” in Fig. 1. In panel (b) theconvergence is shown as a function of the imaginary time τ using Eq. (20). The common energy scale is used inboth panels and the energy obtained from BCS approachand from the exact diagonalization of the pairing Hamil-tonian are shown with horizontal grid lines. The magni-fied energy scale used here allows one to clearly see thedifference between BCS and the exact solution.As appropriate in a variational technique, the BCSground state energy is above the exact one. However,the estimates for the ground state energy, using the lin-ear norm, approach the exact value from below. Thisfeature, as discussed in Sec. III B, is used for providing alower bound for the ground state energy estimate.In order to compare the projection with power func-tion in Eq. (19), and using the exponential in Eq. (20),Fig. 1(b) includes an additional L scale shown at the top.The quantity L is defined as the average number of stepsthat needs to be taken by walkers in order for the se-ries (21) to converge for a given imaginary time τ. While E L no-wf wf(a) FermiBCS 0.5 1 1.5 2 2.5 ExactBCS50 100 150 τ L (b) FIG. 1: (Color online) The half-occupied ladder model with ω = 18 , n = 9 and G = 1 is used to show the convergence of groundstate energy using an approach based on the linear norm. Left panel (a) shows the method based on ground state projectionwith power law Eq. (19), and right panel (b) shows projection using imaginary time evolution Eq. (20). In both panels BCSand exact values of energy are shown with horizontal grid lines; the four curves represent two initial states and two methods:with and without wave functions being built. both panels (a) and (b) look similar, using the exponen-tial as a projector is more computationally expensive asit requires almost three times as many steps.The use of exponent to project a ground state does notprovide any additional numerical stability; fluctuationsat remote times, in cases with no wave function, are seenin both panels of Fig. 1. These fluctuations are removedby reconstructing wave functions at certain steps; thecorresponding curves in Fig. 1 are labeled with “wf”.The origin of these fluctuations and error analysis areaddressed next. Since the exponential projection usingimaginary time is deemed to be less effective we will notdiscuss it any further. B. Error and convergence control
In the CSMC there are generally two kinds of errors.The first one is the statistical error that emerges as aresult of stochastic evaluation, for example, estimatingwave functions using Eq. (8) or evaluation of the linearnorm in Eq. (18). The second error is associated with thealgorithm used to obtain physical quantities of interest;for example, in projection technique this concerns thequality of approximation E ( L ) ≈ E in Eq. (19). In thissubsection we examine both of these errors and methodsof their control.The Central Limit Theorem (CLT) is at the core ofstatistical error control. It is usually expected that asthe number of samples, N, grows the associated stan-dard deviation σ of the ensemble average goes down as σ ∝ / √ N. However, in CSMC the main disadvantage of independent walks is that the variance grows exponen-tially as the path length increases. Therefore, for largenumber of steps the average of bags in Eq. (18) is hardto evaluate because the distribution of bags becomes toobroad. This is the cause of fluctuations seen in Fig. (1)at large L. Let us analyze this problem. Consider an ensemble ofall L -step bags for all possible paths { B L } , let σ { B L } be its variance and B L its mean. According to Eq. (18) B L = || Φ L || . As proved earlier, all bags are positive mak-ing the coefficient of variation Cv { B L } ≡ σ { B L } /B L anappropriate measure of relative error. Indeed CLT im-plies that with N estimates of energy using Eq. (19) therelative error is∆ E ( L ) /E ( L ) ≈ √ N Cv { B L +1 } . (24)The problem with divergent behavior of Cv { B L } as afunction of L arises due to B ( s ) L for each walker s being aproduct of matrix elements weighted by the correspond-ing probability, see Eq. (13). The product of a largenumber of random matrix elements is poorly behaved.Let us assume that c gives the coefficient of variation forall possible matrix elements weighted by chosen probabil-ities, then each term in the product (13) has a coefficientof variation, Cv (cid:26) h n L | V | n L − iP ( n L − → n L ) (cid:27) ≡ c . (25)Then the standard statistical treatment of a productleads to Cv { B L } = q (1 + c ) L −
1; (26)this simple form is obtained under the assumption of uni-form initial distribution in Eq. (12). With the exemptionof some special cases, c >
0; in most situations of inter-est c ≈ . , therefore Cv { B L } grows exponentially with L. Thus, it is practically impossible to compete with theexponentially increasing variance of the distribution byincreasing the number of walkers.In Fig. (2) the behavior of Cv { B L } as a function of L is shown for the half-occupied 18-level ladder model. Theresults for two different starting states are shown. Theexponential divergence in Eq. (26) pertains to the sit-uation involving independent walkers, where the actualwave function is not obtained. The two correspondingcurves, labeled with “no-wf”, both display the same ex-ponential divergence with c ≈ .
42 which corresponds toasymptotic behavior, Cv { B L } ∝ . L . -1 C v L no-wf (a)no-wf (b) wf (a)wf (b) FIG. 2: (Color online) The half-occupied ladder model with n = 9 , ω = 18 , and G = 1 . Coefficient of variation Cv { B L } isshown as a function of L . The figure includes four curves withtwo possible initial states: (a) Fermi state Φ (Fermi)0 Eq. (22)and (b) constant-component state where h n | Φ (const)0 i = 1; andfor calculations with and without wave functions being ob-tained. The limitation on the number of independent steps isrelatively easy to overcome. The exponential growth ofCv { B L } is usually weak, and in most cases, such as theexample in Fig. 1, no problems emerge for L less than 30or 50. Moreover, the choice of probabilities that followsimportance sampling in Eq. (14) would lead to c = 0 . Practically, numerical noise never allows one to reachthis ideal limit but the the growth of variance can bedelayed. In addition to that, with a good initial wavefunction, such as the one from BCS theory, the conver-gence is reached in a few steps, before the onset of sta-tistical problems; see example in Fig. 1.Preventing walkers from taking long independentwalks, by combining them in wave functions after a cer-tain number of steps with Eq. (8), allows one to avoid the problem completely. This is demonstrated in Fig. 1 withcurves labeled “wf”. The intermediate summation at amoment when bags are combined, due to the CLT, pre-vents an exponential increase of the variance. In the im-plementation of the CSMC the statistical error is trackedby controlling the coefficient of variation in the bags as L is increased. This allows one to apply a computation-ally expensive procedure of reconstructing the wave func-tion only when necessary, typically once in every 5–20steps. At a moment when the full wave function is built,the convergence of the projection technique (the secondkind of error) can be assessed using the usual square, L , norm.The second kind of errors, which is convergence E ( L ) → E in these examples, is a part of any itera-tive diagonalization technique, such as Lanczos or David-son algorithms; and it has been well studied in the past.However, the specifics of the pairing problem describedin Sec. II C and the use of the linear norm allow one toplace exact upper and lower limits on the value of energy.The convergence of the projection algorithm is exam-ined in Fig. 3. Here, using the same half-occupied laddermodel with ω = 18 , we show deviation of the predictedenergy from the exact value as a function of the number ofsteps. Three curves, that are essentially indistinct, show E ( L ) − E where E ( L ) was evaluated using the linearnorm L , Eq. (19). The three sets of results are obtainedby evaluating | Φ L i exactly with matrix-vector multiplica-tion (dotted line); using CSMC with wave function beingreconstructed at each step (dashed line); and without thewave function using bag average in Eq. (18) (solid line).The slight difference between exact and CSMC resultsis only due to an intermediate shift by chemical poten-tial. These three curves approach the exact energy frombelow, which is a distinct property of the L norm.The other two, nearly indistinct, curves show the con-vergence of energy evaluated using the traditional squarenorm, labeled as L E ( L ) = h Φ L | H | Φ L ih Φ L | Φ L i . (27)The L norm can be used only when the wave function isavailable, for that reason only the curves for exact (dash-dot) and CSMC with wave function (short dash) appearin Fig. 3. Naturally, the expectation value of the Hamil-tonian in Eq. (27) is subject to the variational principle,and all curves with L norm approach ground state en-ergy from above. Thus, the estimates using L , Eq. (19),and L , Eq. (27) norms give the lower and upper boundsfor the value of ground state energy.To summarize, in our algorithm we rely on computa-tionally inexpensive independent propagation of walkersin configuration space until the coefficient of variationof their bags exceeds some critical value. At that mo-ment the full wave function is reconstructed and is usedto evaluate energy from Eq. (27) and all other operatorsof interest. The combination of energy estimates fromlinear and square norms give lower and upper bounds for -60-50-40-30-20-1001020 10 20 30 40 50 60 70 80 E ( L ) - E L L L EXMC wf MC no-wf EX MC wf
FIG. 3: (Color online) The half-occupied ladder model with n = 9 , ω = 18 , and G = 1 . Deviation of the energy es-timates using linear L and square L norms from Eqs. (19)and (27), respectively, is shown as a function of L. The curvescorrespond to exact, CSMC with and without wave functionreconstruction. the actual value of energy. If the desired convergence isnot reached, the process is continued starting from thecurrent wave function.
C. Weak pairing limit
As mentioned earlier, superconducting paired states insmall systems face a lot of competition from other in-coherent interactions as well as from the single particleshell structure. Thus, relatively weak and fragile super-conducting states is one of the distinct characteristics ofpairing in nuclei. Unfortunately, the BCS theory is notdesigned to work in this limit, and having the CSMC asa computationally inexpensive alternative is one of themain motivations of this work. In Fig. 4 we demonstratethe effectiveness of CSMC in the limit of weak pairingusing our half-filled 18-level ladder model. In the limitwhen the pairing strength G = 0 , the system settles inthe Fermi state with 9 lowest double-degenerate single-particle states being occupied. As soon as G > , pairexcitation promotes particles up, and the occupation ofthe upper 9 levels becomes non-zero. For very strongpairing, G ≫ ǫ, the limit of degenerate model with equaloccupancy of all states is reached. This limit leads tohalf of the 18 particles being on lower 9 levels and halfon the upper ones.In Fig. 4 the net occupation of the upper 9 levels as afunction of G is shown. This plot includes results fromBCS, CSMC (labeled as “MC”) and exact diagonaliza-tion. While on a large scale all results are similar, in the region of low pairing strength, which is shown in in-set, the well known problem with BCS solution, shownin dotted (black) line, is noticeable. At the same time,exact and CSMC results are indistinct; dashed (crimsoncolor) line goes right on top of the solid (sea-green color)line. This test illustrates that the CSMC is well suited O cc up a ti on G BCSMCExact 00.51 0 0.3
FIG. 4: (Color online) The half-occupied ladder model with ω = 18 and G = 1. The net occupancy (total number ofparticles) on the upper 9 orbitals is shown as a function ofthe pairing strength. The weak pairing limit is magnified ininset. The three curves correspond to BCS solution, CSMC(MC) solution, and the exact solution by means of diagonal-ization. The CSMC and exact results are indistinct and thecorresponding curves are overlaid. For the CSMC solution weused N = 7 . × walkers, limiting the number of indepen-dent steps to five. for all limits of pairing strength. Moreover, in the limitof weak pairing, the computational effort in CSMC isreduced as the contribution from rare excursions aboveFermi surface can be easily evaluated with importancesampling. D. Excited states
Obtaining excited states with CSMC is more compu-tationally difficult. One can no longer use a linear normsince all amplitudes cannot be positive definite simulta-neously; that is, statement (vi) in Sec.II C is not validfor excited states. Therefore, the usual quadratic L norm has to be used and the bag values can be neg-ative. Nevertheless, features (i)-(v) in Sec.II C remainvalid and useful. In particular, since the matrix elementsof V are positive definite, the importance sampling is stillan effective strategy and the signs of bags are not alteredby repeated application of the Hamiltonian which cur-tails the typical MC sign problem. Similar, to Lanczostechnique, the CSMC approach requires orthogonaliza-tion, therefore the wave functions have to be built eachtime the orthogonalization is to be performed. The needfor orthogonalization limits the number of independentsteps, which is the main reason for higher computationaldemand.In Fig. 5 we show the CSMC applied to the study ofexcited states in the same half-occupied 18-level laddermodel. The ladder model example is particularly chal-lenging since the density of states above the gap is high.In this model the level spacing between the ground andfirst excited state, which is about twice the BCS gap, is E − E ≈ . ǫ = 1). Atthe same time, the spacing between the following states E − E ≈ .
35 40 45 50 55 60 65 70 10 20 30 40 50 60 70 80 90 100 E L E E E E E E E FIG. 5: (Color online) The half-occupied ladder model with ω = 18 and G = 1. Convergence of CSMC to ground state,to first excited state and to double-degenerate second excitedstate is shown. E. Pairing in Sn isotopes
In order to illustrate the CSMC algorithm in a realisticcase where pairing matrix elements are not all equal toa constant, we consider isotopes of tin. The role of pair-ing in − Sn isotopes has been extensively exploredin the literature [10, 28, 29, 53]. Apart from questionsof scientific interest such as pairing matrix elements andtheir connection to superconducting state in infinite mat-ter, near constancy of the excitation energy of the lowest2 + states, and unexplained behavior of electric quadru-ple transition rates, the tin case emerged as a benchmarkfor computational techniques. In Tab. I we present com-parison of energies and occupation numbers for Sn,
Sn, and
Sn. The model space here includes fivesingle-particle levels (with total ω = 16), their energies and spins are listed in the first two columns of Tab. I,the matrix elements are taken from from the G-matrixcalculation in Ref. [54], the values can be found in Ta-ble 1 of Ref. [10]. The results in Tab. I show expectedlevel of agreement. With increased computational effort,mainly using large number of walkers, any desired levelof precision can be obtained; our goal here was to useminimal effort and to solve the pairing problem with aprecision that exceeds any practical need, which is set tobe 5 keV uncertainty for energy and 0 .
01 for occupationnumbers.
F. Large scale model
As a final illustration of the CSMC algorithm we ex-plore a model of the O nucleus intended to reflect thenature of pairing correlations in a system containing bothbound states and a continuum of scattering states. Ourmain goal is to demonstrate the capabilities of our al-gorithm while addressing the problem of pairing in con-tinuum qualitatively. Quantitative studies require goodknowledge of the effective interaction Hamiltonian; con-struction of this Hamiltonian is outside the scope of thispresentation.For our study we select the Woods-Saxon potentialwith parameters from Ref. [55] to model the mean field ofweakly-bound O nucleus. We discretize this potentialusing a large quantization-box of size 500 fm. This al-lows us to generate a dense continuum of states. We limitscattering states by about 8 MeV of energy, which leadsto ω of about 100. For the pairing interaction betweenneutrons we use a density dependent contact interactionfrom Refs. [56–58] V ( r , r ′ ) = − G (cid:18) − η ρ ( r ) ρ (cid:19) δ ( r − r ′ ) . (28)Here ρ ( r ) /ρ is the nucleonic density expressed relativeto the saturation density. This quantity is assumed to begiven by the Woods-Saxon form factor. The density de-pendence of pairing is controlled by a parameter η whichis selected as η = 0 . . Following Ref. [56, 58], we alsointroduce a momentum cut-off function, that graduallyreduces the pairing matrix elements to zero for scatteringstates at energies above 5 MeV, the diffuseness parameterof the cut-off function is 0.5 MeV, see also Ref. [59].For our example we assume an inert O core whichleaves two bound s / and d / valence single-particlestates. Therefore, the bound states can accommodate n = 4 pairs of valence neutrons in O. The pairingmatrix elements involving these states are known fromthe phenomenological shell model Hamiltonian in Ref.[60]. Following previous studies, we adopt the value of G = 1 GeV · fm for treating pairing interaction involv-ing the continuum of scattered states. This value is alsoconsistent with the pairing strength in phenomenologi-cal Hamiltonians [60]. Here we limit our considerationto s -wave single-particle continuum. Due to centrifugal0 Sn Sn SnExact CSMC Exact CSMC Exact CSMCj E (MeV) -153.766 -153.765 -170.115 -170.113 -185.945 -185.9455 / − .
736 4.99 4.99 5.13 5.13 5.25 5.257 / − .
957 5.88 5.89 6.28 6.27 6.60 6.601 / − .
302 0.67 0.67 0.76 0.76 0.86 0.863 / − .
634 1.18 1.17 1.63 1.63 2.08 2.0911 / − .
544 3.29 3.28 4.21 4.20 5.21 5.21TABLE I: Comparison of exact and CSMC results for selected isotopes of tin. After header, first row shows comparison ofenergies, the remaining five rows show occupation numbers for five singe particle states. The calculation of energies is donewith N = 5 × walkers using linear norm. The final error is about 5 keV. barrier the overlap between bound and unbound d -wavestates is small; this inhibits virtual pair excitations to d -wave states in the continuum.Our goal in this investigation is to estimate the roleof continuum. We do this by comparing full calculationwith the one where the continuum is ignored. In Fig. 6the change in the ground state energy ∆ E is shown as afunction of the single-particle energy ǫ of the s / state.We present two different cases. In case (a) the parame-ters of the pairing Hamiltonian, which includes the singleparticle energies and pairing matrix elements, are firstevaluated with a realistic choice of Woods-Saxon param-eterization for O and then the ǫ which corresponds to s / state, is varied while all other parameters remain un-changed. In the self-consistent case (b) the depth of theWoods-Saxon potential is varied which moves the s / state, and each time a new configuration space Hamilto-nian matrix is calculated and studied.Let us summarize this study. First, the correction frompair excitations into continuum appears to be relativelysmall, here it is of the order of one kilovolt, for all rea-sonable choices of pairing strength G the effect is notexpected to exceed a few tens of kilovolts, see Ref. [59].The smallness of the effect does not seem to contradictobservations. So far there has been no significant near-threshold discontinuity observed in nuclear structure thatcan be attributed to two-body decay or to pair excita-tions. The decay of O is observed to be very slow, seediscussion [61], which through dispersion relations indi-cates weakness of the continuum coupling.Second, as expected, the effect increases sharply as thebound state approaches the continuum threshold. Thisis similar to the results known for single particle states,while the exact near threshold behavior is defined by thephase space volume, see Refs. [62, 63].Third, the difference between the two models high-lights the importance of halo phenomenon and its propertreatment. In model (b) the wave function of thesingle-particle s / state spatially extends as its energyapproaches the threshold. This facilitates pairing inthe continuum and the resulting effect is significantlystronger than that in model (a) where the spatial struc-ture of the single particle wave function was not modified.Finally, we find that this example successfully demon- strates the power of the CSMC method. -0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.10 -3 -2.5 -2 -1.5 -1 -0.5 0 ∆ E ( k e V ) ε (a)(b) FIG. 6: (Color online) The energy correction due to the in-clusion of continuum states for a pairing model with ω = 100and n = 4. The correction, ∆ E , is the difference betweenthe CSMC result with continuum states and an exact answerfor a model including only the two bound states. The ǫ alongthe lower axis is the energy of the s-wave bound state as itis moved closer to the continuum. The Monte Carlo error isnegligible. IV. SUMMARY
In this work we put forward a new Configuration SpaceMonte Carlo method for solving the many-body pairingproblem in finite systems. Unlike previous Monte Carlotechniques that deal with pairing interaction in a waysimilar to MC methods in physics of spin systems, ourapproach does not use evolution in imaginary time, doesnot need Trotter-Suzuki propagator breakup, and doesnot depend on the pairing matrix elements being con-stant. We propose to evaluate Hamiltonian and otherobservables by stochastically evaluating the correspond-ing operators in the Krylov subspace spanned by states1formed as a result of powers of Hamiltonian acting on anarbitrary initial state. States in the Krylov subspace areevaluated using random walks in the many-body config-uration space. Importance sampling is used to effectivelyprobe components of the wave functions. We emphasizeseveral important features of the pairing Hamiltonian,that make the MC approach appealing. In particular,we stress boson-like behavior of nucleon pairs, absenceof the fermion sign problem, potential for probabilisticinterpretation of transitions in configuration space, andprobabilistic interpretation of ground state amplitudes.In addition to traditional quadratic quantum mechanicalnorm, probabilistic interpretation allows us to use a lin-ear norm. We demonstrate that the approach based onthe linear norm is computationally efficient, is perfect forparallelization, and provides effective methods for controlof errors of both stochastic and non-stochastic origins.The workings of the CSMC method are demonstratedwith several examples. With a classic ladder model wedemonstrate convergence using several variations of themethod; we discuss errors and present effective means oftheir control. The effectiveness of CSMC in small sys- tems where pairing can be effectively weak is shown; theCSMC in its most complete form is used for obtainingdegenerate and nearly-degenerate excited states in theladder model.As a realistic example, we use isotopes of tin whichrepresents another well studied classic case of pairingin nuclei. The energies and occupation numbers in thisnon-constant pairing example are consistent with exactresults. Large-scale study of pairing correlations is illus-trated using a model of O that includes a continuum ofscattering states. While our last example is still far fromrealistic, it highlights the effectiveness of CSMC, and sug-gests an arena for future applications of the method.
Acknowledgments
This material is based upon work supported by theU.S. Department of Energy Office of Science, Office ofNuclear Physics under Award Number de-sc0009883. [1] R. A. Broglia and V. Zelevinsky,
Fifty years of nuclearBCS (World Scientific, 2013).[2] A. Spyrou, Z. Kohley, T. Baumann, D. Bazin, B. A.Brown, G. Christian, P. A. DeYoung, J. E. Finck,N. Frank, E. Lunderberg, et al., Phys. Rev. Lett. ,239202 (2012).[3] C. R. Hoffman, T. Baumann, J. Brown, P. A. DeYoung,J. E. Finck, N. Frank, J. D. Hinnefeld, S. Mosby, W. A.Peters, W. F. Rogers, et al., Phys. Rev. C , 031303(2011).[4] Z. Kohley, E. Lunderberg, P. A. DeYoung, A. Volya,T. Baumann, D. Bazin, G. Christian, N. L. Cooper,N. Frank, A. Gade, et al., Phys. Rev. C , 011304(2013).[5] J. Bardeen, L. Cooper, and J. Schrieffer, Phys. Rev. ,1175 (1957).[6] A. Bohr, B. Mottelson, and D. Pines, Phys. Rev. ,936 (1958).[7] S. T. Belyaev, K. Dan. Vidensk. Selsk. Mat. Fys. Medd. , 31 (1959).[8] A. K. Kerman, Ann.Phys. , 300 (1961).[9] J. Dukelsky, G. G. Dussel, J. G. Hirsch, and P. Schuck,Nucl. Phys. A , 63 (2003).[10] V. Zelevinsky and A. Volya, Physics of Atomic Nuclei , 1829 (2003).[11] V. Zelevinsky and A. Volya, Nucl. Phys. A , 299(2004).[12] S. C. Pang and A. Klein, Can. J. Phys. , 655 (1972).[13] J. Y. Zeng and T. S. Cheng, Nucl. Phys. A , 1 (1983).[14] G. D. Dang and A. Klein, Phys. Rev. , 735 (1966).[15] G. D. Dang and A. Klein, Phys. Rev. , 689 (1966).[16] H. J. Lipkin, Ann. Phys. , 272 (1960).[17] Y. Nogami, Phys. Rev. B , B313 (1964).[18] F. Pan, J. P. Draayer, and W. E. Ormand, Phys. Lett. B , 1 (1998). [19] A. Volya and V. Zelevinsky, in , edited byC. Francesco, M. Chadwick, A. Ferrari, T. Kawano,S. Bottoni, and L. Pellegri (CERN Conference Proceed-ings Series, Geneva, 2012).[20] G. Racah, Phys. Rev. , 186 (1942).[21] G. Racah, Phys. Rev. , 438 (1942).[22] G. Racah, Phys. Rev. , 0367 (1943).[23] K. T. Hecht, Phys. Rev. , B794 (1965).[24] J. N. Ginocchio, Nucl. Phys. , 321 (1965).[25] N. Auerbach, Nucl. Phys. , 321 (1966).[26] V. K. B. Kota and J. A. C. Alcaras, Nucl. Phys. A ,181 (2006).[27] J. Dukelsky and G. Ortiz, Int. J. Mod. Phys. E , 324(2006).[28] B. A. Brown, R. R. C. Clement, H. Schatz, A. Volya, andW. A. Richter, Phys. Rev. C , 045802 (2002).[29] W.-C. Chen, J. Piekarewicz, and A. Volya, Phys. Rev. C , 014321 (2014).[30] R. W. Richardson, Phys. Rev. , 1007 (1967).[31] R. Richardson, Phys. Rev. , 874 (1966).[32] R. W. Richardson, Phys. Rev. , 792 (1967).[33] R. W. Richardson and N. Sherman, Nucl. Phys. , 221(1964).[34] R. W. Richardson, J. Math. Phys. , 1034 (1965).[35] L. Cooper, Phys. Rev. , 1189 (1956).[36] J. Dukelsky, S. Pittel, and G. Sierra, Rev. Mod. Phys. , 643 (2004).[37] S. Pittel and J. Dukelsky, Phys. Scripta T125 , 91 (2006).[38] J. Dukelsky and G. Ortiz, Int. J. Mod. Phys. E , 324(2006).[39] J. Dukelsky, C. Esebbag, and S. Pittel, Phys. Rev. Lett. , 062501 (2002).[40] R. R. Whitehead, Nucl. Phys. A , 290 (1972).[41] A. Volya, B. A. Brown, and V. Zelevinsky, Phys. Lett. B , 37 (2001).[42] T. Sumaryada and A. Volya, Phys. Rev. C , 024319(2007).[43] D. P. Landau and K. Binder, A guide to Monte Carlosimulations in statistical physics (Cambridge UniversityPress, Cambridge, 2005), 2nd ed.[44] S. E. Koonin, D. J. Dean, and K. Langanke, Ann. Rev.Nucl. Part. Sc. , 463 (1997).[45] S. E. Koonin, D. J. Dean, and K. Langanke, Phys. Rep. , 2 (1997).[46] S. C. Pieper and R. B. Wiringa, Ann. Rev. Nucl. Part.Sc. , 53 (2001).[47] N. Cerf and O. Martin, Phys. Rev. C , 2610 (1993).[48] N. Cerf, Nucl. Phys. A , 383 (1993).[49] H. F. Trotter, Proc. Amer. Math. Soc. , 545 (1959).[50] M. Suzuki, Phys. Lett. A (1990).[51] A. Mukherjee, Y. Alhassid, and G. F. Bertsch, Phys. Rev.C (2012). [53] D. J. Dean and M. Hjorth-Jensen, Rev. Mod. Phys. ,607 (2003).[54] A. Holt, T. Engeland, M. Hjorth-Jensen, and E. Osnes,Nucl. Phys. A , 41 (1998).[55] N. Schwierz, I. Wiedenhover, and A. Volya,arXiv:0709.3525 (2007).[56] P. Bonche, H. Flocard, P. H. Heenen, S. J. Krieger, andM. S. Weiss, Nucl. Phys. A , 39 (1985).[57] R. R. Chasman, Phys. Rev. C , 1935 (1976).[58] J. Dobaczewski, W. Nazarewicz, T. R. Werner, et al. ,Phys. Rev. C , 2809 (1996).[59] M. Lingle, Ph.D. thesis, Florida State University (2015).[60] B. A. Brown and W. A. Richter, Phys. Rev. C , 034315(2006).[61] A. Volya and V. Zelevinsky, Phys. At. Nucl. , 969(2014).[62] A. Volya, EPJ Web of Conf. , 03003 (2012).[63] A. Volya and V. Zelevinsky, AIP Conf. Proc.1619