Strong boundary and trap potential effects on emergent physics in ultra-cold fermionic gases

The field of quantum simulations in ultra-cold atomic gases has been remarkably successful. In principle it allows for an exact treatment of a variety of highly relevant lattice models and their emergent phases of matter. But so far there is a lack in the theoretical literature concerning the systematic study of the effects of the trap potential as well as the finite size of the systems, as numerical studies of such non periodic, correlated fermionic lattices models are numerically demanding beyond one dimension. We use the recently introduced real-space truncated unity functional renormalization group to study these boundary and trap effects with a focus on their impact on the superconducting phase of the 2D Hubbard model. We find that in the experiments not only lower temperatures need to be reached compared to current capabilities, but also system size and trap potential shape play a crucial role to simulate emergent phases of matter.


Introduction
Ultra-cold atomic gases provide a powerful method for simulations of quantum manybody lattice systems [1][2][3].Many models introduced in the context of solid-state physics, like the Hubbard model [3], the Haldane model [4], the bilayer Hubbard model [5] or even quasicrystalline models [6], can be investigated in laboratories without the restriction of numerical ambiguities or theoretical approximation errors.Beyond this, other exciting many-body phenomena that cannot be easily explored in solids could be realized in optical lattices as well, see [1,[7][8][9] and references therein, making them a very flexible tool to study quantum many-body physics.
Especially in the context of the fermionic Hubbard model in optical lattices, there have been great advances in recent years, including the observation of the Mott insulating phase [10][11][12] and first sign of antiferromagnetic correlations [13,14].Despite these advances, in the context of correlated fermions state-of-the-art experiments still do not reach beyond the capabilities of exact numerical methods, like quantum Monte Carlo (QMC) or exact diagonalization, as all observed phases can also be investigated quantitatively with the available set of theoretical methods.This may change when the experiments move towards the exploration of d-wave superfluidity [15], where, e.g., QMC encounters the notorious sign problem and finite size scaling from exact diagonalization is hopeless.Despite the lack of quantitative studies of trap or finite size effects in this regard, the current hope is that if the experiments can reach lower temperatures, the superfluid phase will emerge naturally.
However, the effects of the trapping potential have always been an issue in fermionic [16][17][18][19][20][21][22][23][24] and bosonic [25][26][27][28] optical lattices, even for the Mott-phase and high temperature properties of one and two dimensional Hubbard models.In this paper we will focus on fermionic setups.It is not expected that the observed effects are the same or even similar for bosons and fermions as their behavior is already very different for non-interacting particles.Additionally, we concentrate on two-dimensional systems only.In continuation of these previous results, we investigate the effects of the trapping potential, as well as the boundary conditions on the Cooper pairing in a Hubbard model.This is of high experimental importance as one of the current goals is to measure the superfluid phase.A general challenge for Cooper pairing in finite-size systems can, e.g., be gleaned from the literature for superconductivity of ultra-small aluminium grains [29,30], with a diameter of the order of ∝ nm.There, it was worked out clearly how the usual superconducting description is valid as long as the energy gap due to pairing is larger than the energy spacing due to the confinement but breaks down if the latter scale is the larger one.A similar crossover is expected to happen also in a Hubbard model.Even from these basic considerations, the question arises whether there is a minimal number of sites required to unveil superconducting hallmarks.What might add in the case of d-wave pairing is the question, if the generation 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 the non-trivial d-wave ordering pattern will be influenced by the trapping potential.In the antiferromagnetic case the trap was not observed to be a significant perturbation [13].In fact the data compared well to Quantum Monte Carlo results obtained for a square lattice with periodic boundary conditions.Here we want to examine whether this holds true also for the d-wave superconducting phase in a square lattice Hubbard model.To extrapolate to the bulk case (e.g., to connect to real materials), one should also understand how closely the finite-size set-up resembles the situation in the thermodynamic limit.Hence, we will also examine at which size the trapped system approximates the bulk case to high fidelity.i , where r i is the radial distance of site i to the center of the lattice and a trap gives the curvature of the trapping potential.The fifth and last setup is the same as the fourth with a different potential: Instead of the simple quadratic one, we reconstruct a potential used in a recent experiment by Mazurenko et.al. [13].There the trap consists of a Gaussian potential created by a digital micro-mirror device (DMD) with an additional circular quadratic potential, whose superposition creates a nearly flat disc at the center.For simplicity we approximate it by a piece-wise definition, see Eq. (14).The square lattice Hamiltonian in second quantization for all cases reads with the operators c ( †) i,σ annihilating (creating) an electron on site i with spin σ.For the hopping amplitudes we set t ij = t if i and j are nearest-neighbors and t ij = t if i and j are next-nearest-neighbors, where we define nearest and next-nearest neighbors in terms of the distance on the underlying lattice as the closest and second closest sites.All other terms are set to zero.For simplicity we choose t = 1 and measure all quantities in units of t.Distances are measured in units of the lattice spacing a which we set to 1.We choose onsite interactions U ij = U δ i,j and fix the 'Van-Hove condition' µ = 4t .The Van Hove condition ensures that our setup is tuned to a van Hove singularity [31] of the thermodynamic limit Hubbard model, at which the gradient of the energy dispersion is flat and therefore the density of states is large.We neglect density dependent hoppings as it scales with U , which is comparably small in our setup, and n i , which is always below half filling, therefore its effect should be negligible [32][33][34].We use an extended Hubbard model as recent numerical studies [35][36][37] suggest that in the relevant U parameter range for cuprates there is no superconductivity in the pure (t = 0) Hubbard model.

Method
We employ the recently introduced 2D real-space functional renormalization group (realspace TUFRG) [38], which can be used to study models with broken translational symmetry.It is based on the standard one-particle-irreducible FRG [39,40], in a level-2 truncation.Additionally, we will neglect the frequency dependence of the vertex, as well as the self-energy.With this we obtain a method, which scales relatively well and is able to predict finite-size precursors to phase transitions and their respective orderings.The main idea of FRG is the introduction of a cutoff in the bare propagator, enabling an interpolation between a solvable theory and the full solution of the model.Variations with respect to the cutoff parameter generate an infinite set of coupled differential flow equations for the Taylor coefficients of the effective action.These equations must be truncated in order to become numerically tractable, which is well defined in the case of U W < 1, with W being the bandwidth of the model, then only the first few terms in the expansion are relevant and a truncation at the level of keeping just the self-energy and the effective interaction is justified.Explicitly the cutoff is introduced in the bare Green's function by with R(Λ) being the cutoff function obeying R(0) = 1 and lim Λ→∞ R(Λ) = 0 in our convention.This reduces the action in the limit λ → ∞ to the non-interacting, exactly solvable system.By integrating the coupled flow equations, an interpolation between this solution of the non-interacting problem and the approximate full solution is calculated.The variations w.r.t the cutoff parameter generate the single scale propagator defined as in our approximation.
The scale-dependent two-particle vertex can be separated in three different channels according to the three different fermionic bilinears or three different diagram types as shown in Fig. 2. Each of those bilinears is directly associated to a bosonic field of an effective Hamiltonian ‡.The particle-particle channel (P Λ ) is associated with the superconducting pairing interactions and therefore a divergence of this channel at total momentum zero indicates a tendency towards superconductivity.The direct particlehole channel (D Λ ) is associated to charge fluctuations and thus indicates, among others, general charge ordering.The crossed particle-hole channel (C Λ ) can be associated with effective spin-spin interactions, therefore its divergence is an indicator of magnetic ordering.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 as soon as an eigenvalue of one of the three channels surpasses a threshold value.The corresponding RG scale is called the critical scale Λ c , which defines an energy scale of the breakdown of the weakly interacting behavior and can be seen as an estimate of an ordering temperature T c .The leading eigenvalue of the diverging channel gives, by its eigenvector, the expected ordering pattern up to an unknown prefactor.A more precise determination of the order would follow from the gap equations for each of the channels.In practice we will normalize the shown ordering pattern from the leading eigenvector such that their maximal absolute value is 1.As we consider finite-size systems here, the interpretation of the flows to strong coupling as phase transitions is not strictly valid, as neither a real discontinuity occurs, nor any correlation length truly diverges.Nevertheless, it can be expected that the finite-size system with a strongly enhanced interaction of a specific type will locally resemble a infinite system with a truly diverging effective interaction of the same type.Thus, we use the same vocabulary as for the infinite system and discuss different phases of the finite system as regimes with qualitatively different flows to strong coupling.
As we do not include self-energy effects, we cannot incorporate possible gap openings, which would allow us to continue the flow into the symmetry broken region [41][42][43][44].Therefore, lower lying divergences in other locations of the trap cannot be ‡ In fact these bilinears are interpreted most easily for the particle-particle (P ) and the crossed particlehole (C=spin) channels, the direct particle-hole channel (D) is more complex as it contains information on the C-channel.The physical relevant quantity is the charge channels given by V ch = 2D − C, which removes the redundant information.
. Diagrammatic representation of the FRG flow equation for the one-particle irreducible interaction vertex in the level-II truncation, adapted from [39].The first term on the right-hand side, annotated with Ṗ Λ , represents the flow of the particleparticle channel P Λ , while the other terms generate the flow of C Λ and D Λ .The dashed internal lines represent a single scale propagator, see Eq. ( 3), and the solid lines represent a propagator, see Eq. ( 2).The arrows indicate the direction of the propagators and fix the order of their arguments.The grey boxes represent the full two-particle interaction.A second class of diagrams with the scale-derivative on the other internal line contributes as well but is now drawn for simplicity.
resolved within our approach.In the case that there are multiple eigenvalues diverging at once, we can access the sub-leading eigenvalues and could perform a post-processing mean-field decoupling in order to obtain the expected gap pattern, thereby, resolving cases in which the dominant eigenvector does not contribute to the order parameter.
The results presented here showed in all cases a clear dominant eigenvector instead of divergences of multiple ones.Therefore the above described procedure was never needed.Next, we discuss the spatial dependence of the effective interactions and the respective channels P , C and D. In principle, each channel can depend on four site indices.Yet, the initial condition in Eq. ( 1) is an onsite interaction, i.e.where all four site 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, 'bilocal' correction that depends on the joint index of the incoming lines i and the outgoing lines j.Like in the P -channel, these RPA-like contributions create a main dependence on different bi-local pairs of indices in each channel.Contributions beyond these bi-local indices are only generated by the feedback in between the channels and, thus, come at a higher order in the interaction.The further apart from the native bi-local indices, the higher is the interaction order [45].Therefore, these terms can be neglected in the spirit of the perturbative motivation of the FRG.The bi-local contributions of each channel on the other hand have to be captured more rigorously, as they are equivalent to the sharp . Generation of new index dependencies during the flow.A density-density interaction U i,i n i,s n i,s receives corrections in the P -channel that can be interpreted as an incoming particle pair c i,s c i,s scattering to an outgoing pair c j,s c j,s .
Formally, this specific structure of the effective interactions can be exploited by an expansion of the subleading dependencies in form factors, which form a complete orthonormal basis set, according to The three channels can then be build up from form factors capturing the internal structure of the bilinears that in turn can mutually interact on longer distances With this at hand we can rewrite the flow equations by an insertion of form-factor unit matrices.The main idea of the truncated unity approach is to restrict the number of form-factors in these unities to much less than the lattice sites [38,40,48] b k ∈L where L is the set of all allowed bonds and analogously In the form-factor space the flow equations can be rewritten in terms of highly efficient block-matrix products.Due to the truncation, the full vertex cannot be recovered exactly but instead projections must be introduced in order to reconstruct the full vertex approximately for the cross-channel feedback.These projections can be found for example in Ref. [38,45,49].Instead of including the full spin dependence we will assume an SU (2) invariant problem.This is by no means necessary but simplifies the 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 with the particle-hole bubble χph and the particle-particle bubble χpp defined as For the discussion of the results we will often refer to the s-wave and d x 2 −y 2wave component of the leading eigenvector, which we define as follows: The leading eigenvector of a channel depends on a single site and a bond index, thus it can be written as v b i i .The s-wave, or onsite component is in this notation defined as v b 0 i , where b 0 refers to the bond with zero length.The definition of the d x 2 −y 2 component follows analogously to momentum space [40], we enumerate next-nearest-neighbor bonds from 1 to 4 in a clockwise fashion starting from the up pointing bond, then we define v

Numerical Implementation
For the integration, an adaptive fifth order Runge-Kutta scheme [50] of the odeint library [51] is used.The matrix products on the right hand side of Eq. ( 9), which are the numerically most demanding parts, are performed on GPU's.To reduce the number of Matsubara frequencies needed for the calculation of the propagator, we use the Padé-approximation of the Fermi-function [52,53].The results are then verified using a heuristic tangent spacing scheme.Both yield fast convergence of the Bubble terms tested against the analytical result.In the calculations we incorporate 300 frequencies, which is equivalent to an error of the bubble calculation of about 10 −8 at the lowest used temperature (T = 10 −3 ) in comparison to the analytic result.We use a mapping of the infinite integration range to a finite one Λ = 1−x x to accelerate the integration [49].The integration is started from a minimum x = 10 −3 as smaller starting values yield no significant difference in the final results.To reduce the load imbalance between the channels we perform a completion of the square for the D-channel.Thereby we only need 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 choose so-called Ω-cutoff, first introduced by Husemann et.al. [54], given by It has the advantage of a real IR-divergence regularization but is numerically more challenging as for example the interaction cutoff [55].Additionally, no numerical instabilities arise if the Hamiltonian has a zero eigenvalue.At zero temperature we employ a sharp cutoff, which drastically simplifies the bubble calculations [45,56].In total, the flow equations scale like O(N 3 • Nb 3 ) (N is the number of orbitals and Nb is the average number of bonds per site) due to the matrix products, the bubble calculation is the number of Matsubara frequencies included in the summation).Upon including the self-energy as well as the frequency dependence of the vertex (in the single channel coupling or ECLA sense [49,57]) the error scales with U 3 , this scaling has been checked for the 1D-Hubbard model comparing to exact diagonalization.Additionally, we verified that our implementation is sufficiently efficient to reach the thermodynamic limit of the 2D-Hubbard model, where the known FRG results are reproduced.

Open boundary conditions -finite size effects
We start by an analysis of the finite size effects in a square lattice Hubbard model with open boundary conditions and without any confining potentials.The Van-Hove condition is fixed by setting µ = 4t and we vary t at different system sizes to find at which size the critical scales, as well as the occurring phases, are converged.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 do vary only very slightly close to half filling.We find a superconducting phase even for the smallest system.The superconducting transition gets shifted to higher values of − t t for larger systems, the largest system investigated shows a phase diagram roughly matching the thermodynamic limit case.It should be noted that the occurring superconducting phase is not a pure d-wave divergence, as one would obtain with periodic boundary conditions, but instead we have mixing between s-and d-wave bond pairing at all sizes at the boundary.Such an s-and d-mixing behavior has been a topic at the [110] surface of square lattice systems (e.g.[58]), but here it is more profound as the onsite s-wave coupling is not put in à priori, but generated by the RG flow.Note that the initial interaction is a repulsive onsite interaction, but the RG flow generates, as a seemingly diverging tendency, an interaction between pair bilinears that include noticeable s-wave onsite terms, next to dominant nearest-neighbor pairs with d-wave pattern.Currently it is not clear which conditions, besides spatial inhomogeneity, are required to obtain this accompanying s-wave component.In our approach, the phase between the two components is fixed to unity as they arise in a single eigenvector.The phase could change in the gap equation as this is a non-linear relation that could alter the onsite or bond pairings independently.Here, we will stick to the notion of d + s mixing without a phase, indicating the dominant d-character.For a proper verification of this finding, the possible gap opening should be investigated along the flow, which is numerically 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 increasing lattice size.In the largest system, the s-wave component has smaller amplitude than the d-wave component.It has additionally a stronger site dependence, i. e. , it shows its highest values at the boundaries.The superconducting phase does not fully resemble the thermodynamic limit case due to the visibility of boundary effects in the whole bulk even for l = 60.In summary, we observe a clear tendency towards an increased relevance of finite size effects with increasing next-nearest neighbor hopping.Additionally, the superconducting and ferromagnetic ordering vectors we found to deviate substantially from the thermodynamic limit.
For the experiments, round trapping potentials are easier to generate than sharp square trapping potentials [5,13,14].Therefore, the second setup under investigation is a square lattice Hubbard model with circular open boundary conditions and without any confining potential, 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 plot in Fig. 4. There, no ambiguity in the resulting phase occurs, an antiferromagnetic divergence is found at all sizes.Even the ordering pattern is not influenced much by varying the system size.This is in agreement with the observations by Mazurenko et.al. [13].Upon increasing t the critical scales are decreasing.The finite size effects become more pronounced with increasing t .At t ≈ −0.26 we find a transition to a pairing divergence for the two largest systems.For these two, the critical scales are comparable to the ones in the thermodynamic limit.We only encounter superconductivity for a radius larger than r = 20, additionally the phase transition between AFM and SC is shifted towards higher values of − t t .This already hints at possible complications for the experimental setups.In the round case we are not able to reach the thermodynamic limit, even for r = 30, as the phase boundaries do not match the ones from the thermodynamic limit even at the largest radius.As for the square lattice case, no pure bulk d x 2 −y 2 superconductor is recovered, as again in the complete bulk a mixing with the s-wave component is present.
We observe a strong size dependence for high − t t values where, in the infinite model, the ferromagnetic phase occurs.This can be explained by considering the orderings, see Fig. 5; Whereas the antiferromagnetic and superconducting phases more or less comply with the circular shape of the lattice (up to boundary effects and the already mentioned mixing), the ferromagnetic phase does not.In fact, the resulting magnetization pattern has a square like shape for the largest size.This of course is not matching to the lattice constraints and therefore ferromagnetic order is suppressed by boundary effects.Before the formation of the square like pattern, it consists of five localized peaks, one in the center and the others at the edges of the square, at increasing sizes they become more and more connected and fuse ultimately to the square-like shape shown.Changing from square to circular open boundary conditions significantly deteriorates the convergence of the phase diagram with respect to system size, but the way the orderings change does not differ significantly.

Effects of the trap potential
We now turn to the third setup, a square lattice Hubbard model with open boundary conditions in which a second box potential with finite height is embedded.This is realized by an on-site potential at the two sites closest to the boundary at the upper and the right edge and at the three sites closest to the boundary at the lower and left edge.By this we can verify that the boundaries existing beyond the box potential do 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 forward open boundaries is that due to the finite height, the electronic wavefunctions are allowed to be non-zero at the boundary.Inside the potential well, the wavefunction decays exponentially.In this case, we are mainly interested in the superconducting phase, thus we choose −t = 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.At the smallest system size, we find a charge divergence.Starting from an effective system length of l = 12 we observe a P -channel divergence.Again the s-wave component is dominant at small system sizes and the divergences still show similar patterns as for the 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 reduced the finite size effects.Thereby, such a potential could reduce the system size necessary to resolve bulk d x 2 −y 2 -superconductivity.Such a square box trap potential is of course hard to realize experimentally, therefore we next proceed with a potential which is closer to experiments.
In experiments, the lattice is generated by a modulated square potential [3], therefore the fourth case investigated is a square lattice Hubbard model with circular open boundary conditions, in which we insert a quadratic trap as on-site potential: The coordinate system is arranged such that its origin coincides with the center of the lattice.The trapping is embedded in a round lattice with r = 30 and the temperature is fixed to T = 10 −3 .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 examine whether superconductivity can still be observed.At low values of − t 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 the critical 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 t values.We observe that the parameter regime in which we expect superconductivity is dominated by spinwave ordering, which is in line with earlier investigations [20].The superconducting phase vanishes due to the confinement of the effective system size in which the Van-Hove condition (VHC) is fulfilled.The emergent spin divergences are pinned to only t .Calculations are performed using U = 3 and including nearest-neighbor correlations at T = 10 −3 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 right plot.The critical scales Λ c are encoded in the color scheme, whereas the type of the diverging phase is given by the shape of the data point.The thermodynamic limit results are given as a reference in the lower/upper parts of the plots.We encounter antiferromagnetism (AFM), general spin density waves (SDW), ferromagnetism (FM), charge density waves (CDW), superconductivity (SC) and no divergence (N) at low scales.
a very small region close to the center of the lattice, the VHC fulfilling region for which µ i + V i ≈ 4t .As we cannot follow the RG flow into the symmetry broken region, we cannot resolve possibly coexisting states at lower energy scales outside of the VHC fulfilling region.All in all, as in experiments a trap cannot be reduced arbitrarily, it is unlikely to observe bulk d-wave superconductivity in this simple setup.More sophisticated methods of trapping the ultracold atomic gas must be applied.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 to a 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.Here we observed first, that a simple trap smoothens the open boundary conditions, reducing the finite size effects and second that with vanishing size of the VHC-fulfilling region the superconducting phase vanishes too.

Trap-shaped superconductivity
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 2 − r 2 V HC ), thereby the VHC fulfilling region is now a ring of radius r = r V HC and not 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 shrinks to zero for large radii.The two effects counteract each other such that an optimal value for r V HC to promote superconductivity exists.The emergence of ring like pairing correlations were observed earlier [19] in a DQMC study, we now want to study its radius dependence and stability w.r.t.variations of the trapping curvature a trap .We choose T = 0 combined with a sharp cutoff in order to enlarge the system size to r = 35.We apply three different trapping curvatures, a trap = 4 • 10 −4 , 1 • 10 −3 and 4 • 10 −3 (note that the smallest value is still one order of magnitude too large to recover the d-wave superconducting phase in the simple setup above at a next-nearest-neighbor hopping t = −0.28,where we expect superconductivity).We vary the VHC fulfilling radius r V HC and track the leading divergences of the channels and their changes.
For the smallest curvature, we find three phase transitions upon increasing the VHC 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 hints at a localized state emerging due to the potential, similar to what has been observed in [24].At r V HC = 4 we find a d x 2 −y 2 -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 we further increase the radius we obtain a superconductor, which at first has mainly a d x 2 −y 2 component, and is circular shaped with no clear hole in the center, similar to the upper plots in Fig. 8.At larger radii, the pairing amplitude is reduced in the center and has main weight on a thin ring at the VHC-radius, this leads to a slight reduction of critical scales.Additionally, the s-wave component is increasing due to the weaker screening of the boundaries, see the lower plots in Fig. 8.At the largest radii, the s-wave component dominates due to the proximity to the boundary, which leads to a slight increase of the critical scale.All these findings underline the importance of the radius and the shape of the VHC fulfilling region.We did not observe a transition to a disconnected ring, which would require a hole in the gap amplitude at the center.Such a 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 the two larger curvatures.In those systems we mostly obtain local divergences bound to the VHC fulfilling region.At r V HC = 10 we find a superconducting divergence in all three setups, where the weight of the s-wave component compared to the d x 2 −y 2 -component is increasing with the trapping curvature.
Remarkably the smooth boundary conditions seem to suppress effects arising due to open boundaries, see for example the small s-wave component for r V HC = 10 in Fig. 8. Thus it could become easier to observe a bulk d x 2 −y 2 -superconductors with a smooth trap than with a sharp trap.Additionally, we were able to show that we can shape the superconducting phase to a circle.Inspired by this, an experiment could try to shape the superconducting phase by applications of specifically designed traps.This could result in new types of topological superconductivity or designed superconducting currents.For this it would be helpful to create a particle depletion inside the ring, which is achievable by the application of a Mexican-Hat potential ar 4 − br 2 (with a, b > 0) for which the minima is set to fulfill the VHC.To this end, further investigation is needed either numerically or experimentally.In general, it will be very interesting to investigate the possible uses of trapping potentials for experimental realisations of exotic phases of matter.

Experimental Trap
The last setup we investigate is a square lattice Hubbard model with circular open boundary conditions, in which a trap reconstructed from experiments is embedded.The information given by Mazurenko et.al. [13] allows for an approximate modelling of their trapping potential.It consists of the square lattice part, which we discussed in detail above, and an additional digital micromirror device (DMD) potential to create a disk of constant number density.We use a piece-wise defined function consisting of a constant disk, a Gaussian which emulates the DMD and a square potential to capture the lattice potential, see Eq. ( 14).The parameters are chosen such that the particle density at half filling as a function of the distance to the center is roughly matching the one given in [13].
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 quadratic prefactor b.The trap is designed continuous but has kinks.In the following we vary the trap radius and choose γ = 85, σ t = 20, ∆ = 5, b = 0.01666.
The VHC is chosen to be fulfilled at all sites within the trap radius.The difference to the radius variation in the open boundary system lies in the smooth boundaries we obtain 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 to what was observed in the setups investigated above.
The finite size effects are again small for weak to intermediate values of t , which is expected in analogy to the results we obtained for the other models as well as earlier studies employing similar trappings [20].The critical scale, and thus the critical temperature, is increasing slightly upon increasing the lattice size.At r = 20 and t = −0.1 we observe that the critical scale is approaching the known FRG critical scale of Λ c = 0.11 at U = 3 [59].For increasing values of t the phases are less stable.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 -channel divergences seems to be a fine-tuned problem as increasing t or the radius can make them vanish and reappear.In the r t = 10 model, there is no superconducting phase transition observed.
Still, the superconducting phases at the lowest system sizes indeed show the expected ordering pattern, namely d+s.At higher radii, the s-component is again more suppressed in the bulk.The ferromagnetic phase is hampered most by the potential and is not fully recovered at the largest system size studied.It does not fill the whole trap region but manifests itself on a square shaped subsection, this behavior is found in all cases where we applied non-periodic boundary conditions.Thus, it should be observable in experiments.We observe that the trapping reduces the effects of the open boundary conditions especially concerning the s-wave component of the superconducting phase.At the largest radius, we obtain a similar distribution of phases in the t domain as for the thermodynamic limit case [40], with the exception that the transition between superconductivity and the ferromagnetic phase is shifted to lower values of − t t .The question arises whether a pure bulk d x 2 −y 2 superconductivity is observable at the largest system size.In Fig. 9 the superconducting order parameters are visualized at t = −0.23 and r t = 20.We observe that the s-wave component of the eigenvector is vanishing in a small region at the center, such that there, we obtain a local bulk d x 2 −y 2 superconductor.With such a setup we thus might be able to recover the thermodynamic limit in the laboratory at reasonable effort.
Comparing the present setup to the simpler trap shapes studied above, we observe that the two ingredients most important for reaching a fast convergence in system size in the studied Hubbard model are large sizes of the VHC fulfilling regions in combination with smooth trapping potentials.An optimal design of the latter will certainly help to reduce the effective system size required in the future.

Conclusions
We studied the finite-size effects in round-and square-shaped fermionic Hubbard models on an underlying square lattice, as well as the influences of three different trapping potentials on the (finite-size) phase diagram of the 2D Hubbard model.We found that the application of open boundary conditions leads to a suppression of the ferromagnetic region of the phase diagram which occurs for larger values of the next-nearest-neighbor hopping amplitude.The finite-size corrections for the weak next-nearest-neighbor hopping are much smaller, which is in line with experiment.The superconducting phase is suppressed at small system sizes and has in general an d + s form induced by 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 specific shape of the trap can tune the divergence encountered, i. e., we were able to create disk like superconductors with main weight at the outer circle by choosing the Van-Hove condition fulfilling radius appropriately.Applying the trapping potential reconstructed from the experiments we observed that the smoothing of the boundaries leads to a smaller suppression of superconductivity.Thus, we observed a stable transition from antiferromagnetism to superconductivity already at r = 15.The s-wave component is observed to be less pronounced compared to standard open boundary conditions, thus it is likely that the equivalence between such an experimental system and periodic boundary conditions model is recovered earlier.
To obtain a more complete description of the trapped Hubbard models, a next step is to incorporate the disorder introduced by the DMD which was neglected here but could lead to further localization effects.The next steps from the theoretical point of view are the inclusion of self-energy feedback as well as a single frequency dependence of the vertex [49,57,60].This would allow for predictions correct up to U 3 .Especially the effects of the self-energy in open boundary systems could play a crucial role.Predictions for the realistic cuprates interaction strength U W ∝ 1 are not reliable with the realspace TUFRG due to its perturbative motivation.This issue could be resolved by the combination of dynamic mean-field theory with real-space TUFRG [61][62][63][64].
Our study illustrates that on the experimental side, systems sizes need to be sufficiently large, in addition to reaching a lower temperature.However, we illustrate that the trapping should not only be understood as a challenge but can also provide an opportunity to freely shape emergent phases to our needs.
We will simulate five different setups, all are different realizations of the square lattice Hubbard model with different open boundary conditions.They differ in the lattice confining geometry, which is either a square box or a circle, and trapping potentials, the different setups are visualized in Fig. 1.The first setup investigated is a square lattice Hubbard model with open boundary conditions and without any confining potentials.The second one is a square lattice Hubbard model with circular open boundary conditions and without any confining potential.These two cases are investigated to distinguish between the effects introduced by the open boundary conditions applied in all cases and the ones introduced by the different lattice shapes.The third case is a square lattice Hubbard model with open boundary conditions in which a second box potential with finite height is embedded.Thereby the open boundary conditions are smoothed.As the fourth setup we consider a square lattice Hubbard model with circular open boundary conditions, in which a quadratic trap is inserted, with V trap i = a trap r 2

Figure 1 .
Figure 1.Schematic visualization of the five different geometries and trapping combinations considered here.In the upper left, the simple square lattice with open boundary conditions is visualized and in the upper middle plot, the circular open boundary conditions case is shown.The upper right plot shows the combination of open boundary conditions with a finite height trapping.The lower plots visualize the two 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 experimental trap.The color scale gives the values of the trapping potential V i .

Figure 4 .
Figure 4.The phase diagram of the square lattice Hubbard model with open boundary conditions and without any confining potentials is shown on the left and the phase diagram of a square lattice Hubbard model with circular open boundary conditions and without any confining potential is shown on the right.Simulations are performed at van Hove filling at varying system sizes and next-nearest-neighbor hopping tt .We choose U = 3, T = 10 −3 and include nearest-neighbor correlations.The critical scales Λ c are encoded in the color scheme, whereas the type of the diverging phase is given by the shape of the data point.The thermodynamic limit results are given as a 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).

Figure 5 .
Figure 5. Examples for the leading eigenvectors for each phase in the square lattice Hubbard model with circular open boundary conditions and without any confining potential at r = 30.Calculations are performed using U = 3, T = 10 −3 and including nearest-neighbor correlations.In the upper row the leading eigenvectors of the Cchannel are shown at t = 0 in the left, which is an antiferromagnet, and t = −0.4 in the right image, forming a bulk square shaped ferromagnet.The lower row shows the s-wave (left) and symmetrized d x 2 −y 2 (right) contribution of the leading eigenvector of the P -channel at t = −0.26.The P -channel eigenvector indicates d + s ordering tendencies.

Figure 6 .
Figure 6.We show the phases of a square lattice Hubbard model with open boundary conditions 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 show the 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 resulting phase.In the insets in the left plot, the s and d x 2 −y 2 component of the leading pairing eigenvector at the largest system size are shown.Calculations are performed using U = 3, T = 0 and including nearest-neighbor correlations at t = −0.28.We encounter antiferromagnetism (AFM), charge density waves (CDW), superconductivity (SC) and no divergence (N) at low scales.

Figure 7 .
Figure 7.The phase diagram of a square lattice Hubbard model with circular open boundary conditions with an embedded quadratic trap is shown on the left and the phase diagram of a square lattice Hubbard model with circular open boundary conditions with embedded experimental trap is shown on the right, both at varying nearest-neighbor hoppings t t .Calculations are performed using U = 3 and including nearest-neighbor correlations at T = 10 −3 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 right plot.The critical scales Λ c are encoded in the color scheme, whereas the type of the diverging phase is given by the shape of the data point.The thermodynamic limit results are given as a reference in the lower/upper parts of the plots.We encounter antiferromagnetism (AFM), general spin density waves (SDW), ferromagnetism (FM), charge density waves (CDW), superconductivity (SC) and no divergence (N) at low scales.

Figure 8 .
Figure 8. Superconducting ordering parameters for two different radii of the V HCfulfilling region at U = 3, T = 0 and t = −0.28.The upper row shows r V HC = 10 at a trap = 10 −3 and the lower row r V HC = 20 at a trap = 4 • 10 −4 , in the left column the s-wave contribution of the largest eigenvector is shown and in the right column d x 2 −y 2 contribution is shown.

Figure 9 .
Figure 9.Leading P -channel eigenvector at t = −0.25 for r t = 5 in the upper, and at t = −0.27for r t = 20 in the lower row.The left pictures show the on-site component of the eigenvectors whereas the right pictures visualize the symmetrized d x 2 −y 2 component of it.Calculations are performed using U = 3, r = 35, T = 0 applying a sharp-cutoff and including nearest-neighbor correlations.