Magnetic phases of spin-1 lattice gases with random interactions

TL;DR

This paper investigates how random spin-dependent interactions in a spin-1 lattice gas lead to complex magnetic phases, including disordered ferromagnetic and dimer phases, using advanced numerical methods.

Contribution

It introduces a detailed analysis of disordered magnetic phases in spin-1 lattice gases with local ferromagnetic or antiferromagnetic interactions, revealing new disordered phases.

Findings

Disordered ferromagnetic phase with spin-glass order parameter

Absence of random singlet order in certain regions

Identification of phase boundaries between ferromagnetic and dimer phases

Abstract

A spin-1 atomic gas in an optical lattice, in the unit-filling Mott Insulator (MI) phase and in the presence of disordered spin-dependent interaction, is considered. In this regime, at zero temperature, the system is well described by a disordered rotationally-invariant spin-1 bilinear-biquadratic model. We study, via the density matrix renormalization group algorithm, a bounded disorder model such that the spin interactions can be locally either ferromagnetic or antiferromagnetic. Random interactions induce the appearance of a disordered ferromagnetic phase characterized by a non-vanishing value of spin-glass order parameter across the boundary between a ferromagnetic phase and a dimer phase exhibiting random singlet order. The study of the distribution of the block entanglement entropy reveals that in this region there is no random singlet order.

Equations

$$H$$

$$-t\sum _{i,\sigma }\left(a_{i,\sigma }^{†}{a}_{i+1,\sigma }+a_{i+1,\sigma }^{†}{a}_{i,\sigma }\right)$$

$$n_i=\sum _{\sigma }{a}_{i,\sigma }^{†}{a}_{i,\sigma }$$

$${\mathbf{S}}_i=\sum _{\sigma ,{\sigma }^{\prime }}{a}_{i,\sigma }^{†}{\mathbf{T}}_{\sigma ,{\sigma }^{\prime }}{a}_{i,{\sigma }^{\prime }}$$

$$H=\sum _{i=1}^{L-1}{H}_i=\sum _{i=1}^{L-1}[J_{1i}\left({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}\right)+J_{2i}{\left({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}\right)}^2],$$

$${\theta }_i=\arctan \left(\frac{J_{2i}}{J_{1i}}\right)-\pi .$$

$$H_{BB}=J\sum _{i=1}^{L-1}[\cos \theta \left({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}\right)+\sin \theta {\left({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}\right)}^2],$$

$$H_{BB}=\sum _{i=1}^{L-1}{J}_i(\cos {\theta }_i({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1})+\sin {\theta }_i({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}{)}^2),$$

$$U_{2i}=U_{2C}+\eta (2{\zeta }_i-1),$$

$$D=[⟨H_i-H_{i+1}⟩{]}_D,$$

$$D_{ϵ}=-\frac{1}{L}\sum _{mn}\sin \left[\frac{\pi }{2}(m+n)\right][G_z(m,n){]}_D,$$

$$G_z(m,n)\equiv ⟨S_{z_m}{S}_{z_n}⟩-⟨S_{z_m}⟩⟨S_{z_n}⟩.$$

$$m_A=\frac{1}{L}\sum _i[⟨S_{zi}⟩{]}_D.$$

$$q=\frac{1}{L}\sum _i[⟨S_{zi}{⟩}^2{]}_D-[⟨S_{zi}⟩{]}_D^2.$$

$$E=-\mathrm{Tr}{\rho }_{\ell }{\log }_2{\rho }_{\ell }$$

$${\rho }_{\ell }={\mathrm{Tr}}_{L-\ell }|{\psi }_G⟩⟨{\psi }_G|,$$

$$H_0=\frac{U_0}{2}\sum _i{n}_i\left(n_i-1\right)+\frac{U_2}{2}\sum _i\left({\mathbf{S}}_i^2-2n_i\right)-\mu \sum _i{n}_i$$

$$|\left\{m\right\}⟩=\prod _i|n_i=n;s_i=s;m_i⟩,$$

$$H_0|\left\{m\right\}⟩=E_0|\left\{m\right\}⟩.$$

$$T=-t\sum _{i,\sigma }\left(a_{i,\sigma }^{†}{a}_{i+1,\sigma }+a_{i+1,\sigma }^{†}{a}_{i,\sigma }\right),$$

$$H_{\text{eff}}=H_0+\sum _{\left\{m\right\},{\left\{m\right\}}^{\prime }}|{\left\{m\right\}}^{\prime }⟩⟨\left\{m\right\}|\sum _{\alpha }\frac{⟨\left\{m\right\}|T|\alpha ⟩⟨\alpha |T|{\left\{m\right\}}^{\prime }⟩}{{\Delta }_{\alpha }},$$

$$H_0|\alpha ⟩=E_{\alpha }|\alpha ⟩$$

$${\Delta }_{\alpha }=E_0-E_{\alpha }.$$

$$⟨m_1+1,m_1-1⟩H_{\text{eff}}{m}_1,m_2$$

$$⟨m_1+2,m_1-2⟩H_{\text{eff}}{m}_1,m_2$$

$$D_1(m_1,m_2)$$

$$D_2(m_1,m_2)$$

$$W(m_1,m_2)$$

$$W_1(m_1,m_2)$$

$$⟨1,-1⟩H_{\text{eff}}0,0$$

$$⟨1,-1⟩H_{\text{eff}}-1,1$$

$$a_{\sigma }^{†}|n;s;m⟩$$

$$a_{\sigma }|n;s;m⟩$$

Full text

I Introduction

The interest for quantum magnetism has grown recently together with the possibility to simulate a huge variety of spin models with ultracold atoms trapped in optical lattices [1; 2]. The magnetic properties of the gas derive from their low-energy hyperfine level structure which provides an extra spin degree of freedom ${F=I+J}$, where $I$ is the nuclear spin and $J$ is the total electronic angular momentum. Several experiments have been realized with spinor atomic gases [3; 4], especially with alkali atoms, having electronic angular momentum $1/2$ and two hyperfine levels with ${F=\left|I\pm 1/2\right|}$. Here, we are interested in bosonic particles with lower energy hyperfine manifold corresponding to $F=1$; this is the case, for example, of ${}^{23}\textrm{Na}$, ${}^{87}\textrm{Rb}$, ${}^{41}\textrm{K}$, and ${}^{7}\textrm{Li}$ [4].

A system of bosons trapped in a deep optical lattice potential is well described by a Bose-Hubbard model, consisting of two competing terms: the hopping between lattice sites and the repulsive interaction produced by local scattering. If the spin orientation is fixed by an external magnetic field, as happens when the gas is confined in a magnetic trap, a scalar model is sufficient to describe the system [5; 6]. Conversely, if the spin orientation is not externally constrained, as in the case of optical trapping, the spinor character of the gas has to be taken into account.

Bosonic interactions, due to two-body $s$-wave collisions, are sensitive to the spin degree of freedom and contribute to the ordering properties at zero temperature. For spin-1 particles, two possible scattering channels open, one with total spin ${s=0}$ and one with total spin ${s=2}$, corresponding respectively to a scattering length $a_{0}$ and $a_{2}$. This leads to a spin-dependent interaction [7; 8; 9] which, depending on the scattering lengths $a_{0}$ and $a_{2}$, can be ferromagnetic or antiferromagnetic. Spin-1 condensates have been observed in Refs. [10; 11] and the two-body interaction has been studied in Refs. [12; 13; 14].

At zero temperature, the phase diagram of interacting spin-1 bosons on a lattice consists of the Mott insulator (MI) lobes as in the scalar case while for smaller interactions the system is superfluid (SF) [9]. In the MI phase, density fluctuations are suppressed while the spin degrees of freedom give rise to a variety of possible magnetic phases. In order to explore the magnetic behavior inside the MI phase, it is possible to map the problem into a spin Hamiltonian by a perturbative expansion in the hopping parameter. For unit filling, the MI phase has one boson per site, so the Bose-Hubbard spin-1 chain is mapped onto the bilinear-biquadratic spin-1 Hamiltonian [15; 16; 17; 18; 19; 20; 21], preserving the SO(3) symmetry of the hyperfine atomic splitting.

Among all the opportunities ultracold atoms can offer, there is the possibility to produce a disordered potential in a controlled way by extra speckle potentials [22; 23], optical superlattices of incommensurate frequencies [24; 25; 26; 27; 28; 29; 30; 31], holographic masks [32; 33], or adding different atomic species randomly trapped in sites distributed across the sample and acting as impurities [34; 35; 36]. In the so-called dirty bosons problems, the competition between the local disorder and the interactions gives rise to a Bose glass (BG) between the MI and SF which exhibits localization but, unlike the MI, is compressible and gapless. The phase diagram in the presence of this type of disorder has been studied both for spinless [5; 37; 38; 26; 39; 40] and for spinful bosons [41; 42; 43; 44; 45; 4; 46].

On the other hand, disorder can also be present in the interaction terms of the Hamiltonian. Such a possibility is more difficult to control and reproduce experimentally since it requires, for example, varying the scattering length of the $s$-wave collision at random. Recently, the Bose-Hubbard model in the presence of random interactions has been analyzed numerically in the context of many-body localization [47]. Magnetic Feshbach resonances can be exploited for this purpose [48; 49]; nevertheless, for spinor bosons, an external magnetic field would freeze the spin degrees of freedom and it would map the system into a spinless one. However, from a theoretical point of view, disorder in the interaction is especially interesting for spin-1 bosons, where different spin phases emerge inside the MI phase, contrary to the spinless case, where the BG separates MI to SF. Therefore, in order to implement random local interactions without freezing the spin degrees of freedom, different techniques are in order. Adopting optical Feshbach resonances could be a possible route. This technique employs a laser tuned near a photoassociation transition with an excited molecular level. The method was developed theoretically [50; 51; 52; 53] and has been proposed to quantum simulate frustrated magnetism with spinor Bose gases in 2D lattices[54]. Experimental observations of optical Feshbach resonances have recently been carried out with alkali-metal atoms [55; 56] and alkali-earth-metal-like atoms [57; 58; 59; 60], where atom losses have been shown to be reduced as suggested in Ref. [61]. Nevertheless, atomic losses are still limiting the lifetime of the gas below the millisecond regime and novel ways to address Feshbach resonances for spinor condensates are highly desirable.

In this paper, we study the case of a spin-1 system in the MI phase with one boson per site, when disorder mixes locally ferromagnetic and antiferromagnetic interactions. By mapping the spin-1 Bose-Hubbard system into the corresponding spin model (bilinear-biquadratic), we find numerically the ground state of the system for each disorder realization using the density matrix renormalization group (DMRG) [62; 63]. To identify the possible phases of the model, we analyze the ferromagnetic and the dimer order parameters, together with a normalized Edwards-Anderson (EA) order parameter [64]. Our results indicate that the phase diagram consists of three phases: a ferromagnetic phase, a random singlet phase with a nonzero dimer order parameter, and an intermediate phase between the first two phases. This intermediate phase exhibits nonvanishing EA parameter and we argue that it can be identified as a disordered “large-spin” regime [65]. By studying the entanglement entropy between two half chains, we exclude the existence of random singlet ordering in the large-spin phase. A linear scaling with the system size of the ferromagnetic domain walls suggests that the intermediate phase could be a locally disordered ferromagnet containing microscopic magnetized droplets.

The paper is organized as follows: In Sec. II, we describe the spin-1 Bose-Hubbard model and the corresponding bilinear-biquadratic spin-1 model for the MI phase and unit density. Starting from disordered two-body atomic interactions, we derive the corresponding random coefficients of the spin model. In Sec. III, we present numerical results for local order parameters, showing the existence, together with the ferromagnetic and the dimer phases, of a large-spin phase. In Sec. III.4, we analyze the entanglement properties of the system by studying the block entanglement entropy. In Sec. III.5, we study the scaling of the domain walls with the system size, and finally in Sec. IV, we summarize and present open questions.

II Model

We consider a one-dimensional chain of spin-1 bosons in an optical lattice. The effective Bose-Hubbard Hamiltonian for spin-1 bosons is [7; 9]

[TABLE]

with ${a}_{i,\sigma}^{\dagger}$ and ${a}_{i,\sigma}$ being the creation and annihilation operators for site $i$, spin component $\sigma=0,\pm 1$, respectively,

[TABLE]

is the total numbers of particles on site $i$, and

[TABLE]

the total spin on site $i$ with ${\mathbf{T}}_{\sigma,\sigma^{\prime}}$ the spin-1 angular momentum matrix elements. The on-site interactions are described by the first two terms of Eq. (1). The local density coupling is assumed to be site independent with coupling constant $U_{0}$. The second interaction term takes into account the local spin and we assume its coupling constant $U_{2i}$ to be disordered and site dependent. The third term is an on-site energy, with $\mu$ being the chemical potential, and finally the fourth term of (1) is the tunneling term describing hopping between nearest-neighbor sites with tunneling amplitude $t$. Our focus is on the strong interaction regime ${(t\ll U_{0})}$; further, we set $U_{0}$ as the energy scale unit ${(U_{0}=1)}$. Our analysis is restricted to the first Mott lobe with unit filling factor that extends, as it is easy to see in the atomic limit, in the interval ${U_{2}\in(-1;0.5)}$ [41; 42].

By considering a spin-1 chain of $L$ sites with open boundary conditions, deep in the Mott insulator, with strong interactions and unit filling per site, ${{n}_{i}=1}$, the Hamiltonian of Eq. (1) can be mapped into a spin-1 bilinear-biquadratic Hamiltonian:

[TABLE]

where ${\mathbf{S}_{i}=(S_{xi},S_{yi},S_{zi})}$ are the $i$th site angular momentum operators and ${J_{1i}=f\left(U_{0},U_{2i},U_{2(i+1)}\right)}$ and ${J_{2i}=g\left(U_{0},U_{2i},U_{2(i+1)}\right)}$ [see Appendix A for the derivation of this mapping] are the linear and quadratic exchange coefficients between spin sites $i$ and $i+1$.

For fixed $U_{0}$, distinct values of $U_{2i}$ and $U_{2(i+1)}$ in the region of interest correspond to a unique value of the ratio $J_{2i}$ to $J_{1i}$, which is conveniently parametrized by the angle $\theta_{i}$ defined as

[TABLE]

II.1 Homogeneous coupling constants

In the case of homogeneous exchange constants ${(U_{2i}=U_{2},\,\forall\,i)}$, the Hamiltonian of Eq. (4) becomes the general bilinear-biquadratic Hamiltonian [15; 16; 17; 18; 20; 19]:

[TABLE]

where $J_{1}$ and $J_{2}$ have been replaced using Eq. (5) and ${J=\sqrt{J_{1}^{2}+J_{2}^{2}}}$.

The ground-state phase diagram for the model with $H_{BB}$ is as follows: For ${\pi/2<\theta\leq\pi}$ and ${-\pi\leq\theta<-3\pi/4}$, the system is in the ferromagnetic phase characterized by a net spontaneous magnetization along a direction which breaks spontaneously space isotropy. For ${-3\pi/4<\theta<-\pi/4}$, the system is in the dimerized phase, characterized by a nonvanishing dimer order parameter. In this phase, translational invariance is broken and imperfect singlets of neighboring spins appear. For ${-\pi/4<\theta<\pi/4}$, the system is in the gapped Haldane phase exhibiting absence of any local order, a nonzero string order parameter and entanglement spectrum with even number of degeneracies. Finally, for ${\pi/4<\theta<\pi/2}$, the system is in the so-called critical phase characterized by quasi-long-range quadrupolar order.

The region of interest for ultracold atoms, ${U_{2}\in(-1;0.5)}$, corresponds thus to the region ${\theta\in(-\pi+\arctan(1/3);-\pi/2)}$ which encompasses part of the ferromagnetic phase, part of the dimer phase and the transition between them, occurring at $U_{2}=0$. The full phase diagram could be realized by preparing the highest energy state of 1D optical lattices or by using spin ladders [66].

The phase transition that occurs at ${U_{2}=0}$, which corresponds to ${\theta=-3\pi/4}$, is the focus of this paper. Chubukov conjectured [67] that the ferromagnetic and dimerized phases were separated by a possible spin nematic phase, but due to an abnormal divergence in correlation length this transition was extremely hard to categorize. The existence of the nematic phase was widely studied numerically [16; 18; 20; 9; 68] and it is now accepted that the extension of the nematic phase is extremely small, suggesting a direct first-order transition between the ferromagnetic and dimer phase.

II.2 Disordered coupling constants

We implement disorder in the spin system by considering random exchange coefficients $J_{1i}$ and $J_{2i}$ obtained by random interaction constants $U_{2i}$. The resulting Hamiltonian has the form

[TABLE]

with ${J_{i}=\sqrt{J_{1i}^{2}+J_{2i}^{2}}}$. In this work, disorder is confined to a range around a central point, $U_{2C}$, of width $\eta$. Thus $U_{2i}$ is generated for $i=1,2,\dots,L$ via

[TABLE]

where $\zeta_{i}$ is a uniformly distributed random variable between [math] and $1$, leading to $U_{2i}$ being uniformly distributed between ${[U_{2C}-\eta,U_{2C}+\eta]}$. To enforce locally the constraint ${U_{2}\in(-1;0.5)}$, we choose $U_{2C}\in(-1+\eta;0.5-\eta)$. Throughout this paper, the value of $\eta$ is fixed at $\eta=0.1$. A situation where disorder appears mostly on $U_{2}$ can be realized by considering anticorrelated fluctuations of the scattering lengths. Since the coupling constants depend on the $s$-wave scattering lengths $a_{s}$ $(s=0,2)$ as $U_{0}\propto a_{0}+2a_{2}$ and $U_{2}\propto a_{2}-a_{0}$ (details are given in Refs. [12; 13; 14]), it is sufficient to consider a fluctuation $\delta_{i}$ in the scattering lengths such that ${a_{2}^{(i)}=a_{2}+\delta_{i}}$ and ${a_{0}^{(i)}=a_{0}-2\delta_{i}}$.

III Numerical Results

Results were obtained using the finite-size DMRG with open boundary conditions; 500 random configurations of exchange coefficients were used. In order to reduce stability issues related to the degeneracy of the ground state, a weak magnetic field $(-10^{-5}S_{z1})$ was added to the first spin in the chain in order to promote one of the almost degenerate ground states. It is important to stress that even if the total angular momentum commutes with the random bilinear-biquadratic Hamiltonian Eq. (4), the minimum energy state can assume values, between $-L$ and $L$ of the projection of the total angular momentum ${M_{z}=\sum_{i}S_{zi}}$, depending on the disorder configuration. In the DMRG simulations we used five complete sweeps and kept 80 states during the renormalization procedure corresponding to a maximum discarded weight of $10^{-5}$. For all figures presented, the results obtained from a Hamiltonian with homogeneous exchange coefficients are labeled H and those with disordered exchange coefficients are labeled D.

III.1 Local magnetization and nearest-neighbor correlations

We start our analysis by computing the local magnetization $\langle S_{zi}\rangle$ and the nearest-neighbor spin-spin correlations $\langle S_{zi}S_{zi+1}\rangle$, shown in Fig. 1 for a single configuration of $\{\zeta_{i}\}$. For ${U_{2C}=-0.2}$, the system is deep in the ferromagnetic phase and the local magnetization is everywhere equal to $+1$ (thanks to the negative-symmetry-breaking field mentioned above), while the correlations are exactly $1$, showing perfect spin alignment. For ${U_{2C}=0.2}$, the system has zero magnetization and an overall negative, thus antiferromagnetic, nearest-neighbor correlations. The plot also shows an even-odd effect that can be explained, as we show in the next section, assuming that the ground state is dimerized. Finally for ${U_{2C}=0.02}$ there is an intermediate behavior: The system is not fully magnetized but exhibits islands of nonzero magnetization with large nearest-neighbor correlations surrounded by regions of very small magnetization with overall negative nearest-neighbor correlations. This is our first evidence of a disordered intermediate region between the dimer and ferromagnetic phases.

III.2 Dimerization versus ferromagnetism

To characterize the dimer phase, a common choice is the local dimer order parameter [16; 18]:

[TABLE]

where $\langle\,\rangle$ is the quantum ground state average and $[\,]_{D}$ is the disorder average. Because of the disordered nature of the chain we consider, we prefer to use a similar dimer order parameter $D_{\epsilon}$ defined in Ref. [69]

[TABLE]

where $G_{z}(m,n)$ is the two-point correlation function:

[TABLE]

While for the homogeneous model $D$ and $D_{\epsilon}$ have been shown to be equivalent [69], for the model with disorder, $D_{\epsilon}$ converges faster thanks to a self-averaging effect due to the summation on all sites. Moreover, $D_{\epsilon}$ can be observed experimentally with a quantum nondemolition measurement based on the Faraday effect [69].

As is customary, the average magnetization along the $z$ axis is defined as

[TABLE]

Figure 2(a) shows the average dimer order parameter and Fig. 2(b) shows the average magnetization of the spin-1 chain throughout the region of interest for both the homogeneous and disordered cases. For short lengths, we observe nonzero dimerization below the first-order transition point for the homogeneous model $(U_{2C}=0)$. Disorder slightly suppresses the dimer order parameter, which has a fairly large finite-size effect. On the other hand, magnetization does not show strong finite-size effects, but it is considerably affected by the presence of disorder which rounds the sudden jump observed in the homogeneous case. In Fig. 2(c), we show the results for the average dimer order parameter and magnetization after a finite-size extrapolation to infinite lengths (see details in Appendix B). In contrast to the homogeneous case, the results clearly show that the magnetization is continuous at $U_{2C}=0$.

III.3 Spin glass order parameter in the intermediate phase

To ascertain the nature of the intermediate phase, we compute the Edwards-Anderson order parameter, defined as [70]

[TABLE]

The term subtracted in Eq. (13) takes into account the possibility of finite magnetization in the disordered magnetic phase. It is clear that this parameter is zero in the homogeneous case and takes into account correlations of the local magnetization with the disorder variables. This parameter is shown in Fig. 3 for different chain lengths. The spin glass order parameter is not affected by finite-size effects and its behavior shows that the maximum of correlations between local magnetization and disorder occurs for $U_{2C}>0$. From the Edwards-Anderson order, parameter we infer a range for the intermediate phase of $-0.051\lesssim U_{2C}\lesssim 0.073$.

The question thus arises as to whether the intermediate phase that spans from the vanishing of the ferromagnetic phase to the onset of the dimer phase is a new phase that can be characterized, for example, by random singlet order or by disorder induced magnetic domains.

In a previous work, Quito et al. [65], using the strong-disorder renormalization group (SDRG) [71], found a large-spin phase separating the ferromagnetic phase and the random singlet phase. The term large-spin is in reference to the renormalization procedure to find the ground state of the spin chain.

Here we summarize the renormalization group procedure of Quito et al. [65]. The SDRG uses a decimation procedure similar to the procedure originally pioneered by Dasgupta and Ma [72; 73]. Each decimation replaces the spin pair corresponding to the largest energy gap between the ground state and the first excited state of the local Hamiltonian, $H_{i}$, with a new effective spin. For a spin-1 chain, the Hamiltonian $H_{i}$ corresponds to a spin pair located at sites $i$ and $i+1$. If their coupling is ferromagnetic, i.e., $U_{2i},\,U_{2(i+1)}\in(-1;0)$, the pair of spins is replaced by a new effective spin of size spin 2 and new effective couplings are produced between the new spin and the spins on sites $i-1$ and $i+2$. For $U_{2i},\,U_{2(i+1)}\in(0;0.5)$, corresponding in the homogeneous model to the dimer phase, the two spins form a singlet. In the decimation procedure, a spin pair of this type is removed and a new effective coupling between the spins on sites $i-1$ and $i+2$ is found. Repeating this decimation procedure will eventually reduce a chain in which all $U_{2i}$ lie in the range $(-1;0)$ to a single effective spin of size $L$, whereas a chain in which all $U_{2i}$ lie in the range $(0;0.5)$ will reduce to an effective spin of size [math]. This large-spin phase describes a region in which a chain may have a final effective spin which lies somewhere between [math] and $L$, indicating that it is neither fully ferromagnetic nor fully dimerized.

Our DMRG calculations also show that in the intermediate phase different instances of disorder exhibit a wide distribution of total magnetization, precisely as predicted in Ref. [65]. This indicates that the intermediate phase could coincide with the large-spin phase found in Ref. [65]. It should be noted, however, that Quito et al. [65] modeled their disorder in such a way that the random coefficients $J_{1i}$ and $J_{2i}$ were independent of each other, while in this work $J_{1i}$ and $J_{i2}$ are both dependent on $U_{2i}$ and $U_{2(i+1)}$.

III.4 Entanglement entropy

Further to the three order parameters discussed previously, we consider the block entanglement entropy. The advantage of this bipartite entanglement measure is that as well as being discontinuous at a first-order transition [74; 75] its spatial distribution in a disordered model can give insight into the bonds existing between blocks of spins. In the so-called random singlet phase, found for instance in the random coupling $XX$ model and in the random transverse-field Ising model for spin-1/2 [76; 77], the probability distribution of the entanglement entropy is characterized by pronounced peaks at integer values, indicating the occurrence of maximally entangled states of pairs of spins belonging each to one side of the partition. Very recently the scaling of the entanglement entropy distribution away from integer points has been studied in Refs. [78]. Finding the probability distribution of the entanglement entropy is particularly interesting for our problem because we can confirm the existence or absence of a random singlet phase in this model.

The block entanglement entropy measures the entanglement between two blocks of lengths $\ell$ and $L-\ell$ and can be calculated as the von Neumann entropy of the reduced density matrix of one of the two blocks [79]:

[TABLE]

with

[TABLE]

where $|\psi_{G}\rangle$ is the ground state of Hamiltonian (4) and we set $\ell=L/2$.

In Fig. 4 we show the results for the average entanglement entropy in the disordered model and the entanglement entropy for the homogeneous case. In agreement with the magnetization reported in Fig. 2, $E$ exhibits a strong discontinuity at ${U_{2C}=0}$ in the homogeneous case while varying continuously in the presence of disorder. It is interesting to note that while the disordered system is unable to reach the entropy values obtained in the weak dimerized regime ${0<U_{2C}\ll 1}$, the disordered system has on average a nonzero entropy in the ferromagnetic phase. This confirms the results Fig. 2 suggests: For negative $U_{2C}$ the chain is not fully magnetized and spins show strong correlations signaled by a nonzero entanglement entropy.

In Fig. 5, we show histograms of the entanglement entropy for ${U_{2C}=0.05,0.1,0.2,0.39}$. These values correspond in the homogeneous system to the dimer phase. Since we are interested in detecting the random singlet phase, we plot the histogram so that the horizontal axis displays the ratio $E/E_{S}$, where ${E_{S}=\log_{2}(3)\simeq 1.585}$ is the entanglement of a spin-1 singlet (which is maximally entangled state as for spin-1/2 particles). We find that the entanglement distribution exhibits a broad peak for ${E>E_{S}}$, whose midpoint decreases for increasing $U_{2C}$ and increases for increasing $L$. For ${U_{2C}=0.39}$, the distribution is much wider. In the cases considered, we observe the major peak tending towards a value of $E_{S}$ as predicted for the random singlet phase for the model presented in Ref. [65]. However, probably because of the finite lengths considered in our simulations, we do not observe peaks occurring at other integer multiples of $E_{S}$.

III.5 Domain walls

Thus far we have shown that, through the introduction of disorder, an intermediate phase can be observed between the ferromagnetic and dimer phases. This phase is characterized by a nonzero EA order parameter coexisting with a nonzero magnetization and a possible exponentially small dimer order parameter. Let us concentrate on the local magnetization properties. From our analysis in Sec. III.1, we know that in the ferromagnetic phase neighboring spins are aligned, thus having ${\langle S_{zi}S_{zi+1}\rangle>0}$. Conversely, in the dimer phase, spins tend to form imperfect singlets and thus have a tendency to antialign, giving rise to a negative value of $\langle S_{zi}S_{zi+1}\rangle$. Thus in both of these phases, the sign of $\langle S_{zi}S_{zi+1}\rangle$ is constant and by observing sign changes in its value we can detect a domain wall.

Figure 6 shows the average density of domain walls, $\bar{\rho}_{d}$, occurring in a given chain of length $L$ for different $L$. A thermodynamically large number of domain walls, corresponding to small droplets of magnetized spins, is found in the intermediate region with the maximum density reached around $U_{2C}\simeq 0.017$, indicating that the dimer phase is more susceptible to disorder than the ferromagnetic phase. This could be ascribed to the much smaller gap in the dimer phase as $U_{2C}\to 0^{+}$. By comparing Fig. 6 and Fig. 3, we see that the density of domain walls and the EA order parameter follow the same qualitative behavior. In other words, they are indeed markers of the intermediate phase. It is easy to show that within small disorder expansion, the EA order parameter is proportional to the local susceptibility with respect to changes in the local spin-spin coupling. On the other hand, a sizable amount of domain walls measure how the ground state locally changes from antiferromagnetic to ferromagnetic short-range correlations due to different values of the spin couplings, which qualitatively explains the similar behavior of the two quantities.

IV Conclusions

In conclusion, we analyzed the zero-temperature phase diagram of spin-1 bosons in one-dimensional optical lattices with disordered interactions generated with uniformly distributed random scattering lengths. We found that between the ferromagnetic and dimer phases occurring in the clean case, there exists an intermediate phase showing features of a disordered ferromagnet and characterized by a finite EA order parameter. It displays disorder-induced magnetized microscopic droplets whose concentration, as well as the EA parameter, does not show a strong dependence on the system size. The strength of the local magnetization in each droplet is locally correlated with the disorder interactions for each site. We cannot exclude that a vanishingly small dimer order is still present in this disordered phase. If this were the case, there exists the possibility that a weak local dimer order parameter is present in the regions where the magnetization is smaller. Thus we cannot exclude the existence of a Griffiths-type phase between the intermediate phase and the random singlet phase.

We thus conjecture that this intermediate phase coincides with the large-spin phase predicted by Quito et al. [65] by observing that the average magnetization per particle has an intermediate value between zero and one.

The intermediate disordered phase could potentially be close to the nematic region conjectured by Chubukov [67]. It would therefore be interesting to analyze the stability of the intermediate phase with respect to the appearance of the nematic phase when a uniaxial field is added [19]. This study is beyond the scope of the current paper and we leave it for future investigations.

In analogy with the classical version of the bilinear-biquadratic spin-1 model [80], which for infinite range interactions exhibits a spin glass phase, one can ask if also in the present case the disordered ferromagnet displays some nonequilibrium properties typical of spin glasses. It remains an open problem to determine if the replica symmetry is actually broken and in order to give a precise answer a further study of the long-time dynamics would be in order. We think this is an interesting subject for a subsequent project.

Acknowledgements.

We acknowledge helpful discussions with C. Cartwright. K.D.M. wishes to thank the Department for Employment and Learning (DEL) and the Professor Caldwell Travel Studentship for support. S.P. is supported by a Rita Levi-Montalcini fellowship of MIUR. A.S. acknowledges financial support from the Spanish MINECO Projects FIS2013-40627-P and FIS2016-80681-P and the Generalitat de Catalunya CIRIT (2014-SGR-966).

Appendix A DERIVATION OF THE EFFECTIVE SPIN-HAMILTONIAN

In this Appendix, we explain how to derive the bilinear-biquadratic Hamiltonian Eq. (4) from the spin-1 Bose-Hubbard Hamiltonian Eq. (1) in the MI phase by second-order perturbation theory with respect to the tunneling amplitude.

In the MI phase, at zero temperature and without any disorder, the zero-order ground state of (1) is given by the degenerate ground state of the unperturbed interaction Hamiltonian

[TABLE]

given by the product state

[TABLE]

where $n$ and $s$ are fixed and the same for each site and $m_{i}$ is free to change. The configuration of all the $m_{i}$s in the lattice is denoted with $\left\{m\right\}$ and

[TABLE]

We focus on the case with filling $n=1$ and $s=1$: A small amount of hopping does not affect the on-site density so we can construct an effective Hamiltonian $H_{\text{eff}}$ acting only in this subspace. To determine $H_{\text{eff}}$, one can perform a second-order perturbation expansion in the hopping term

[TABLE]

following the procedure described in Ref. [81]

[TABLE]

where $\ket{\alpha}$ are intermediate eigenstates of $H_{0}$ with

[TABLE]

and

[TABLE]

The effective Hamiltonian is expected to maintain the same SO(3) symmetry of the original model, so it must have the form of the bilinear-biquadratic Hamiltonian $H_{BB}$ (4) with $J_{1i}$ and $J_{2i}$ depending on $t$, $U_{0}$, $U_{2i}$, and $U_{2(i+1)}$. Following Ref. [7], the coefficients $J_{1i}$ and $J_{2i}$ are determined by comparing $H_{BB}$ with (20).

In our case, the spin interaction $U_{2i}$ is a site-dependent disordered variable, thus giving rise to coefficients $J_{1i}$ and $J_{2i}$ dependent on the bond ${(i;i+1)}$. To determine them explicitly, we consider a two-site system (only one bond) with a vector basis $\ket{\left\{m\right\}}=\ket{m_{1},m_{2}}$ and calculate the matrix elements

[TABLE]

with

[TABLE]

and

[TABLE]

To our purpose, it is sufficient to consider only two elements:

[TABLE]

For the calculation of the matrix elements, we use (20), taking into account that the creation and annihilation operators act as

[TABLE]

with coefficients $A$, $B$, $\bar{A}$, and $\bar{B}$ given in Ref. [8]. There are four involved intermediate states:

[TABLE]

with

[TABLE]

from which one obtains

[TABLE]

The same calculation can be done for every bond and, putting together (A), (28), and (29), one obtains

[TABLE]

The general relations for the coefficients with fixed odd $n$ and $s=1$ are

[TABLE]

[TABLE]

Appendix B SUPPLEMENTARY FIGURES

In this section of the Appendix, we show an example of the extrapolation analysis for the dimer order parameter, mentioned in Sec. III.2. The piecewise cubic Hermite interpolating polynomial (pchip) extrapolation method in MATLAB is used to extrapolate the values of the order parameters obtained for $1/L=1/12,1/24$, and $1/48$ to the case where $1/L=0$, therefore obtaining the corresponding value for the order parameter at infinite chain length. Figure 7 shows how the extrapolated values of the dimer order parameter were obtained for both the homogeneous and disordered cases at $U_{2C}=0.1$.

Figure 8 shows the results of this method when applied to the dimer order results from Fig. 2. These results, as previously discussed, show a significant reduction when compared to the values obtained for the dimer order parameter at finite lengths.