Phase separations induced by a trapping potential in one-dimensional fermionic systems as a source of core-shell structures

TL;DR

This paper explores how trapping potentials induce phase separation in one-dimensional ultracold fermionic gases, leading to diverse core-shell structures influenced by system parameters and phase diagram relations.

Contribution

It investigates the formation and properties of core-shell structures in trapped fermionic systems, highlighting the influence of trap shape, spin imbalance, and interaction strength on phase separation.

Findings

Core-shell structures can include BCS and FFLO states within harmonic traps.

Spatial separation of BCS and BEC regimes observed in the system.

External parameters significantly affect the phase configurations.

Abstract

Ultracold fermionic gases in optical lattices give a great opportunity for creating different types of novel states. One of them is phase separation induced by a trapping potential between different types of superfluid phases. The core-shell structures, occurring in systems with a trapping potential, are a good example of such separations. The types and the sequences of phases which emerge in such structures can depend on spin-imbalance, shape of the trap and on-site interaction strength. In this work, we investigate the properties of such structures within an attractive Fermi gas loaded in the optical lattice, in the presence of the trapping potential and their relations to the phase diagram of the homogeneous system. Moreover, we show how external and internal parameters of the system and parameters of the trap influence their properties. In particular, we show a possible occurrence of the core-shell structure in a system with a harmonic trap, containing the BCS and FFLO states. Additionally, we find a spatial separation of two superfuild states in the system, one in the BCS limit as well as the other one in the tightly bound local pairs (BEC) regime.

Equations

$$\mathcal{H}=\sum _{⟨i,j⟩\sigma }\left[-t-({\mu }_i+\sigma h){\delta }_{ij}\right]c_{i\sigma }^{†}{c}_{j\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$$n_{i\uparrow }{n}_{i\downarrow }={\Delta }_i^{\ast }{c}_{i\downarrow }{c}_{i\uparrow }+{\Delta }_i{c}_{i\uparrow }^{†}{c}_{i\downarrow }^{†}-|{\Delta }_i{|}^2,$$

$${\mathcal{H}}^{MF}=\sum _{⟨i,j⟩\sigma }\left[-t-({\mu }_i+\sigma h){\delta }_{ij}\right]c_{i\sigma }^{†}{c}_{j\sigma }+U\sum _i\left({\Delta }_i^{\ast }{c}_{i\downarrow }{c}_{i\uparrow }+H.c.\right)-U\sum _i|{\Delta }_i{|}^2.$$

$$c_{i\sigma }=\sum _n\left(u_{in\sigma }{\gamma }_{n\sigma }-\sigma v_{in\sigma }^{\ast }{\gamma }_{n\bar{\sigma }}^{†}\right),$$

$$\displaystyle\sum_{j}\left(\begin{array}[]{cc}H_{ij\sigma}&\Delta_{ij}\\ \Delta_{ij}^{\ast}&-H_{ij\bar{\sigma}}^{\ast}\end{array}\right)\left(\begin{array}[]{c}u_{jn\sigma}\\ v_{jn\bar{\sigma}}\end{array}\right)=\mathcal{E}_{n\sigma}\left(\begin{array}[]{c}u_{in\sigma}\\ v_{in\bar{\sigma}}\end{array}\right)$$

$${\Delta }_i=⟨c_{i\downarrow }{c}_{i\uparrow }⟩=\sum _n\left[u_{in\uparrow }{v}_{in\downarrow }^{\ast }f({\mathcal{E}}_{n\uparrow })-u_{in\downarrow }{v}_{in\uparrow }^{\ast }f(-{\mathcal{E}}_{n\downarrow })\right].$$

$$\Omega \equiv -k_{B}T\sum _{n\sigma }\ln \left[1+\exp \left(\frac{-{\mathcal{E}}_{n\uparrow }}{k_{B}T}\right)\right]-\sum _i\left({\mu }_i+h+U|{\Delta }_i{|}^2\right).$$

$${\Delta }_i=\frac{1}{N_x}\sum _{\mathbf{q}}{\Delta }_{\mathbf{q}}\exp \left(i{\mathbf{R}}_i\cdot \mathbf{q}\right),$$

Full text

Introduction

Recently, an extensive progress in experimental techniques in ultracold quantum gases in optical lattices occurs. These systems can be give as the greatest examples of practical realization of “quantum simulators” [1, 2, 3, 4, 5]. Because of the fact that they are systems with fully controllable parameters, they open a new avenue to study fundamental phenomena of condensed matter physics. Moreover, such simulators can provide information about properties of physical systems in a context of different effects and mechanisms, which are difficult to observe in solid state materials due to their complexity. An experimental realization of such systems with bosonic and fermionic gases in optical lattices have already been conducted [6, 7, 8, 9, 10]. It creates a very attractive field for further development. Additionally, in these systems trapping potential and lattice geometry can be modified to study new exotic phases. Hence, ultracold gases in optical lattices allow to simulate the well-known Hubbard model in various regimes of parameters [11].

The lattice geometry is a fundamental characteristic of many-body systems and has large influence on their physical properties [12]. Ultra-cold atomic gases give the opportunity of realization of systems with different geometries. Experimentally, the geometry of the lattice can be changed by different spatial arrangement of laser beams [1]. Recently, the occurrence of antiferromagnetic spin correlations in the repulsive fermionic gas, in different lattice geometries of varying dimensionality, also including crossover configurations between different geometries, was investigated [13, 14, 15, 16, 17, 18, 19].

Equally important is the better comprehension of unconventional phases that can appear in fermionic superfluids with population imbalance, the latter being the effect of a magnetic (Zeeman) field or the result of preparing a mixture with the desired composition. In the weak coupling limit, states with nontrivial Cooper pairs can exist at large population imbalance [20, 21, 22, 23, 24]. One of examples of such pairing is the Fulde–Ferrell–Larkin–Ovchinnikov (FFLO) state [25, 26], in which the Cooper pairs have non-zero total momentum as a result of pairing across the spin-split Fermi surface. The properties of this state attracted a lot of theoretical and experimental attention [27, 28, 29, 30, 31]. Another example of unconventional coherent state is the homogeneous spin-polarized superconductivity with a gapless spectrum for the majority spin species [32]. For this phase, the coexistence of the normal and the superfluid states in the isotropic phase is characteristic. This phase was firstly proposed by Sarma [33] who studied the case of a superconductor in an external magnetic field within the BCS theory, neglecting orbital effects. He showed that self-consistent mean field solutions with gapless spectrum [$\Delta(h)$] are energetically unstable at $T=0$, in contrary to the fully gapped BCS solutions. On the other hand, a spin-polarized superconducting state can be stabilized by non-zero temperature.

Experimental studies of the spin–imbalanced fermionic gases give new possibilities for research in the field of condensed matter systems with strong correlations. There are experimental evidences for the occurrence of core-shell structure — an unpolarized superfluid core in the center of the trap surrounded by a polarized normal state [34, 35, 36]. The structure has been observed in the density profiles of trapped spin–imbalanced fermionic mixtures. In this system a phase separation between these two states appears [37, 38, 39, 40, 41, 42, 43, 2, 44, 45]. Additionally, in the case of two-component fermionic gases in one-dimensional optical lattice, the exact thermodynamic Bethe ansatz solution indicates the occurrence of a mixed phase with two-shell structure [3]. In such a state, a partially polarized superfluid core (the FFLO phase) is surrounded by a fully paired (BCS-like) or fully polarized (normal) phases [46, 47]. Similar observation has been also performed in the case of a one-component trapped gas [48]. It is impotant to emphasize that the one-dimensional system is a good candidate for observing the FFLO phase because of a nesting effect, which makes the state much more robust than in the three-dimensional case [11, 12].

In this paper we show that the occurrence of core-shell structures is a consequence of the presence of inhomogeneity in the system, in particular, of the changes of chemical potential or magnetic field (or equivalently effective spin-imbalance) depending on the trapping potential and the short-range interactions between atoms. Moreover, we provide an analysis according to which, depending on the attractive interaction, multiple core-shell structures can appear in the system, including different phases, such as the spatially homogeneous spin-unpolarized superconducting state, i.e., the BCS state, as well as the spatially inhomogeneous superconducting FFLO state. We show that it is possible to prepare the system in such a way that one can observe the two different phases separately in space and we study the influence of the trapping potential on such spatially separated phases.

Model and method

In this paper, we study a one-dimensional system with the $s$-wave superconductivity. The attractive Hubbard Hamiltonian (i.e., with on-site pairing, $U<0$) in presence of the magnetic field ($h$) has a following form:

[TABLE]

where $c_{i\sigma}^{\dagger}$ ($c_{i\sigma}$) denotes the creation (annihilation) operator of the particles at site $i$ and spin $\sigma=\{\uparrow,\downarrow\}$. $t$ is the hopping integral between nearest-neighbor sites, $U<0$ is the on-site pairing interaction, $h$ is the external magnetic Zeeman field, while $\mu_{i}$ is an effective on-site chemical potential. The interaction term is treated within the mean-field broken-symmetry approximation:

[TABLE]

where $\Delta_{i}=\langle c_{i\downarrow}c_{i\uparrow}\rangle$ is the s-wave superconducting order parameter (SOP). Hence, the mean-field Hamiltonian in real space has the form:

[TABLE]

For a general case (i.e., for any distribution of $\mu_{i}$), Hamiltonian (1) can be exactly diagonalized within the Bogoliubov–Valatin transformation:

[TABLE]

where $\gamma_{n\sigma}$ and $\gamma_{n\sigma}^{\dagger}$ are the new quasi-particle fermionic operators, whereas ${\bm{u}}$ and ${\bm{v}}$ are the Bogoliubov–de Gennes (BdG) eigenvectors. One can obtain the self-consistent BdG equations in real space:

[TABLE]

where $H_{ij\sigma}=-t\delta_{\langle i,j\rangle}-(\mu_{i}+\sigma h)\delta_{ij}$ is the single-particle Hamiltonian and $\Delta_{ij}=U\Delta_{i}\delta_{ij}$ are the on-site SOPs. Using transformation (4), the SOPs can be found as:

[TABLE]

Equations (11) can be solved self-consistently with respect to the distribution of $\Delta_{i}$. In this case, one can find the grand canonical potential for a given state as:

[TABLE]

From several solutions of the BdG equations, only those with a minimal value of grand canonical potential $\Omega$ (at fixed $\mu$ and $h$) indicate a thermodynamically stable state of the system. In the absence of a trap (i.e., $\mu_{i}=\mu$ at each site), the distribution of $\Delta_{i}$ can be rewritten in momentum space, using the Fourier transform:

[TABLE]

where ${\bm{q}}$, which are restricted to the first Brillouin zone, are total momenta of the Cooper pairs.

Numerical results

In this section, the obtained numerical results are presented. First, the homogeneous system (without the trap) is investigated. The magnetic field $h$ versus chemical potential $\mu$ phase diagram is determined. Next, we show that phase separation occurs in the presence of a harmonic trap. Such a phase separation induced by a trapping potential we will call (artificially) enforced phase separation further in the text. For chosen parameters (i.e., $\mu$ and trapping potential), the phase separation between different phases (e.g., a state with the BCS core and the FFLO shell as well as the FFLO core and the BCS shell) can be realized in the system.

Phase diagram: homogeneous system

Numerical calculations presented in this section have been performed for a one-dimensional chain with the periodic boundary conditions and $N=200$ sites. For the homogeneous system, i.e., $\mu_{i}=\mu=\text{const.}$, the magnetic field $h$ versus chemical potential $\mu$ phase diagram is shown in Fig. 1. The results for different values of $U/t$ are presented (cf. also Ref. 24 and 49). In each case, they consist of three regions. For low values of magnetic field, the conventional BCS phase is stable (with $\Delta_{i}=\text{const}$ and $\Delta_{i}\neq 0$). With increasing $h$, one finds a discontinuous phase transition from the BCS to the FFLO phase (with $\Delta_{i}$ changing from site to site). In the third region the normal (non-ordered, NO) phase is stable (with $\Delta_{i}=0$). The transition from the BCS to the NO phase is continuous (for changing $h$ and fixed $\mu$) or discontinuous (at the vertical boundary in the phase diagram). Notice also that for small $\mu/t$ and large $h/t$, the so-called $\eta$ phase can exist [50] (the FFLO phase with maximal $q=\pi$, at least for the FF ansatz [24, 49]). With increasing $U/t$, the region of the BCS phase extends, whereas the continuous FFLO–NO boundary weakly depends on $U/t$. In Fig. 1(a), the solid lines correspond to the case in which only the FF state is considered, i.e., all Cooper pairs have the same momentum ${\bm{q}}$ [25]. Including the fact that generally in the “full” FFLO state, Cooper pairs with different momenta can contribute (i.e., given by Eq. (12), one gets that the region of the FFLO phase occurrence is slightly extended in comparison to the one obtained only for the FF ansatz. This is an expected and well known relation between superconducting phases with different numbers of allowed ${\bm{q}}$, i.e., a phase with larger number of ${\bm{q}}$’s is more stable than the FF phase with only one ${\bm{q}}$ [27, 51].

It should be also mentioned that for low filling (i.e., for $|\mu|\sim 2t$), the BCS–BEC crossover can be realized. The detailed discussion of this issue has been presented in Ref. 24. In the context of the present work (particularly for the problem of the phase separation region in the system), the following property of the phase diagram is also very important. Namely, depending on values of the magnetic field $h$ and on-site attraction $U$, with increasing $|\mu|$ the transitions from the BCS phase to the FFLO phase or from the FFLO phase to the BCS phase can occur [red and blue stars, respectively, in Fig. 1(b)].

However, it is important to keep in mind that the use of the mean-field approximation is generally restricted to the weak coupling limit and ground state properties. The limitations of the mean-field method affect the one-dimensional system the most because the pair fluctuations become very important in this case [52]. For such system geometry, the nature of phase transition between the BCS and FFLO phases can be faultily predicted. However, the mean-field approximation can give some useful description in the weak and intermediate couplings, which are comparable with the Bethe ansatz results [52]. While the mean-field FF-type calculations do not predict the correct type of the phase transition from the BCS phase to the FFLO phase, the self-consistent Bogoliubov–de Gennes results are in a good agreement with those obtained from the Bethe ansatz. In the continuum model, the first order phase transition is simply an artifact of the FF-type calculation [52]. This suggests a similar problem in the case of a lattice in the thermodynamic limit, when average particle concentration and lattice constant go to zero (i.e., $n\rightarrow 0$ and $a\rightarrow 0$ — the dilute limit).

Phase separation and its types

Generally, the phase separation is a state of the system where two or more uniform phases (e.g., those which have been defined previously) occur in different parts (so-called domains) of the system. In this section, we distinguish two different types of phase separation, which can emerge in the system, i.e.: (a) a spontaneous phase separation, which can occur in a homogeneous system and (b) an artificially enforced phase separation, which can emerge in inhomogeneous system. In the present work these cases correspond to the system with $\mu_{i}=\text{const}$ at each site and to the system with inhomogeneous spatial distribution of $\mu_{i}$, respectively. Below, we characterize briefly these two types of phase separations.

Spontaneous (macroscopic) phase separation

The discontinuous transitions between two (homogeneous) phases in the diagram, as a function of $\mu$ are usually related to the discontinuous change of particle concentration from value $n_{+}$ to value $n_{-}$ ($n_{+}>n_{-}$). Such a transition can be associated with the occurrence of the (macroscopic) phase separation in a defined range $n_{-}<n<n_{+}$ of particle concentration [53]. In such a phase-separated state two domains with different particle concentrations $n_{-}$ and $n_{+}$ coexist (there can be also regions differing in the magnitude of the order parameter as well as thermodynamic phases). In this approach, because of neglecting the interface energy at the boundaries of the domains, such states can exist only in the thermodynamic limit (i.e., when $N\rightarrow\infty$) [54, 55, 56, 57, 58]. In a finite system, the interface energy can lead to an occurrence of states with other textures [59, 60, 61, 62, 63], besides the homogeneous states and the phase separated states discussed above. Due to the fact that the BCS–FFLO, BCS–NO, and FFLO–NO boundaries in Fig. 1(a) can be discontinuous in the diagram as a function of $n$, one expects an occurrence of the following (macroscopic) phase separation regions: BCS/FFLO, BCS/NO, and FFLO/NO [24, 49]. Notice that only the BCS–FFLO boundary is discontinuous for all model parameters. To distinguish the phase separation from other states discussed above, we will call it a spontaneous one, because it can occur spontaneously in the homogeneous system.

Artificially enforced phase separation: a model example

From the analysis of the phase diagram and shapes of boundaries between the phases, it is clear that for $U/t=-1.5$ and $h/t=0.125$ (set $\mathcal{A}$ of the model parameters), the phase with the lowest energy is (a) the FFLO phase at $\mu=\mu_{1}=0.75t$ and (b) the BCS phase at $\mu=\mu_{2}=1.75t$ [the red stars in Fig. 1(b)]. The SOP at each site for these two solutions is presented in Fig. 2(a) and Fig. 2(b). Analogously, we also choose the following parameters for $U/t=-3.0$ and $h/t=0.55$ [set $\mathcal{B}$ of the parameters, the blue stars in Fig. 1(b)]: $\mu_{1}=0.75t$ (the BCS phase) and $\mu_{2}=1.5t$ (the FFLO phase), and present the local dependence of the SOP in Fig. 3(a) and Fig. 3(b). One can notice that in the FFLO phase $\Delta_{i}$ changes from site to site periodically, whereas for the BCS solution $\Delta_{i}$ is homogeneous in the whole system.

Now, we investigate the system with a particular distribution of the chemical potential: $\mu_{i}=\mu_{1}$ in one half of system (favoring the FFLO phase or the BCS phase for set $\mathcal{A}$ or $\mathcal{B}$, respectively) and $\mu_{i}=\mu_{2}$ in the other half (favoring the BCS phase or the FFLO phase for set $\mathcal{A}$ or $\mathcal{B}$, respectively). After solving the self-consistent set of the BdG equations in real space for such a system with $N=200$ sites (and with open boundary conditions), one gets that solutions obtained in parts of the system resemble the ones for the homogeneous system [Fig. 2(c) and Fig. 3(c)]. In such a case, we have a coexistence of two homogeneous solutions, which are spatially separated. Notice that small changes (in comparison to the solutions for the homogeneous system in each part of the system) are distinguishable only in the neighborhood of interfaces between two such defined domains. Thus, one can conclude that the effects associated to the interfaces are almost irrelevant for the state under consideration. Moreover, the spin polarization (defined as the difference of concentrations of particles with spin up and spin down, $m_{i}=n_{i\uparrow}-n_{i\downarrow}$) as a function of site $i$ in these inhomogeneous states is shown in Fig. 2(d) and Fig. 3(d), for $\mathcal{A}$ and $\mathcal{B}$, respectively. As can be expected, the spin polarization has a modulation, which is two times faster than that of the SOP. Additionally, maximal values of the spin polarization are located at the nodal points of the SOP [64, 65]. This is a consequence of the interplay between the unpaired polarized particles and those in the superconducting state (i.e., the Cooper pairs). A small number of pairs (at sites with $\Delta_{i}\approx 0$) supports the occurrence of polarized states (i.e., a large value of $m_{i}$). Note also that there is no coexistence of the superconducting and magnetic ordering in the BCS phase.

The nature of this inhomogeneous state is, however, different than the origin of the spontaneous (macroscopic) phase separation discussed in the previous point. In the present setup, the system parameters are inhomogeneous (different values of chemical potential in two parts of the system). We will call such separated state an (artificially) enforced phase separation in contrary to the spontaneous phase separation, which could occur in the homogeneous system.

Note, that even if $\mu_{1}=\mu_{2}=\mu_{c}$ (where $\mu_{c}$ is a critical value of the chemical potential at the FFLO–BCS boundary for particular $U/t$ and $h/t$), due to the fact that the interface has small but still finite energy, the finite system cannot exhibit spontaneous phase separation. In such a case, the whole system is in the FFLO phase or in the BCS phase (both states have equal energy and both solutions correspond to local minima of the grand canonical potential).

System with a trap

After studying the model systems, let us investigate more realistic situations, in which the on-site potential $\mu_{i}$ changes from site to site. An example of experimental realization of such systems are ultracold atomic gases in optical lattices. In the following, we consider two types of traps: (i) a linear trap with $\mu_{i}=V_{0}|r_{i}-r_{c}|$ and (ii) a harmonic trap with $\mu_{i}=V_{0}(r_{i}-r_{c})^{2}$, where $r_{i}$ is a location of $i$-th site, whereas $r_{c}$ is a location of the center of the trap (i.e., $i=N/2=500$). In this part, we performed the calculations for $N=1000$ sites with open boundary conditions. $V_{0}$ was chosen in such a way to ensure a value much larger than the chemical potential corresponding to the BCS-NO transition in the homogenous system at $h=0$ and fixed $U$.

The self-consistent solutions of the BdG equations for the case of a linear trap are presented in Fig. 4. In this case, we show SOP $U\Delta_{i}$, magnetization $m_{i}$, particle concentration $n_{i}$, and the trapping potential $\mu_{i}\equiv V(i)=V_{0}|r_{i}-r_{c}|$ (with $V_{0}=6.0t/N$) as a function of the lattice site $i$ (blue, green, red and dashed black lines, respectively). Due to the fact that $\mu_{i}$ changes from site to site, $n_{i}$ changes at each site and one does not observe the solutions with, e.g., constant value of $\Delta_{i}$, corresponding to the BCS phase (but the region without oscillations of the SOP can be identified as corresponding to this phase, see below). However, one can indicate the regions along the chain which correspond to the phases indicated in Fig. 1(a). Namely, starting from the left side of Fig. 4(a) (for $U/t=-1.5$ and $h=0.125t$, set $\mathcal{A}$), one can find a region of the NO (empty) phase, where $\Delta_{i}=0$ (and $n_{i}=0$). Next, there is the BCS region, where $\Delta_{i}\neq 0$ exhibits no oscillations (and $m_{i}=0$). The oscillations of $\Delta_{i}$ and nonzero magnetization $m_{i}\neq 0$ are clearly visible in the center of the system. Such a behaviour of $\Delta_{i}$ indicates the occurrence of the FFLO phase in this part of the system. Hence, in the linear trapping potential, the FFLO phase in the core (i.e., in the center of the system) is surrounded by the shell of the BCS states. For $U/t=-3.0$ and $h=0.55t$ (set $\mathcal{B}$), the situation is analogous, but going to the center of the trap, the following sequence of regions is presented: the NO phase [with $n_{i}=0$ (vacuum state), as well as filled state with $n_{i}\neq 0$], the FFLO phase, and the BCS phase (as a core). The NO phase with $n_{i}\neq 0$ is fully polarized (i.e., $m_{i}/n_{i}=1$).

Note, that the location $i$ of the boundaries between the regions along the chain corresponds, approximately, to the values of $\mu_{i}$ for which the phase transition occurs in the homogeneous system [cf. Fig. 1(b)]. It is attributed to the fact that the interactions in the system are short-ranged and the interface energy is relatively small. Moreover, for set $\mathcal{A}$, only the NO phase with $n=0$ can exist, whereas for set $\mathcal{B}$ also the NO phase with $n\neq 0$ is possible (cf. Ref. 24).

The results for a harmonic trap, presented in Fig. 5, do not differ qualitatively from those obtained for a linear trap. We take $\mu_{i}\equiv V(i)=V_{0}(r_{i}-r_{c})^{2}$, with $V_{0}=4.0t/(N/2)^{2}$. Such chosen values correspond to a relatively flat trap. The sequences of the regions for both sets of the model parameters ($\mathcal{A}$ and $\mathcal{B}$) are NO–BCS–FFLO and NO–FFLO–BCS, respectively, with spatially extended regions in the center (cores) and spatially shrunken external regions (shells).

The results presented above show that in the system with a trap at each site, the state of a homogeneous system corresponding to $\mu=\mu_{i}$ is approximately realized. The trap shapes considered above are chosen in such a way to change $\mu_{i}\equiv V(i)$ from a range that all of possible phases occurring in the system are present, at fixed $h$ cf. Fig. 1 (from the NO phase at large $|\mu|$ at the boundaries of the system, i.e., at $i=0,N-1$, to the half-filled phases at $\mu=0$ in the center, i.e., at $i=N/2$). Thus, in other words, all phases for a homogeneous system at fixed $\mu$ can be realized in the setup with the trap. However, by modification of parameters of the trap for instance, one can realize a system which corresponds only to a part of the phase diagram (i.e., some limited part of $\mu$ range). A sample of the results is shown in Fig. 6. They are obtained for the system with a harmonic trap in the form of $\mu_{i}=V_{0}(r_{i}-r_{c})^{2}+\bar{V}$, with $V_{0}=2.0t/(N/2)^{2}$ and $\bar{V}=1.25t$. In this case, the minimum value of the trapping potential is shifted to some finite value, i.e., $\mu_{N/2}=\bar{V}$ (black doted line). The chosen parameters allow to realize a part of the phase diagram in which two superconducting domains are separated by the NO state [cf. Fig. 1(b)]. Depending on a fixed value of magnetic field $h$, the states with BCS or BCS-FFLO core surrounded by the BCS shell, separated by the NO sub-shell-part can be realized in the system.

The results shown in Fig. 6 correspond to a part of the phase diagram [cf. Fig. 1(b)], where the BCS–BEC crossover is realized. One should notice that the system under consideration is in a trap, which can change the state at the site with $\mu_{i}$ in comparison to the homogeneous solution with $\mu=\mu_{i}$. Indeed, numerical results show that this state can be realized outside the BCS core state [Fig. 6(a)], for $h$ smaller than the one at which the FFLO phase occurs in a homogeneous system. The realization of the two BCS-like domains, separated by the NO state, is possible by the decrease of $h$ [Fig. 6(b)]. At the boundary of the internal BCS core, one observes the FFLO-type gap oscillation with changing of its sign [66, 67, 68]. It is important to emphasize that, as a consequence of low filling at the external BCS shell state, we expect realization of the BEC condensate (the tightly bound local pairs region) of the Cooper pairs. According to the Leggett criterion [69], the BEC (i.e., Bose–Einstein condensate) begins when the effective chemical potential is smaller than the lower band edge.

Summary

In conclusion, we have shown that properties of the system in the presence of a trap are strongly associated to the phase diagram of the homogeneous system (in the absence of a trap). It is clearly seen from the spatial dependence of different quantities in the system with linear and harmonic trap (Figs. 4 and 5, respectively) and provides a confirmation about the validity of so-called local density approximation, that is very often employed in theoretical calculations for trapped system starting from previous results for homogeneous systems [70, 48, 22]. We have shown that in such a system the core-shell structures can be created due to the trapping potential. For chosen parameters, the order of states realized in such structures corresponds to the sequence of phases occurring in the phase diagram for the homogeneous system (Fig. 1). The states occurring in a particular sequence depend mainly on values of model parameters (mainly $h$ and $U$). Generally, the shape of the trapping potential does not change the structure of phases occurrence. Moreover, even in the same type of trap, for different values of the model parameters (i.e., interaction $U$ or magnetic field $h$), the structure with the BCS core and the FFLO shell as well as with the FFLO core and the BCS shell can be realized (depending on $U$ and $h$).

An increase of the on-site interaction at low filling can lead to the BCS–BEC crossover in the system. By tuning of the trap parameters, we can realize the part of the phase diagram in which this crossover occurs. Thus, one can obtain the more complex core-shell structures. In particular, one can obtain two BCS states separated by the NO state (e.g., Fig. 6). Due to the low filling, at the outer BCS shell, the BEC should emerge in the system.

In this work, we have shown that different core-shell structures can occur in a trapped system. These structures are examples of so-called (artificially) enforced phase separation, occurring in spatially inhomogeneous systems, whose origin is different than that of the macroscopic phase separation in homogeneous systems. It is important to emphasize that such experimental setup with a trap allows to investigate phase diagrams of homogeneous systems with short-range interactions. The theoretical prediction presented in this work should be realizable experimentally in a relatively simple way.

Acknowledgements

We thank Krzysztof Cichy for careful reading of the manuscript, valuable comments and discussions. A.P. is grateful to Laboratoire de Physique des Solides (CNRS, Université Paris-Sud) for hospitality during a part of the work. This work was supported by the National Science Centre (Narodowe Centrum Nauki, NCN, Poland) under grants nos.: UMO-2017/24/C/ST3/00357 (A.C.), UMO-2017/24/C/ST3/00276 (K.J.K.), and UMO-2017/25/B/ST3/02586 (A.P.).

Author contributions statement

A.P. initialized and coordinated the project. A.P. performed numerical calculation. All authors discussed the results. A.C. and K.J.K. prepared the first version of the manuscript. All authors contributed to the final version of the manuscript and all of them reviewed and accepted it.

Competing interests

The authors declare no competing interests.