Strong Boundary and Trap Potential Effects on Emergent Physics in Ultra-Cold Fermionic Gases
SStrong Boundary and Trap Potential Effects onEmergent Physics in Ultra-Cold Gases
J.B. Hauck , C. Honerkamp , and D.M. Kennes , , Institut f¨ur Theorie der Statistischen Physik, RWTH Aachen, 52056 Aachen,Germany and JARA - Fundamentals of Future Information Technology Institute for Theoretical Solid State Physics, RWTH Aachen University, 52074Aachen, Germany JARA-FIT, J¨ulich Aachen Research Alliance - Fundamentals of Future InformationTechnology, Germany Max Planck Institute for the Structure and Dynamics of Matter and Center forFree Electron Laser Science, 22761 Hamburg, GermanyE-mail: [email protected]
February 2021
Abstract.
The field of quantum simulations in ultra-cold atomic gases has beenremarkably successful. In principle it allows for an exact treatment of a variety ofhighly relevant lattice models and their emergent phases of matter. But so far thereis a lack in the theoretical literature concerning the systematic study of the effects ofthe trap potential as well as the finite size of the systems, as numerical studies of suchnon periodic, correlated fermionic lattices models are numerically demanding beyondone dimension. We use the recently introduced real-space truncated unity functionalrenormalization group to study these boundary and trap effects with a focus on theirimpact on the superconducting phase of the 2D Hubbard model. We find that inthe experiments not only lower temperatures need to be reached compared to currentcapabilities, but also system size and trap potential shape play a crucial role to simulateemergent phases of matter.
Submitted to:
New J. Phys.
1. Introduction
Ultra-cold atomic gases provide a powerful method for simulations of quantum many-body lattice systems [1, 2]. Many models introduced in the context of solid-state physics,like the Hubbard model [2], the Haldane model [3], the bilayer Hubbard model [4] or evenquasicrystalline models [5], can be investigated in laboratories without the restriction ofnumerical ambiguities or theoretical approximation errors. Beyond this, other excitingmany-body phenomena that cannot be easily explored in solids could be realized inoptical lattices as well [6–10], making them a very flexible tool to study quantum many-body physics. a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b trong Boundary and Trap Potential Effects d -wave superfluidity [16], where, e.g., QMCencounters the notorious sign problem and finite size scaling from exact diagonalizationis hopeless. Despite the lack of quantitative studies of trap or finite size effects in thisregard, the current hope is that if the experiments can reach lower temperatures, thesuperfluid phase will emerge naturally.However, the effects of the trapping potential have always been an issue in opticallattices [17–20], even for the Mott-phase and high temperature properties of one andtwo dimensional Hubbard models. In extension of these previous results, we investigatethe effects of the trapping potential, as well as the boundary conditions on the Cooperpairing in a Hubbard model. This is of high experimental importance as one of thecurrent goals is to measure the superfluid phase. A general challenge for Cooper pairingin finite-size systems can, e.g., be gleaned from the literature for superconductivity ofultra-small aluminium grains [21, 22], with a diameter of the order of ∝ nm. There, itwas worked out clearly how the usual superconducting description is valid as long asthe energy gap due to pairing is larger than the energy spacing due to the confinementbut breaks down if the latter scale is the larger one. A similar crossover is expected tohappen also in a Hubbard model. Even from these basic considerations, the questionarises whether there is a minimal number of sites required to unveil superconductinghallmarks. What might add in the case of d -wave pairing is the question, if thegeneration of the attractive interaction, at weaker coupling mainly by spin fluctuations,is somehow altered by the trap confinement. Finally, it is an open question how thenon-trivial d -wave ordering pattern will be influenced by the trapping potential. In theantiferromagnetic case the trap was not observed to be a significant perturbation [14].In fact the data compared well to Quantum Monte Carlo results obtained for a squarelattice with periodic boundary conditions. Here we want to examine whether thisholds true also for the d-wave superconducting phase in a square lattice Hubbardmodel. To extrapolate to the bulk case (e.g., to connect to real materials), oneshould also understand how closely the finite-size set-up resembles the situation in thethermodynamic limit. Hence, we will also examine at which size the trapped systemapproximates the bulk case to high fidelity.
2. Setup
We will simulate five different setups, all are different realizations of the square latticeHubbard model with different open boundary conditions. They differ in the lattice trong Boundary and Trap Potential Effects V trapi = a trap r i , where r i is the radial distance of site i to the center of the latticeand a trap gives the curvature of the trapping potential. The fifth and last setup is thesame as the fourth with a different potential: Instead of the simple quadratic one, wereconstruct a potential used in a recent experiment by Mazurenko et. al. [14]. There thetrap consists of a Gaussian potential created by a digital micro-mirror device (DMD)with an additional circular quadratic potential, whose superposition creates a nearlyflat disc at the center. For simplicity we approximate it by a piece-wise definition, seeEq. (14). The Hamiltonian in second quantization for all cases reads H = − (cid:88) i,j,σ [ t ij + ( µ + V trapi ) δ i,j ] c † i,σ c j,σ + 12 (cid:88) i,j,σ,σ (cid:48) U i,j n i,σ n j,σ (cid:48) , (1)with the operators c ( † ) i,σ annihilating (creating) an electron on site i with spin σ . For thehopping amplitudes we set t ij = t if i and j are nearest-neighbors and t ij = t (cid:48) if i and j arenext-nearest-neighbors. All other terms are set to zero. For simplicity we choose t = 1and measure all quantities in units of t . Distances are measured in units of the latticespacing a which we set to 1. We choose onsite interactions U ij = U δ i,j and fix the ’Van-Hove condition’ µ = 4 t (cid:48) . In experiments so far, the model mostly used was a standardHubbard model with nearest-neighbor hoppings t only. We use an extended Hubbardmodel as recent numerical studies [23–25] suggest that in the relevant U parameter rangefor cuprates there is no superconductivity in the pure ( t (cid:48) = 0) Hubbard model.
3. Method
We employ the recently introduced 2D real-space functional renormalization group (real-space TUFRG) [26], which can be used to study models with broken translationalsymmetry. It is based on the standard one-particle-irreducible FRG [27, 28], in a level-2truncation. Additionally, we will neglect the frequency dependence of the vertex, as wellas the self-energy. With this we obtain a method, which scales relatively well and is ableto predict finite-size precursors to phase transitions and their respective orderings. Themain idea of FRG is the introduction of a cutoff in the bare propagator, enabling aninterpolation between a solvable theory and the full solution of the model. Variations trong Boundary and Trap Potential Effects −
20 0 20 x − − y V i −
20 0 20 x − − y V i −
20 0 20 x − − y V i −
20 0 20 x − − y V i −
20 0 20 x − − y V i Figure 1.
Schematic visualization of the five different geometries and trappingcombinations considered here. In the upper left, the simple square lattice with openboundary conditions is visualized and in the upper middle plot, the circular openboundary conditions case is shown. The upper right plot shows the combination ofopen boundary conditions with a finite height trapping. The lower plots visualize thetwo different choices of potentials applied in the circular open boundary conditions case,where the left is the simple quadratic trap, and the right approximates the experimentaltrap. The color scale gives the values of the trapping potential V i . with respect to the cutoff parameter generate an infinite set of coupled differential flowequations for the Taylor coefficients of the effective action. These equations must betruncated in order to become numerically tractable, which is well defined in the case of UW small, then only the first few terms in the expansion are relevant and a truncation at thelevel of keeping just the self-energy and the effective interaction is justified. Explicitlythe cutoff is introduced in the bare Green’s function by G ( ω n ) = 1 iω n − H −→ G Λ0 ( ω n ) = R (Λ) iω n − H (2)with R (Λ) being the cutoff function obeying R (0) = 1 and lim Λ →∞ R (Λ) = 0 in ourconvention. This reduces the action in the limit λ → ∞ to the non-interacting,exactly solvable system. By integrating the coupled flow equations, an interpolationbetween this solution of the non-interacting problem and the approximate full solution iscalculated. The variations w.r.t the cutoff parameter generate the single scale propagatordefined as S Λ ( ω n ) = ˙ R (Λ) iω n − H (3)in our approximation.The scale-dependent two-particle vertex can be separated in three different channelsaccording to the three different fermionic bilinears or three different diagram types as trong Boundary and Trap Potential Effects dd Λ = ˙ P Λ + ˙ C Λ + + +˙ D Λ Figure 2.
Diagrammatic representation of the FRG flow equation for the one-particleirreducible interaction vertex in the level-II truncation, adapted from [27]. The firstterm on the right-hand side, annotated with ˙ P Λ , represents the flow of the particle-particle channel P Λ , while the other terms generate the flow of C Λ and D Λ . Thedashed internal lines represent a single scale propagator, see Eq. (3), and the solidlines represent a propagator, see Eq. (2). The arrows indicate the direction of thepropagators and fix the order of their arguments. The grey boxes represent the fulltwo-particle interaction. A second class of diagrams with the scale-derivative on theother internal line contributes as well but is now drawn for simplicity. shown in Fig. 2. Each of those bilinears is directly associated to a bosonic field ofan effective Hamiltonian ‡ . The particle-particle channel ( P Λ ) is associated with thesuperconducting pairing interactions and therefore a divergence of this channel at totalmomentum zero indicates a tendency towards superconductivity. The direct particle-hole channel ( D Λ ) is associated to charge fluctuations and thus indicates, among others,general charge ordering. The crossed particle-hole channel ( C Λ ) can be associatedwith effective spin-spin interactions, therefore its divergence is an indicator of magneticordering. In typical flows, the effective interaction undergoes a flow to strong coupling,i.e. diverges at a certain finite RG scale. When this occurs, we stop the flow assoon as an eigenvalue of one of the three channels surpasses a threshold value. Thecorresponding RG scale is called the critical scale Λ c , which defines an energy scale ofthe breakdown of the weakly interacting behavior and can be seen as an estimate ofan ordering temperature T c . The leading eigenvalue of the diverging channel gives, byits eigenvector, the expected ordering pattern up to an unknown prefactor. A moreprecise determination of the order would follow from the gap equations for each of thechannels. In practice we will normalize the shown ordering pattern from the leadingeigenvector such that their maximal absolute value is 1. As we consider finite-sizesystems here, the interpretation of the flows to strong coupling as phase transitions isnot strictly valid, as neither a real discontinuity occurs, nor any correlation length truly ‡ In fact these bilinears are interpreted most easily for the particle-particle ( P ) and the crossed particle-hole ( C =spin) channels, the direct particle-hole channel ( D ) is more complex as it contains informationon the C -channel. The physical relevant quantity is the charge channels given by V ch = 2 D − C , whichremoves the redundant information. trong Boundary and Trap Potential Effects U i,i = ii → P (2) = ii jj Figure 3.
Generation of new index dependencies during the flow. A density-densityinteraction U i,i n i,s n i,s (cid:48) receives corrections in the P -channel that can be interpretedas an incoming particle pair c i,s c i,s (cid:48) scattering to an outgoing pair c j,s c j,s (cid:48) . diverges. Nevertheless, it can be expected that the finite-size system with a stronglyenhanced interaction of a specific type will locally resemble a infinite system with a trulydiverging effective interaction of the same type. Thus, we use the same vocabulary asfor the infinite system and discuss different phases of the finite system as regimes withqualitatively different flows to strong coupling.Next, we discuss the spatial dependence of the effective interactions and therespective channels P , C and D . In principle, each channel can depend on four siteindices. Yet, the initial condition in Eq. (1) is an onsite interaction, i.e. where all foursite indices have to be the same § . Looking at the second-order correction e. g. in the P -channel as in Fig. 3, we observe that this adds a non-local, potentially long-ranged, ’bi-local’ correction that depends on the joint index of the incoming lines i and the outgoinglines j . Like in the P -channel, these RPA-like contributions create a main dependenceon different bi-local pairs of indices in each channel. Contributions beyond these bi-localindices are only generated by the feedback in between the channels and, thus, come ata higher order in the interaction. The further apart from the native bi-local indices, thehigher is the interaction order [29]. Therefore, these terms can be neglected in the spiritof the perturbative motivation of the FRG. The bi-local contributions of each channel onthe other hand have to be captured more rigorously, as they are equivalent to the sharpmomentum structures that builds up in the translationally invariant case [28, 30, 31].Formally, this specific structure of the effective interactions can be exploited byan expansion of the subleading dependencies in form factors, which form a completeorthonormal basis set, according to (cid:88) i f b k ( i ) f ∗ b k (cid:48) ( i ) = δ b k ,b k (cid:48) , (cid:88) b k f b k ( i ) f ∗ b k ( i (cid:48) ) = δ i,i (cid:48) . (4)The three channels can then be build up from form factors capturing the internal § This argumentation generalises to non-onsite interactions but we focus on the simplest case for thesake of simplicity trong Boundary and Trap Potential Effects P [Γ] b i ,b j i,j = (cid:88) k,l Γ( i, k ; j, l ) f b i ( k ) f ∗ b j ( l ) , (5)ˆ C [Γ] b i ,b j i,j = (cid:88) k,l Γ( i, k ; j, l ) f b i ( l ) f ∗ b j ( k ) , (6)ˆ D [Γ] b i ,b l i,l = (cid:88) k,l Γ( i, k ; j, l ) f b i ( j ) f ∗ b l ( k ) . (7)With this at hand we can rewrite the flow equations by an insertion of form-factor unitmatrices. The main idea of the truncated unity approach is to restrict the number ofform-factors in these unities to much less than the lattice sites [26, 28, 32] (cid:88) b k ∈ L f b k ( i ) f ∗ b k ( i ) = (cid:88) b k δ Lb k f b k ( i ) f ∗ b k ( i ) ≈ , (8)where L is the set of all allowed bonds and analogously δ Lb k = 1 if b k ∈ L with bondlength up to b max .In the form-factor space the flow equations can be rewritten in terms of highlyefficient block-matrix products. Due to the truncation, the full vertex cannot berecovered exactly but instead projections must be introduced in order to reconstructthe full vertex approximately for the cross-channel feedback. These projections can befound for example in Ref. [26, 29, 33]. Instead of including the full spin dependence wewill assume an SU (2) invariant problem. This is by no means necessary but simplifiesthe presentation. The diagrammatic flow equation for the interaction vertex, for such a SU (2) invariant model, are displayed in Fig. 2.These diagrams translate to dd Λ P Λ = − ˆ P [Γ] Λ · ˙ χ pp · ˆ P [Γ] Λ ,dd Λ C Λ = − ˆ C [Γ] Λ · ˙ χ ph · ˆ C [Γ] Λ ,dd Λ D Λ = 2 ˆ D [Γ] Λ · ˙ χ ph · ˆ D [Γ] Λ (9) − ˆ C [Γ] Λ · ˙ χ ph · ˆ D [Γ] Λ − ˆ D [Γ] Λ · ˙ χ ph · ˆ C [Γ] Λ , with the particle-hole bubble ˙ χ ph and the particle-particle bubble ˙ χ pp defined as˙ χ b i ,b j ph ( i,j ) = 2 (cid:88) ω> (cid:60) (cid:0) G Λ ( ω ) i,j S Λ ( ω ) j + b j ,i + b i + G ↔ S (cid:1) (10)˙ χ b i ,b j pp ( i,j ) = 2 (cid:88) ω> (cid:60) (cid:0) G Λ ( ω ) i,j S Λ ( − ω ) i + b i ,j + b j + G ↔ S (cid:1) . (11)For the discussion of the results we will often refer to the s -wave and d x − y -wave component of the leading eigenvector, which we define as follows: The leadingeigenvector of a channel I b i ,b j i,j depends on a single site and a bond index, thus it can trong Boundary and Trap Potential Effects v b i i . The s -wave, or onsite component is in this notation defined as v b i ,where b refers to the bond with zero length. The definition of the d x − y componentfollows analogously to momentum space [28], we enumerate next-nearest-neighbor bondsfrom 1 to 4 in a clockwise fashion starting from the up pointing bond, then we define v d x − y i = v b i − v b i + v b i − v b i . For the integration, an adaptive fifth order Runge-Kutta scheme [34] of the odeintlibrary [35] is used. The matrix products on the right hand side of Eq. (9), whichare the numerically most demanding parts, are performed on GPU’s. To reduce thenumber of Matsubara frequencies needed for the calculation of the propagator, we use thePad´e-approximation of the Fermi-function [36, 37]. The results are then verified usinga heuristic tangent spacing scheme. Both yield fast convergence of the Bubble termstested against the analytical result. In the calculations we incorporate 300 frequencies,which is equivalent to an error of the bubble calculation of about 10 − at the lowestused temperature ( T = 10 − ) in comparison to the analytic result. We use a mapping ofthe infinite integration range to a finite one Λ = − xx to accelerate the integration [33].The integration is started from a minimum x = 10 − as smaller starting values yieldno significant difference in the final results. To reduce the load imbalance between thechannels we perform a completion of the square for the D -channel. Thereby we onlyneed to calculate a single matrix product, as the missing part for the completion is dC Λ d Λ which needs to be calculated anyway. As a regulator at finite temperature we chooseso-called Ω-cutoff, first introduced by Husemann et. al. [38], given by R (Λ , ω ) = ω ω + Λ ⇒ ˙ R (Λ , ω ) = − ω ( ω + Λ ) . (12)It has the advantage of a real IR-divergence regularization but is numerically morechallenging as for example the interaction cutoff [39]. Additionally, no numericalinstabilities arise if the Hamiltonian has a zero eigenvalue. At zero temperature weemploy a sharp cutoff, which drastically simplifies the bubble calculations [29, 40]. Intotal, the flow equations scale like O ( N · ¯ N b ) ( N is the number of orbitals and ¯ N b is theaverage number of bonds per site) due to the matrix products, the bubble calculationscales like O ( N f · N · ¯ N b ) ( N f is the number of Matsubara frequencies included inthe summation). Upon including the self-energy as well as the frequency dependenceof the vertex (in the single channel coupling or ECLA sense [33, 41]) the error scaleswith U , this scaling has been checked for the 1D-Hubbard model comparing to exactdiagonalization. Additionally, we verified that our implementation is sufficiently efficientto reach the thermodynamic limit of the 2D-Hubbard model, where the known FRGresults are reproduced. trong Boundary and Trap Potential Effects
4. Results
We start by an analysis of the finite size effects in a square lattice Hubbard modelwith open boundary conditions and without any confining potentials. The Van-Hovecondition is fixed by setting µ = 4 t (cid:48) and we vary t (cid:48) at different system sizes to find atwhich size the critical scales, as well as the occurring phases, are converged. ∞ AFMSC CDWFM SDW0 . . . . − t t l − − − Λ c ∞ AFMSC CDWN FMSDW0 . . . . . − t t r − − − − Λ c Figure 4.
The phase diagram of the square lattice Hubbard model with open boundaryconditions and without any confining potentials is shown on the left and the phasediagram of a square lattice Hubbard model with circular open boundary conditionsand without any confining potential is shown on the right. Simulations are performedat van Hove filling at varying system sizes and next-nearest-neighbor hopping t (cid:48) t . Wechoose U = 3, T = 10 − and include nearest-neighbor correlations. The criticalscales Λ c are encoded in the color scheme, whereas the type of the diverging phase isgiven by the shape of the data point. The thermodynamic limit results are given asa reference in the upper part of the plot. We encounter antiferromagnetism (AFM),charge density waves (CDW), general spin density waves (SDW), superconductivity(SC), ferromagnetism (FM) and no divergence at low scales (N). The finite size effects turn out to be rather smooth except for the smallest linear size l = 10 corresponding to 100 lattice sites, see the left plot in Fig. 4. The critical scales dovary only very slightly close to half filling. We find a superconducting phase even for thesmallest system. The superconducting transition gets shifted to higher values of − t (cid:48) t forlarger systems, the largest system investigated shows a phase diagram roughly matchingthe thermodynamic limit case. It should be noted that the occurring superconductingphase is not a pure d -wave divergence, as one would obtain with periodic boundaryconditions, but instead we have mixing between s - and d -wave bond pairing at all sizesat the boundary. Such an s - and d -mixing behavior has been a topic at the [110] surfaceof square lattice systems (e.g. [42]), but here it is more profound as the onsite s -wavecoupling is not put in `a priori, but generated by the RG flow. Note that the initial trong Boundary and Trap Potential Effects s -waveonsite terms, next to dominant nearest-neighbor pairs with d -wave pattern. Currentlyit is not clear which conditions, besides spatial inhomogeneity, are required to obtainthis accompanying s -wave component. In our approach, the phase between the twocomponents is fixed to unity as they arise in a single eigenvector. The phase couldchange in the gap equation as this is a non-linear relation that could alter the onsiteor bond pairings independently. Here, we will stick to the notion of d + s mixingwithout a phase, indicating the dominant d -character. For a proper verification ofthis finding, the possible gap opening should be investigated along the flow, which isnumerically very demanding and thus is not done in this paper. The magnitude of the s -wave component relative to the one of the d -wave component decreases with increasinglattice size. In the largest system, the s -wave component has smaller amplitude thanthe d -wave component. It has additionally a stronger site dependence, i. e. , it shows itshighest values at the boundaries. The superconducting phase does not fully resemblethe thermodynamic limit case due to the visibility of boundary effects in the whole bulkeven for l = 60.For the experiments, round trapping potentials are easier to generate than squaretrapping potentials. Therefore, the second setup under investigation is a square latticeHubbard model with circular open boundary conditions and without any confiningpotential, which can be seen as a prototypical model for a sharp and very high trap.Again finite size effects close to perfect nesting in the infinite system, or half filling,are more or less negligible with only slightly varying critical scales, see the right plotin Fig. 4. There, no ambiguity in the resulting phase occurs, an antiferromagneticdivergence is found at all sizes. Even the ordering pattern is not influenced much byvarying the system size. This is in agreement with the observations by Mazurenkoet. al. [14]. Upon increasing t (cid:48) the critical scales are decreasing. The finite size effectsbecome more pronounced with increasing t (cid:48) . At t (cid:48) ≈ − .
26 we find a transitionto a pairing divergence for the two largest systems. For these two, the criticalscales are comparable to the ones in the thermodynamic limit. We only encountersuperconductivity for a radius larger than r = 20, additionally the phase transitionbetween AFM and SC is shifted towards higher values of − t (cid:48) t . This already hints atpossible complications for the experimental setups. In the round case we are not able toreach the thermodynamic limit, even for r = 30, as the phase boundaries do not matchthe ones from the thermodynamic limit even at the largest radius. As for the squarelattice case, no pure bulk d x − y superconductor is recovered, as again in the completebulk a mixing with the s -wave component is present.We observe a strong size dependence for high − t (cid:48) t values where, in the infinite model,the ferromagnetic phase occurs. This can be explained by considering the orderings, seeFig. 5; Whereas the antiferromagnetic and superconducting phases more or less complywith the circular shape of the lattice (up to boundary effects and the already mentionedmixing), the ferromagnetic phase does not. In fact, the resulting magnetization pattern trong Boundary and Trap Potential Effects −
20 0 20 x − − − y − . − . . . . S D W p a tt e r n −
20 0 20 x − − − y − . − . . . . S D W p a tt e r n −
20 0 20 x − − − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − − − y − . − . . . . P a i r i n g p a tt e r n Figure 5.
Examples for the leading eigenvectors for each phase in the square latticeHubbard model with circular open boundary conditions and without any confiningpotential at r = 30. Calculations are performed using U = 3, T = 10 − and includingnearest-neighbor correlations. In the upper row the leading eigenvectors of the C -channel are shown at t (cid:48) = 0 in the left, which is an antiferromagnet, and t (cid:48) = − . s -wave (left) and symmetrized d x − y (right) contribution of the leading eigenvectorof the P -channel at t (cid:48) = − .
26. The P -channel eigenvector indicates d + s orderingtendencies. trong Boundary and Trap Potential Effects We now turn to the third setup, a square lattice Hubbard model with open boundaryconditions in which a second box potential with finite height is embedded. This isrealized by an on-site potential at the two sites closest to the boundary at the upperand the right edge and at the three sites closest to the boundary at the lower andleft edge. By this we can verify that the boundaries existing beyond the box potentialdo not influence our finding. The on-site potential is set higher than the bandwidth,i. e. V trap = 5 .
0, suppressing particles in this region. The difference to straight forwardopen boundaries is that due to the finite height, the electronic wavefunctions are allowedto be non-zero at the boundary. Inside the potential well, the wavefunction decaysexponentially. In this case, we are mainly interested in the superconducting phase, thuswe choose − t (cid:48) = 0 .
28. For simplicity we set T = 0 and apply a sharp frequency cutoff.The results are summarized in the left plot in Fig. 6.
10 20 30 40 l − − Λ c SC CDW xd x − y y s − P a i r i n g p a tt e r n r V HC − − − Λ c a trap =4E-04 a trap =1E-03 a trap =4E-03AFM SC N Figure 6.
We show the phases of a square lattice Hubbard model with open boundaryconditions in which a second box potential with finite height is embedded on the left.Here, the side length l of the inner boxes are given on the x -axis. On the right we showthe phases of a square lattice Hubbard model with circular open boundary conditions,in which we insert a quadratic trap and vary r V HC at different values of a trap . The y -axis gives the critical scales, and the shape of the data points encode the resultingphase. In the insets in the left plot, the s and d x − y component of the leading pairingeigenvector at the largest system size are shown. Calculations are performed using U = 3, T = 0 and including nearest-neighbor correlations at t (cid:48) = − .
28. We encounterantiferromagnetism (AFM), charge density waves (CDW), superconductivity (SC) andno divergence (N) at low scales.
At the smallest system size, we find a charge divergence. Starting from an effectivesystem length of l = 12 we observe a P -channel divergence. Again the s -wave componentis dominant at small system sizes and the divergences still show similar patterns as forthe open boundary case. But compared to the case without the embedded potential,the s -components are less pronounced in the bulk; thus, the confining potential reducedthe finite size effects. Thereby, such a potential could reduce the system size necessaryto resolve bulk d x − y -superconductivity. Such a square box trap potential is of course trong Boundary and Trap Potential Effects H = H + (cid:88) i,σ a trap | (cid:126)r i | c † i,σ c i,σ . (13)The coordinate system is arranged such that its origin coincides with the center of thelattice. The trapping is embedded in a round lattice with r = 30 and the temperatureis fixed to T = 10 − . The potential is at first fixed at the center, meaning that the Van-Hove condition is only fulfilled there. We now vary the trap curvature a trap to examinewhether superconductivity can still be observed. . . . a t r a p AFMSC NFM SDW0 . . . . . − t t − − − − Λ c ∞ AFMSC CDWFM SDW0 . . . − t t r t − − − Λ c Figure 7.
The phase diagram of a square lattice Hubbard model with circularopen boundary conditions with an embedded quadratic trap is shown on the left andthe phase diagram of a square lattice Hubbard model with circular open boundaryconditions with embedded experimental trap is shown on the right, both at varyingnearest-neighbor hoppings t (cid:48) t . Calculations are performed using U = 3 and includingnearest-neighbor correlations at T = 10 − and r = 30 for the left plot and T = 0 and r = 35 for the right plot. The y-axis gives the curvature of the trapping potential a trap in the left and the radius of the constant particle number disc in the rightplot. The critical scales Λ c are encoded in the color scheme, whereas the type ofthe diverging phase is given by the shape of the data point. The thermodynamic limitresults are given as a reference in the lower/upper parts of the plots. We encounterantiferromagnetism (AFM), general spin density waves (SDW), ferromagnetism (FM),charge density waves (CDW), superconductivity (SC) and no divergence (N) at lowscales. At low values of − t (cid:48) t the critical scales do not differ much, see the left plot in Fig. 7.Additionally, the emergent phase is matching the expectation. With lowering a trap thecritical scale converges to the same values observed in Fig. 6 for the same system size.The phase diagram is cut-off by the temperature for high − t (cid:48) t values. We observe that trong Boundary and Trap Potential Effects µ i + V i ≈ t (cid:48) . As we cannot follow the RG flow into thesymmetry broken region, we cannot resolve possibly coexisting states at lower energyscales outside of the VHC fulfilling region. All in all, as in experiments a trap cannotbe reduced arbitrarily, it is unlikely to observe bulk d -wave superconductivity in thissimple setup. More sophisticated methods of trapping the ultracold atomic gas must beapplied. We did not recover superconductivity for any reasonable trapping curvature.In contrast, enlarging the VHC fulfilling region by shifting the minima to a ring, leads toa recovery of the superconducting phase which is then bound to this region. Therefore,this offers a possible way to manipulate the spatial dependence of the order parameter,e.g. restricting superconducting order to a ring. To examine the possibility of shaping the ordering by changing the trapping potential,we change the radius of the region fulfilling the VHC by changing the potential to a trap ( r − r V HC ), thereby the VHC fulfilling region is now a ring of radius r = r V HC andnot a point at r = 0 anymore. The slope at the zero crossing is increasing with increasing r V HC , such that the width of the ring in which the VHC approximately holds shrinksto zero for large radii. The two effects counteract each other such that an optimal valuefor r V HC to promote superconductivity exists. We choose T = 0 combined with a sharpcutoff in order to enlarge the system size to r = 35. We apply three different trappingcurvatures, a trap = 4 · − , · − and 4 · − (note that the smallest value is stillone order of magnitude too large to recover the d -wave superconducting phase in thesimple setup above at a next-nearest-neighbor hopping t (cid:48) = − .
28, where we expectsuperconductivity). We vary the VHC fulfilling radius r V HC and track the leadingdivergences of the channels and their changes.For the smallest curvature, we find three phase transitions upon increasing theVHC radius r V HC , see right plot in Fig. 6. The first occurs between r V HC = 0 and r V HC = 4. At r V HC = 0 we obtain a localized spin-divergence of the vertex, which hintsat a localized state emerging due to the potential, similar to what has been observedin [20]. At r V HC = 4 we find a d x − y -divergence which vanishes again upon increasing r V HC further. Between r V HC = 4 and r V HC = 10 we find an antiferromagnet. If wefurther increase the radius we obtain a superconductor, which at first has mainly a d x − y component, and is circular shaped with no clear hole in the center, similar to theupper plots in Fig. 8. At larger radii, the pairing amplitude is reduced in the centerand has main weight on a thin ring at the VHC-radius, this leads to a slight reductionof critical scales. Additionally, the s-wave component is increasing due to the weakerscreening of the boundaries, see the lower plots in Fig. 8. At the largest radii, the trong Boundary and Trap Potential Effects −
20 0 20 x − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − y − . − . . . . P a i r i n g p a tt e r n Figure 8.
Superconducting ordering parameters for two different radii of the
V HC -fulfilling region at U = 3, T = 0 and t (cid:48) = − .
28. The upper row shows r V HC = 10 at a trap = 10 − and the lower row r V HC = 20 at a trap = 4 · − , in the left column the s -wave contribution of the largest eigenvector is shown and in the right column d x − y contribution is shown. s-wave component dominates due to the proximity to the boundary, which leads to aslight increase of the critical scale. All these findings underline the importance of theradius and the shape of the VHC fulfilling region. We did not observe a transition to adisconnected ring, which would require a hole in the gap amplitude at the center. Sucha pure ring has a different topology and might host exotic topological quantum states.Here, the dependence of the critical scale on r V HC is rather smooth, in contrast to thetwo larger curvatures. In those systems we mostly obtain local divergences bound to theVHC fulfilling region. At r V HC = 10 we find a superconducting divergence in all threesetups, where the weight of the s -wave component compared to the d x − y -componentis increasing with the trapping curvature.Remarkably the smooth boundary conditions seem to suppress effects arising dueto open boundaries, see for example the small s -wave component for r V HC = 10 in trong Boundary and Trap Potential Effects d x − y -superconductors with asmooth trap than with a sharp trap. Additionally, we were able to show that we canshape the superconducting phase to a circle. Inspired by this, an experiment could tryto shape the superconducting phase by applications of specifically designed traps. Thiscould result in new types of topological superconductivity or designed superconductingcurrents. For this it would be helpful to create a particle depletion inside the ring, whichis achievable by the application of a Mexican-Hat potential ar − br (with a, b >
0) forwhich the minima is set to fulfill the VHC. To this end, further investigation is neededeither numerically or experimentally.
The last setup we investigate is a square lattice Hubbard model with circular openboundary conditions, in which a trap reconstructed from experiments is embedded.The information given by Mazurenko et. al. [14] allows for an approximate modellingof their trapping potential. It consists of the square lattice part, which we discussed indetail above, and an additional digital micromirror device (DMD) potential to create adisk of constant number density. We use a piece-wise defined function consisting of aconstant disk, a Gaussian which emulates the DMD and a square potential to capturethe lattice potential, see Eq. (14). The parameters are chosen such that the particledensity at half filling as a function of the distance to the center is roughly matching theone given in [14]. a ( x ) = r < r t γ (cid:32) e ( r − rt − ∆)22 σ t − e ∆22 σ t (cid:33) if r t < r < r t + ∆ γ (cid:18) − e ∆22 σ t (cid:19) − b · ( r − r t − ∆) if r > r t + ∆ (14)In total we have five free parameters, the gaussian prefactor γ , the radius of the disc r t , the gaussian regime parameter ∆, the width of the gaussian σ t , and the quadraticprefactor b . The trap is designed continuous but has kinks. In the following we varythe trap radius and choose γ = 85 , σ t = 20 , ∆ = 5 , b = 0 . . The VHC is chosen to be fulfilled at all sites within the trap radius. The differenceto the radius variation in the open boundary system lies in the smooth boundaries weobtain due to the Gaussian. Thus, we expect different finite size effects. For example,we expect a less pronounced s-wave component in the superconducting phase due towhat was observed in the setups investigated above.The finite size effects are again small for weak to intermediate values of t (cid:48) , which isexpected in analogy to the results we obtained for the other models. The critical scale,and thus the critical temperature, is increasing slightly upon increasing the lattice size.At r = 20 and t (cid:48) = 0 we observe that the critical scale is approaching the known FRG trong Boundary and Trap Potential Effects c = 0 .
11 at U = 3 [43]. For increasing values of t (cid:48) the phases are lessstable. We encounter a superconducting divergence at the two smallest system sizes,but they have support on only a few sites, see Fig. 9. The occurrence of these P -channeldivergences seems to be a fine-tuned problem as increasing t (cid:48) or the radius can makethem vanish and reappear. In the r t = 10 model, there is no superconducting phasetransition observed. −
20 0 20 x − − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − − − y − . − . . . . P a i r i n g p a tt e r n −
20 0 20 x − − − y − . − . . . . P a i r i n g p a tt e r n Figure 9.
Leading P -channel eigenvector at t (cid:48) = − .
25 for r t = 5 in the upper,and at t (cid:48) = − .
27 for r t = 20 in the lower row. The left pictures show the on-sitecomponent of the eigenvectors whereas the right pictures visualize the symmetrized d x − y component of it. Calculations are performed using U = 3, r = 35, T = 0applying a sharp-cutoff and including nearest-neighbor correlations. Still, the superconducting phases at the lowest system sizes indeed show theexpected ordering pattern, namely d + s . At higher radii, the s -component is again moresuppressed in the bulk. The ferromagnetic phase is hampered most by the potential andis not fully recovered at the largest system size studied. It does not fill the whole trapregion but manifests itself on a square shaped subsection, this behavior is found in allcases where we applied non-periodic boundary conditions. Thus, it should be observable trong Boundary and Trap Potential Effects s -wave component of the superconducting phase.At the largest radius, we obtain a similar distribution of phases in the t (cid:48) domain asfor the thermodynamic limit case [28], with the exception that the transition betweensuperconductivity and the ferromagnetic phase is shifted to lower values of − t (cid:48) t .The question arises whether a pure bulk d x − y superconductivity is observable atthe largest system size. In Fig. 9 the superconducting order parameters are visualizedat t (cid:48) = − .
23 and r t = 20. We observe that the s -wave component of the eigenvector isvanishing in a small region at the center, such that there, we obtain a local bulk d x − y superconductor. With such a setup we thus might be able to recover the thermodynamiclimit in the laboratory at reasonable effort.
5. Conclusions
We studied the finite-size effects in round- and square-shaped fermionic Hubbard modelson an underlying square lattice, as well as the influences of three different trappingpotentials on the (finite-size) phase diagram of the 2D Hubbard model. We found thatthe application of open boundary conditions leads to a suppression of the ferromagneticregion of the phase diagram which occurs for larger values of the next-nearest-neighborhopping amplitude. The finite-size corrections for the weak next-nearest-neighborhopping are much smaller, which is in line with experiment. The superconductingphase is suppressed at small system sizes and has in general an d + s form inducedby the open boundary conditions, with the somewhat unexpected admixture of onsite s -wave paring to dominant d -wave nearest-neighbor pairing. For the square lattices,the superconducting phase is less suppressed. We observed that designing a specificshape of the trap can tune the divergence encountered, i. e., we were able to create disklike superconductors with main weight at the outer circle by choosing the Van-Hovecondition fulfilling radius appropriately. Applying the trapping potential reconstructedfrom the experiments we observed that the smoothing of the boundaries leads to asmaller suppression of superconductivity. Thus, we observed a stable transition fromantiferromagnetism to superconductivity already at r = 15. The s -wave componentis observed to be less pronounced compared to standard open boundary conditions,thus it is likely that the equivalence between such an experimental system and periodicboundary conditions model is recovered earlier.To obtain a more complete description of the trapped Hubbard models, a next stepis to incorporate the disorder introduced by the DMD which was neglected here butcould lead to further localization effects. The next steps from the theoretical point ofview are the inclusion of self-energy feedback as well as a single frequency dependence ofthe vertex [33, 41, 44]. This would allow for predictions correct up to U . Especially theeffects of the self-energy in open boundary systems could play a crucial role. Predictionsfor the realistic cuprates interaction strength UW ∝ trong Boundary and Trap Potential Effects Acknowledgments
The authors acknowledge the input on how to write and optimize a high-performancecode by E. di Napoli, D. Rohe and S. Achilles. We thank J. Ehrlich, L. Klebl, N. Caci andJ. Beyer for fruitful discussions. The Deutsche Forschungsgemeinschaft (DFG, GermanResearch Foundation) is acknowledged for support through RTG 1995 and underGermany’s Excellence Strategy - Cluster of Excellence Matter and Light for QuantumComputing (ML4Q) EXC 2004/1 - 390534769. We acknowledge support from the MaxPlanck-New York City Center for Non-Equilibrium Quantum Phenomena. Simulationswere performed with computing resources granted by RWTH Aachen University underproject rwth0514.
References [1] Bloch I, Dalibard J and Zwerger W 2008
Reviews of Modern Physics https://link.aps.org/doi/10.1103/RevModPhys.80.885 [2] Tarruell L and Sanchez-Palencia L 2018 Comptes Rendus Physique https://linkinghub.elsevier.com/retrieve/pii/S1631070518300926 [3] Jotzu G, Messer M, Desbuquois R, Lebrat M, Uehlinger T, Greif D and Esslinger T 2014 Nature http://arxiv.org/abs/1406.7874 [4] Gall M, Wurz N, Samland J, Chan C F and K¨ohl M 2021
Nature https://doi.org/10.1038/s41586-020-03058-x [5] Mac´e N, Jagannathan A and Duneau M 2016
Crystals
124 ISSN 2073-4352 arXiv: 1609.08509URL http://arxiv.org/abs/1609.08509 [6] Honerkamp C and Hofstetter W 2004
Physical Review Letters https://link.aps.org/doi/10.1103/PhysRevLett.92.170403 [7] Rapp A, Zar´and G, Honerkamp C and Hofstetter W 2007 Phys. Rev. Lett. (16) 160405 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.160405 [8] Rosch A, Rasch D, Binz B and Vojta M 2008 Physical Review Letters
URL https://doi.org/10.1103/physrevlett.101.265301 [9] Sun K, Liu W V, Hemmerich A and Das Sarma S 2012
Nature Physics [10] B¨uhler A, Lang N, Kraus C, M¨oller G, Huber S and B¨uchler H 2014 Nature Communications [11] J¨ordens R, Strohmaier N, G¨unter K, Moritz H and Esslinger T 2008 Nature [12] Schneider U, Hackermuller L, Will S, Best T, Bloch I, Costi T A, Helmes R W, Rasch D and RoschA 2008
Science https://doi.org/10.1126/science.1165449 [13] Cheuk L W, Nichols M A, Lawrence K R, Okan M, Zhang H and Zwierlein M W 2016
Physical trong Boundary and Trap Potential Effects Review Letters https://link.aps.org/doi/10.1103/PhysRevLett.116.235301 [14] Mazurenko A, Chiu C S, Ji G, Parsons M F, Kan´asz-Nagy M, Schmidt R, Grusdt F, DemlerE, Greif D and Greiner M 2017
Nature [15] Hart R A, Duarte P M, Yang T L, Liu X, Paiva T, Khatami E, Scalettar R T, TrivediN, Huse D A and Hulet R G 2015
Nature [16] Hofstetter W, Cirac J I, Zoller P, Demler E and Lukin M D 2002
Physical Review Letters https://link.aps.org/doi/10.1103/PhysRevLett.89.220407 [17] Rigol M, Muramatsu A, Batrouni G G and Scalettar R T 2003 Physical Review Letters http://arxiv.org/abs/cond-mat/0304028 [18] Scarola V W, Pollet L, Oitmaa J and Troyer M 2009 Physical Review Letters https://link.aps.org/doi/10.1103/PhysRevLett.102.135302 [19] Mendes-Santos T, Paiva T and dos Santos R R 2015
Physical Review A https://link.aps.org/doi/10.1103/PhysRevA.91.023632 [20] Chanda T, Yao R and Zakrzewski J 2020 Phys. Rev. Research (3) 032039 URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.032039 [21] von Delft J 2001 Annalen der Physik https://onlinelibrary.wiley.com/doi/abs/10.1002/1521-3889%28200103%2910%3A3%3C219%3A%3AAID-ANDP219%3E3.0.CO%3B2-I [22] von Delft J, Golubev D S, Tichy W and Zaikin A D 1996 Physical Review Letters http://arxiv.org/abs/cond-mat/9604072 [23] Qin M, Chung C M, Shi H, Vitali E, Hubig C, Schollw¨ock U, White S R, Zhang S and SimonsCollaboration on the Many-Electron Problem 2020 Physical Review X https://link.aps.org/doi/10.1103/PhysRevX.10.031016 [24] Zheng B X, Chung C M, Corboz P, Ehlers G, Qin M P, Noack R M, Shi H, White S R,Zhang S and Chan G K L 2017 Science [25] LeBlanc J, Antipov A E, Becca F, Bulik I W, Chan G K L, Chung C M, Deng Y, Ferrero M,Henderson T M, Jim´enez-Hoyos C A, Kozik E, Liu X W, Millis A J, Prokof’ev N, Qin M,Scuseria G E, Shi H, Svistunov B, Tocchio L F, Tupitsyn I, White S R, Zhang S, Zheng B X,Zhu Z, Gull E and Simons Collaboration on the Many-Electron Problem 2015
Physical ReviewX https://link.aps.org/doi/10.1103/PhysRevX.5.041041 [26] Hauck J B, Honerkamp C, Achilles S and Kennes D M 2020 ArXiv: 2008.13667 URL http://arxiv.org/abs/2008.13667 [27] Metzner W, Salmhofer M, Honerkamp C, Meden V and Schoenhammer K 2012 Reviews of ModernPhysics http://arxiv.org/abs/1105.5289 [28] Lichtenstein J, Pe˜na D S d l, Rohe D, Napoli E D, Honerkamp C and Maier S A 2017 ComputerPhysics Communications
100 – 110 ISSN 0010-4655 URL [29] Markhof L, Sbierski B, Meden V and Karrasch C 2018
Physical Review B http://arxiv.org/abs/1803.00272 [30] Husemann C and Salmhofer M 2009 Phys. Rev. B (19) 195125 URL https://link.aps.org/doi/10.1103/PhysRevB.79.195125 [31] Wang W S, Xiang Y Y, Wang Q H, Wang F, Yang F and Lee D H 2012 Phys. Rev. B (3) 035414URL https://link.aps.org/doi/10.1103/PhysRevB.85.035414 [32] Eckhardt C J, Honerkamp C, Held K and Kauch A 2020 Physical Review B http://arxiv.org/abs/1912.07469 trong Boundary and Trap Potential Effects [33] Weidinger L, Bauer F and von Delft J 2017 Physical Review B http://arxiv.org/abs/1609.07423 [34] Dormand J and Prince P 1980 Journal of Computational and Applied Mathematics https://linkinghub.elsevier.com/retrieve/pii/0771050X80900133 [35] Ahnert K and Mulansky M 2011 arXiv:1110.3397 [physics] http://arxiv.org/abs/1110.3397 [36] Han X J, Liao H J, Xie H D, Huang R Z, Meng Z Y and Xiang T 2017 Chinese Physics Letters http://arxiv.org/abs/1705.10016 [37] Ozaki T 2007 Physical Review B https://link.aps.org/doi/10.1103/PhysRevB.75.035123 [38] Husemann C and Salmhofer M 2009 Physical Review B http://arxiv.org/abs/0812.3824 [39] Honerkamp C, Rohe D, Andergassen S and Enss T 2004 Physical Review B http://arxiv.org/abs/cond-mat/0403633 [40] Klebl L, Kennes D M and Honerkamp C 2020 Physical Review B
ISSN 2469-9969 URL http://dx.doi.org/10.1103/PhysRevB.102.085109 [41] Reckling T and Honerkamp C 2018
Physical Review B http://arxiv.org/abs/1803.08431 [42] Honerkamp C, Wakabayashi K and Sigrist M 2000 Europhysics Letters (EPL) http://dx.doi.org/10.1209/epl/i2000-00280-2 [43] Lichtenstein J 2018 Functional Renormalization Group Studies on Competing Orders in the SquareLattice
Ph.D. thesis RWTH Aachen Aachen[44] Bauer F, Heyder J and von Delft J 2014
Physical Review B https://link.aps.org/doi/10.1103/PhysRevB.89.045128 [45] Wentzell N, Taranto C, Katanin A, Toschi A and Andergassen S 2015 Physical Review B https://link.aps.org/doi/10.1103/PhysRevB.91.045120 [46] Vilardi D, Taranto C and Metzner W 2019 Physical Review B https://link.aps.org/doi/10.1103/PhysRevB.99.104501 [47] Katanin A A 2019 Physical Review B https://link.aps.org/doi/10.1103/PhysRevB.99.115112 [48] Taranto C, Andergassen S, Bauer J, Held K, Katanin A, Metzner W, Rohringer G and Toschi A2014 Physical Review Letters196402 ISSN 0031-9007, 1079-7114 URL