Assembly of 2N entangled fermions into multipartite composite bosons

TL;DR

This paper investigates how groups of 2N fermions can form composite bosons, proposing an entanglement-based method to evaluate their bosonic behavior and studying the ground state transitions in a Hubbard model.

Contribution

It introduces a novel entanglement-based approach to determine the bosonic nature of fermionic arrangements and applies it to a Hubbard model analysis.

Findings

Entanglement measures can assess bosonic quality of fermionic composites.

Ground state properties change with interaction strength in the Hubbard model.

Multiple fermionic arrangements influence the emergent bosonic behavior.

Abstract

An even number of fermions can behave in a bosonic way. The simplest scenario involves two fermions which can form a single boson. But four fermions can either behave as two bipartite bosons or further assemble into a single four-partite bosonic molecule. In general, for 2N fermions there are many possible arrangements into composite bosons. The question is: what determines which fermionic arrangement is going to be realized in a given situation and can such arrangement be considered truly bosonic? This work aims to find the answer to the above question. We propose an entanglement-based method to assess bosonic quality of fermionic arrangements and apply it to study how the ground state of the extended one-dimensional Hubbard model changes as the strength of intra-particle interactions increases.

Equations

$$|\psi ⟩=\sum _{i_1,i_2,\dots ,i_{2N}}{\alpha }_{i_1,i_2,\dots ,i_{2N}}{a}_{i_1}^{†}{a}_{i_2}^{†}\dots a_{i_{2N}}^{†}|0⟩,$$

$$|\psi ⟩=\sum _{i,j}{\alpha }_{i,j}{a}_i^{†}{a}_j^{†}|0⟩,$$

$$|\phi ⟩=\sum _{i,j}{\alpha }_{i,j}U{a}_i^{†}U{a}_j^{†}|0⟩=\sum _k\sqrt{{\lambda }_k}{a}_{2k}^{†}{a}_{2k+1}^{†}.$$

$$|\phi ⟩=\sum _{k=0}^{d-1}\sqrt{{\lambda }_k}{a}_k^{†}{b}_k^{†}|0⟩.$$

$$P=\sum _k{\lambda }_k^2.$$

$$\sum _{k=0}^{d-1}\sqrt{{\lambda }_k}{a}_k^{†}{b}_k^{†}|0⟩\equiv c^{†}|0⟩=|1⟩,$$

$$⟨0|c^N{c}^{†N}|0⟩=N!.$$

$$⟨0|c^N{c}^{†N}|0⟩={\chi }_{N}N!,$$

$$|N⟩\equiv \frac{c^{†N}}{\sqrt{{\chi }_{N}N!}}|0⟩$$

$$c^{†}|N-1⟩$$

$$c|N⟩$$

$${\alpha }_N=\sqrt{\frac{{\chi }_N}{{\chi }_{N-1}}}$$

$$⟨{\epsilon }_N|{\epsilon }_N⟩=1-N\frac{{\chi }_N}{{\chi }_{N-1}}+(N-1)\frac{{\chi }_{N+1}}{{\chi }_N}.$$

$$1-NP\le \frac{{\chi }_N}{{\chi }_{N-1}}\le 1-P,$$

$${\chi }_2=2\sum _{i<j}{\lambda }_i{\lambda }_j=\sum _{i,j}{\lambda }_i{\lambda }_j-\sum _i{\lambda }_i^2=1-P.$$

$$c^{†}|0⟩=\frac{1}{\sqrt{d}}\sum _{k=0}^{d-1}{a}_k^{†}{b}_k^{†}|0⟩,$$

$${\chi }_N=\frac{d!}{d^N(d-N)!}$$

$${\alpha }_N=\sqrt{\frac{d-N+1}{d}},⟨{\epsilon }_N|{\epsilon }_N⟩=0.$$

$$c_{s,r}^{†}|0⟩=\frac{1}{\sqrt{d}}\sum _{k=0}^{d-1}{e}^{i\frac{2\pi }{d}kr}{a}_k^{†}{b}_{k+s}^{†}|0⟩,$$

$$⟨0|c_{s,r}{c}_{s^{\prime },r^{\prime }}^{†}|0⟩={\delta }_{s,s^{\prime }}{\delta }_{r,r^{\prime }},$$

$$\sum _{s,r=0}^{d-1}{\alpha }_{s,r}{c}_{s,r}^{†}|0⟩$$

$$\mathcal{H}=J{\mathcal{H}}_0+U{\mathcal{H}}_p$$

$${\mathcal{H}}_0=-\sum _{i=0}^{d-1}\left(a_i^{†}{a}_{i+1}+b_i^{†}{b}_{i+1}+h.c.\right)$$

$${\mathcal{H}}_p=-\sum _{i=0}^{d-1}{a}_i^{†}{a}_i{b}_i^{†}{b}_i$$

$$\mathcal{H}|s,r⟩=$$

$$P=\underset{d\to \infty }{\lim }\frac{1}{d}=0,$$

$$|{\psi }_0⟩=\sum _{s=-\infty }^{\infty }{\alpha }_s|s,0⟩,$$

$$\mathcal{H}|{\psi }_0⟩=\epsilon |{\psi }_0⟩,$$

$$-\frac{\epsilon }{2J}{\alpha }_s={\alpha }_{s+1}+{\alpha }_{s-1},$$

$$-\frac{(U+\epsilon )}{2J}{\alpha }_0={\alpha }_1+{\alpha }_{-1}.$$

$${\alpha }_s=A{\left(\frac{\sqrt{U^2+16J}-U}{4J}\right)}^{|s|},$$

Full text

I Introduction

Most bosons studied in laboratories are in fact composed of elementary fermions. It is therefore important to understand what makes them behave in a bosonic way. One may think that it is the binding force that keeps them together, however recently an intriguing idea was proposed that the bosonic behaviour stems from the intra-fermionic entanglement. More precisely, the more entangled two fermions are, the more bosonic they behave Law ; Wootters . This result established a new field of research, the quantum information oriented studies on composite bosons pawel2011 ; Bouvrie2012a ; Bouvrie2012b ; Bouvrie2013 ; Bouvrie2014 ; Bouvrie2016 ; Bouvrie2017 ; Bouvrie2018 ; Pawel2012 ; Pawel2013 ; Pawel2015 ; zak2017 ; zak2018 . Up to now scientists focussed on composite bosons made of two elementary components, whose quality can be quantified by bipartite entanglement. Here, we propose a generalization to composite bosons made of $2N$ fermions, whose quality is described by genuine multipartite entanglement.

The main problem of this work can be formulated in the following way. Suppose that $2N$ fermions are in the state

[TABLE]

where $a^{\dagger}_{i_{k}}$ creates a fermion in the mode $i_{k}$ and $\{\alpha_{i_{1},i_{2},\ldots,i_{2N}}\}$ is the set of antisymmetric coefficients. We ask: does the state $|\psi\rangle$ describe a single composite bosonic particle made of $2N$ fermions? Or perhaps it describes two bosonic particles, each made of $N$ fermions? Or maybe it describes $N$ bipartite bosonic particles? In fact, for $2N$ fermions there is a number of possible bosonic assemblies. How to decide to which assembly the state $|\psi\rangle$ corresponds to and how to quantify its bosonic quality? We are going to address the above questions. In particular, we will focus on fermionic states that are ground states of the one-dimensional Hubbard model. We will study how assemblies of $2N$ spin-$1/2$ particles on a lattice, and their corresponding bosonic qualities, depend on the strength of intra-particle interactions.

Apart from contributing to the new field of research, there are two additional motivations behind our studies. The first one is fundamental and is related to general investigations on complexity in the quantum domain. In particular, we want to understand what makes a complex quantum system to behave as a single entity. The second one is more pragmatic. We would like to understand how to engineer and control the creation of composite particles. This is related to the following problem – in some situations spontaneous emergence of composite structures can affect the property one wants to observe. A very illustrative example is the problem of the Bose-Einstein condensation of atomic Hydrogen Fried1998 . Hydrogen atoms naturally try to recombine into $H_{2}$ molecules and one needs to find out a method to prevent it (for example by spin polarization in high magnetic fields Silvera1980 ). Here, we consider a system of $2N$ interacting fermions on a one-dimensional lattice and show that the formation of composite structures can be controlled by a proper tuning of the nearest neighbour interaction.

II Bipartite composite bosons

In this section we recall properties of composite bosons made of two fermions and show that such composite bosons naturally describe ground states of the Hubbard model for two fermionic particles.

II.1 Bipartite composite bosons and the role of correlations

Consider a system made of two fermions whose general state is given by

[TABLE]

where the matrix of coefficients $\alpha_{i,j}$ satisfies $\alpha_{i,j}=-\alpha_{j,i}$, $\alpha_{i,i}=0$ and $\sum_{i,j}|\alpha_{i,j}|^{2}=1$. The above bi-fermionic state admits a Schmidt-like (Slater) representation Slater1 ; Slater2 ; Slater3 ; Eckert2002 , i.e., there exist a unitary transformation $Ua^{\dagger}_{i}=\sum_{k}\beta_{i,k}a^{\dagger}_{k}$, such that

[TABLE]

The above representation allows to divide the fermions into two different groups, the ones occupying even modes and the ones occupying odd modes. We therefore set $a^{\dagger}_{2k}\equiv a^{\dagger}_{k}$ and $a^{\dagger}_{2k+1}\equiv b^{\dagger}_{k}$ and assume $k=0,1,\ldots,d-1$, which allows us to write

[TABLE]

In many situations the division into the two groups of fermions, call them $A$ and $B$, is quite natural from the physical point of view. For example, $a^{\dagger}_{k}$ can correspond to creation of a particle with spin up and $b^{\dagger}_{k}$ to creation of a particle with spin down. In the remaining part of this work we assume that we deal with these two types of fermions and we explicitly use the operators $a_{k}^{\dagger}$ and $b_{k}^{\dagger}$.

The real non-negative coefficients $\{\lambda_{k}\}$ determine the correlations between the two fermions. More precisely, one can define the purity $\frac{1}{d}\leq P\leq 1$ of the system as

[TABLE]

If $P=1$ the system is separable and in any other case the system is entangled, the smaller the purity the more entangled it is. For a pair of $d$-level systems the smallest purity, corresponding to a maximally entangled state, is $1/d$.

While it is true that an even number of fermions can behave like a boson, the exact conditions under which these fermions can be treated as a single bosonic particle has been studied for a long time Monique2001 ; Monique2003a ; Monique2003b ; Monique2009 ; Monique2010 ; Monique2011a ; Monique2011b ; Thilagam2013 ; Thilagam2015 ; Monique2015 ; Monique2016 , for a review cf. Monique2008 . An interesting contribution was done recently by Law , who showed that the bi-fermionic state (4) has properties of a single bosonic particle in the limit $P\rightarrow 0$. This result established a connection between the theory of entanglement and the studies on composite particles. The keystone is the idea that the bosonic behaviour of composite system does not stem exactly from interactions, but from entanglement. The above paradigm shift, from interactions to entanglement, is quite subtle. Nevertheless, although entanglement cannot be created without interactions, in principle it does not require interactions to last, once it is created. This allows to study the concept of composite particles that are bound solely by entanglement zak2017 ; zak2018 .

The idea of Law is based on the ladder structure of bosonic operators. If the state (4) is to be treated like a state of a single boson

[TABLE]

then the creation operator $c^{\dagger}$ should obey

[TABLE]

However, due to the fact that $c^{\dagger}$ is not a perfect bosonic operator one gets

[TABLE]

where $\chi_{N}$ is a factor describing a departure from perfect bosonic behaviour. In addition

[TABLE]

and

[TABLE]

where

[TABLE]

and $|\varepsilon_{N}\rangle$ is a state of $N-1$ fermionic A-B pairs. This state is orthogonal to the Fock state $|N-1\rangle$. Its norm is

[TABLE]

Law observed that in the limit $P\rightarrow 0$ one has $\frac{\chi_{N}}{\chi_{N-1}}\rightarrow 1$ for all $N$. In this case $\alpha_{N}\rightarrow 1$ and $\langle\varepsilon_{N}|\varepsilon_{N}\rangle\rightarrow 0$, therefore the composite boson operator $c^{\dagger}$ becomes perfect bosonic. Moreover, Law showed that $\frac{\chi_{N}}{\chi_{N-1}}\leq 1-P$. This result was further strengthen by Chudzicki, Oke and Wootters Wootters , who provided the following bounds

[TABLE]

and by Tichy, Bouvrie and Mølmer who provided the tights bounds for a given $P$ Bouvrie2012a . As a result, in the limit $P\ll 1/N$ the composite creation operators become bosonic.

Interestingly, $P$ is directly related to $\chi_{2}$

[TABLE]

Therefore, because of bounds (13), the properties of two composite bosons tell us about properties of more than two such particles.

II.2 Maximally entangled bi-fermions

In our recent work zak2017 we introduced the concept of a maximally entangled composite boson

[TABLE]

i.e., $\lambda_{k}=1/d$ for all $k$. The above state has a few nice properties. First of all, the internal structure is described by only one integer $d$. In addition, it is easy to evaluate

[TABLE]

and as a result one gets

[TABLE]

Next, let us define the following bi-fermionic states

[TABLE]

where $r,s=0,1,\ldots,d-1$. Each of these states corresponds to a good composite boson for sufficiently large $d$, since for each of these states $P=1/d$. This is the smallest possible $P$ that can be achieved for a bi-partite system with $d$ modes, therefore the creation operators $c^{\dagger}_{s,r}$ provide the best possible bosonic quality. The $d^{2}$ states form an orthonormal basis because

[TABLE]

where $\delta_{x,y}$ is the Kronceker delta. Therefore, any bi-fermionic state can be represented as a linear combination of these states. Interestingly, the above orthogonality relation does not imply orthogonality of states corresponding to more than one composite boson. In general $c^{\dagger N}_{s,r}|0\rangle$ and $c^{\dagger N}_{s^{\prime},r^{\prime}}|0\rangle$ are not orthogonal. The most extreme case corresponds to $N=d$, for which $c^{\dagger d}_{s,r}|0\rangle\equiv c^{\dagger d}_{s^{\prime},r^{\prime}}|0\rangle$ for any $s$, $s^{\prime}$, $r$ and $r^{\prime}$. In simple words, for $N=d$ all available modes are filled with fermions and there is exactly one state describing this possibility. The fact that higher power of different composite boson operators do not preserve orthogonality has been already noticed in Monique2008 .

Finally, notice that although a perfect boson remains a boson even if it is in a superposition, the superposition of composite boson states

[TABLE]

may not correspond to a good composite boson. In the worst case scenario this superposition can correspond to a product state $a^{\dagger}_{k}b^{\dagger}_{k^{\prime}}$. However, if the number of terms in superposition is large, the system should still manifest good bosonic behaviour. Still, good bosonic behaviour is always a question of how many composite bosons occupy the same state.

II.3 Bosonic quality of the Hubbard ground state

In the limit of strong interactions one can use the hard-core boson approximation (see the end of this section) to prove that the ground state of the Hubbard model is a composite boson state of the form (4) with $\lambda_{k}=1/d$. In general, the one dimensional Hubbard model can be solved exactly using Bethe ansatz Lieb1968 . Here, we show that the bipartite ground state can be found using the above maximally entangled states (18).

The Hubbard Hamiltonian for our system is given by

[TABLE]

where the parameters $J,U\geq 0$,

[TABLE]

is the kinetic energy (hopping) term and

[TABLE]

describes the attractive point interaction between fermions A and B. In the above formula we assume the periodic boundary conditions $d\equiv 0$.

Next, consider the action of the Hamiltonian (21) on $|s,r\rangle\equiv c^{\dagger}_{s,r}|0\rangle$

[TABLE]

The Hamiltonian does not change the parameter $r$, therefore for each $r$ one can consider a separate decoupled set of equations. To find the ground state, we need to choose $r$ which minimizes the energy. The hopping amplitude is $-J(1+e^{\pm i\frac{2\pi}{d}r})$ and the greatest negative contribution occurs for $r=0$. Therefore, we fix $r=0$ so the hopping term becomes $-2J$.

In addition, we assume for the moment that the Hamiltonian describes particles hopping on infinite discrete line. Therefore, we get $-\infty\leq s\leq\infty$. This corresponds to $d\rightarrow\infty$. As a result, the states (18) have an infinite number of terms and the corresponding purity is

[TABLE]

therefore they can be considered perfect bosonic for any number of composite bosons $N$.

We represent the candidate ground state as

[TABLE]

for which

[TABLE]

where $\varepsilon$ is the ground state energy. The goal is to find the coefficients $\alpha_{s}$ and the energy $\varepsilon$. The corresponding set of recurrence equations consists of typical equations

[TABLE]

which apply to cases $s\neq 0$, and an atypical equation

[TABLE]

The solution to the above equations (provided in Appendix A) yields

[TABLE]

where $A$ can be determined from normalization. The corresponding energy is

[TABLE]

In general, the ground state (26) is a superposition of many bi-fermionic maximally entangled states (18), therefore it does not need to describe a perfect composite boson. However, in the limit $U\gg J$ the energy becomes $\varepsilon\rightarrow-U$, $\alpha_{s}\rightarrow 0$ for $s\neq 0$ and $\alpha_{0}\rightarrow 1$. Therefore, in this case the ground state of $\mathcal{H}$ is dominated by $c^{\dagger}_{0,0}|0\rangle$, i.e., it can be considered a perfect composite boson.

The system allows us to address the difference between entanglement and interactions. If the bosonic behaviour of two fermions were determined solely by interaction, one could choose $J=0$ and the Hamiltonian would consist only of the interaction part. The corresponding ground state would be degenerated and would be of the form $\eta^{\dagger}_{k}|0\rangle\equiv a^{\dagger}_{k}b^{\dagger}_{k}|0\rangle$, i.e., any pair of fermions A and B occupying the same mode would be considered a ground state. The operator $\eta_{k}^{\dagger}$ obeys the following commutation relations

[TABLE]

These relations are not bosonic. The A-B pairs generated by such operators are sometimes called hard-core bosons. Note, that $\eta_{k}^{\dagger}$ creates two fermions in the product state. This is a clear manifestation of the fact that bosonic behaviour needs fermionic entanglement. This entanglement is provided by the introduction of the kinetic energy term which lifts of the degeneracy.

III Four-partite composite bosons

In this section we consider composite bosons made of four fermions. As shown above, the Hubbard Hamiltonian with point interaction can lead to formation of bipartite composite bosons of the A-B type. Now, we would like to observe formation of a single four-partite bosonic molecule (see Fig. 1). Due to Pauli exclusion the bound pair described by $\eta^{\dagger}_{k}\equiv a^{\dagger}_{k}b^{\dagger}_{k}$ cannot interact via point interaction with any other such pair. Therefore, in order to observe formation of larger compounds we need to introduce a nearest neighbour interaction.

In the following subsections we introduce a four-partite molecular state and compare it with a state of two bi-partite composite bosons. Next, we analyse which of these states dominate the ground state of the extended Hubbard Hamiltonian as the strength of the nearest neighbour interaction changes. Finally, we analyse bosonic properties of the four-partite molecular state.

III.1 Two bipartite composite bosons

The bipartite ground state of (21) in the limit $U\gg J$ can be approximated as $c^{\dagger}_{0,0}|0\rangle$. Let us represent it in terms of operators $\eta^{\dagger}_{k}$

[TABLE]

From now on we use the operators $\eta^{\dagger}_{k}$ whenever it is helpful to simplify our notation.

Next, let us consider Fock state representing two such composite bosons

[TABLE]

where $\chi_{2}=1-P=(d-1)/d$. We are going to show that for sufficiently weak nearest neighbour interactions and for $U\gg J$ the above state dominates the ground state of the extended one-dimensional Hubbard model.

III.2 Four-partite entangled states

Consider the following states of four fermions represented via two $\eta^{\dagger}_{k}$ operators

[TABLE]

where $r=0,1,\ldots,d-1$ and $s=1,\ldots,d/2$. The latter results from $\eta_{k}^{\dagger}\eta_{k+s}^{\dagger}=\eta_{k+s}^{\dagger}\eta_{k}^{\dagger}$ and we assume that $d$ is even.

Just like (18), the above states form an orthonormal basis

[TABLE]

therefore they can be used to represent any state of two A-B pairs created by $\eta_{k}^{\dagger}$. In particular, the state (34) can be represented as

[TABLE]

Note, that the structure of (35) differs from the one of (34). Due to correlations between A-B pairs the number of terms in the first one is quadratically smaller than in the second one. Writing these states explicitly using operators $a^{\dagger}_{k}$ and $b^{\dagger}_{k}$ we see that terms in (35) are of the form $a^{\dagger}_{k}b^{\dagger}_{k}a^{\dagger}_{k+s}b^{\dagger}_{k+s}$, which indicates four-partite correlations since $s$ is a constant. On the other hand, terms in (34) are of the form $a^{\dagger}_{k}b^{\dagger}_{k}a^{\dagger}_{k^{\prime}}b^{\dagger}_{k^{\prime}}$, which indicates bipartite correlations since $k$ and $k^{\prime}$ are independent indices. We already know that bipartite entangled fermionic states exhibit bosonic properties. At the end of this section we are going to examine bosonic properties of four-partite entangled states.

III.3 Extended Hubbard model

In order to study the formation of four-partite composite bosons we extend the Hamiltonian (21) to include the nearest neighbour interaction. The new Hamiltonian is of the form

[TABLE]

where $J,U,\gamma\geq 0$ and the new term is

[TABLE]

It represents an attractive nearest neighbour interaction between fermions A and B. This is a simplified model, since we could also introduce an attractive nearest neighbour interaction between fermions A and A or B and B. However, because in our model A interacts with B via point interaction, the nearest neighbour interaction between the fermions of the same type is going to be mediated by the fermions of the other type. For example, the fermion A in mode $k$ is going to interact with fermion B in mode $k$ that interacts with another fermion A in mode $k+1$. This leads to indirect interactions between A in mode $k$ and A in mode $k+1$. Moreover, we are going to show below that in the particularly interesting limit $U\gg J\gg\gamma$ we can focus on the effective Hamiltonian. The form of this Hamiltonian does not depend on whether we choose interaction between A and B, A and A, or B and B.

In general we need to consider a four-partite problem, however in the limit $U\gg J\gg\gamma$ the above Hamiltonian can be represented in an effective form using operators $\eta^{\dagger}_{k}$ (for details see Appendix B)

[TABLE]

for which our four-partite problem reduces to a bipartite problem.

In order to find the ground state, let us consider the action of (40) on $|s,r\rangle_{q}\equiv q^{\dagger}_{s,r}|0\rangle$

[TABLE]

where $\bar{J}=2\frac{J^{2}}{U}$ and $\bar{\gamma}=2(\gamma-\bar{J})$. This resembles the previous case. As before, the parameter $r$ is not affected by (40), therefore we choose $r=0$ to minimize the kinetic energy. Moreover, we consider $d\rightarrow\infty$, therefore there is no upper bound on $s$.

We assume that the ground state is of the form

[TABLE]

and $\mathcal{H}_{eff}|\psi_{0}\rangle=\bar{\varepsilon}|\psi_{0}\rangle$ hence we obtain the following set of typical recurrence equations

[TABLE]

which apply to cases $s>1$, and an atypical equation

[TABLE]

Using the same methods as before (see Appendix C) we find that

[TABLE]

where $B$ is the normalization constant. The above solution works for $U\gg\gamma>2\bar{J}$, since in this case $\lim_{s\rightarrow\infty}\beta_{s}=0$. This is the prerequisite for the bound state and is the key assumption behind the solution (see Appendix C). For $\gamma\leq 2\bar{J}$ there is no bound state.

Let us analyse the properties of the ground state. When $\gamma\gg 2\bar{J}$ the term $\beta_{1}$ becomes much larger than any other $\beta_{s}$ and the ground state tends to a four-partite state describable by $q_{1,0}^{\dagger}|0\rangle$. On the other hand, for $\gamma\rightarrow 2\bar{J}$ all coefficients become equal, i.e., $\beta_{s}\rightarrow B$. This corresponds to two bipartite composite bosons represented by the state (34) or (37).

It is somehow surprising that the second case occurs for non-zero $\gamma$. In order to understand this effect let us recall the form of (40) and the relations (32). In order to have (34) as a ground state each configuration of two pairs $\eta_{k}^{\dagger}\eta_{k^{\prime}}^{\dagger}$ needs to contribute with the same amount of energy. More precisely, we want the ground state to be of the form

[TABLE]

where $\mathcal{N}$ is the normalization constant. After the application of (40) the state changes to

[TABLE]

where $\varepsilon_{k,k^{\prime}}$ are different energy contributions. However, the above state should be of the same form as before. Only the normalization constant may change. But this requires that $\varepsilon_{k,k^{\prime}}$ is the same for all $k$ and $k^{\prime}$. If the interaction term were zero, the kinetic energy would result in $\varepsilon_{k,k^{\prime}}=-4\bar{J}$ for $k^{\prime}\neq k+1$ and $\varepsilon_{k,k+1}=-2\bar{J}$. This is due to Pauli exclusion principle, which imposes (32). To compensate this, we need to add the nearest neighbour interaction of the strength $-\gamma=-2\bar{J}$.

III.4 Numerical simulations

In order to confirm our predictions, we studied numerically the ground state in the limit $U\gg J\gg\gamma$ and for $d=8$ (assuming periodic boundary conditions). We found that the results of numerical simulations using (38) coincide with the ones obtained using (40). For the fixed $J$ and $U$ we studied how the ground state $|\psi(\gamma)\rangle$ changes as the parameter $\gamma$ increases. We calculated the fidelities $|\langle\psi(\gamma)|q_{1,0}^{\dagger}|0\rangle|^{2}$ and $|\langle\psi(\gamma)|\frac{c^{\dagger 2}_{0,0}}{\sqrt{2\chi_{2}}}|0\rangle|^{2}$. The results are plotted in Fig. 2. The parameter $U$ was chosen to be of the order $10^{5}$, $J$ of the order $10^{2}$ and $\gamma$ of the order of $10^{0}$. Interestingly, if one fidelity approaches one, the other fidelity approaches $1/4$. This is because of the finiteness of space in our simulations. Note that $|\langle 0|q_{1,0}\frac{c^{\dagger 2}_{0,0}}{\sqrt{2\chi_{2}}}|0\rangle|^{2}=d\binom{d}{2}^{-1}$ (see Eqs. (37) and (35)). In the limit $d\rightarrow\infty$ the two states become orthogonal.

III.5 Bosonic properties of the four-partite entangled state

Finally, let us investigate if the state $q^{\dagger}_{1,0}|0\rangle$ can be interpreted as a single four-partite bosonic molecule. We do this by analysing the bosonic quality of the operator $q^{\dagger}_{1,0}$ via the same methods as in Law . In particular, we analyse the parameter $\chi^{(2)}_{N}$ defined by

[TABLE]

The operator $q^{\dagger}_{1,0}$ generates the bosonic-like ladder structure analogous to (II.1). If the ratio $\frac{\chi^{(2)}_{N}}{\chi^{(2)}_{N-1}}\rightarrow 1$ for all $N$, then $q^{\dagger}$ can be considered perfect bosonic. From the results presented in Appendix D we get that

[TABLE]

and

[TABLE]

The above is upper bounded by 1 and lower bounded by

[TABLE]

which in the limit $d\gg N$ approaches 1. Therefore, for $d\gg N$ the action of operator $q^{\dagger}_{1,0}$ can be considered as a creation of a single four-partite bosonic particle.

IV Multipartite composite bosons

Here we consider composite boson made of $2N$ fermions. Such composite particles are expected to appear in the effective extended Hubbard model (40) if the particle attraction $\gamma$ is stronger than the effective kinetic energy contribution $\bar{J}$. On the other hand, if the effective kinetic energy contribution is two times stronger than the attraction between the particles, the system should be describable by $N$ A-B independent pairs. The transition between these two types of behaviour should occur when effective kinetic energy and attraction are of the same strength. A schematic representation of the above multifermionic assemblies are shown in Fig. 3.

IV.1 Multipartite entangled states

It is natural to expect that for N fermions of type A, N fermions of type B, and sufficiently strong nearest neighbour interaction the Hamiltonian (40) has the following ground state

[TABLE]

Note, that for $M=1$ the above state becomes the bi-partite state (18) corresponding to $c^{\dagger}_{0,0}|0\rangle$ and for $M=2$ it becomes the four-partite state (35) corresponding to $q^{\dagger}_{1,0}|0\rangle$. Therefore, the above expectation is in accordance with the few particle cases studied above.

IV.2 Bosonic quality of multipartite entangled states

The state $q^{\dagger}_{(M)}|0\rangle$ can be considered a single composite boson made of $2M$ fermions if $d\gg M$. We provide a detailed proof in the Appendix D. As a result, the state

[TABLE]

represents $N$ such composite particles, provided $d\gg NM$, where

[TABLE]

The above bosonic behaviour stems from the fact that state $q_{(M)}^{\dagger}|0\rangle$ is multipartite entangled. In fact, multipartite entanglement is a necessary condition for a bosonic behaviour of 2M fermions. Firstly, note that 2M fermions in a fully separable state, i.e., having Slater rank one Eckert2002 ,

[TABLE]

cannot be considered a composite boson. The simplest proof is that for such states there is no two-boson state. We get $c_{sep}^{\dagger 2}|0\rangle=0$, which is a consequence of Pauli exclusion. Next, note that the same happens if at least one fermion is in a well defined mode, for example

[TABLE]

Therefore, each fermion needs to be entangled to prevent $q_{(M)}^{\dagger 2}=0$. This implies multipartite entanglement. However, multipartite entanglement is not enough. In Appendix E we show that a multifermionic state representing a single composite boson needs to be genuinely multipartite entangled. This is the only way to recover the ladder structure (II.1) of composite bosonic operators.

IV.3 Composite bosons of various sizes

Let us define the following states of 2N fermions representing $k$ composite bosons of various sizes

[TABLE]

where we assume the decreasing order, i.e., $M_{1}\geq\ldots\geq M_{k}$. Moreover, $M_{1}+\ldots+M_{k}=N$ and $\mathcal{N}$ is a normalization factor. This normalization factor is necessary due to two reasons. Firstly, it may happen that in the above state there are $m$ composite bosons of the same size, i.e., $M_{i}=\ldots=M_{i+m}=M$. In this case, if there were no normalization factor the norm of the state would be proportional to $\chi_{m}^{(M)}m!$. Secondly, even if all composite bosons were of a different size, the state $q^{\dagger}_{(M_{1})}\ldots q^{\dagger}_{(M_{k})}|0\rangle$ would not be normalized, despite the fact that each of the states $q^{\dagger}_{(M_{1})}|0\rangle,\ldots,q^{\dagger}_{(M_{k})}|0\rangle$ is of norm one. This is because of indistinguishability of fermions, Pauli exclusion and finiteness of $d$. For example, consider a state

[TABLE]

It is easy to find (using the already applied methods) that $\mathcal{N}^{2}=\frac{d^{2}}{d^{2}-4d}$, which in the limit $d\rightarrow\infty$ tends to one.

Our next goal is to investigate to which state $|M_{1}+\ldots+M_{k}\rangle$ the ground state of (40) corresponds to. For weak interaction the ground state should correspond to $|1+\ldots+1\rangle$, whereas for strong interaction it should correspond to $|N\rangle$. In addition, we are interested in how the one state changes into the other as the interaction strength $\gamma$ increases.

IV.4 Numerical simulations for N=3 and N=4

Here, we discuss numerical results for $N=3$ and $N=4$, i.e., composite bosons made of six and eight fermions. As before, we assume the strong point interaction limit, therefore we consider three and four A-B pairs, respectively, to which we apply the effective Hamiltonian (40).

We numerically found the ground state $|\psi(\gamma)\rangle$ and evaluated the fidelities $|\langle\psi(\gamma)|M_{1}+\ldots+M_{k}\rangle|^{2}$. Due to high computational complexity we considered $d=10$, which does not imply a perfect bosonic quality, but still allows us to see some important qualitative behaviour of the model. The corresponding fidelities are plotted in Fig. 4.

In case of $N=3$ we observe that for low value of the ratio $\gamma U/J^{2}$ the ground state is dominated by the state $|1+1+1\rangle$. The fidelity $|\langle\psi(\gamma)|1+1+1\rangle|^{2}$ reaches one for $\gamma U/J^{2}=4$. This maximum occurs for $\gamma\neq 0$ due to the same reason as in the case $N=2$. For $\gamma U/J^{2}>4$ there is a small region in which the ground state is dominated by $|2+1\rangle$, although the fidelity of this state never reaches one. After this region the ground state is dominated by $|3\rangle$ and as $\gamma U/J^{2}\gg 4$ the fidelity $|\langle\psi(\gamma)|3\rangle|^{2}\rightarrow 1$.

In case of $N=4$ the transition from $|1+1+1+1\rangle$ to $|4\rangle$ seems to be more complex. As before, the ground state is dominated by $|1+1+1+1\rangle$ for low values of $\gamma U/J^{2}$ and the corresponding fidelity reaches maximum of one for $\gamma U/J^{2}=4$. As the ratio increases, the state $|2+1+1\rangle$ and then $|3+1\rangle$ starts to take over, but soon the state $|4\rangle$ becomes to dominate and its fidelity approaches one as the ratio becomes large. Interestingly, the state $|2+2\rangle$ seems to play no significant role in the above transition. The fact that it appears in the Fig. 4 is rather due to its overlap with other states, which is caused by the relatively small value of $d$. Note, that for high values of $\gamma U/J^{2}$ the ground state is describable by $|4\rangle$, but the fidelities corresponding to other states are still high. In particular, the fidelity $|\langle\psi(\gamma)|2+2\rangle|^{2}$ does not change much after it reaches its maximum of $\approx 0.3$.

The above observations allow us to speculate that the transition from N A-B pairs to a single N-partite bosonic molecule follows the pattern,

[TABLE]

This means that bi-partite bosonic particles assemble into a large bosonic molecule by adding particle by particle into a single large compound. In this situation one does not observe creation of two, or more, larger compounds. In the next subsection we support this hypothesis by showing that the transition (60) is most energetically favourable. It is worth to add here that during transition process (when strength of $\gamma$ is increased) we assume that system follows its ground state.

IV.5 Transition from N bipartite composite bosons to a single bosonic molecule

Before we go into details, let us first discuss an important property of states $|M_{1}+\ldots+M_{k}\rangle$. They are superpositions of terms

[TABLE]

The sequence $\eta^{\dagger}_{j_{i}}\ldots\eta^{\dagger}_{j_{i}+M_{i}-1}$ represents the i’th compound made of $M_{i}$ A-B pairs. There is a possibility that in this superposition some two compounds are next to each other, i.e., for $\eta^{\dagger}_{j_{i}}\ldots\eta^{\dagger}_{j_{i}+M_{i}-1}$ and $\eta^{\dagger}_{j_{i+1}}\ldots\eta^{\dagger}_{j_{i+1}+M_{i+1}-1}$ we have $j_{i+1}=j_{i}+M_{i}$. Nevertheless, the total number of terms in the superposition is of the order $d^{k}$, whereas the number of terms with adjacent compounds scales as $d^{k-1}$. Therefore, for $d\gg N$ one can assume that almost all terms in the superposition correspond to non-adjacent compounds (non-adjacency assumption). This leads to

[TABLE]

for two different configurations $\{M_{1},\ldots,M_{k}\}$ and $\{M^{\prime}_{1},\ldots,M^{\prime}_{k^{\prime}}\}$, although the total number of particles in both configurations is the same.

Let us once more consider the Hamiltonian (40) and estimate its expectation value for a state $|M_{1}+\ldots+M_{k}\rangle$

[TABLE]

where we explicitly split the kinetic and the potential energy parts. First, we consider the kinetic part. Note that under the action of $\mathcal{H}_{k}$ the compounds made of more than a single A-B pair split into smaller compounds, i.e.,

[TABLE]

where $\bar{J}=2J^{2}/U$. Therefore, under the action of $\mathcal{H}_{k}$ the state $|M_{1}+\ldots+M_{k}\rangle$ changes into superposition of states $|M^{\prime}_{1}+\ldots+M^{\prime}_{k^{\prime}}\rangle$. Interestingly, if in the original state there were no single A-B pairs ($M_{k}>1$), then none of the states $|M^{\prime}_{1}+\ldots+M^{\prime}_{k^{\prime}}\rangle$ in the effective superposition is equal to the original one. Therefore, the non-adjacency assumption implies

[TABLE]

Next, assume that in $|M_{1}+\ldots+M_{k}\rangle$ there are $r$ single A-B pairs. The kinetic energy term moves these pairs one step to the right and one step to the left. Recall that the above state consists of a superposition of all possible (non-adjacent) configurations of such pairs. Therefore,

[TABLE]

This is because each term in $|M_{1}+\ldots+M_{k}\rangle$ can be obtained from $2r$ other therms by shifting some A-B pair either one step to the right or one step to the left.

The expectation value of the potential energy part is much easier to evaluate. It is straightforward to show that under the non-adjacency assumption

[TABLE]

where $\bar{\gamma}=2(\gamma-\bar{J})$. Interestingly, this value depends only on the total number of compounds $k$, not on the way the A-B pairs are distributed between these compounds $\{M_{1},\ldots,M_{k}\}$. As a result

[TABLE]

Now, consider two states with the same number of compounds $k$, but different number of single A-B pairs, $r$ and $r^{\prime}<r$ (e.g. $|3+1+1\rangle$ and $|2+2+1\rangle$). The corresponding average energies are $-2r\bar{J}-(N-k)\bar{\gamma}$ and $-2r^{\prime}\bar{J}-(N-k)\bar{\gamma}$. It is clear that the average energy is lower for the state containing more single A-B pairs. Therefore, only the states of the form $|M+1+\ldots+1\rangle$ need to be taken into account during the transition between $|1+\ldots+1\rangle$ into $|N\rangle$, which confirms our previous hypothesis. For such states the average energy equals

[TABLE]

where $\delta_{M,1}$ is the Kronecker delta.

In Fig. 5 we plot few of these values for $N=10$. The lowest energy corresponds either to $|1+1+\ldots+1\rangle$ or to $|10\rangle$. This observation implies the following conjecture: for $d\gg N$ (no-adjacency assumption) the state $|1+1+\ldots+1\rangle$ is directly transformed into $|N\rangle$ as $\gamma$ increases. Referring to Fig. 3, the above means that the assembly (b) goes directly to (d) and the assembly (c) never occurs. Using the formula (69) one can show that the transition should occur for $\bar{\gamma}/\bar{J}=2N/(N-1)$, which corresponds to $\gamma U/J^{2}=2+2N/(N-1)$. This agrees with plots in Fig. 4. From the point of view of entanglement analysis, the conjecture implies that the system of 2N bipartite entangled fermions is directly transformed into genuinely multipartite entangled state.

IV.6 Entanglement and correlation function between fermion pairs

In general, in the ground state $|\psi(\gamma)\rangle$ the $N$ A-B pairs are correlated. For $\gamma=0$ an effective repulsive interaction between them, inherent in one-dimensional fermion systems, leads to particle anti-bunching (see red dot in Fig. 6). This repulsive interaction is compensated by a nearest neighbour interaction of strength $\gamma=4J^{2}/U$ (orange squares) – in this case the pairs independently occupy the sites of the lattice. For large interaction strength $\gamma\gg 4J^{2}/U$ (blue diamonds) the pairs become maximally multipartite entangled.

Entanglement has been demonstrated to be useful for the identification of quantum phase transitions in the extended 1D Hubbard model ShiJian2004 . Here we will show that the entanglement between fermion pairs can be used to identify the transition between bunching and anti-bunching.

Since fermion pairs are identical, the amount of entanglement can be characterized by the purity

[TABLE]

of the single-pair reduced density matrix

[TABLE]

Note that, as in case of indistinguishable fermions BouvrieValdesetall2016 , the normalization factor of the reduced density matrix $\rho_{ij}^{(1)}$ is $1/N$. We show in Fig. 7 the entanglement $1-P_{1}(\gamma)$ between one pair and the remaining $N-1$ pairs. The correlations are minimized for the interaction $\gamma U/J^{2}=4$. In this interaction regime the purity has the analytical solution

[TABLE]

In the limit $d\to\infty$ the ground state corresponds to $N$ independent pairs ($P_{1}(4J^{2}/U)=1$), where each pair is maximally entangled. On the other hand, the maximal correlations are reached for $\gamma\gg 4J^{2}/U$, for which

[TABLE]

The entanglement takes its maximum value in the limit $d\to\infty$, where $P_{1}(\infty)=0$. In this case all particles are maximally multipartite entangled.

The above has observable consequences in the second order correlation function

[TABLE]

of the fermions pairs located at positions $i$ and $j$. The mean occupations of sites are constant $\langle\psi(\gamma)|\eta_{i}^{\dagger}\eta_{i}|\psi(\gamma)\rangle=N/d$. As can be seen in Fig. 8, small nearest neighbour interaction yields anti-bunching between nearest neighbours $g^{(2)}_{\gamma}(i,i+1)<1$. For $\gamma=4J^{2}/U$ the anti-bunching is only due to the Pauli exclusion principle between identical pairs $g^{(2)}_{4J^{2}/U}(i,j)=\frac{d(N-1)}{N(d-1)}~{}~{}~{}\forall i,j$ and in the limit $d\to\infty$ the usual boson correlation function $g_{\text{bosons}}^{(2)}=(N-1)/N$ is observed. Finally, in the limit of large interaction $g^{(2)}_{\infty}(i,j)=d(N-|i-j|)/N$ if $|i-j|\leq N$ and $g^{(2)}_{\infty}(i,j)=0$ otherwise. Therefore, bunching dominates the particle statistics, even for long range particle correlations $g^{(2)}_{\infty}(i,j)>1$ for $(|i-j|)<N-N/d$. Comparing figures 7 and 8, the minimal correlations clearly indicate the transition between bunching and anti-bunching.

V Summary

We studied multifermionic states and showed that in order to form a composite boson made of 2N fermions the system needs to exhibit a genuine 2N-multipartite fermionic entanglement. Next, we tested our results in a specific physical model. We examined emergence of multipartite composite bosons in the ground state of the extended one-dimensional Hubbard model. In particular, we focussed on which different bosonic assemblies emerge as the intra-particle interaction increases. We solved the problem analytically for two and four fermions. In the second case we assumed the strong point interaction limit $U\gg J$. For higher number of fermions we also considered the strong point interaction limit and performed numerical simulations. They confirmed that if the nearest neighbour interactions are strong the system forms a single multipartite bosonic molecule, whereas if they are weak the system consists of many bipartite bosonic particles. Then, we made an analytical estimation of the average energy for different fermionic assemblies. This lead us to the following conjecture: as the nearest neighbour interaction increases it is energetically favourable that the fermions transform from the state representing many bipartite composite bosons into the state corresponding to a single multipartite composite boson, without going through other possible assemblies. This conjecture should be directly verifiable using the analytical methods of Lieb and Wu based on Bethe ansatz Lieb1968 .

Acknowledgements. We would like to thank Ana Majtey, Cecilia Cormick and Paula Cespedes for sharing relevant information concerning the effective Hamiltonian (Eq. (40)), Marcin Karczewski and Andrzej Grudka for stimulating discussions and suggestions. Z.K., A.S. and P.K. were supported by the National Science Center in Poland through NCN Grant No. 2014/14/E/ST2/00585. P.A.B. gratefully acknowledge support by the Conselho Nacional de Desenvolvimento Científico e Tecnológico do Brasil and by the Spanish MINECO project FIS2014-59311-P (co-financed by FEDER).

Appendix A

Here we provide solution to the set of recurrence equations (28) and (29). The typical equations (28) have general solution of the form

[TABLE]

where $r_{0}$ and $r_{0}^{-1}$ are roots of the following quadratic equation

[TABLE]

and $A$ and $B$ are constants. In addition

[TABLE]

Note, that the ground state should be a bound state. Such states are normalizable, therefore we have the commonly used constraint $\lim_{s\rightarrow\pm\infty}\alpha_{s}=0$. Since $r_{0}^{-1}$ is the inverse of $r_{0}$, the modulus of one of them is less than one, whereas the modulus of the other one is greater than one. Say, $|r_{0}|\leq 1$. Therefore, due to the normalization constraint we get

[TABLE]

To determine $\alpha_{0}$ we plug the above and (77) to

[TABLE]

and obtain

[TABLE]

which imply $A=B=\alpha_{0}$.

Finally, we plug the above results to the atypical equation (29) and obtain

[TABLE]

where $\bar{U}=U/2J$, which leads to the following quadratic equation for $r_{0}$

[TABLE]

We get

[TABLE]

Since we assumed $|r_{0}|\leq 1$ we choose the solution with minus sign and get the final form of the ground state

[TABLE]

and the corresponding energy

[TABLE]

Appendix B

Here we derive the effective Hamiltonian (40). We start with the original one (38), which is of the form

[TABLE]

First, lets consider the Hamiltonian $\mathcal{H}^{(1)}=J\mathcal{H}_{0}+U\mathcal{H}_{p}$. The projector on the ground state of $\mathcal{H}_{p}$ reads

[TABLE]

where $\left|j_{n}j_{n}\right\rangle=a_{n}^{\dagger}b_{n}^{\dagger}|0\rangle$, and the unitary evolution operator associated to $\mathcal{H}_{p}$

[TABLE]

where $P_{e}$ is the projector on the first exited state which consists of $N-1$ bounded pairs with energy $(N-1)U$ and one separated pair with zero energy. In the Master Thesis of P. Cespedes cespedes it has been shown that at first order of perturbation

[TABLE]

since $P_{g}H_{0}P_{g}=0$, and that at second order

[TABLE]

where

[TABLE]

The effective Hamiltonian $\mathcal{H}^{(1)}$ for $U\gg J$ is given by

[TABLE]

where we have used that $P_{g}H_{p}P_{g}=-NU$.

Now let us consider the full Hamiltonian $\mathcal{H}$. At first order we have

[TABLE]

Contributions at second order of $\mathcal{H}_{nn}$ can be neglected because they depend on $\gamma^{2}/U\ll J^{2}/U$. Therefore the effective Hamiltonian for $U\gg J\gg\gamma$ is

[TABLE]

Since for a fixed number of particles the first term is constant, we drop it in (40).

Appendix C

Recurrence equations (43) and (44) can be solved in basically the same way as (28) and (29). Firstly, note that (43) is of the same form as (28), therefore

[TABLE]

where $A$ and $B$ are constants and

[TABLE]

We assume $r_{0}\leq 1$ and since $\lim_{s\rightarrow\infty}\beta_{s}=0$ we get $A=0$. Therefore,

[TABLE]

We plug the above to (44) and obtain

[TABLE]

which leads to

[TABLE]

Appendix D

Consider the following state of $M$ fermionic A-B pairs represented by operators $\eta_{k}^{\dagger}=a_{k}^{\dagger}b_{k}^{\dagger}$

[TABLE]

where

[TABLE]

We are looking for the parameter

[TABLE]

Firstly, we write

[TABLE]

The above operators obey

[TABLE]

We get

[TABLE]

where $*$ over the sum denotes that we take into account the relation (108). This leads to

[TABLE]

The binomial coefficient represents the combination with repetitions. The reason we get this value is the following. We consider $N$ compounds, each composed of $M$ A-B pairs. A single compound occupies $M$ positions. In total there are $d$ positions, therefore after putting $N$ compounds there are $d-NM$ unoccupied spaces left. Each unoccupied space can be placed in one of $N+1$ possible places: before first compound, between the first and the second compounds, $\ldots$, after the $N$-th compound. Hence, we choose one of $d-NM$ positions from $N+1$ possibilities and the choices can repeat. In general, a choice of $x$ positions out of $y$, including repetitions, is given by ${x+y-1}\choose{x}$.

The formula (109) can be written as

[TABLE]

Finally, we estimate the ratio $\frac{\chi_{N+1}^{(M)}}{\chi_{N}^{(M)}}$

[TABLE]

The above value is upper-bounded by 1 and lower-bounded by

[TABLE]

In the limit $d\gg NM$ this value approaches 1. In Fig. 9 we plot the above lower bound for $d=10000$.

Appendix E

Let us show that a multifermionic state representing a single composite boson needs to be genuinely multipartite entangled. We analyse what happens if the state is not genuinely multipartite entangled. This means that it is separable with respect to some partition. In case of fermions separability implies that we can write it as

[TABLE]

where $w_{i_{1},i_{2},\ldots}$ and $v_{j_{1},j_{2},\ldots}$ are normalized antisymmetric coefficients corresponding to two partitions, such that

[TABLE]

For the moment we assume that these partitions can have an arbitrary number of fermions (here, without the loss of generality, we drop the division of fermions into A and B).

The main point of our argumentation is the following: if $c^{\dagger}|0\rangle$ represents a single composite boson, then $c^{\dagger 2}|0\rangle$ represents two composite bosons and its norm is approximately $2$. This is a necessary condition for bosonic creation and annihilation operators to obey the ladder structure (II.1). We write explicitly

[TABLE]

For non-overlapping terms (the ones containing different fermionic creation operators) we have

[TABLE]

(same holds for the terms in the other partition). In the above $x=+1$ ($x=-1$) if the number of fermions in each partition is even (odd). If $x=-1$ all the terms cancel and $c^{\dagger 2}=0$. Therefore, we allow only for partitions with even number of fermions. We can write

[TABLE]

The factor of 4 stems from (115), i.e., the terms corresponding to products within each partition are counted twice. If there were $s$ partitions, the overall factor would be $2^{s}$. In addition, the symbol $\ast$ above the sums indicates that we do not sum over the terms with overlapping indices. For example, due to Pauli exclusion we do not sum over the terms corresponding to $w_{1,4,7,\ldots}w_{2,4,9,\ldots}$ etc.

We obtain

[TABLE]

The above can be rewritten as

[TABLE]

where

[TABLE]

and the sum is taken only over terms with overlapping indices. If the number of modes that can be occupied by fermions is arbitrarily large and the states of each partition are highly entangled (e.g., they are the ones discussed in Appendix D), then the parameter $\omega(\ast)$ can be arbitrarily small. Therefore, the norm of $c^{\dagger 2}|0\rangle$ is approximately four, not two. If instead of two partitions we considered $s$ partitions, the norm would be approximately $2^{s}$. For a proper bosonic behaviour we need $s=1$, which implies genuinely multipartite entangled state.

Finally, we need to make one remark. The norm of $c^{\dagger 2}|0\rangle$ can be in principle equal to two, even for $s=2$. This requires $\omega(\ast)=1/2$, i.e., a large number of terms with overlapping indices (occurring for not that strongly entangled states). However, the operator $c^{\dagger}$ also needs to recover the ladder structure for higher powers than two. If $\omega(\ast)=1/2$, then $\langle 0|c^{N}c^{\dagger N}|0\rangle\neq N!$ for $N>2$. In fact, for a large number of overlapping terms we get $\langle 0|c^{N}c^{\dagger N}|0\rangle\rightarrow 0$ for relatively small N. This is the consequence of Pauli exclusion.