Sanjoy Datta
email: sanjoy.datta@grenoble.cnrs.fr

Viveka Nand Singh

Pinaki Majumdar
email: pinaki@hri.res.in

Radio frequency spectroscopy of the attractive Hubbard model in a trap

Attractive interaction between fermions can lead to pairing and superfluidity in an optical lattice. In contrast to the 'continuum', on a lattice the trap induced density variation can generate a non monotonic profile of the pairing amplitude, and completely modify the spectral signatures of any possible pseudogap phase. Using a tool that fully captures the inhomogeneity and strong thermal fluctuations, we demonstrate how the crucial radio frequency signatures of pairing are 'inverted' in a trapped attractive fermion lattice compared to the traditional continuum case. These features would be central in interpreting any spectroscopic hint of fermion pairing and superfluidity.

Optical lattices allow controllable cold atom realisation [1][2][3][4] of interacting quantum lattice models. The achievements include the observation of a Fermi surface [5] and Mott insulating phase [6,7] for repulsive fermions, and the evidence of superfluidity (SF) [8] and anomalous expansion [9] in the attractive case. While the canonical antiferromagnetic state [10,11] of repulsive fermions and superfluidity in the attractive Hubbard model [12] (AHM) remain inaccessible, the observation of precursors to these states would already be a major advance.

Even if a pairing induced gapped, or pseudogap (PG), phase is thermally accessible, the spectroscopic signatures would be hard to interpret. The well developed theory of pairing in the 'flat' AHM [13][14][15] provides no obvious guidance on the angle resolved spectrum of the trapped lattice. The complication has a simple origin. Trapping potentials lead to a monotonic increase in density, as one moves from the edge to the center of the trap, but the pairing amplitude variation becomes non monotonic once the central density crosses unity. The non monotonicity affects the spatial character of excitations, and generates a spectroscopic response differing drastically from the famed 'backbending' that one observes in the flat lattice or the trapped continuum gas [18,19].

We completely solve this problem, using a Monte Carlo (MC) method that handles both the inhomogeneity and thermal fluctuation on large lattices. We predict the following: (i) Increasing confinement leads to rapid decrease in the overall spectral gap, pushing weight to low frequency, and quick suppression of the coherence feature at the gap edge. (ii) Radio frequency spectroscopy (RFS), the cold atom analog of angle resolved photoemission spectroscopy (ARPES), shows 'backbending', the traditional signature of a pairing gap, only for weak trapping and low temperature, with the momentum dependent gap smallest near k ∼ k F ∼ {π/2, π/2}. For stronger confinement, however, this inverts to 'forward bending' with the gap largest near k ∼ k F , despite the presence of strong pairing. (iii) This 'inversion' is generic, and arises when the density at the trap center exceeds 1. It survives beyond T c , but vanishes for T T c .

We provide an analysis in terms of the quasiparticle states in the trap, and demonstrate an approximate "local density" approach that captures most of the MC based features and can yield reliable RF spectra on very large, experimentally relevant, lattices.

Model and method: We study the two dimensional (2D) attractive Hubbard model in the presence of a harmonic potential:

H = H 0 -|U | i n i↑ n i↓ , where H 0 = -t ij σ c † iσ c jσ + iσ (V i -µ)n iσ .
The first term denotes the nearest neighbour tunneling amplitude of atoms on the optical lattice, the confining potential has form

V i = V 0 (x 2 i + y 2 i )
, µ is the chemical potential, and U > 0 is the strength of attractive on-site interaction. x i and y i are measured in units of lattice spacing a 0 . On a L × L lattice, the corner value Beyond the small window at weak Vc the system has only SF order at low temperature, with an increasing BI core for Vc/t > ∼ 4. N refers to the normal state. The Tc (on left axis) falls monotonically with Vc as does ∆0 (right axis), the T = 0 gap in the spectrum. The Tc at Vc = 0 would vanish in the infinite volume limit, the results here are for a 24 × 24 system. arXiv:1312.5761v2 [cond-mat.str-el] 24 Dec 2013

V c = V {L/2,L/2} = V 0 * 2 * (L/2) 2 . We use L = 24.
The spatial variation in mean value, and the thermal fluctuation about the mean pairing amplitude are crucial in describing the physics of this system. Unbiased calculations in the homogeneous limit employ determinantal quantum Monte Carlo (DQMC) [13][14][15] to access finite temperature properties. While there are a few recent calculations using large system size [11,16,17], they are focused on thermodynamic properties and have not touched upon the spectral functions of the AHM.

We use a strategy used earlier on moderately sized systems [20,21], augmented by a cluster Monte Carlo technique [22] that readily allows access to system size ∼ 30×30. We first derive an effective Hamiltonian by decoupling the interaction term simultaneously [23] in the pairing and density channels via a Hubbard-Stratonovich (HS) transformation. The exact transformation puts a constraint on the coupling constants in these two channels [24]. We choose both couplings to be unity, and neglect the time dependence of the auxiliary fields, to reproduce Hartree-Fock-Bogoliubov-de Gennes (HFBdG) theory at T = 0. Our model is:

H ef f = H 0 +H coup +H stif f , where H coup = i (∆ i c † i↑ c † i↓ + ∆ i c i↓ c i↑ ) -i φ i n i , and
H stif f = 1 U i (|∆ i | 2 + φ 2 i ). ∆ i = |∆ i |e iθi is
a complex scalar and φ i is a real scalar field. The inclusion of φ i is essential to capture the Hartree shift in the inhomogeneous system. The T = 0 state corresponds to solving δE/δ∆ i = 0 and δE/δφ i = 0, where E is the energy in the {∆, φ} background, and reproduces mean field theory [25]. Finite temperature configurations {∆ i , φ i } follow the distribution P {∆ i , φ i } ∝ T r c,c † e -βH ef f and may fluctuate significantly from the mean field state.

We use the Metropolis algorithm to update the |∆|, θ and φ variables. This involves solution of the HFBdG equation [25,[START_REF] De Gennes | Superconductivity of metals and alloys[END_REF] for each attempted update, to compute the fermion trace. For determining the acceptance of a move we solve the HFBdG equation on a 8 × 8 cluster around the update site. Global properties like pairing field correlation, density of states, etc, are computed via solution of the HFBdG equation on the full 24×24 system in equilibrium {∆ i , φ i } configurations. We have checked (see Supplement) that our T c matches the DQMC estimate [15] over a wide U/t window.

The parameter space of the trap problem involves U/t, V c /t, average density n av , and temperature T /t. To keep the effort manageable we set U/t = 6, where the T c in the flat system is maximum. We have explored the variation from weak to strong confinement over a wide density window but will show detailed results mainly at n av = 1.

For V 0 = 0 the model is known [13,14] to have a SF ground state for 0 < n < 2, except at n = 1 where there is coexistence of SF and DW correlations. For n = 1 the SF has Bardeen-Cooper-Schrieffer (BCS) character at U/t

n av = N f /L 2 ,
where N f is the number of fermions, and corner potential V c . At finite V there is a small window near n av = 1 where DW correlations survive, upto V c /t ∼ 0.8 [START_REF] Datta | [END_REF]. Beyond this window the system has only SF order. However, the spatial extent of the SF shrinks with increasing V c or n av since the central part of the trap becomes doubly occupied (n i = 2) suppressing ∆ i . There is a narrow SF+DW window at small V c , beyond which there is only SF order, with the T c (left axis) decreasing quickly with increasing confinement. The T = 0 spectral gap ∆ 0 (right axis) falls even more sharply, dropping from ∼ 4.6t at

V c = 0 to ∼ 1.5t at V c = 6t.
Fig. 2 shows the radial variation of the thermal average of

n i = σ c † iσ c iσ (left), |∆ i | (center), and Φ i = |∆ i ||∆ i+δ |cos(θ i -θ i+δ ) (right). The coordinate i is r =
x 2 + y 2 /(L/2), varying along the diagonal. Φ i tracks nearest neighbour correlation in that direction. We have set V c = 3t and n av ∼ 1 and T = 0, 0.08t, 0.3t. The full 2D spatial maps are shown in the Supplement.

Fig. 2(a) shows the expected monotonic fall in n r at all T . The cloud at T = 0.3t is slightly broader than at T = 0. The pairing field amplitude in 2(b) is more interesting. It is non-monotonic at all T , a peculiarity of the lattice where it grows with n till n = 1 and falls beyond. The T = 0 result for |∆ r | is what is expected from mean field HFBdG theory, with a clear peak in the region where n r ∼ 1. At T = 0.08t the amplitude profile looks similar to T = 0, but with a large growth in the corner where it was zero at T = 0! The trend amplifies at T = 0.3t where |∆ r | is much less inhomogeneous than at T = 0. This is due to the low amplitude stiffness in regions with low |∆ i | at T = 0. We provide a connection to the flat system physics in the Supplement.

Fig. 2(c) is meant to highlight the suppression of phase correlation with temperature. At T = 0 the phases are locked, so

Φ i = |∆ i ||∆ i+δ | ≈ |∆ i | 2 . At T = 0.08t ∼ 0.
7T c while the amplitudes are not very different from T = 0 the phase correlation is weakened. By T = 0.3t while amplitudes have grown, NN phase correlations have weakened to about 20% of the T = 0 value. Long range phase correlation is of course lost at T c . The spatial characteristics are not directly accessible so we move to the spectral signatures that RF spectroscopy can probe.

V c =0.1t V c =1.0t V c =3.0t (a) -4 -2 0 2 4 ω 0 0.1 0.2 V c =0.1t V c =1.0t V c =3.0t (b) -4 -2 0 2 4 ω 0 0.1 V c =0.1t V c =1.0t V c =3.0t
Fig. 3 shows the single particle density of states (DOS). Fig. 3.(a) shows V c dependence at T = 0. There are two primary effects of trapping: (i) the effective gap reduces with increasing V c due to appearance of low frequency spectral weight, and (ii) the 'coherence peak' and sharp gap edge are blurred. The decrease in the gap arises from the smaller pairing amplitude in regions which have density n i → 0 or n i → 2. We have explicitly checked this from the local density of states (LDOS). In fact at n ∼ 1.9 the pairing gap in the flat system is 0.8t, not very different from the threshold that we observe. The n i ∼ 1 region contributes to spectral weight at |ω| > ∼ 2.5t, consistent with results from the flat system. In a flat system the threshold, ω gap , and the coherence peak location, ω coh , coincide.

At T ∼ 0.08t, Fig. 3.(b), the DOS for V c = 0.1t and V c = t look very similar, with a reduction of ω gap from the T = 0 value and suppression of the coherence peak. The V c = 3t case also shows reduction of ω gap with respect to T = 0, but continues to be distinct compared to the weaker V c cases. Since |∆ i | has not changed significantly with respect to T = 0 (Fig. 2.(b)) these changes are attributable to phase disorder.

By the time T = 0.3t, Fig. 3.(c), the DOS in the three cases are essentially similar, since the |∆ i | homogenises even in the trap (Fig. 2). The density does continue to be inhomogeneous, affecting φ i , but |∆ i | is more important for the low frequency spectrum.

The momentum resolved spectral function, Fig. 4, is more dramatically affected by trapping. The 3 × 3 panel shows the spectrum A(k, ω). The formal definition in terms of HFBdG eigenstates is given in the Supplement. In each panel, the x-axis corresponds to the k scan from {0, 0} to {π, π}, the y-axis is the frequency ω, and A(k, ω) is colour coded as indicated. The columns are for V c = 0.1t, t, 3t (left to right), the rows are T = 0, 0.08t, 0.3t (top to bottom). The size dependence of our results is shown in the Supplement.

The left column at V c = 0.1t shows the thermal evolu- where, for us, k ≈ µ. (ii) Middle: at T = 0.08t coherent particle-hole mixing is almost lost. For k ∼ {{0, 0} → {π/2, π/2}} the spectrum is mainly 'particle-like', while for k ∼ {{π/2, π/2} → {π, π}} it is 'hole-like'. There is significant mixing only near k ∼ {π/2, π/2}. There is a faint surviving trace of the mean field, ±E k , dispersion, the +E k branch for k ∼ {0, 0} and the -E k branch for k ∼ {π, π}. Effectively there are three branches in A(k, ω) at each k. (iii) Bottom: at T = 0.3t there is no trace of the mean field E k , the spectrum is an incoherent combination of upper and lower band features at all k.

For V c = t, middle column, the moderate confinement already shows signatures in A(k, ω). The third column shows results at V c = 3t where the trap center density is n i ≈ 2. The ARPES differs qualitatively from the flat case. (i) Top: at T = 0 the A(k, ω) is very broad since a large number of trap eigenstates overlap with |k . The interband separation now has a maximum for k ∼ {π/2, π/2} and is minimum for k → {0, 0} or {π, π}. This is a case of 'forward bending' rather than backbending. If RF spectroscopy probes the edge of the lower band it would obtain a concave pattern, rather than the convex result that traditionally indicates a pairing gap. The gap, as is obvious from the full A(k, ω) is nevertheless present. (ii) Middle: at T = 0.08t all gaps are smaller compared to T = 0 but the unusual k dependence persists. (iii) Bottom: at 0.3t there is only the hint of the k dependent gap observed at lower T . How do we relate these results to spatial structure?

The overall DOS is N (ω) = -(1/π)Im i G ii (ω), i.e, a sum of the local DOS over the system where G ii (ω) is the local projection of the spin averaged fermion Green's function. If the density n i were slowly varying then as a starting approximation we could use G trap ii (ω) ≈ G f lat (ω, n = n i ). We have checked that this works reasonably even on our 24 × 24 system. The overall DOS is then given by N (ω) ≈ dnP (n)N f lat (ω, n), where N f lat (ω, n) is the flat system DOS at density n and the density distribution

P (n) = 1 N i δ(n -n i ) can be com- puted from the MC density profile.
This immediately creates a connection between the density (and auxiliary field) variation in the trap and the features observed in the DOS. The ARPES, however, involves the overlap k|m of a plane wave state with a BdG eigenstate ψ m . If all ψ m were extended over the system, and overlap all |k , the strange gap modulation with k would not arise.

We find that the BdG states are radially localised to a remarkable degree, the Supplement shows typical real space and momentum space patterns. The lowest energy excitation at T = 0, at E m ∼ 0.9t, is localised near the corners, where n i → 0. This has fourier modes only near k = {0, 0}. For E m > ∼ 1.3t the excitations shift to the center of the trap, and involve modes near k ∼ {π, π}.

Only for E m > ∼ 2.5t, where the BdG states have large weight on the n i ≈ 1 annulus do we see contribution at k ∼ {π/2, π/2}.

Although our system size is larger than accessible in typical DQMC studies, it is well below the ∼ 100 × 100 lattices used in experiments. This is where the local density approximation (LDA) to P (n) becomes useful. LDA prescribes that n trap i ≈ n f lat (µ i ), where µ i = µ -V i and n f lat (µ) can be computed from DQMC or analytic approximations. In the Supplement we compare the MC based ARPES data with results obtained using LDA on the same size. The agreement is remarkable. We extended this to a huge ∼ 200 × 200 system, and all the qualitative features of our original result survive.

Conclusions:

We provide the first solution to the angle resolved spectral properties of an attractive fermion lattice in the presence of confinement, crucial for any cold atom experiment. Even a moderate trapping potential creates a 'core' with low pairing amplitude and generates spectral features that are widely different from the well studied 'continuum' model and the 'flat' Hubbard lattice. We point out a novel 'forward bending' feature that would be the RF spectroscopy signature of a pairing gap, clarify the spatial origin of this feature, and illustrate a scheme that allows access to the spectrum on very large experimentally realised lattices. Large scale determinantal Monte Carlo (DQMC) results are not available for the trapped problem so we compared the results of our method to DQMC data in the 'flat' problem. DQMC results for the superfluid transition temperature (T c ) are available at density n = 0.7 on a 10 × 10 lattice for U/t varying from 2 -8. Fig. 5 compares our results to this benchmark. We capture the non monotonic character, the correct peak location, and our T c estimate is within 10% of the DQMC result at all U/t. This is far superior to mean field theory which would have generated a T c growing monotonically with U/t, with an order of magnitude overestimate already at U/t = 6.

Supplementary material for

II. Homogenisation of |∆ i | with growing temperature

We have seen in Fig. 2 of the main text that the mean value of |∆ i | tends to become independent of position r i with growing T , even though the density n i remains inhomogeneous. We found that this is related to the lower amplitude stiffness of regions with low |∆| at T = 0, and has a direct correspondence with the behaviour in flat systems. We studied the mean value ∆(n, T ) = (1/N ) i |∆ i | n,T in the flat system and discovered that although ∆ is strongly n dependent at T = 0, with a 70% variation as n changes from 1.0 to 0.1, at T = 0.3t that variation is only ∼ 30%. This flat system effect shows up in the trap as a local amplitude stiffness that depends on the T = 0 magnitude of |∆| in that region.

III. Behaviour of Bogoliubov-de Gennes (BdG) wavefunctions

We analysed the BdG wavefunctions in real space and in terms of their momentum content, and show a few illustrative examples in Fig. 7 for V c = 3t. BdG wavefunctions in real space have been represented by u i , v i and in momentum space as U k , V k . The center of the spatial maps is r i = (0, 0). For the momentum maps the center is k = (0, 0) and the corners are (±π, ±π). A. The 1st row of Fig. 7 corresponds to the lowest energy excitation. One can see that (i) the state has large amplitude in the low density region at the corners and, (ii) U k and V k are large near k = 0. The low gap in A( k = 0, ω) arises due to overlap with this excitation. B. Second row, E = 1.326t, higher up in the spectrum. This state is (i) mainly localised at the center of the trap, i.e, the highest density region and (ii) is the first state with significant k = (π, π) content. The |∆ i | here is small, but larger than in the corner region. C. Rows 3 and 4 show states with contribution at k = (π/2, π/2). These are spread over the system but have significant weight in the n ≈ 1 annulus, where |∆| at T = 0 is largest. The states are at significantly higher energy than the states in rows 1 and 2.

IV. Computation of the momentum resolved spectral function.

The spectral function A( k, ω) for a given configuration of ∆ i , φ i has been calculated via the following expression

A( k, ω) = n,En≥0 u n ( k) 2 δ(ω -E n ) + v n ( k) 2 δ(ω + E n ) , (1)
where

u n ( k) = 1 N 1/2 N i=1 e i k• ri u n ( r i ) v n ( k) = 1 N 1/2 N i=1 e i k• ri v n ( r i ).
To get the final A( k, ω) it has been averaged over many equilibrium configuration of ∆ i , φ i .

V. Comparison of exact and local density approximation (LDA) based spectral functions

Pushing the 'local density' approach to the momentum resolved spectral function we checked the accuracy of this approach in capturing A( k, ω) in the trap. We computed the 'local density' based spectral function A L trap ( k, ω) as follows:

A L trap ( k, ω) = P (n)A f lat (n, k, ω)dn, (2)
This prescription is incomplete without specifying P (n). The first approximation is to use the P M C (n) that emerges from the MC itself. This approach, although it does not require BdG based information, still requires MC generated data, and is impractical on large sizes, ∼ 100 × 100, that are likely to be used in experiments. For that P (n) itself needs to be approximated.

We tested the standard prescription that, for a slowly varying density field, one can relate n i to a local chemical potential µ i = µ -V i , where n i and µ i are related by the same equation of state as in the homogeneous system. That relation we infer from numerical results on the flat system. The µ LDA itself is fixed by requiring 1 N i n i (µ i ) = n av . From n i one can generate the 'local density approximation' result P LDA (n). This can be computed easily on any size, and we generated it on 24 × 24 and 192 × 192 lattices.

Fig. 8 compares the 'exact' spectral function at V c = 3t with three approximations (along the row) and three temperatures (down the columns). The first column shows the HFBdG based result for A( k, ω), while the second column shows A L trap ( k, ω) based on P M C (n) integration. The third column shows A L trap ( k, ω) based on P LDA (n) on a 24 × 24 lattice, the fourth column shows the result on a 192×192 lattice. For the larger lattice the corner potential is kept at V c = 3t, as in the small system, so that the larger and smaller systems are roughly equivalent. All the main features of the HFBdG based calculation in column one survive in the P (n) based result, provided an accurate reference is used for the flat system.

V. Size dependence of momentum resolved spectral function.

In this section we present the size dependence of A( k, ω).

We demonstrate that the system size chosen in the main text is sufficient to capture the physics. We have computed (see Fig 9 of this reply) the spectral function on four sizes, 8 × 8, 16 × 16, the original 24 × 24, and 32 × 32. We can go to even larger sizes, but these will already make the point. To compare results at different system sizes in this inhomogeneous system we have kept the 'corner potential' fixed at V c = 3t. The rows, from top to bottom, are for L = 8, 16, 24, 32. The columns, from left to right, are for temperatures, T /t = 0, 0.08, 0.30 as in the paper.

The intent is to identify an unusual 'forward bending' feature in A(k, ω), in the confined system, as the momentum crosses k F .

Looking at the first column, T = 0, it is impossible to infer anything reliable about backbending or 'forward bending' from the L = 8 result. At L = 16 the hints are already clear. At L = 24 the features observed at L = 16 are more prominent, and the wider gap at k ∼ {π/2, π/2}, compared to k = {0, 0}, is convincing. L = 32 is indistinguishable from L = 24.

The persistence of this feature to intermediate temperatures is again clear in L = 24 and L = 32, somewhat ambiguous in L = 16, and non existent in L = 8. So, (i) the system size matters: a small size calculation would not observe the effect, and (ii) our sizes are large enough for the result to be scaled up to experimental lattices (with proper trap scaling).

FIG. 1 .

 1 FIG. 1. Colour online: (a) Ground state of the 2D AHM at U/t = 6 for varying average density nav and trapping strength Vc/t. The tiny region near nav = 1 corresponds to strong density wave (DW) correlation coexisting with superfluidity. The band insulator (BI) refers to the presence of a ni = 2 core. The B and I regions, separated by the dotted line, correspond to 'backbended' and 'inverted' RFS. (b) Finite temperature phase diagram at U/t = 6 and nav = 1.Beyond the small window at weak Vc the system has only SF order at low temperature, with an increasing BI core for Vc/t > ∼ 4. N refers to the normal state. The Tc (on left axis) falls monotonically with Vc as does ∆0 (right axis), the T = 0 gap in the spectrum. The Tc at Vc = 0 would vanish in the infinite volume limit, the results here are for a 24 × 24 system.

FIG. 2 .

 2 FIG. 2. Colour online: Spatial variation and temperature dependence at U = 6t, Vc = 3t. (a) density ni , (b) pairing field magnitude |∆i| , (c) nearest neighbour pairing field correlation. All patterns are thermally averaged.

Fig. 1 .

 1 Fig.1.(a) shows the ground state for varying n av = N f /L 2, where N f is the number of fermions, and corner potential V c . At finite V there is a small window near n av = 1 where DW correlations survive, upto V c /t ∼ 0.8[START_REF] Datta | [END_REF]. Beyond this window the system has only SF order. However, the spatial extent of the SF shrinks with increasing V c or n av since the central part of the trap becomes doubly occupied (n i = 2) suppressing ∆ i .Fig.1.(b) shows the V c -T phase diagram at n av = 1. There is a narrow SF+DW window at small V c , beyond which there is only SF order, with the T c (left axis) decreasing quickly with increasing confinement. The T = 0 spectral gap ∆ 0 (right axis) falls even more sharply, dropping from ∼ 4.6t atV c = 0 to ∼ 1.5t at V c = 6t.Fig.2shows the radial variation of the thermal average ofn i = σ c † iσ c iσ (left), |∆ i | (center), and Φ i = |∆ i ||∆ i+δ |cos(θ i -θ i+δ ) (right). The coordinate i is r =x 2 + y 2 /(L/2), varying along the diagonal. Φ i tracks nearest neighbour correlation in that direction. We have set V c = 3t and n av ∼ 1 and T = 0, 0.08t, 0.3t. The full 2D spatial maps are shown in the Supplement.Fig.2(a) shows the expected monotonic fall in n r at all T . The cloud at T = 0.3t is slightly broader than at T = 0. The pairing field amplitude in 2(b) is more interesting. It is non-monotonic at all T , a peculiarity of the lattice where it grows with n till n = 1 and falls beyond. The T = 0 result for |∆ r | is what is expected from mean field HFBdG theory, with a clear peak in the region where n r ∼ 1. At T = 0.08t the amplitude profile looks similar to T = 0, but with a large growth in the corner where it was zero at T = 0! The trend amplifies at T = 0.3t where |∆ r | is much less inhomogeneous than at T = 0. This is due to the low amplitude stiffness in regions with low |∆ i | at T = 0. We provide a connection to the flat system physics in the Supplement.Fig.2(c) is meant to highlight the suppression of phase correlation with temperature. At T = 0 the phases are locked, so Φ i = |∆ i ||∆ i+δ | ≈ |∆ i | 2 . At T = 0.08t ∼ 0.7T c while the amplitudes are not very different from T = 0 the phase correlation is weakened. By T = 0.3t while amplitudes have grown, NN phase correlations have

Fig. 1 .

 1 Fig.1.(a) shows the ground state for varying n av = N f /L 2, where N f is the number of fermions, and corner potential V c . At finite V there is a small window near n av = 1 where DW correlations survive, upto V c /t ∼ 0.8[START_REF] Datta | [END_REF]. Beyond this window the system has only SF order. However, the spatial extent of the SF shrinks with increasing V c or n av since the central part of the trap becomes doubly occupied (n i = 2) suppressing ∆ i .Fig.1.(b) shows the V c -T phase diagram at n av = 1. There is a narrow SF+DW window at small V c , beyond which there is only SF order, with the T c (left axis) decreasing quickly with increasing confinement. The T = 0 spectral gap ∆ 0 (right axis) falls even more sharply, dropping from ∼ 4.6t atV c = 0 to ∼ 1.5t at V c = 6t.Fig.2shows the radial variation of the thermal average ofn i = σ c † iσ c iσ (left), |∆ i | (center), and Φ i = |∆ i ||∆ i+δ |cos(θ i -θ i+δ ) (right). The coordinate i is r =x 2 + y 2 /(L/2), varying along the diagonal. Φ i tracks nearest neighbour correlation in that direction. We have set V c = 3t and n av ∼ 1 and T = 0, 0.08t, 0.3t. The full 2D spatial maps are shown in the Supplement.Fig.2(a) shows the expected monotonic fall in n r at all T . The cloud at T = 0.3t is slightly broader than at T = 0. The pairing field amplitude in 2(b) is more interesting. It is non-monotonic at all T , a peculiarity of the lattice where it grows with n till n = 1 and falls beyond. The T = 0 result for |∆ r | is what is expected from mean field HFBdG theory, with a clear peak in the region where n r ∼ 1. At T = 0.08t the amplitude profile looks similar to T = 0, but with a large growth in the corner where it was zero at T = 0! The trend amplifies at T = 0.3t where |∆ r | is much less inhomogeneous than at T = 0. This is due to the low amplitude stiffness in regions with low |∆ i | at T = 0. We provide a connection to the flat system physics in the Supplement.Fig.2(c) is meant to highlight the suppression of phase correlation with temperature. At T = 0 the phases are locked, so Φ i = |∆ i ||∆ i+δ | ≈ |∆ i | 2 . At T = 0.08t ∼ 0.7T c while the amplitudes are not very different from T = 0 the phase correlation is weakened. By T = 0.3t while amplitudes have grown, NN phase correlations have

FIG. 3 .

 3 FIG. 3. Colour online: DOS for increasing degree of confinement at three temperatures (a) T = 0, (b) T = 0.08t, and (c) T = 0.3t.

FIG. 4 .

 4 FIG. 4. Colour online:The spectral function A(k, ω) for a 'diagonal scan' from k = {0, 0} → {π, π}. Along the row, Vc varies from 0.1t, t, 3t (left to right). Down the column T varies from 0, 0.08t, 0.3t. The Tc of the unconfined system is ∼ 0.14t, at Vc = 3t it is ∼ 0.1t.

  (i) Top panel: the T = 0 spectral functions are broad since k states overlap with multiple trap eigenstates. Low k states have large (and broad) weight in the lower band while k ∼ {π, π} involves broad weight in the upper band. The gap between the upper and lower bands is still smallest at k ∼ {π/2, π/2} and the backbending feature has not vanished. (ii) At T = 0.08t and T = 0.3t the results are similar to what we saw for the flat case, with some extra (trap induced) broadening noted above.

"FIG. 5 .

 5 FIG. 5. Comparison of superfluid Tc obtained within our static auxiliary field (SAF) scheme, with DQMC (Ref. 15 of main text). The results are for a flat system with density n = 0.7 and lattice size 10 × 10.

FIG. 6 .

 6 FIG.6. The thermal average of |∆| on a 'flat' system for varying density (n) and temperature (T ). The T = 0.001t result corresponds to mean field theory, the finite T results involve fluctuations within the SAF scheme.

FIG. 7 .

 7 FIG. 7. BdG eigenfunctions for 4 different excitation energies. Along the row u m i , U m k , v m i , V m k . Along the column different m, starting with the lowest Em.

and a Bose-Einstein condensed (BEC) form at U/t 1. What is the effect of confinement?

Acknowledgments: We acknowledge use of the Beowulf Cluster at HRI. PM acknowledges support from a DAE-SRC Outstanding Research Investigator Award.