Phase transitions of the coherently coupled two-component Bose gas in a square optical lattice

TL;DR

This paper explores how coherent coupling influences phase transitions in a two-component Bose gas in an optical lattice, revealing shifts in the Mott insulator to superfluid transition and mapping to an effective spin model.

Contribution

It introduces a detailed analysis of the effects of coherent light-matter interactions on quantum phases using advanced theoretical methods.

Findings

Coherent coupling shifts the Mott lobe tip to higher tunneling values.

Strong inter-species repulsion affects the superfluid-Mott transition.

The Bose-Hubbard model maps onto an effective spin Hamiltonian.

Abstract

We investigate properties of an ultracold, two-component bosonic gas in a square optical lattice at unit filling. In addition to density-density interactions, the atoms are subject to coherent light-matter interactions that couple different internal states. We examine the influence of this coherent coupling on the system and its quantum phases by using Gutzwiller mean field theory as well as bosonic dynamical mean field theory. We find that the interplay of strong inter-species repulsion and coherent coupling affects the Mott insulator to superfluid transition and shifts the tip of the Mott lobe toward higher values of the tunneling amplitude. In the strongly interacting Mott regime, the resulting Bose-Hubbard model can be mapped onto an effective spin Hamiltonian that offers additional insights into the observed phenomena.

Equations

$${\hat{H}}_{BH}=$$

$$+\frac{1}{2}\sum _i^L\left(U_a{\hat{n}}_{ia}({\hat{n}}_{ia}-1)+U_b{\hat{n}}_{ib}({\hat{n}}_{ib}-1)\right)$$

$$+\sum _i^L\left(U_{ab}{\hat{n}}_{ia}{\hat{n}}_{ib}-\Omega {\hat{a}}_i^{†}{\hat{b}}_i^{}+h.c.\right),$$

$${\tilde{n}}_i=\frac{⟨{\hat{n}}_{ia}⟩-⟨{\hat{n}}_{ib}⟩}{⟨{\hat{n}}_{ia}⟩+⟨{\hat{n}}_{ib}⟩}.$$

$${\hat{a}}_i^{†}{\hat{a}}_j^{}\approx {\varphi }_i^{\ast }{\hat{a}}_j^{}+{\varphi }_j{\hat{a}}_i^{†}-{\varphi }_i^{\ast }{\varphi }_j,$$

$$\left(\begin{array}[]{cc}-\mu&-\Omega\\ -\Omega&-\mu\end{array}\right).$$

$${\hat{H}}_{\text{eff}}=$$

$$-J_z\sum _i{\hat{S}}_i^z-2\Omega \sum _i{\hat{S}}_i^x,$$

$$\frac{{\Omega }_c}{t}\propto \left(1-\frac{1}{\gamma }\right)\frac{t}{U}.$$

$$J_{zz}$$

$$J_z$$

$$J_{\perp }$$

$$\chi =J_{zz}+|J_z|-J_{\perp }-\Omega .$$

Full text

I Introduction

Due to their highly controllable properties, systems of ultra-cold atoms are promising platforms for quantum simulations. One of the early successes in this direction was the observation of a superfluid-Mott insulator transition in a lattice Bose gas Greiner et al. (2002); Bloch et al. (2008), as a prototype of a quantum phase transition. With recent advances in experimental techniques, present-day cold atom experiments feature finite-range interactions, for example in Rydberg dressed systems Löw et al. (2012); Schauß et al. (2015), as well as artificial gauge potentials that mimic magnetic fields Dalibard et al. (2011); Goldman et al. (2014). All these achievements bring these setups closer to simulating complex condensed matter systems.

Multi-component systems, such as mixtures of different atoms or different hyperfine states of the same atomic species, introduce additional degrees of freedom that can be treated as pseudo spin. In the weakly interacting limit, depending on the ratio of intra-component and inter-component interactions, two-component mixtures may exhibit phase separation Pethick and Smith (2002). The process of phase separation is substantially altered by introducing a coherent coupling term that enables a conversion of one internal atomic state into the other Blakie et al. (1999); Search and Berman (2001); Tommasini et al. (2003); Lee et al. (2004); Merhasin et al. (2005); Abad and Recati (2013); Butera et al. (2017); Chen and Malomed (2017); Shchedrin et al. (2016). In advanced experiments this conversion is commonly implemented by an external laser that couples two internal atomic states as, for example, in Rydberg systems Löw et al. (2012).

Effects of the coherent coupling are much less explored in the regime of strong repulsive interactions. In this paper, we address this topic in detail for the case of a two-component Bose gas on a square lattice. Strong interactions drive a transition from a condensate into a Mott insulator state, and in combination with coherent coupling, the system gives rise to a rich phase diagram that we investigate below. The subject has received some attention recently: DMRG studies for one-dimensional lattices have been performed Barbiero et al. (2016); Zhan et al. (2014) and the possibility of an effective $xy$-antiferromagnetic state in two dimensions has been explored Graß et al. (2014). In the recent papers de Forges de Parny et al. (2015); de Forges de Parny and Rousseau (2017) it was shown that the coherent coupling between atoms and molecules can stabilize an additional quantum phase.

The remainder of the paper is organized as follows: in the next section we introduce the main model of our study and briefly outline the methods that we use throughout the paper. In Sec. III we address the case of finite bosonic coherences. In particular we examine the neutral to polarized phase transition on top of the underlying condensate. The results for the case of stronger local interactions that lead to the superfluid-Mott transition with coherent coupling are presented in Sec. IV. The case of imbalanced hopping amplitudes is the subject of Sec. V. Finally, our main conclusions are summarized in Sec. VI.

II Model and Methods

We investigate the phase diagram of an extended Hubbard model describing a two-component Bose gas in a square optical lattice with an additional coherent coupling term. In second quantization the model explicitly reads

[TABLE]

where $\hat{a}_{i}^{(\dagger)}$ and $\hat{b}_{i}^{(\dagger)}$ are the annihilation (creation) operators for bosonic species $a$ and $b$ respectively, and $\hat{n}_{ia}$ and $\hat{n}_{ib}$ are their number operators. The tight binding hopping amplitudes of the respective species are denoted as $t_{a}$ and $t_{b}$ and $U_{a}$ and $U_{b}$ are their on-site interaction. Interactions are generally assumed to be local and repulsive in this study. An on-site interaction between particles of different species is designated as $U_{ab}$ and the term proportional to $\Omega$ allows for the conversion of one bosonic species into another on the same site. This last term is called the coherent coupling term. Local terms in the Hubbard Hamiltonian are summed over lattice sites $1,\ldots,L$, and in the first term $<i,j>$ stands for the sum over all nearest neighbor sites $i$ and $j$ of the square lattice. Throughout the paper we will consider the case of $U_{a}=U_{b}=U$ and we set all scales by fixing $U=1$. Our aim in this paper is to investigate possible ground states of the model (1).

It is well known that by reducing the ratio of the tunneling amplitude $t_{a,b}$ over the repulsive interaction strength $U$ at commensurate lattice filling, bosons exhibit a superfluid to Mott-insulator transition Fisher et al. (1989); Bloch et al. (2008). On top of this, an effective magnetic ordering emerges on the Mott side in two-component bosonic mixtures Altman et al. (2003); Kuklov and Svistunov (2003). These results are typically obtained for $\Omega=0$ in the regime of $\gamma\leq 1$, where we introduce the ratio between inter- and intra-species interaction $\gamma=U_{ab}/U$. In the opposite case of $\gamma\geq 1$, it is energetically favorable for the two bosonic species to be spatially separated Pethick and Smith (2002). However, a finite $\Omega$ will suppress this tendency by enabling a conversion between the two species. In this way, it allows us to address the regime $\gamma>1$. Note that the model (1) conserves only the total number of particles and not the particle numbers of each species separately. In the limiting case $\gamma\gg 1$, only atoms of one species are present as the system avoids the high energy cost of $U_{ab}$.

This reasoning already suggests that in addition to the well understood superfluid to Mott phase transition, for the model (1) we are able to distinguish another phase transition characterized by the polarization order parameter

[TABLE]

A strong interspecies interaction $U_{ab}$ favors a polarized phase with finite $\tilde{n}_{i}$. In contrast, the $\Omega$ term favors strong local coherence $\langle a_{i}^{\dagger}b_{i}+b_{i}^{\dagger}a_{i}\rangle$ that corresponds to a neutral phase with $\tilde{n}_{i}=0$. To establish boundaries between different phases as a function of the physical parameters of the Hamiltonian (1), we will use two approximate methods that we briefly outline here.

Features of the lattice Bose gas can be explored conveniently by means of the Gutzwiller mean-field theory Rokhsar and Kotliar (1991); Krauth et al. (1992), which amounts to decoupling non-local terms as

[TABLE]

where $\phi_{i}^{(*)}=\langle\hat{a}_{i}^{(\dagger)}\rangle$ is a condensate order parameter that is obtained in a self-consistent way. The approximation becomes an exact description in several limits: in the limit of infinite lattice coordination number $z$, in the atomic limit ($t=0$) and in the weakly interacting limit, in the superfluid phase. However, for a vanishing condensate order parameter the lattice sites are completely decoupled whithin Gutzwiller mean-field theory and the description of the Mott domain is oversimplified. In order to go beyond this limitation, we will use bosonic dynamical mean field theory (BDMFT) Hubener et al. (2009); Anders et al. (2010, 2011); Li et al. (2011); Snoek and Hofstetter (2013).

Formally, BDMFT is derived as a second order expansion of the full model (1) in terms of the inverse lattice coordination number. In comparison to Gutzwiller mean-field theory, we increase the order of the expansion by one Freericks et al. (2009). The approximate effective problem obtained in this way is given by a bosonic Anderson impurity model. Parameters of the effective model are set by imposing a self-consistency in terms of the condensate order parameter and the local Green’s function. The main approximation is the assumption of the locality of self-energies, which is consistent with the second order expansion in the inverse of the coordination number Freericks et al. (2009); Metzner and Vollhardt (1989). In order to take into account states that break translational symmetries, in this paper we use real-space BDMFT Li et al. (2011) and address a lattice size of $6\times 6$, with a special focus on a possible two-sublattice ordering.

Both the Gutzwiller mean-field theory and BDMFT are implemented in the grand canonical ensemble. To this end, we introduce a single chemical potential $\mu$, $\hat{H}_{BH}\rightarrow\hat{H}_{BH}-\mu\sum_{i}^{L}\left(\hat{a}^{\dagger}_{i}\hat{a}^{\phantom{\dagger}}_{i}+\hat{b}^{\dagger}_{i}\hat{b}^{\phantom{\dagger}}_{i}\right)$, as the model (1) conserves only the total number of particles. In the next sections we present and discuss results obtained by using these approaches for describing possible ground states of the model (1).

III The Superfluid Regime

In this section we present the Gutzwiller analysis of the coherently coupled spinor Bose gas in the superfluid phase. In order to address states with a finite condensate fraction, we choose relatively high tunneling amplitudes $t_{a}/U=t_{b}/U=1/4$. As the neutral to polarized phase transition is driven by the strong inter-component interactions $U_{ab}$, we plot the polarization order parameter $\tilde{n}_{i}$, defined in Eq. (2), over the ratio $\gamma=U_{ab}/U$ in Fig. 1a).

Results for the model (1) without coherent coupling are plotted for reference (red solid line) and the transition from $\tilde{n}_{i}=0$ to $|\tilde{n}_{i}|=1$ is found at $\gamma_{c}=1$. We understand from previous works, that it is energetically favorable to have components of both species on each lattice site only in the case of weak inter-species repulsion. Strong interspecies repulsion leads to the polarized phase, where only particles of one species can be found on a single lattice site. For fixed atom densities the system will thus undergo phase separation. We notice that positive and negative values of $\langle\hat{n}_{a,i}\rangle-\langle\hat{n}_{b,i}\rangle$ ($\tilde{n}_{i}=\pm 1$) appear equally as results of numerical calculations with different initial conditions. This indicates two degenerate ground states in the polarized phase, as we will confirm in the following. In these calculations, initial parameters of the self-consistent loop or root search routines determine which of the ground states is selected, whereas in actual experiments the occurrence probabilities for both ground states are equal.

We now turn to the effects of a finite coherent coupling term and set $\Omega/U=0.1$. While the polarization order parameter changes abruptly for vanishing coherent coupling (red solid line in Fig. 1a), it exhibits a continuous change for finite coherent coupling (blue dashed line in Fig. 1a). Moreover, we notice that the coherent coupling shifts the transition point $\gamma_{c}$ to higher values of the inter-species interaction. The same qualitative behavior was reported in Refs. Abad and Recati (2013); Barbiero et al. (2016) in the quasi one-dimensional geometry for the regime of weak interactions, both with and without a lattice. Within Gross-Pitaevskii theory, it was found analytically that the polarized phase sets in at $U_{ab}^{c}=U_{a,b}+2\Omega/n$. For the parameters given in Fig. 1, this relation yields a transition point at $\gamma_{c}=1.2$. However, since the mentioned derivation is strictly valid only in the weakly interacting limit, where all bosons are condensed, the phase transition in Fig. 1 is expected to appear at a slightly different value of $\gamma$. In particular, we find that as the ratio $t/U$ is lowered further, the transition point between the polarized and neutral phase is shifted in favor of the neutral phase. The region of the neutral phase extends toward higher values of $\gamma$ and deviations with respect to the result obtained in the weakly interacting limit become more pronounced. This effect will be further explored in the next section.

The observation of discontinuities in the polarization order parameter (Fig. 1a)) draws our attention to the order of the observed phase transitions. We analyze this in Fig. 1c) and d), where the grand potential and its first derivative are plotted as functions of $\gamma$. As we use the grand canonical description at zero temperature and explicitly include the chemical potential term, the grand potential is given by the expectation value of our mean-field Hamiltonian. For $\Omega=0$, we find that the first derivative of the grand potential is discontinuous at $\gamma_{c}$. This leads us to the conclusion, that the neutral to polarized phase transition is of first order for $\Omega=0$. In contrast, we observe a cusp in the same derivative for finite $\Omega$, implying that the phase transition is of second order.

To explore this in more detail, we plot the lowest energy eigenvalue of the mean field Hamiltonian as a function of the condensate order parameters in Figs. 2 and 3 (left series). In both cases, for vanishing and for finite $\Omega$, we find a single energy minimum for the neutral phase at $\gamma<\gamma_{c}$, see Figs. 2a1) and 3a1). The condensate order parameters $\phi_{a}$ and $\phi_{b}$ corresponding to this minimum are equal. In the polarized regime for $\gamma>\gamma_{c}$, however, two degenerate energy minima are present (Figs. 2a4) and 3a4)). The degeneracy stems from the fact that the ground state breaks the symmetry between the two species, while the Hamiltonian is symmetric with respect to the interchange of these two species. In the polarized phase with $\Omega=0$, the condensate order parameter of one of the species is strictly zero, while the other one has a finite value. In contrast, both order parameters are finite but non-equal at finite values of $\Omega$ as also shown in Fig. 1b).

Having established the properties of the neutral and polarized phase, we now discuss the transition between them.

For $\Omega=0$ and $\gamma=\gamma_{c}$, we find an infinitely degenerate energy minimum for constant $\sqrt{\phi_{a}^{2}+\phi_{b}^{2}}$. In the vicinity of the transition point, a single minimum at $\phi_{a}=\phi_{b}$ for $\gamma<\gamma_{c}$ abruptly transforms into two minima found at ”distant” positions in the $\phi$-space for $\gamma>\gamma_{c}$. The lowest eigenvalue of the Gutzwiller mean-field Hamiltonian plotted along the unit density line over $\alpha=\arctan(\phi_{a}/\phi_{b})$ (right series of Fig. 2) provides the reason for that. The energy has a parabolic shape around $\alpha=45^{\circ}$. The coefficient of this parabola changes from being positive in the neutral phase, which results into a minimum at $\alpha=45^{\circ}$, to being negative in the polarized phase, which results into two well separated minima at $\alpha=0^{\circ},90^{\circ}$.

At finite $\Omega$, we do not find an abrupt change in the minima at the transition point. The neutral minimum at $\phi_{a}=\phi_{b}$ splits into two degenerate minima that evolve towards their final values with increasing $\gamma$. Accordingly, the energy plot along the unity density line (right series, Fig. 3) is not a parabola that simply flips the sign of its coefficient, but one that develops a bump at $\alpha=45^{\circ}$ and thus gradually shifts its minima away from that point.

IV Superfluid-Mott Transition

Strong on-site interactions suppress density fluctuations and deplete the condensate. At commensurate densities, this mechanism drives a transition from a superfluid into a Mott insulator state. In this section we map out the phase diagram of the model defined in Eq. (1) as a function of the tunneling amplitude $t/U$ and coherent coupling $\Omega/t$ by using BDMFT at zero temperature.

In the limit of vanishing coherent coupling $\Omega\rightarrow 0$, for $\gamma>1$ (inter-species interaction stronger than intra-species interaction) our BDMFT simulations, implemented in the grand-canonical ensemble, recover the well-known results for the Mott-insulator to superfluid transition for a single bosonic species on a square lattice Capogrosso-Sansone et al. (2008); Anders et al. (2010). At finite $\Omega$ and $t=0$ (the atomic limit) and at a total filling $\langle n_{ai}+n_{bi}\rangle=1$, the ground state has no polarization. This can easily be seen by considering a single-site Hamiltonian in the subspace spanned by $|1,0\rangle=a_{i}^{\dagger}|0\rangle$ and $|0,1\rangle=b_{i}^{\dagger}|0\rangle$, which is given by

[TABLE]

The ground state is $|\text{GS}\rangle=(|0,1\rangle+|1,0\rangle)/\sqrt{2}$ and it is neutral, since $\langle\text{GS}|n_{ai}-n_{bi}|\text{GS}\rangle=0$. In the following we investigate a range of tunneling amplitudes $t_{a}/U=t_{b}/U=t/U\in[0.001,0.1]$ and coherent couplings $\Omega/t\in[0,1]$. We scan phase diagrams spanned by $t/U$ and $\mu/U$ at fixed $\Omega/t$ to access points with total unit filling as shown in Fig. 4 for $\gamma=8/5$.

We begin our analysis with a small fixed ratio $\Omega/t$. In Fig. 4a) we present the absolute value of the polarization order parameter $|\tilde{n}|$ (Eq. (2)) for $\Omega/t=0.12$. Increasing the tunneling $t/U$ from a starting point near the atomic limit, we encounter a transition from the neutral into the polarized state, before reaching the tip of the first Mott lobe. With further increase of $t/U$, we find a second transition from the polarized Mott state into the polarized condensate state (see Fig. 4a)). At stronger $\Omega/t\geq 0.28$, we find that the whole Mott lobe is neutral and around the tip of the lobe, we have a transition from the neutral Mott directly into the polarized superfluid. An example of this behavior is shown in Fig. 4b). Finally, for $\Omega/t>0.44$, there is a transition from the neutral Mott insulator into the neutral superfluid, followed by a second transition from a neutral into a polarized superfluid state, see Fig. 4c. The change in the polarization on the superfluid side of the diagrams as a function of the chemical potential $\mu$ (Figs. 4a-c, a vertical cut) can be understood as follows: by increasing the chemical potential $\mu$, the total density increases (not explicitly shown in figures) and the enhanced contribution of repulsive interactions can overcome the effect of the coupling $\Omega$ that favors a neutral state.

A complete phase diagram as a function of $\Omega/t$ and $t/U$ at unit filling obtained from the previous type of calculation is presented in Fig. 5a). The color plot gives the polarization $|\tilde{n}|$, the black squares form a transition line between the polarized and unpolarized states, and gray circles show the Mott insulator (left part) to superfluid (right part) transition line. The two transition lines, marking the Mott insulator to superfluid transition and, neutral to polarized transition, coincide in the intermediate region. The Mott region extends up to $t/U\approx 0.06$ for a fully polarized Mott to polarized superfluid transition (weak $\Omega/t$), which is very close to the result for the Mott-lobe tip in a single-component case. At strong $\Omega/t$ we have a neutral Mott to neutral superfluid transition that takes place at a higher value of $t/U\approx 0.074$, as strong $U_{ab}$ plays a more pronounced role in this case. In the horizontal cut shown in Fig. 5b) we explicitly show that at weak $\Omega/t$ there are two second-order phase transitions. These transitions merge into a single transition point with a jump in polarization $|\tilde{n}|$ at intermediate $\Omega/t$, as shown in Fig. 5c). At even stronger $\Omega/t$, we find two separate transitions of the second order again.

In order to explain the neutral-polarized transition on the Mott side, we complement numerical BDMFT results with an insight obtained from an effective spin Hamiltonian, which is valid in the limit of strong interactions. The spin model is derived via second order perturbation theory in the hopping amplitude and for unit filling. Starting from model (1), it is expressed in the pseudo-spin basis $|\uparrow\rangle=|n_{a}=1,n_{b}=0\rangle$ and $|\downarrow\rangle=|n_{a}=0,n_{b}=1\rangle$ Kuklov and Svistunov (2003); Altman et al. (2003); Petrescu and Le Hur (2013). Equivalently, one can also use the Schrieffer-Wolff transformation Schrieffer and Wolff (1966) to obtain the same effective Hamiltonian:

[TABLE]

where we introduce \hat{S}_{i}^{l}=(1/2)(\hat{a}_{i}^{\dagger},\hat{b}_{i}^{\dagger})\sigma^{l}\left(\begin{array}[]{c}\hat{a}_{i}\\ \hat{b}_{i}\end{array}\right) using the Pauli matrices $\sigma^{l}$ with $l=x,y,z$. For the parameters considered in this section ($t_{a}/U=t_{b}/U=t/U$), the spin coupling constants simplify to $J_{zz}=4t^{2}\left(2\gamma-1\right)/\gamma U$, $J_{\perp}=4t^{2}/\gamma U$ and $J_{z}=0.$ Based on the model (5), in the strongly interacting limit of the Bose-Hubbard model with coherent coupling from Eq. (1), we expect to find different phases depending on the magnitude of the coefficients $J_{z}$, $J_{zz}$, $J_{\perp}$ and $\Omega$. In particular, the spin ordering along $z$-direction is equivalent to the finite polarization order parameter from Eq. (2), while spin alignment along $x$-direction corresponds to the neutral phase. At the transition line of the neutral and polarized phase we expect the spin couplings $J_{\perp}$ and $\Omega$, that favor spin alignment in the $x$ direction, to be comparable to the $z$-ordering term $J_{zz}$. This reasoning leads to an approximate condition for the transition line

[TABLE]

We fit the numerical data according to this argumentation to $\Omega_{c}/t\propto\alpha t/U$, where $\alpha$ is the fitting parameter Zhan et al. (2014). For small $t/U$, we find $\Omega_{c}/t\approx 3.2t/U$ for $\gamma=8/5$ (see Fig. 5) and $\Omega_{c}/t\approx 2.15t/U$ for $\gamma=4/3$. These fitting constants explicitly fulfill the $(1-1/\gamma)$ dependence in Eq. (6).

V Imbalanced hopping amplitudes

Up to now we considered two fully equivalent bosonic components described by $t_{a}=t_{b}$ and $U_{a}=U_{b}$. As a consequence, two degenerate solutions with $\pm\tilde{n}$ were found in the polarized regime. In this section we address a more general case of imbalanced hopping amplitudes $t_{a}\neq t_{b}$. In particular, we investigate how this imbalance affects the neutral to polarized transition within the Mott regime at unit filling. We consider both positive and negative tunneling amplitudes. The latter case is less common, but it is experimentally accessible in shaken optical lattices Eckardt et al. (2010).

In Fig. 6a) we present the absolute value of the polarization order parameter $|\tilde{n}|$ as a function of tunneling rates $t_{a}$ and $t_{b}$. Other parameters of the Hamiltonian (1) are set to $\gamma=U_{ab}/U=8/5$, $\Omega/U=5\times 10^{-3}$ and $U_{a}=U_{b}=U$. We limit ourselves to small absolute values of the ratios $|t_{a}|/U$ and $|t_{b}|/U$ and unit filling $\langle n_{a}+n_{b}\rangle=1$, so that the system is in the Mott phase. Low values of the polarization are found for very weak tunneling amplitudes, in the region given by $|t_{a,b}|<t_{c}$, where we have $t_{c}/U\approx 2.5\times 10^{-3}$ for the data shown in Fig. 6. Along the diagonals $t_{a}=\pm t_{b}$, the neutral phase with $\tilde{n}=0$ extends up to the largest values of $t_{a}$ and $t_{b}$. We notice a clear difference in the extension of the neutral phase for the case of $t_{a}=t_{b}$ in comparison with the case of $t_{a}=-t_{b}$.

In order to explain the features observed in Fig. 6a), we use the effective spin model defined in Eq. (5), that captures low energy properties of the Hamiltonian (1) deep in the Mott domain. For the parameters considered in this section, the coefficients of the spin model (5) simplify to

[TABLE]

where we typically consider $1<\gamma<2$. According to Eq. (5), the coefficient $J_{z}$ plays the role of an effective chemical potential that selects which of the two species is preferred in the polarized phase. We notice that $J_{zz}$ increases quadratically both with $t_{a}$ and $t_{b}$, which leads to the largest increment in this spin coupling term isotropically around $t_{a}=t_{b}=0$. In contrast, $J_{z}$ depends on the difference of the square of both hopping amplitudes and exhibits the strongest increase perpendicular to the diagonals $t_{a}=\pm t_{b}$. In our analysis, the value of $\Omega$ is kept constant. Thus, it is the dominant quantity in the region around $t_{a,b}\approx 0$, where all other spin couplings (7)-(9) are weak and where we find the neutral phase accordingly. The asymmetry between the negative and positive side of $t_{a}$ in the plot shown in Fig. 6a) arises due to the $J_{\perp}$ coupling. For hopping amplitudes of the same sign, the $J_{\perp}$ coupling is positive. As such, it lowers the energy of the neutral phase and thus shifts the phase transition to higher values of $t_{a,b}$. The opposite is true for hopping amplitudes of different sign.

To sum up the implications of the spin model and compare these to our numerical results, we examine the interplay of spin coupling amplitudes with the help of an auxiliary function $\chi$. All spin couplings that favor the polarized phase are marked with a positive sign, while the spin couplings favoring the neutral phase carry a negative sign

[TABLE]

The roots of this function give an estimate for the neutral to polarized phase transition line and the resulting plot in Fig. 6b) qualitatively recovers the structure of the numerical phase diagram.

We now investigate how the transitions from Fig. 6 are affected by the change in the interaction ratio $\gamma$ and in the coherent coupling $\Omega$. It turns out that the corresponding $t_{a}-t_{b}$ plots look qualitatively similar to the plot Fig. 6a). In order to make a quantitative comparison, we plot the absolute value of the polarization $|\tilde{n}|$ as a function of $t_{a}$ for $t_{a}=t_{b}$ (Fig. 7, left column) and for $t_{b}=-t_{a}$ (Fig. 7, right column). In the plots shown in Fig. 7a) and b), we set $\gamma=8/5$ and vary $\Omega$. Our results show good agreement with the expectation $t_{c}/U\sim\sqrt{\Omega/U}$ from Eq. (6) for the transition point $t_{c}$. For the two cases considered ($t_{a}=t_{b}$ and $t_{a}=-t_{b}$) only the proportionality constants are different.

In the plots presented in Fig. 7c) and d) we keep the value $\Omega/U=2\times 10^{-3}$ fixed and change $\gamma$. For $t_{a}=t_{b}$ our results are well approximated by $t_{c}\sim\left(1-1/\gamma\right)^{-1/2}$, in agreement with Eq. (6). In contrast, for $t_{a}=-t_{b}$ we find that our results do not depend on the value of $\gamma$. This can be explained by looking at the spin coupling constants from Eqs. (7)-(10). As mentioned earlier, the $J_{\perp}$ coupling opposes the impact of the $\Omega$ term for hopping amplitudes of different sign. The relevant contribution to the auxiliary function $\chi$ (defined in Eq. (10)) that leads to the aforementioned $\gamma$ invariance, however, is given by the sum $J_{zz}+|J_{\perp}|$. The $\gamma$ parts of both terms cancel, leaving $\chi$ independent of $\gamma$.

VI Summary and Discussion

In this work, we focused on a two-component mixture of coherently coupled bosons in a square optical lattice. We analyzed the phase transition between the polarized (finite $\tilde{n}_{i}$, Eq. (2)) and the neutral phase, driven by an interplay of the coherent coupling $\Omega$ and the interspecies repulsion $U_{ab}$ at unit filling $\langle n_{ia}+n_{ib}\rangle=1$.

By comparing Gutzwiller results with more demanding BDMFT calculations, we found that the former provide a reasonable description of the system in the superfluid regime. We investigated the energy landscape of the mean-field Hamiltonian as a function of the two condensate order parameters and established that the coherent coupling leads to a second order phase transition between the polarized and the neutral phase. Furthermore, we found that the neutral phase is suppressed as the ratio of inter- and intra-species interactions $\gamma$ increases.

On the Mott side of the phase diagram, where BDMFT calculations provide a necessary extension of the simpler mean-field theory, we found the polarized phase to be favored only at very low values of the coherent coupling (Fig. 5a). From this, we concluded that the long-range order of the condensate seems to favour the polarized phase. To better understand our numerical results in the Mott phase, we derived an effective spin model using second order perturbation theory in the hopping amplitudes. From the coupling constants of the effective model we inferred that the polarized to neutral transition line is approximately given by $\frac{\Omega^{c}}{t}\propto\left(1-1/\gamma\right)\frac{t}{U}$ in very good agreement with our numerical results (Fig. 5). The dominance of the neutral phase in the deep Mott regime at unit filling agrees well with the other findings in this paper, as well as in the literature, especially with recent DMRG calculations Zhan et al. (2014). Furthermore, our BDMFT results indicate three possible transitions with increasing $t/U$, depending on the value of coherent coupling $\Omega$ and interaction ratio $\gamma=U_{ab}/U$: polarized Mott states turn into a polarized superfluid, the polarized Mott phase turns directly into a neutral superfluid, and from the neutral Mott phase there is a transition to a neutral superfluid. The tip of the Mott lobe is positioned at the smallest value of $t/U$ in the first case, while the lobe extends up to the largest value of the tunneling amplitude for the third case.

Finally, we explored the effects of imbalanced hopping amplitudes for the two species for both (common) positive and negative tunneling terms. The latter were realized only recently in periodically driven optical lattices Eckardt et al. (2010). For now, we considered the case of $\gamma>1$ and found that $\Omega$, enforcing the neutral phase, has a strong influence only for very small hopping amplitudes or for $t_{a}=\pm t_{b}$. An interesting asymmetry, that shows up in the two latter cases, was traced back to the sign change of one of the coupling constants in the effective spin model. In future work, we plan to consider the case of $\gamma=1/2$ and $t_{a}=-t_{b}$, where recent calculations Graß et al. (2014) suggested an occurrence of the $xy$-antiferromagnetic phase.

VII Acknowledgments

Support by the Deutsche Forschungsgemeinschaft via DFG SFB/TR 49, DFG FOR 801, DFG SPP 1929 GiRyd and the high-performance computing center LOEWE-CSC is gratefully acknowledged. This work was supported in part by DAAD (German Academic and Exchange Service) under project BKMH. I. V. acknowledges support by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under projects ON171017 and BKMH, and by the European Commission under H2020 project VI-SEEM, Grant No. 675121. Numerical simulations were partly run on the PARADOX supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade. U. B. acknowledges support by the CQT PhD program. The Centre for Quantum Technologies (CQT) is a Research Centre of Excellence founded by the Ministry of Education, Singapore and the Singapore National Research Foundation.