# Strongly repulsive anyons in one dimension

###### Abstract

To explore the static properties of the one-dimensional anyon-Hubbard model for a mean density of one particle per site, we apply perturbation theory with respect to the ratio between kinetic energy and interaction energy in the Mott insulating phase. The strong-coupling results for the ground-state energy, the single-particle excitation energies, and the momentum distribution functions up to 6th order in hopping are benchmarked against the numerically exact (infinite) density-matrix renormalization group technique. Since these analytic expressions are valid for any fractional phase of anyons, they will be of great value for a sufficiently reliable analysis of future experiments, avoiding extensive and costly numerical simulations.

## I Introduction

Particles are usually classified as either bosons or fermions, depending on whether the wave function is symmetric or antisymmetric with respect to the exchange of two identical particles.
Some systems, however, may realize quasiparticles with fractional
statistics, called anyons, that acquire a complex phase factor
with under exchange Leinaas and Myrheim (1977); Wilczek (1982).
Most notably, anyons have been used in the description of the fractional quantum Hall effect Tsui *et al.* (1982); Laughlin (1983).
While anyons are usually restricted to two-dimensional systems, fractional statistics can in principle be defined in arbitrary dimensions Haldane (1991).

One dimensional (1D) anyon models can be expressed in terms of bosonic operators by using a generalized Jordan–Wigner transformation.
There are several proposals to utilize this equivalence to implement an anyon-Hubbard model (AHM) by loading ultracold atoms in optical lattices. The fractional exchange statistics is thereby translated into an occupation-dependent hopping phase that experimentally may be implemented by assisted Raman tunneling Keilmann *et al.* (2011); Greschner and Santos (2015) or lattice-shaking-assisted tunneling against potential offsets Sträter *et al.* (2016). One of the advantages of any optical lattice setup is the high degree of control of system parameters including the statistical angle . As yet, however, an experimental realization of anyons in optical lattices has not been archived.

Since the introduction of the AHM Keilmann *et al.* (2011), several theoretical and numerical studies have been carried out, inter alia, exploring the effect of fractional statistics on momentum distributions Tang *et al.* (2015)
and the position of the quantum phase transition between the Mott insulator (MI) and superfluid (SF) Keilmann *et al.* (2011); Arcila-Forero *et al.* (2016) as well as revealing additional phases such as an exotic two-component partially-paired phase Greschner and Santos (2015).
It has been shown that the various superfluid phases of the AHM can be qualitatively understood within a generalized Gutzwiller mean-field ansatz Tang *et al.* (2015); Zhang *et al.* .
Here, we instead focus on the MI phase, using strong-coupling perturbation theory as it has been applied to the Bose-Hubbard model (BHM) Ejima *et al.* (2012a); Damski and Zakrzewski (2006); Elstner and Monien .
In addition to the perturbative analysis, we study the model numerically with the density-matrix renormalization group (DMRG) White (1992); McCulloch ; Schollwöck (2011) and a variational matrix-product state (MPS) ansatz for dispersion relations Haegeman *et al.* (2012, 2013).

The outline of this paper is as follows: In Sec. II we introduce the 1D AHM and apply the anyon-boson mapping by the fractional version of the Jordan–Wigner transformation to the AHM in order to rewrite the Hamiltonian with the bosonic operators. In Sec. III we describe the strong-coupling analysis for the ground-state energy, the momentum-dependent single-hole and single-particle excitation energies, and the momentum distribution functions. To evaluate the validity of the proposed strong-coupling approach we perform an extensive comparison with unbiased data obtained by the MPS-based (infinite) DMRG (iDMRG) technique. Finally, Sec. IV summarizes our results and gives a brief outlook.

## Ii Anyonic Hubbard model

On a linear chain of sites with periodic boundary conditions (PBCs), the Hamiltonian of the 1D AHM is defined as , with

(1) |

and

(2) |

describing the nearest-neighbor anyon transfer ()
and the on-site anyon repulsion (), respectively.
Here, , ,
and are the anyon creation, annihilation,
and particle number operators on site , respectively, which fulfill the generalized commutation relations Keilmann *et al.* (2011)

(3) | |||||

(4) |

Since , regular bosonic commutation relations apply for particles on the same site. Anyons with the fractional angle represent the so-called “pseudofermions,” namely, they behave as ordinary fermions off-site, while being bosons on-site.

Carrying out a fractional Jordan–Wigner transformation Keilmann *et al.* (2011),

(5) |

of Eq. (1) can be rewritten with boson creation () and annihilation () operators as

(6) |

To be more precise, when an anyon hops to the left from site to site , an occupation dependent phase is picked up in the bosonic operator description. Note that , so that the on-site repulsion term is form-invariant under the anyon-boson mapping (5). In order to study the model deep in the Mott-insulating regime, we apply an strong-coupling expansion to . Throughout this work, we restrict ourselves to unit filling.

## Iii Strong-coupling expansions

### iii.1 Ground state

At integer filling , the AHM has a unique ground state,

(7) |

in the limit . The state can be used as a starting point for a perturbative calculation of the ground state in the MI phase.

Executing the unitary Harris-Lange transformation Harris and Lange (1967),
the strong-coupling Hamiltonian of the AHM is derived
in a similar way as for the (BHM) Ejima *et al.* (2012a):

(8) | |||||

(9) |

Practically, we keep a finite order in the expansion of . Retaining for denotes the “th-order expansion.” The operators are defined by requiring that in the th order for , the transformed Hamiltonian conserves the number of double occupancies to th order, that is, for . Higher-order terms in the expansion of are neglected, so is an eigenstate of the strong-coupling Hamiltonian.

Following this recipe, the leading-order terms for and are obtained as

(10) | |||||

(12) | |||||

(13) |

where is the projection operator onto the subspace of eigenstate with interactions, . In the above sums it is implicitly suggested that all indices are different from each other. Higher orders are generated recursively as described in Ref. van Dongen (1994), where the necessary bookkeeping can be done by a computer algebra program. The resulting expansion differs from the one for the BHM only by the hopping operator .

Within the strong-coupling expansion the ground state and ground-state energy of the original Hamiltonian are

(14) |

where is the ground state of , see Eq. (7). Since the Harris–Lange transformation is unitary the operators and ground-state expectation values are translated with , , and .

Calculating the various observables in the strong-coupling expansion then amounts to evaluating chains of hopping operators in the unperturbed ground state which are weighted depending on how they change the number of double occupancies at each step. In doing so, one only has to sum over connected hopping processes that can be evaluated using finite clusters. The difference between the AHM and the BHM enters the strong-coupling expansion through the phase factors picked up by the hopping processes or an explicit -dependence of the observables.

The ground-state energy is simply given by

(15) |

Up to 6th order in , we obtain for the rescaled ground-state energy per site

(16) | |||||

in agreement with Refs. Damski and Zakrzewski (2006); Ejima *et al.* (2012a) in the BHM limit
.

Figure 1 compares the strong-coupling perturbation theory
with iDMRG results for various .
Similar to the case in the BHM Ejima *et al.* (2012a), for small
[e.g., Fig. 1(a) for ], the strong-coupling
series expansion is
in reasonable accordance
with the numerically exact result,
as indicated clearly by the relative errors in Fig. 1.
For 6th order in , the deviation starts in the intermediate-coupling regime
at .
As expected, the quality of the perturbation analysis improves
as higher-order corrections are taken into account. This is valid
for all as demonstrated in Fig. 1.

Figure 1 also shows that the range of validity of the strong-coupling theory becomes worse with increasing . The deviation starts already at in the case of , see panel (c).

### iii.2 Excitation energies

Similar to the ground state , Eq. (7), the energy levels of a single-hole excitation, , and of a single-particle excitation, , can be extracted from the strong-coupling expansion to high order in , since the perturbation analysis for these energy levels also starts from nondegenerate states, i.e., in the case of ,

(17) | |||||

(18) |

Therefore, the single-hole and single-particle excitation energies can be obtained from

(19) | |||||

(20) |

Carrying out the above perturbation analysis up to and including 6th order in , we obtain

(21) | |||||

and

(22) | |||||

In the BHM limit (), we obtain
Eqs. (39) and (40),
which agree with Eqs. (24) and (25) of Ref. Ejima *et al.* (2012a),
respectively, when correcting some misprints; see App. A.

To calculate the dispersion relations of the particle and hole
excitations numerically, we use the variational MPS ansatz introduced
in Refs. Haegeman *et al.* (2012, 2013) that works directly
in the thermodynamic limit. In the following, we give a rough
description of the method. Starting point is an infinite MPS (iMPS)
approximation of the ground state

(23) |

where , the indices label the states of the local Hilbert spaces, are site independent complex matrices, and and are -dimensional vectors. The boundary vectors and will not affect the bulk properties and can therefore be ignored. It is assumed that the transfer matrix has one eigenvalue and that its other eigenvalues are smaller in magnitude. To calculate , we use the iDMRG. The ansatz for the elementary excitations is a momentum superposition of local perturbations which are introduced by replacing the matrices at a single site with matrices :

(24) |

This includes all excitations that are induced by one-site operators but can also describe, to some degree, those corresponding to operators with larger support. Increasing the bond dimension of results, in addition to a better approximation for the ground state energy, in a more general ansatz for the excitations. One can define matrices and such that

(25) | ||||

(26) |

where is the (infinite) ground-state energy
and the matrices
have been combined and reshaped into a vector.
The approximate excitation energies for any momentum
can then be obtained by solving the generalized eigenvalue problem
for the effective Hamiltonian and
the normalization matrix .
As described in detail in Ref. Haegeman *et al.* (2013), must be appropriately
parameterized to exclude zero modes which would result
in and to impose orthogonality
to the ground state. A linear parametrization fulfilling these requirements
can be chosen such that the normalization matrix becomes
the identity and only a regular eigenvalue problem needs to be solved.
Since the number of particles is a good quantum number,
we can separately target particle and hole excitations
to obtain both and .

In Fig. 2 we compare the strong-coupling results
up to 6th order in with the lowest excitation energy obtained
by the above mentioned iMPS technique. First of all, are
clearly symmetric about in the BHM limit ,
although they become asymmetric for reflecting the
influence of the fractional angle . By considering the
strong-coupling expansions up to 1st order only, this asymmetry
of the excitation energies can be understood well: the minimum
of excitation energy, (),
shifts from to , consistent with
the positive sign of in cosine terms up to
1st order. Quantitatively, the 6th-order expansions agree perfectly
with iMPS data up to . The deviation between both
results starts about especially in .
The single-particle excitation from the perturbation theory
is clearly higher than the lowest excitation energy by iMPS,
e.g. for and with .
Most probably, the lowest excitation by iMPS stems from a
many-particle excitation such as two particles and one hole that are forced into an artificial bound state by the iMPS ansatz Haegeman *et al.* (2013).
Moreover, plotting the higher excitation energies, it is obvious that a continuum of excitations starts to
arise in this regime.
Figure 3 demonstrates a typical example
for with the parameter sets of Fig. 2.
Because of the finite bond dimension , the continuous part of the spectrum is approximated by a finite number of discrete energy levels.
With increasing further, the results of
strong-coupling expansions start to oscillate, see Fig. 2
for .

From the single-hole and single-particle dispersions we can obtain the phase boundaries between MI and SF in the grand-canonical ensemble. The chemical potentials at which the phase transitions for fixed take place are determined by the minimum energies for adding a particle or hole to the MI ground state: and . In general, the minima of the strong-coupling expressions (21) and (22) have to be found numerically. However, in the BHM limit () we have and and thus the gap is given by . In this way, we can reproduce the single-particle gap in the BHM

(27) |

in agreement with Ref. Elstner and Monien .

As in the case of the BHM Kühner *et al.* (2000); Ejima *et al.* (2011), in the AHM
can be also determined numerically by DMRG using the following definitions
of the chemical potentials for finite system sizes

(28) |

where and denotes the corresponding ground-state energies.

Figure 4 shows the ground-state phase diagram of the 1D AHM, exhibiting MI and SF regions as a function of the chemical potential and the anyon transfer amplitude . The strong-coupling expansions of the chemical potentials via Eqs. (21) and (22) up to 6th order are compared with DMRG results. For small , both methods essentially agree up to [see Figs. 4(a) and (b)], while in the case of the pseudofermions [ in Fig. 4(c)] even 6th-order results start to deviate around . In the intermediate-coupling regime () sudden changes will appear in the perturbation results, especially for (not shown). This is because the perturbation expansions fall into the wrong minima as will be discussed in App. B.

### iii.3 Momentum distribution function

Anyons might be characterized most significantly by momentum
distribution functions as has been demonstrated for both hardcore Hao *et al.* (2008)
and softcore Tang *et al.* (2015) anyons.
For the current model, we can define two different types of single-particle correlation functions

(29) | |||||

(30) |

corresponding to boson or anyon representations. Results for the boson correlation function should be relevant for the proposed realization of the model in optical lattices. The anyon correlation function can be expressed in terms of the boson operators as follows:

(31) |

Within the strong-coupling expansion, the above correlators are translated according to

(32) | |||||

(33) |

In App. C we compare the perturbation results for the two-point correlation functions up to 4th order in with the iDMRG data.

The Fourier-transformed single-particle density matrices give the momentum distribution functions for bosons and anyons as

(34) | |||

(35) |

Up to and including 4th order in , we obtain for the momentum distribution functions of bosons

and for anyons

Taking the limit in Eqs. (LABEL:eq-nkb) or (LABEL:eq-nka), we obtain the momentum distribution function in the BHM

(38) | |||||

in agreement with the former studies of the strong-coupling
expansions in the BHM up to and including third order
in Freericks *et al.* (2009).

Using DMRG with PBCs, the momentum distribution functions of anyons and bosons can be extracted from Eqs. (34) and (35) after calculating the two-point correlation functions, as demonstrated in Fig. 5 by the comparison with the strong-coupling expansions (LABEL:eq-nkb) and (LABEL:eq-nka). While for analytical and numerical methods agree [Fig. 5(a)], small deviations appear for [Fig. 5(b)]. For [Fig. 5(c)] the oscillations become significant in the 4th-order strong-coupling expansions which are clearly an artifact.

Analogous to the momentum-dependent excitation energies
in Sec. III.2,
the characteristic asymmetry in the momentum
distribution functions can be understood by considering the
main -dependent contributions in the strong-coupling expansion
of .
In the BHM limit (), is always symmetric
about , where the position of the maximum is located.
These peak positions of [] shift to the negative
[positive] momentum with increasing for ,
which is consistent with the sign of in the cosine term of the
main -dependent contribution of ,
i.e., the positive [negative]
sign of in [] of
Eq. (LABEL:eq-nkb) [Eq. (LABEL:eq-nka)]. Moreover, the peak positions
of depend more strongly on than those of ,
this is because the -dependent main contribution in
shows up first in the 2nd-order expansion, while in the case of
it can already be seen in the 1st-order expansion.
To 1st order, the peak of is located at , that is, its position depends linearly on the fractional angle.
When increasing from to , the boson momentum distribution becomes flatter because of cancellations in the 1st-order terms. Close to the pseudofermion limit, a second peak appears for that can be attributed mainly to 2nd-order contributions. At , the cancellation of 1st-order terms becomes exact for all and one ends up with .
Our results for the boson momentum distribution function in the MI are in contrast to the results of Ref. Tang *et al.* (2015) for the SF where, depending on the filling, a single peak at either or has been found.

## Iv Summary and outlook

We studied the MI phase of the anyon-Hubbard model at filling factor one for arbitrary fractional angle using strong-coupling perturbation theory. Explicit expressions for the -dependence of the ground state energy per site, the single-particle and single-hole excitation energies as well as the momentum distribution functions were obtained for up to 6th order in (hopping/interaction).

In the BHM both single-particle and single-hole dispersions have their minimum at . For finite , the minimum of the dispersion is shifted differently for single-particle and single-hole excitations, that is, there is an indirect gap for particle-hole excitations.

The momentum distribution functions become asymmetric for with the peak shifted to negative (positive) momentum in the boson (anyon) description of the model. A stronger -dependence is found for the boson momentum distribution than for the anyon one. In particular, the boson momentum distribution function becomes almost flat in the pseudofermion limit .

While the series generated by the strong-coupling expansion might be asymptotic, the results for finite order agree well with numerically exact MPS calculations for small hopping . Increasing , however, the accuracy starts to deteriorate and the perturbative description is no longer sensible. At fixed order, the region of validity of the perturbative expansion seemingly decreases when the fractional angle is increased even though the MI region becomes larger.

Obviously, there are several other directions to extend our work.
The perturbation results in this paper are limited to the Mott
insulator for . First natural extension might be the
strong-coupling study in higher integer fillings. Another possibility is the inclusion of a nearest-neighbor interaction which leads to additional Haldane-insulator and density-wave phases Lange *et al.* (2017), the latter being susceptible to a perturbative treatment.
Furthermore, the strong-coupling approach could be applied to the AHM in higher dimensions. In this case perturbation theory should be particularly useful since no quasi-exact results from MPS-based methods are available.
Finally, for the comparison with future experiments it would be desirable
to investigate the dynamical quantities, such as the single-particle
spectral functions, the dynamical structure factor, and
the dynamical current and kinetic-energy correlation functions
as demonstrated in the BHM Ejima *et al.* (2012b, a); zu Münster *et al.* (2014).

## Acknowledgments

We thank F. Gebhard for valuable discussions. DMRG simulations were performed using the ITensor library ITe . This work was supported by Deutsche Forschungsgemeinschaft (Germany) through SFB 652.

## Appendix A Excitation energies in the Bose-Hubbard limit

Taking the limit in Eqs. (21) and (22) we obtain the single-hole and single-particle excitation energies in the momentum space for the BHM as

(39) | |||||

and

(40) | |||||

Direct comparison of Eqs. (39) and (40)
with iMPS results is demonstrated in Fig. 2.
Note that in Eqs. (24) and (25) of Ref. Ejima *et al.* (2012a)
the coefficient of the term in ,
as well as the coefficients of the and the terms in
are flawed; these errors were corrected in Eqs. (39) and (40).

## Appendix B Minima of excitation energies

In Fig. 4, chemical potentials , obtained from strong-coupling expansions, show a sudden increase for large , especially for [see, e.g., the results of Fig. 7(a), which were not included in Fig. 4.]. In this section we explain the origin of this shortcoming in detail.

Figure 7(b) shows the single-hole excitation energies by the 6th-order strong-coupling expansion (21) for and compared with iMPS results. In this intermediate-coupling region the perturbation results oscillate strongly, so that the position of the minimum for to estimate shift from negative (star symbol) to positive (cross symbol) momentum, while iMPS data still indicate that the minimum should be located at the negative momentum as in the case of . This sudden change of the location of minima leads to the artificial upturns of the strong-coupling expansions in the intermediate-coupling region of Fig. 4.

## Appendix C Correlation function

In this section we will give the strong-coupling expressions for the boson and anyon correlation functions. Note that in the anyonic case Eq. (31) should be taken into account. As explained in the main text, we employ the Harris-Lange transformation and obtain for the distance to 4 up to 4th order in as

(41) | |||||

(43) | |||||

and

(45) | |||||

(47) | |||||