Impurities in systems of noninteracting trapped fermions

We study the properties of spin-less noninteracting fermions trapped in a confining potential but in the presence of one or more impurities which are modelled by delta function potentials. We use a method based on the single particle Green's function. For a single impurity placed in the bulk, we compute the density of the Fermi gas near the impurity. Our results, in addition to recovering the Friedel oscillations at large distance from the impurity, allow the exact computation of the density at short distances. We also show how the density of the Fermi gas is modified when the impurity is placed near the edge of the trap in the region where the unperturbed system is described by the Airy gas. Our method also allows us to compute the effective potential felt by the impurity both in the bulk and at the edge. In the bulk this effective potential is shown to be a universal function only of the local Fermi wave vector, or equivalently of the local fermion density. When the impurity is placed near the edge of the Fermi gas, the effective potential can be expressed in terms of Airy functions. For an attractive impurity placed far outside the support of the fermion density, we show that an interesting transition occurs where a single fermion is pulled out of the Fermi sea and forms a bound state with the impurity. This is a quantum analogue of the well-known Baik-Ben Arous-P\'ech\'e (BBP) transition, known in the theory of spiked random matrices. The density at the location of the impurity plays the role of an order parameter. We also consider the case of two impurities in the bulk and compute exactly the effective force between them mediated by the background Fermi gas.


Introduction
Noninteracting fermions in a confining trap is a topic of much current interest, especially in the context of cold atoms. In the presence of a trap, the density of the Fermi gas in the ground state is confined in a finite region of space. Indeed, the density vanishes outside a finite interval in one-dimension. Inside this interval, usually referred to as the "bulk", the fermion density can be estimated, for a large number of fermions N , using a semi-classical approximation, or equivalently the so-called local density approximation (LDA) [1][2][3]. Near the edge where the density vanishes, the quantum fluctuations play a dominant role and the local properties of the fermions are very different from that of the bulk [4][5][6][7]. This edge region is called the "Airy gas" because the Airy functions play an important role in describing the quantum correlations. It is well known that the LDA is a very good approximation when the confining potential is smooth. However, if this potential has singularities, such as a step or delta-function, this method fails. In a recent paper [8], we have examined the effect of a step potential, using exact methods based on the determinantal properties of the Fermi gas.
Here, we consider instead the case when the smooth confining potential is modulated by introducing one or more delta functions. This situation naturally arises when one introduces one or more immobile impurities in the Fermi gas, where the impurities are modelled by deltafunction potentials (attractive or repulsive). In the absence of the trap, i.e., for a free Fermi gas, the effects of such impurities have been well studied in the literature. For example, when a single impurity is introduced in the Fermi gas, the density of the Fermi gas is modulated near the impurity. At distances far from the impurity, the density exhibits decaying oscillations, famously known as "Friedel oscillations" [9][10][11]. How this density gets modulated close to the impurity has been studied in [12] in one dimension for delta function potentials. In addition, the effective Casimir interaction between two impurities (mediated by the background Fermi gas) has also been studied [13,14], however, again, the results obtained are only valid at large distances.
In this paper we employ a method based on the single particle Green's function that allows us to obtain exact results for the trapped Fermi gas. We first consider the case when a single impurity is added to the trapped Fermi gas in three different locations: (i) in the bulk, (ii) near the edge and (iii) outside the edge. In the bulk (case (i)), where the trapped Fermi gas behaves locally as a free Fermi gas, we obtain the explicit form of the density near the impurity at all scales, not necessarily large. At large distances, we recover Friedel oscillations, our short distance results agree with those of [12] when the impurity is repulsive but we show that a correction is needed for attractive impurities . However in cases (ii) and (iii), the presence of the trap considerably modifies the bulk results and we obtain new results for the density close to the impurity. In addition, in all the three cases, we compute the effective potential felt by the impurity due to the background Fermi gas. In case (iii), for an attractive impurity, we show that an interesting transition occurs where a single fermion is pulled out of the Fermi sea and forms a bound state with the impurity. This is a quantum analogue of the classical Baik-Ben Arous-Péché (BBP) transition [15,16], known in the theory of spiked random matrices, where an eigenvalue detaches from the bounded support of the eigenvalues, due to a rankone perturbation. This rank-one perturbation is the analogue of the delta-function potential induced by the impurity in the Fermi gas and the eigenvalue is the analogue of the fermion position. We then go beyond the case of a single impurity and study the effects of adding two impurities. In this case we compute the effective Casimir-like interaction between the two impurities mediated via the trapped Fermi gas. For two impurities, we restrict our analysis to case (i) where both the impurities are placed in the bulk. In this case, we obtain exact results for this effective interaction at arbitrary separation between the impurities. At large distances, our results are in agreement with the one found in Refs. [13,14] obtained by a different method.
We restrict ourselves here to the case of immobile impurities. The case of mobile impurities has been extensively studied for fermionic systems both theoretically [17][18][19][20][21] and experimentally [22]. The set up of immobile impurities that we study in this paper is more difficult to access experimentally, however it has been suggested that impurities could be introduced by superimposing an optical lattice on an overall trapping potential [23].
2 Model and summary of main results

The model
We consider a gas of identical spin-less fermions of mass m in a trap generated by a potential V (x) at zero temperature. We then add n delta-impurities of strengths g i . The single particle Hamiltonian is then given by where H 0 is the Hamiltonian associated with the trap and ∆H corresponds to the impurities. The eigenfunctions and eigenvalues of H 0 are denoted by ψ 0 k (x) and 0 k . Similarly, ψ k (x) and k denote the eigenfunctions and eigenvalues of H. Consider first the case with no impurity (i.e. g i = 0 for all i = 1, · · · , n). At zero temperature, the system is in the many-body ground-state such that all single-particle states of H 0 below the Fermi level, denoted by µ, are occupied, each by a single fermion. The ground-state energy is given by where θ(x) is the Heaviside function where one usually uses the definition θ(0) = 1. We now switch on the g i 's, i.e., we introduce the impurities in the system. The single-particle Hamiltonian then changes from H 0 to H. This will change the single-particle eigenfunctions and eigenvalues. Consequently the ground state energy will also change. Here we work in the grand-canonical ensemble where the Fermi level µ remains fixed, while the number of fermions is not fixed and the system is in contact with a reservoir of particles. Then the new ground state energy, in the presence of the impurities, is given by Similarly one can define the number of particles N 0 (µ) and N (µ) below the Fermi level µ as Since we are working in the grand-canonical setting, N 0 (µ) and N (µ) can be different and the quantity which will play a crucial role is the grand-potential at zero temperature At zero temperature, thanks to the Wick's theorem, all the correlation functions are given by determinants constructed from the so-called kernel, which reads for H 0 and H respectively [5] Setting x = y in the Eq. (6) we find the fermion density

Outline and main results
In Section 3 A, we first recall the method of the Green's function introduced in [8] and present in Section 3 B the explicit expressions for the Green's function in the absence of impurities. This is then used as the building block to obtain an exact expression for the Green's function in the presence of the impurities in Section 3 C. In Section 4 we examine the effect of impurities in the bulk of the system. In the absence of the impurities, the density in the bulk is given by Here k F (x) is just the local Fermi wave vector. The density vanishes at the edge x e where V (x e ) = µ. The bulk of the system is thus defined such that V (x) µ (see below in Eq. (29) for a more precise definition). We now add a single impurity in the bulk of the system, say at x = 0, and investigate how the kernel and the density change near the impurity. We reparametrize the impurity strength in terms of an inverse length scale λ defined as We show that the change in the kernel upon adding this impurity is given by where Im denotes the imaginary part and E 1 (z) = ∞ z dt e −t /t denotes the exponential integral [24]. Here k F = k F (0) is the Fermi wave vector at the impurity position. Putting x = y in (10) we obtain the change in the density due to the impurity (see also Figs. 4 and 5) We show how this formula recovers the density change computed in [12] for a free Fermi gas when the impurity λ is repulsive. The result given here and in [12] is exact for a homogeneous free fermion system at all impurity strengths and distances and we note that formulas given for these Friedel oscillations [9][10][11] are often given in the regime of linear response or at long distances. However our formula in (11) holds for any x and λ. In particular, we obtain an explicit formula for the density at the position of the impurity We then compute the effective potential V eff (x 1 ) felt by the impurity, where x 1 denotes the position of the impurity in the bulk, due to its interaction with the Fermi gas. This is obtained by computing the change in the ground state energy of the many-body system due to the addition of an impurity. This effective potential can be expressed in the scaling form where Ω 0 (µ) = E 0 (µ) − µN 0 (µ) and the scaling function is given by Figure 1: The scaling function for the effective potential W (γ) defined in Eq. (14). Note that although W is discontinuous at γ = 0 the potential V eff is continuous. Note that although W is discontinuous at γ = 0, the potential V eff (x 1 ) is continuous as a function of λ at λ = 0.
The function W (γ) is shown in Fig. 1 and has the asymptotic properties In Section 5 we investigate what happens when the impurity is placed at x 1 to the right of the edge x e , such that V (x 1 ) > V (x e ) = µ. When the impurity is attractive we show that a phase transition occurs as the reduced strength λ A = −λ > 0 of the impurity is increased beyond a critical value. We call this transition a filling transition. Let us recall from elementary quantum mechanics of a single particle that a delta potential introduced at x 1 , in addition to a flat potential V (x) = V 0 , introduces a single bound state with wave we find that there are two different phases: (i) weak impurity, where E b > µ implying that this bound state is unoccupied and consequently ρ µ (x 1 ) 0; (ii) strong impurity, where E b < µ implying that this bound state is occupied and consequently ρ µ (x 1 ) λ A . The transition occurs exactly at E b = µ, which corresponds to λ A = κ µ (x 1 ) = 2(V (x 1 ) − µ). This filling transition has some similarities with the BBP transition in random matrix theory, where a rank one perturbation to a random matrix can displace the maximal eigenvalue of the matrix [15,16].
In Section 6 we study the effect of adding two impurities in the bulk, say at x 1 and x 2 , close to x = 0. We assume that x 1 and x 2 are such that |V ( is the Fermi wave vector at x = 0. This condition ensures that the potential remains effectively constant on the scale of the separation r = |x 1 −x 2 | between the impurities. We show that the effective Casimir-like interaction between these two impurities, mediated by the background Fermi gas, is given by where ζ = k F r and γ i = λ i /k F are the scaled impurity strengths. Our result is valid for all ζ. Interestingly, using free fermionic field theory, an expression for V int (r, k F , γ 1 , γ 2 ) was derived in Refs. [13,14], which reads where Li 2 (z) = ∞ n=1 z n /n 2 is the di-logarithm function and Re denotes the real part. Our formula (16) can be shown to reduce to (17) when ζ 1. However, this form (17) is an approximate form that holds only for large ζ. As ζ → 0, V int (r, k F , γ 1 , γ 2 ) in Eq. (17) diverges, which is not physical. Instead, our exact result (16), which holds for all ζ, approaches to a constant as ζ → 0 (see Fig. 6).
In Section 7 we analyse the effect due to an impurity close to the edge of the Fermi gas, i.e. x 1 ≈ x e . The Friedel oscillations [9][10][11] around impurities in the bulk are strongly suppressed near the edge. However weak oscillations, that are already present at the edge without any impurities, still persist in the presence of impurities. The main effect of the impurity is to alter the phase of these oscillations. For an attractive impurity, we show that the filling transition discussed above becomes a smooth crossover on the scale of the inter-particle distance at the edge and the local density profile is described by a universal scaling function. Finally we obtain an analytic expression for the effective potential acting on the impurity placed in the edge region.
In Section 8 present our general conclusions and perspectives for future studies.

Basic formalism and set up
We now describe how the single particle Green's function can be used to extract the change in the kernel due to the addition of an impurity at a fixed position in the system. The results given below are derived in detail in a recent paper [8] and we refer the reader there for detailed derivations.

Kernels via Green's functions
The Green's function G µ (x, y) associated to the Hamiltonian H in (1) is defined for an arbitrary running Fermi energy µ as The Green's function has poles at µ = k + i0 + , i.e., infinitesimally above the real axis in the complex plane. In operator notation we also have the equivalent resolvent representation from which we see that G µ (x, y) is solution to the equation The kernel can be obtained from the Green's function from the following formula where the imaginary part can be taken outside the integral as the integration contour is real. Noting that when µ → ∞ the kernel becomes the sum over a complete set of states, we can derive an alternative representation In this paper we will be interested in the change in the kernel due to an added impurity. If we denote the Green's function in the absence of impurities by G 0µ (x, y), then we can write the Green's function as G µ (x, y) = G 0µ (x, y) + ∆G µ (x, y) where ∆G µ (x, y) is the change in the Green's function due to the impurities. It is then easy to see that the change in kernel if one uses Eq. (21), or alternatively if one uses Eq. (22). The fact that these two representations (23) and (24) are equivalent can also be seen from the following argument. These integrals can be interpreted as contour integrals in the complex µ plane along the real axis. In Fig. 2 these two contours are represented as Γ 1 = (µ, ∞) for (24) and Γ 2 = (−∞, µ) for (23), along with the position of the poles which are infinitesimally above the real axis and are shown by crosses. In terms of the contours Γ 2 and Γ 1 shown on the figure we have ∆K µ (x, y) = 1 π Im Γ 1 dµ ∆G µ (x, y) = − 1 π Im Γ 2 dµ ∆G µ (x, y). The equivalence of the two representations can also be demonstrated as follows. First, using Cauchy's theorem which, as there are no poles in the lower half of the complex plane, gives Γ 1 ∪Γ 2 ∪Γ 3 dµ ∆G µ (x, y) = 0, where Γ 3 is taken to be an infinite semicircle, with center at the origin, in the lower half of the complex plane. One can then show that Γ 3 dµ ∆G µ (x, y) = 0 to obtain the desired result. Finally we should point out that by rotating the contour Γ 1 about z = µ by −π/2 (so it is parallel to the imaginary axis) gives an integral representation that corresponds to the sum over Matsubara frequencies in the fermionic field theory setting [13,14].
It is important to note here that the representation given in Eq. (24) has a number of advantages over that in Eq. (23). First, it involves an integral over µ > µ therefore we do not need to know the Green's function ∆G µ (x, y) for small µ < µ. Since µ is large, we just need to know the Green's function for large µ which can be conveniently computed using the semi-classical approximation. Furthermore, we will see that the representation in Eq. (24) is more suitable to asymptotic analysis of certain formulas. . The contour Γ 3 is used to close the contour Γ 1 ∪ Γ 2 and is taken to be a semi-circle in the lower half of the complex plane whose radius is taken to ∞.

Bulk and edge Green's function
In this section, we recall the results obtained in Ref. [8] for the kernel both in the bulk as well as at the edges, using the Green's function method, for a smoothly varying trapping potential V (x).
In the bulk. We start with the bulk and consider the potential around a point x 0 where we assume that V (x) ≈ V (x 0 ) and solve Eq. (20) with a constant potential V (x 0 ). This gives where is the local Fermi wave vector given the Fermi energy µ . The positive value of the square root is taken and it is understood to have a negative infinitesimal imaginary part as indicated above. By inserting this expression (25) in Eq. (21), setting x = y and performing the integral over µ one finds the density in the bulk These results are obviously exact for a flat potential V (x) = V 0 . However they are also accurate as long as the relative variations of k µ (x) on microscopic scales, of order O(1/k µ ) are negligible, i.e., Using (27), this condition (28) translates to [8] Note that this argument naturally introduces a length scale which sets the size of the region over which this assumption that V (x) is constant holds. The condition in (29) clearly gets violated in two cases: (i) when the potential is not smooth, for instance when the potential exhibits a step structure as studied in [8] or when there is a delta function contribution to the potential (ii) when the analysis is carried out at the edge of the trap where the density in Eq. (27) vanishes.
At the edge. The edge x e of the Fermi gas occurs where the density vanishes, i.e. when We see from Eq. (29) that R diverges as x 0 → x e . The physics near the edge region (the so called Airy gas), can be studied by linearizing the potential V (x) near x = x e , i.e.
In this region, by solving Eq. (20) with a linear potential (32), the Green's function can be written in the scaling form [8] where w e and α e are respectively the length scale associated with 1/ρ e , where ρ e is the fermion density at the edge [5] and the energy scale associated with the edge given by The function g e (ζ, ζ ) is given by For later purposes, we note that Outside the bulk. In the classically forbidden region, far outside the edge where the condition in (29) holds, we will again assume that V (x) is slowly varying around a point x 0 . In this case, the solution of Eq. (20) reads where is positive. This result is similar to the one found in the bulk in Eq. (25) with a local Fermi wavevector which is imaginary. One can check, using the asymptotic properties of the Airy functions Ai(z) and Bi(z), that the edge result for the Green's function (33) matches (i) on the left with the bulk result (20) and (ii) on the right with the result far outside the bulk (38).

Treating delta function potentials
In this section we show how one can use the Green's function method to study systems where the potential has a delta-function part. The delta function potential has been extensively studied in quantum systems to model impurities [25]. Here we show that the Green's function method is particularly well suited to study this problem.
We start with a single particle Hamiltonian H 0 and denote its Green's function by We now add n impurities so that the total Hamiltonian is The coordinates x i 's are the positions of the impurities and g i 's denote their interaction strengths with the fermions. Note that the sign of g i can be positive (repulsive) or negative (attractive). Many methods have been found [25] to extract the Green's function G µ = (µ − i0 + − H) −1 . A simple way to do this is to observe that Note that the A i 's depend implicitly on both the x i 's and y. The solution of this linear equation (44) can be expressed as where the n × n matrix R has components R ij given by This leads to the change in the Green's function For general x and y the above expression is quite complicated. However, if we consider the Green's function at points where there are impurities and define the n × n matrices G 0 and G such that G 0ij = G 0µ (x i , x j ) and G ij = G µ (x i , x j ) things simplify a bit. In this case, using Eq. (47) and adding to G 0µ , one gets in the matrix form where the matrix Λ g has components Λ gij = g i δ ij . Noting that Eq. (46) implies R + Λ g G 0 = I and multiplying both sides by where we further used Eq. (46) for R. Thanks to this relation, one can compute G if one knows G 0 . Using, furthermore, the formula for the derivative of the determinant of a matrix with respect to a parameter t, i.e.
we can rewrite Eq. (49) as We will use this representation later. Let us consider the simplest case of a single impurity located at x 1 with an amplitude g 1 = g (this corresponds to n = 1 in Eq. (42)). In this case, the full Green's function can be computed from Eq. (43), (45) and (46) which has a simple "Schwinger-Dyson" form. Expanding the denominator in powers of g, this formula (52) has a simple physical interpretation: it adds up contributions to the Green's function arising from no scattering, one scattering, two scatterings, etc, from the impurity. Exactly at the position of the impurity it reads We will see below that this formula is particularly useful to deduce the fermion density at the impurity as well as the energy change induced by the introduction of an impurity. In the following, we will use the formulae derived in this section in the case where H 0 corresponds to a smooth trapping potential in the absence of a delta-potential.

Impurity in the bulk
Here we consider the effects of delta-function impurities in the bulk. The overall smooth trapping potential V (x) described by H 0 (see Eq. (1)), as discussed before, can be taken to be locally constant and without loss of generality we set V (x) = 0. Thus the Fermi wave vector at Fermi energy µ is given by k F = k µ = √ 2mµ/ . The local density of the system without impurity is given by ρ 0µ = k F π . The effect of impurities has been well studied in the literature, notably the density around impurities exhibits the well known Friedel oscillations [9][10][11]. However, in most previous studies only the behavior of the density at distances greater than the inter-particle distance 0 = 1/ρ 0µ or in the linear response regime (i.e. to first order in λ). In [12] the density change induced by a delta function impurity was studied exactly. Here, by using a more versatile method, we show how this result can be rigorously extended to inhomogeneous bulk systems. We also show how the result given in [12] is, simply, modified to take properly into account the appearance of a bound state in the case of an attractive impurity.
It is important to know the exact behavior of this density at the location of the impurity since, as we will see, it determines the effective interaction between the impurity and the surrounding fermions. We show that this density at the impurity is finite and depends on the local Fermi-wave vector. We also compute the effective interaction and show that it is given by Eq. (14).

Friedel oscillations
We start by computing the kernel in the presence of a single impurity of strength g 1 = g placed at x 1 = 0. Here the change in the Green's function, using Eqs. (24), (25) and (52) for the bulk Green's function yields a change in the kernel ∆K µ (x, y) = K µ (x, y) − K 0µ (x, y) which is given by Now we change variables µ = 2 k/2m (so dµ = 2 kdk/m) and we assume that µ > 0 so that the integration over k is along the real axis. The way in which all the contours in Fig. 2 transform under the transformation k = 2m 2 µ is shown in Fig. 3 along with the position of the original poles in the Green's function shown again as crosses. This gives is an inverse length scale associated with the impurity, and k F = k F (0) = √ 2mµ/ is the Fermi wave vector at the position of the impurity.
The integral in Eq. (55) corresponds to the contour Γ 1 . The contour Γ 1 can be deformed onto the contour Γ 4 as the integral over the contour Γ 3 is zero and there are no poles crossed during this deformation [8], to give Figure 3: Contour integrals used in the integral representations of the kernel in terms of the variable k = 2mµ / 2 . Crosses indicate the poles of the Green's function in the complex plane of k. The contours Γ i for i = 1, 2, 3 correspond to the contours Γ i shown in Fig. 2 when mapped into the k plane. Note that the contour Γ 2 has two components: one vertical, and one horizontal going from k = 0 to k = k F . The contour Γ 3 denotes a portion on Γ 3 , shown in bold, which is useful for later purpose. The contour Γ 4 = (k F , k F − i∞), which is useful for asymptotic analysis, is also shown.
The above can also be written in terms of standard functions via the change of variables where we recall that is the exponential integral function [24]. Taking x = y we find that the change in the local density ∆ρ µ (x) around the position of the impurity at x = 0 is given by This is an important result of this paper, since it is valid for all values of λ and at all distances x. The comparison between this result (60) and the result of Ref. [12] is performed in Appendix A. We now study this form of the density profile in Eq. (60) at distances respectively very close and very far from the impurity.
Density close the impurity. Here we analyse the density in Eq. (60) very close to the impurity, i.e. the limit x → 0. For small x we can use the asymptotic expansion [24]. In particular at x = 0, taking the imaginary part of the logarithm, we get the total local density (with the bulk term where tan −1 above is the principle branch such that tan −1 (0) = 0. The change in the density is equal to zero when λ = 0 as because k F ≡ k F − i0 + one has that arg(ik F ) = π/2. The change in the local density is significant if λ ∼ k F . Repulsive case. In the limit of strong repulsion if λ = λ R with λ R > 0 and λ R k F we find while for λ R k F one has Attractive case. In this case, setting λ = −λ A with λ A > 0, we find, for strong attraction Note that the dominant term λ A in Eq. (64) is the contribution from the single bound state associated to the attractive delta function. Indeed such a bound state has a wave function which gives rise to a density |ψ bs (0)| 2 = λ A . On the other hand, for weak attraction Density far from the impurity and Friedel Oscillations. In this case the asymptotic expansion [24] for |ζ| 1 and | arg(ζ)| < 3π/2, can be used but one can also use the representation given in Eq. (57) where the long distance behavior comes from an expansion about κ = 0. We find from (60) that the density for large x decays as which can be rewritten as At large |x|, both the repulsive and attractive cases are described by the formula (69) with λ = λ R > 0 in the repulsive case while λ = −λ A < 0 in the attractive case. In the limit of small λ, one finds which is the standard formula for large distance one-dimensional Friedel oscillations in the regime of linear response [12].
In Fig. 4 we plot the relative perturbation of the exact density as a function of ζ = k F x where we recall that ρ 0µ is the local density in the absence of the impurity and where we have written λ = γk F [26]. In Fig. 4 we plot n(ζ, γ) in the repulsive case with γ = 1. The asymptotic expansion Eq. (69) is also shown as an orange dashed line. For ζ > 3, this asymptotic form describes accurately the exact result. However, when extrapolated to small values of ζ, this asymptotic form diverges while the exact result approaches a finite value as ζ → 0 (see Eq. (61)). In Fig. 5 we plot n(ζ, γ) for an attractive impurity with γ = −1. We see the signature of the localized wave function about the impurity which causes an increase in the local density. Again we see that the asymptotic approximation for Friedel oscillations becomes accurate only for ζ ∼ 3.
Going beyond the density we can also analyse the kernel K µ (x, y) in Eq. (58) for large |x| and |y| by using the same asymptotics for the Exponential integral. We find at large |x| and |y|, both for the repulsive and attractive cases,

The effective potential acting on an impurity
Here we examine how the total energy of the fermion system is changed by adding of a single impurity at fixed Fermi energy µ. The Hamiltonian H(λ) depends explicitly on the parameter λ = gm/ 2 where g is the impurity strength. We denote by k (λ) the k th energy level, as a function of λ. We define the total energy E(µ, λ) and the total number of fermions N (µ, λ) as The goal is to compute the effective potential felt by the impurity at position x 1 which can be identified as the change in the grand-potential (since we are working at fixed Fermi energy µ) To perform this computation, we use the Hellmann-Feynman theorem which states that where ψ k (x, λ) is the eigenstate associated to the energy level k (λ). In the present case of the delta-impurity this theorem (75) gives From this, one sees that every energy level is moved up for repulsive impurities and down for attractive ones. Thus at fixed µ, the derivatives of the energy E(µ, λ) and of the number of particles N (µ, λ) in Eq. (73) with respect to λ read Therefore, using Eq. (74), the derivative of Ω(µ, λ) with respect to λ, using (77) and (78) together with (75) is given by where we have made explicit the dependence of the fermion density on λ. The effective interaction of the impurity with the fermion system is thus given by We can now use the expression for the fermion density given in Eq. (61) to obtain where the function W (γ) is given in Eq. (14).
In the limit of a weak impurity strength |λ| k F (x 1 ) we find from (81) Therefore for λ > 0 the impurity is pushed away from a dense region, while for λ < 0 it is attracted by dense regions. On the other hand, for a strong impurity strength, |λ| k F (x 1 ), we find In the case where λ < 0 we see that V eff (x 1 ) contains a term corresponding to the bound state energy, E b = − 2 λ 2 /2m of the state localized around the impurity. Again we see that repulsive impurities are repelled from dense regions and attractive impurities are attracted by dense regions, which clearly agrees with physical intuition.
Link with the problem of mobile impurities. The above results on the effective interaction potential felt by an impurity in an inhomogeneous system is to our knowledge new, the problem for a homogeneous system was discussed briefly in [27]. However a problem with a similar flavor, involving a mobile impurity, has been studied. McGuire [17,18] considered the problem of N identical spin-less fermions with no mutual interactions. In this system, one introduces an additional particle, with coordinate x 0 , which interacts with each of the N fermions via a delta function potential. In its most general form we can consider the N + 1 body Hamiltonian given by where M is the mass of the impurity particle and V(x) the effective potential it feels due to the trap. The problem examined by McGuire corresponds to identical masses, i.e. M = m and to homogeneous system with V (x) = 0, hence k F (and thus the density) being constant. The change in the energy due to the additional particle is found to be Interestingly, we note that this formula (85) is strikingly similar to the expression obtained here in the case of an immobile impurity and one can write with the same scaling function W (γ) given in Eq. (14). Note that the problem considered in this paper corresponds to the limit M → ∞ where the position x 0 is fixed. More recently the problem introduced by McGuire was revisited in the presence of an external harmonic potential V (x) = mω 2 x 2 /2, which is the same on both the fermions and the impurity particle. The effect of inhomogeneity was treated by combining McGuire's result with an LDA-like approximation, which turns out to be remarkably accurate even for systems with a small number of fermions [20].

Impurity far from the bulk and the filling transition
We consider an impurity of strength g 1 = g placed at x 1 far from the bulk of the Fermi gas, where the condition in Eq. (29) holds and the density vanishes. One may ask the question whether an attractive impurity can pull fermions out of the bulk. In this region the kernel in the region of the point x 1 is given by where in the above integral µ < µ V (x 1 ) (and thus the contour Γ 2 in Fig. 2 is the appropriate one to use). More precisely the condition in Eq. (29) holds in the above integral and so we can use Eq. (38) for the Green's function. When an impurity is placed at the point x 1 the induced change in the Green's function is where κ µ (x 1 ) is given by Eq. (39). Now using dµ = − 2 κdκ/m we find that the change in the kernel is Here and below we denote κ µ (x 1 ) = κ F (x 1 ). We now use the standard identity which gives, for λ > 0, as κ F (x 1 ) > 0. Hence a repulsive impurity far from the bulk has no effect on the Fermi gas. However for an attractive impurity, writing λ = −λ A with λ A > 0 we find The above result is easily interpreted physically. It can be written as where ψ b (z) is the bound state wave-function for a delta potential at z = 0 and in the absence of any other potential given in Eq. (65). The kernel well outside the bulk is thus generated by a single particle bound state. The fermion density at the position of the impurity is thus given by It exhibits a transition as a function of λ A . It vanishes when λ A < κ F (x 1 ) and is nonzero for λ A > κ F (x 1 ). This corresponds to a filling transition of the bound state where the density exhibits a "jump" by κ F (x 1 ). The transition is sharp for large κ F (x 1 ). At smaller values of κ F (x 1 ), i.e. close to the edge, it is replaced by a smooth crossover, which is analysed in Section 7. This transition can be interpreted by the following energy argument. The energy of this bound state in the presence of a local potential is given by Hence this state is occupied if E * b < µ, which corresponds to λ A > κ F (x 1 ). In contrast, when E * b > µ this bound state energy level exceeds the Fermi energy and hence it remains unoccupied at zero temperature. We note that this type of transition is quite generic, and not specific to a delta-function impurity. For instance, it can occur in a more general context when there is a second additional potential well (i.e. a second minimum of the trapping potential). When the Fermi energy increases above this second minimum, a new disjoint interval arises in the support of the density. However the case of the delta potential yields a particularly simple tractable example, since it corresponds to a rank one perturbation.
A natural question to ask is whether there is an effective potential felt by the particle when it is outside the bulk. This potential can, again, be derived using the Hellmann-Feynman theorem (80) together with the expression for the density outside the bulk given in Eq. (95). We find Furthermore we note that using Eq. (39) we can write Hence we see that for λ A > κ F (x 1 ), the x 1 -dependence of the effective potential V eff (x 1 ) is the same as the trapping potential V (x 1 ), this is due to the fermion which forms a state bound about the impurity.
An analogy in random matrix theory. This filling transition is reminiscent of the Baik-Ben Arous-Péché (BBP) transition in random matrix theory [15,16].
, with Gaussian fluctuations). The analogy with our quantum problem is suggested by the fact that (i) the joint PDF of the eigenvalues of M 0 is identical to the joint PDF of the fermion positions in the ground state of the harmonic oscillator V (x) = 1 2 x 2 (ii) the perturbation is of rank one in each problem. It is then tempting to establish an analogy between the spike/outlier from the Wigner sea, and the fermion bound to the delta impurity outside the Fermi sea. A difference is that the analog of the typical value of λ max would be x 1 , which in our problem is given. However, in both cases the order parameter of the strong coupling phase (i.e. γ > 1 in RMT or λ A > κ F (x 1 ) in the fermion problem) is the overlap of the state of the system with the perturbation. Note that the BBP transition has a non-trivial critical regime when γ − 1 = O(N −1/3 ), where λ max gradually leaves the edge of the spectrum. The critical region in our problem corresponds to the case where the impurity is placed near the edge studied below in Section 7.

Interaction between two impurities in the bulk
In this section we compute the effective interaction between two impurities in the bulk separated by a distance r. We place impurity 1 at x 1 and impurity 2 at x 2 . We assume that |x 2 − x 1 | ξ where the length ξ is defined in Eq. (30) such that the trapping potential can be considered as constant. For two impurities we can still apply the Hellmann-Feymnam theorem, for example differentiating with respect to λ 2 , which measures the interaction strength of impurity 2. Writing explicitly the dependence on the coupling constants λ 1 and λ 2 , we obtain the analogue of Eq. (79) valid for two impurities where Ω(µ, λ 1 , λ 2 ) is the grand-potential of the system in the presence of the two impurities (it depends on both x 1 and x 2 but we omit the explicit dependence for notational simplicity). In Eq. (99), ρ µ (x 2 , λ 1 , λ 2 ) denotes the density of the Fermi gas at the location of the second impurity. Let us define the effective interaction between the two particles as which depends only on the distance r = |x 2 − x 1 | since the system is translationally invariant on scale of the order O(ξ). The interaction potential V int (r) in Eq. (100) can be written as Now using the representation of ρ µ (x 2 , λ 1 , λ 2 ) in terms of the Green's function in Eq. (21) by We now use Eq. (51) which can be written as thus facilitating the integration with respect to λ 2 in Eq. (101). This then yields Using the fact that G 0µ (x, y) = G 0µ (x − y), we get where we recall that r = |x 2 − x 1 |. Now using the expression for the Green's function in Eq. (25), valid in the bulk, and again changing the integration variable to k where µ = 2 k 2 /2m, we find where the contour Γ 2 is shown in Fig. 3. As λ 1 and λ 2 are real, there are no poles inside the region enclosed by the contours Γ 2 , Γ 4 and Γ 3 shown in Fig. 3. Therefore Cauchy's theorem tells us that the contour integral around this region is identically zero. Using further the fact that the integrand vanishes on the contour Γ 3 as its radius is extended to ∞ we obtain Making the substitution k = k F − iκ where k F = k F (x 1 ) ≈ k F (x 2 ), we find the exact result Writing λ i = k F γ i and making the change of variable κ = k F u we find the interaction energy, with its dependence on the physical parameters r, k F , γ 1 and γ 2 made explicit, is given by For k F r 1 one can expand the integrand about u = 0. The corrections due to terms of order u are of order 1/k F r (see Eq. (112) below). We thus find for r large This is exactly the result obtained in [13,14] via a field theoretic method based on the summation of Matsubara frequencies, which becomes a continuous integral at zero temperature.
In this large distance approximation the integral can be evaluated to give where Li 2 is the di-logarithm function [24].
The interaction potential can then be written in terms of the Fermi energy E F = 2 k 2 F 2m and the variable ζ = k F r to give and with the asymptotic form for large ζ given by In Fig. 6a we show the interaction energy in units of E F , U (ζ, γ 1 , γ 2 ) = V int (ζ, k F , γ 1 , γ 2 , )/E F , for the exact result Eq. (112) for two repulsive impurities with (γ 1 , γ 2 ) = (1, 1), as a function of ζ, along with the corresponding large distance approximation of Eq. (113). As was the case for Friedel oscillations the asymptotic result is accurate for ζ > 3 but diverges towards −∞ as ζ → 0. The potential oscillates in a manner reminiscent of the Friedel oscillations, exhibiting local minima, but is attractive for ζ < 1. In Fig. 6b we show the corresponding result for two attractive impurities (γ 1 , γ 2 ) = (−1, −1), again the potential oscillates and presents local minima. A sharp barrier appears at ζ ∼ 1/2, but a deep minimum is formed for small ζ.
Here, in contrast to the case of impurities of the same type, the short distance behavior of the potential is repulsive.

Impurity near the edge
In this section we investigate the effect of adding a delta function impurity near the edge at x = x e such that µ = V (x e ). In this region, in the absence of the impurity the Green's function satisfies the scaling form given in Eq. (33). The width of this region is given by w e and the corresponding energy scale is denoted by α e , both displayed in Eq. (34). Consider now an impurity located at x = x 1 = x e + z 0 where z 0 will be of the same order as w e . Substituting the scaling form from (33) into (52), we obtain the change in the Green's function due to the impurity as where g e (ζ, ζ ) is given in (35) and is a dimensionless measure of the impurity strength in the edge region. The kernel K 0µ at the edge in the absence of impurity is given by, in terms of the Airy kernel K Ai . The change in the kernel, ∆K µ = K µ − K 0µ , is obtained from Eq. (24) by integrating over µ between µ and +∞, making the change of variables z 0 we + µ−µ αe = −u and setting z 0 = cw e then gives where function g e is given in (35) and the dimensionless number c given by measures the relative position of the impurity compared to the edge. The above integral has an integrand which oscillates and decays like 1/u for large u and it cannot, at least in any obvious sense, be evaluated analytically. However if we use Eq. (23) we find the alternative expression which converges quickly for u → +∞ allowing an efficient numerical integration (see below).

Density at the edge
Of particular interest is how the density is modified by the presence of a delta function at the edge. The average density of the Fermi gas around the impurity, can now be obtained by setting coinciding points in the total kernel which, in terms of the scaled position a measured from the position of the delta impurity, leads to where n(a, c, λ * ) = ∞ c du Ai(u + a) 2 + λ * π Im g e (a + u, u) 2 and c is the scaled relative position of the delta impurity with respect to the edge. In the above, when λ * = 0 we recover the usual edge density of the Airy gas which, in random matrix theory, corresponds to the eigenvalue density for the Gaussian Unitary Ensemble at the edge where the Wigner semi-circle law vanishes [28]. Using (35) the integrand in the second term can be written more explicitly as Let us recall the asymptotic behavior of the Airy functions. For u → +∞ one has which implies that Ai(u)Bi(u) 1 2π √ u , and so we see that for u → +∞ the r.h.s. of (123) behaves as −2πλ * Ai(u) 3 Bi(u) and thus decays very quickly. On the other side, for u → −∞ both Ai(u) and Bi(u) decay as 1/|u| 1/4 with oscillating prefactors. Hence it seems that the integral can be easily evaluated numerically (apart from a subtlety arising from the denominator for large negative λ * see below). Note also that as a → ±∞ the change in the density due to the impurity decays to zero.
In Fig. 7 we have plotted n(a, 0, λ * ) as a function of a for the values λ * = 0, the case where there is no impurity at the edge, and for value λ * = 1 (repulsive impurity) and the value λ * = −1 (attractive impurity) for an impurity placed exactly at the edge c = 0. We see that in all cases, the density oscillates in the region to the left of the edge (as one moves towards the bulk) but decays monotonically to the right as one moves away from the bulk. The presence of a repulsive impurity decreases the density, as expected, and induces a small phase shift in the oscillations to the left. However the attractive impurity increases the density at the original edge and leads to a larger density of fermions to the right, again with monotonic decay. In addition, the oscillations in the density experience a substantial phase shift with respect the the case of no impurity and a repulsive impurity. The presence of an impurity introduces a discontinuity in the derivative of the density at x 1 , i.e. at a = 0. Similar effects are seen when the impurity is not placed exactly at the edge x 1 = x e . For an attractive impurity placed on the right of the edge, the density increases around the impurity.
Filling transition of an attractive impurity located far from the edge. Here we examine, within the edge region, the filling transition already discussed in Section 5 in the context of the bulk. Clearly if the delta impurity is placed far to the right of the edge, i.e. c = x 1 −xe we 1, the local potential at the position of the impurity, V (x 1 ), is large compared to the value V (x e ) at the edge. Being well above the Fermi sea it should play no role in the Fermi gas, unless the amplitude of the delta impurity, λ = −λ A < 0, is tuned to be sufficiently attractive. Indeed, in that case we can revisit the qualitative argument given in Section 5. An attractive delta impurity in a uniform potential produces a bound state with a binding energy We can now surmise that when this energy is lowered below the Fermi energy µ = V (x e ), this bound state should be filled and be part of the ground state of the Fermi gas. If one equates this binding energy with the energy shift V ( We see that the real part of the denominator vanishes at u = u c = λ * 2 A /4 when λ * = −λ * A < 0, i.e. for an attractive impurity. To study the transition, from (125) we should consider λ A large, hence u c 1. In this case since the imaginary part is very small we can make the approximation where P denotes the Cauchy principle part. This means that the local density of states has a sharp resonance at u = u c . Inserting into (123) and (121) we see that the term which converges to a delta function gives a contribution Using Re[g e (a + u c , u c )] 1 2 √ uc e −2a √ uc for large u c , we find that the corresponding contribution to the local density reads This contribution was obtained under the assumption that u c 1. We see that it is non zero if u c > c, that is for λ * A > 2 √ c, exactly the same condition as (125). In that case the above contribution (129) corresponds precisely to a total of one particle, since its integral over a is equal to 1. Hence for λ * A > 2 √ c there is a local density peak corresponding to a total of one fermion. When λ * A < 2 √ c this extra fermion is no longer present. We see that this transition is very sharp for c 1 and thus coincides with the transition discussed in section 5. One can perform a slightly more precise estimate of the above formula (121) and (123) for the density near the impurity and obtain for aλ * A = O(1) and c 1, the Lorentzian dependence in the impurity position c near the transition at c = u c n(a, c, −λ * A ) where Hence the width η is exponentially small, i.e. η e − It is important to also compute the fermion density at the position of the impurity, e.g. to derive the effective potential in the next section. It is obtained by setting x 1 = x and thus a = 0 in (132). Recalling that the first term is the imaginary part of 1 π g e (u, u) we see that the formula simplifies into In the limit where c = (x 1 − x e )/w e → −∞, i.e. when the position of the impurity enters the bulk, one can easily check, using the explicit expression of w e in Eq. (34) and of λ * in Eq.
(115), that ρ µ (x 1 ) k F (x 1 )/π, independently of the sign of λ, i.e. both for a repulsive and an attractive impurity. This behavior matches perfectly with the behavior found in the bulk in Eqs. (63) and (66) in the limit |λ| k F . This is expected since that the edge scaling form in Eq. (132) holds for finite λ * , which implies λ 1/w e [see Eq. (115)], and thus |λ| k F (since w e 1/k F for large µ).

Effective potential at the edge
We now calculate the effective potential felt by an impurity in the edge region, as defined in (74). We can again use the Hellmann-Feynman theorem as in (80) which requires the density at the location of the impurity as given in (132). Since this density is expressed in terms of λ * it is convenient to write the Hellmann-Feynman formula in terms of λ * using (115). It reads This equation can then easily be integrated with respect to λ * using Eq. (132) to obtain with W edge (c, λ * ) = 1 π ∞ c du tan −1 λ * πAi(u) 2 1 + λ * πAi(u)Bi (u) . The function W edge (c, λ * ) can be evaluated numerically and is plotted in Fig. 8 as a function of c for several values of λ * . At large negative values of λ * the numerical evaluation becomes difficult, presumably due to the formation of a bound state about the impurity to the left of the edge. A more detailed analysis of this regime would be interesting to pursue and we leave this for future work. One can also verify, numerically, that Eq. (135) matches with the bulk form given in Eq. (13) as it should (see the discussion below Eq. (132)), although an analytic demonstration is not obvious given the highly oscillatory nature of the integrand.

Discussion
In this paper we studied non interacting fermions in a trap at zero temperature, in the presence of a singular potential created by delta function impurities. The presence of these impurities changes the density of the Fermi gas around the impurities. For a single impurity the change in the density profile has been studied using a number of techniques from condensed matter physics. These methods have allowed the characterization of the density at distances far from the impurity, which shows the celebrated Friedel oscillations. In this paper, using a Green's function method developed in our previous work we have computed the exact form of the density at all distances from the impurity. Furthermore our method goes beyond the one point function and also allowed to obtain the quantum correlations by computing the central object known as the kernel. In addition this allowed us to compute the effective potential felt by the impurity. We have shown how the behavior of the density and of the effective potential changes as one moves the impurity from the bulk of the Fermi gas to the edge created by the confining potential. We also unveiled an interesting "filling transition" which occurs when the impurity is moved outside of the support of the density of the Fermi gas. All these results are exact and non-perturbative in the strength of the impurity.
In addition when a pair of impurities is placed in the bulk of the Fermi gas at a distance r from each other, the fermion background gives rise to an effective interaction V int (r) between them, much like the Casimir effect in quantum electrodynamics. We have calculated exactly this effective interaction V int (r) at all distances, and our formula agrees with previous results known only for large distances.
In this paper all calculations in the presence of impurities are performed in the ensemble where the Fermi energy µ is fixed, and the system is in contact with a reservoir, so that the number of fermions can vary. This corresponds to the grand-canonical ensemble (here at zero temperature). In the Appendix B we briefly discuss the possible differences which may appear if instead one works in the canonical ensemble where the number of fermions is fixed (isolated system) as the impurity strength and position may vary.
We have focused here on the zero temperature limit, however it is important to derive the results at finite temperature since experiments are usually conducted at finite temperature. Indeed, the results derived here can be extended to finite temperature T in a straightforward way. As shown in [5,29] the kernel at finite temperature in the grand canonical ensemble at chemical potentialμ can be obtained from the zero temperature kernel, a relation which in the present framework can be written as Kμ(x, y) = 1 π dµ 1 1 + e β(µ −μ) Im G µ (x, y) , with β = 1/(k B T ). Using this expression, integral formulas can be obtained for all of the quantities studied in this paper. It would be challenging to analyse these formula in the future.
Another line of investigation would be to study the Wigner function [32,33] in the neighborhood of impurities both in the bulk and at the edge [34]. As well as being interesting in its own right, this might be a first step to understand the dynamics of systems in the presence of impurities as the Wigner function turns out to be a useful tool in the context of dynamics [35,36].
Finally, another interesting problem for further investigation is the question about a mobile impurity in a Fermi gas. This problem has been studied in several works, notably by McGuire [17,18]. It would be interesting to see if one could develop a general theory which extrapolates between the static impurity case studied here and the mobile impurity problem, which could explain the similarities between the two cases which we unveiled in Eq. (86).
A Comparison with the results of Ref. [12] In this appendix, we compare our exact result for the density in the presence of a deltafunction impurity in the bulk, given in Eq. (60), to the formula obtained in Ref. [12] by a quite different method. For this purpose, it is convenient to start from the formula given in Eq. (57). This formula can also be represented using the contour Γ 2 (see Fig. 3) which yields Note however that this representation is only valid for uniform systems, as it assumed that the bulk approximation for the Green's function is valid for small k, which in general is not.
The case λ > 0: in this case there is no bound state and there is no contribution to ∆K µ (x, y) in Eq. (137) coming from the part of Γ 2 along the negative imaginary axis. Assuming that the system is homogeneous, taking the imaginary part in Eq. (137) Using this representation when one sets x = y we obtain the formula of [12] for the change in the density -however we note that there is a factor of 2 difference as in [12] spin 1/2 fermions were treated.
The case λ < 0: When λ < 0 the integral over the contour Γ 2 in (137) picks up a half pole contribution at k = iλ, we thus find ∆K µ (x, y) = − λ π Im where the last term comes from the negative imaginary axis and corresponds to a bound state. The contribution from this bound state was in fact overlooked in Ref. [12]. In fact the formula Eq. (138) was computed in [12] via a direct summation of eigenfunctions, however when λ < 0 this formula misses the bound state which is introduced by an attractive impurity. Note however that the omission of this bound state does not affect the behavior at large distance from the impurity of the Friedel oscillations since the contribution of the bound state to the density decays exponentially at large distance.
The small λ limit: We note that taking the small λ limit in Eq. (60) gives, to order O(λ) where is the sine integral [24]. This matches perfectly with the linear response formula derived in [12].

B Canonical ensemble
We discuss here how one would approach the problem of adding impurities in the canonical ensemble where the number of particles if fixed and equal to N . Consider for example adding one impurity of strength λ. The number of single particle energy levels below the Fermi energy µ (i.e. the integrated density of states) is given by which is a function of λ. In order that N be fixed µ must be a function of λ, µ(λ) such that Let us define the change, due to the introduction of the impurity, in the integrated density of states as ∆N (µ, λ) = N (µ, λ) − N (µ, 0).
Given that the energy change is of order 1 and that the perturbation in the density is local it is clear that ∆N (µ, λ) is also of order 1 [30]. Denoting ∆µ = µ(λ) − µ(0) with µ(0) = µ we can rewrite (143) as N (µ + ∆µ, 0) + ∆N (µ + ∆µ, λ) = N = N (µ, 0), where ∆µ is the shift in the Fermi energy due to the impurity. To analyse what happens in the canonical ensemble one must carry out the computations in this paper at chemical potential µ + ∆µ, so the total fermion number is fixed upon adding the impurity. However, if ∆µ is zero, then the results in this paper can simply be applied to the canonical ensemble. A first example of where the canonical and grand canonical ensembles are equivalent is in a bulk system of volume V where one has a total particle number N = N (µ, 0) = k F (µ)V π . Now using Eq. (26) for a bulk system, k F (µ) = √ 2mµ/ , we see that µ = 2 π 2 N 2 2mV 2 and so ∆µ − 2 π 2 N ∆N (µ+∆µ,λ) mV 2 . From this we see that ∆µ → 0, since ∆N (µ + ∆µ, λ) is of order one, in the thermodynamic limit where N → ∞ and with N/V fixed. Note that a bulk system can have a varying periodic potential, and so the results given here are not just valid for constant potentials.
For a generic trap, we assume that for large µ one has N (µ, 0) = ( µ µ 0 ) z , where µ 0 is an intrinsic energy scale, and as N (µ, 0) must increase with µ we must have z > 0. Indeed, using the LDA in the bulk to compute N (µ, 0) as a function of µ for potentials of the form V (x) ∼ x p we find where B denotes the bulk region where µ − V (x) is real. Writing V (x) = v|x| p then gives