Symmetry between repulsive and attractive interactions in driven-dissipative Bose-Hubbard systems

TL;DR

This paper reveals a symmetry in driven-dissipative Bose-Hubbard systems, showing that their steady states and phenomena like fermionization and superbunching are invariant under simultaneous sign changes of key parameters, enabling equivalent behaviors with attractive or repulsive interactions.

Contribution

The study uncovers a fundamental symmetry in the Lindblad master equation for driven-dissipative Bose-Hubbard models, linking attractive and repulsive interaction regimes.

Findings

Steady states are invariant under simultaneous sign change of parameters.

Local boson number distribution and correlations are unchanged.

Phenomena like fermionization and superbunching occur in both interaction regimes.

Abstract

The driven-dissipative Bose-Hubbard model can be experimentally realized with either negative or positive onsite detunings, inter-site hopping energies, and onsite interaction energies. Here we use one-dimensional matrix product density operators to perform a fully quantum investigation of the dependence of the non-equilibrium steady states of this model on the signs of these parameters. Due to a symmetry in the Lindblad master equation, we find that simultaneously changing the sign of the interaction energies, hopping energies, and chemical potentials leaves the local boson number distribution and inter-site number correlations invariant, and the steady-state complex conjugated. This shows that all driven-dissipative phenomena of interacting bosons described by the Lindblad master equation, such as "fermionization" and "superbunching", can equivalently occur with attractive or repulsive interactions.

Tables

Table 1. Table 1: Simulation parameters for the DDBHM trimer. Case 1 corresponds to a number-conserving transformation in which the signs of all the parameters are flipped simultaneously. Case 2 corresponds to a transformation in which the sign of the hopping energy ( J 𝐽 J ) is kept fixed while the sign of the interaction strength ( U 𝑈 U ) and detuning ( Δ Δ \Delta ) are changed. γ = 1 𝛾 1 \gamma=1 for all cases. (a) Uniform trimer, (b) Disordered trimer with non-uniform parameters.

	(a) Uniform			(b) Disordered
Case	$J$	$U$	$\Delta$	$J_{1,2}$, $J_{2,3}$	$U_1$, $U_2$, $U_3$	${\Delta }_1$, ${\Delta }_2$, ${\Delta }_3$
1	$\pm 1$	$\pm 10$	$\pm 1$	$\pm 1,\mp 3$	$\pm 8,0,\mp 10$	$0,\mp 10,\pm 1$
2	$1$	$\pm 10$	$\pm 1$	$1,-3$	$\pm 8,0,\mp 10$	$0,\mp 10,\pm 1$

Equations

$$\frac{d}{dt}\rho$$

$$H$$

$$|n⟩⟨n|=(b^{†}{)}^n|0⟩⟨0|(b{)}^n\to (-b^{†}{)}^n|0⟩⟨0|(-b{)}^n=|n⟩⟨n|.$$

$$H=$$

$$+\sum _l\frac{U_l}{2}{b}_l^{†}{b}_l^{†}{b}_l{b}_l+\sum _l(\Omega b_l^{†}{e}^{-i{\omega }_{d}t}+{\Omega }^{\ast }{b}_l{e}^{+i{\omega }_{d}t}).$$

$$\frac{d}{dt}\tilde{\rho }=$$

$$=$$

$$\tilde{H}=$$

$$+\sum _l\frac{U_l}{2}{b}_l^{†}{b}_l^{†}{b}_l{b}_l+\Omega \sum _l(b_l^{†}+b_l),$$

Full text

Introduction

The non-equilibrium behaviour of Bose-Hubbard systems has received considerable theoretical attention recently Carusotto et al. (2009); Hartmann (2010); Grujic et al. (2013); Le Boité et al. (2014); Biella et al. (2015); Weimer (2015); Foss-Feig et al. (2016); Debnath et al. (2017). However, to our knowledge the dependence of the non-equilibrium physics of the Bose-Hubbard model (BHM) on the signs of the hopping and interaction energies has yet to be explored. In superconducting circuits, which are a natural setting for studying the non-equilibrium physics of driven-dissipative many-body systems Houck et al. (2012); Schmidt and Koch (2013); Hur et al. (2016), strong interactions are more accessible with attractive interaction energies than with repulsive interaction energies Bourassa et al. (2012); Hacohen-Gourgy et al. (2015). On the other hand, theoretical studies of the driven-dissipative BHM (DDBHM) have focused on the case of repulsive interactions. Finding a theoretical link between the attractive and repulsive interaction regimes of the DDBHM would therefore be of practical experimental benefit.

In this work we point out a symmetry in the Lindbladian equation of motion for the DDBHM that implies that the driven-dissipative physics of repulsive interactions can be replicated with attractive interactions, irrespective of the magnitude of the interaction strength. To illustrate this, we employ a fully quantum (i.e. non-mean-field) numerical treatment of a DDBHM trimer. We show that simultaneously changing the signs of the interaction, hopping, and detuning while keeping their magnitudes fixed changes the NESS but does not change the three-site number correlator nor the statistics of the on-site boson number. We also demonstrate that this observable symmetry persists even in the presence of strong disorder in all of the sign-flipped parameters. This symmetry can be experimentally tested with existing superconducting circuit technology, which has the potential to realize the BHM such that the chemical potential, on-site interaction energy, and inter-site hopping energy are all tunable in situ (within a limited range) in both magnitude and sign Hacohen-Gourgy et al. (2015); Deng et al. (2015, 2016).

Model

We investigate the open boundary dissipative Bose-Hubbard chain under homogeneous coherent driving in a frame rotating at the drive frequency. With on-site dissipation to a Markovian bath, the effective equation of motion (EOM) (see Appendix for derivation) is given by the following Lindblad master equation ($\hbar=1$):

[TABLE]

where $J_{l,l+1}$ denotes the hopping amplitude between the $l$th and ($l+1$)th site, $U_{l}$ denotes the boson interaction energy on the $l$th site, $\gamma$ is the local dissipation rate, $\Omega$ denotes the drive amplitude (assumed real), and $\Delta_{l}=\omega_{l}-\omega_{d}$, which plays the role of a chemical potential, is the site-dependent drive detuning when $\omega_{l}$ is the bare frequency of the $l$th site and $\omega_{d}$ is the drive frequency.

The NESS of the DDBHM, denoted ${\rho}_{\infty}$, is defined as the fixed point of the evolution given by equation (1), $\frac{d}{dt}{\rho}_{\infty}=0$. We observe that the EOM for ${\rho}$ given by equation (1) is the same as the EOM for ${\rho}^{*}$ if the Hamiltonian is negated (${H}\rightarrow-{H}$). Therefore the NESS attained by evolving with ${H}$ is equal to the complex conjugate of the NESS attained by evolving with $-{H}$. However, the transformation ${\rho}\rightarrow{\rho}^{*}$ does not change the observable statistics of the state. The observables of the NESS are therefore invariant under negation of the Hamiltonian. We note that this symmetry applies not just to the DDBHM, but to any model described by the Lindblad master equation where the dissipation operators are invariant under complex conjugation.

For the DDBHM there is a further simplification of the symmetry. The transformation ${H}\rightarrow-{H}$ entails $\Omega\rightarrow-\Omega$, which is equivalent to $b_{l}\rightarrow-b_{l}$. However, $b_{l}\rightarrow-b_{l}$ itself does not change the boson number statistics. To see this, note that if $b\rightarrow-b$, then

[TABLE]

We therefore conclude that $\Omega\rightarrow-\Omega$ is unnecessary to preserve the boson number statistics in the NESS; the invariance only requires $U_{l}\rightarrow-U_{l}$, $J_{l,l+1}\rightarrow-J_{l,l+1}$, and $\Delta_{l}\rightarrow-\Delta_{l}$.

Numerical Simulation

The numerical simulation is performed by employing a matrix product density operator (MPDO) representation of ${\rho}$ Zwolak and Vidal (2004); Verstraete et al. (2004), which amounts to a quantum mechanical treatment characterized by a refinement parameter $\chi$ that designates the maximum size of the tensors that represent each site, and therefore the maximum amount of total correlations (classical plus quantum) between bipartitions of the chain that can be captured by the MPDO. Linking each site tensor with its neighbor in the MPDO is a diagonal matrix of $\chi$ “singular values” that represents these correlations.

In the MPDO picture the system density matrix ${\rho}$ becomes a vector, denoted $|{\rho}\rangle$, and the superoperator $\mathcal{{L}}$ becomes a regular operator $\mathcal{{L}}_{\sharp}$ such that $\langle{\rho}|\mathcal{{L}}_{\sharp}|{\rho}\rangle=0$ at the NESS. To obtain an approximation for ${\rho}_{\infty}$ under a given set of system parameters $U_{l}$, $J_{l}$, $\Delta_{l}$, $\Omega$, and $\gamma$, we first use the hybrid evolution method of Ref. Gangat et al. (2016) to evolve the MPDO representation of a random initial state ${\rho}$ under a desired choice of parameters until convergence in achieved. We then sweep the value of $\Omega$ in increments, converging the MPDO with real time evolution at each increment. Convergence is considered complete when $\langle\mathcal{{L}}_{\sharp}\rangle\lesssim 10^{-3}$ and the singular values between the first two sites of the MPDO are converged on a logarithmic scale. We find that $\chi=15$ and a timestep size of $10^{-1}$ is sufficient to achieve this for all of the cases that we consider. We verify uniqueness of the NESS by performing the sweep of $\Omega$ in both directions. We truncate the Hilbert space on each site at four quanta, and always choose $\gamma=1$.

Results

.1 Uniform trimer

To test the arguments set forth above, we perform numerical investigations on a DDBHM trimer system. We first test the boson number symmetry when the parameters are uniform across the trimer. We specifically look at two cases: Case 1 examines the change in the NESS under the number-conserving transformation argued above (the hopping energy $J$, the interaction strength $U$, and the detuning $\Delta$ all change signs simultaneously); Case 2 examines the change in the NESS under a transformation that is different from the number-conserving transformation discussed in the previous section: the sign of $J$ is kept fixed while the sign of $U$ and $\Delta$ are changed. The simulation parameters are summarized in Table 1. Both cases are examined at several different values of the drive strength $\Omega$. At each value of $\Omega$ we denote the NESS for upper and lower sign choices by $|{\rho}^{(+)}_{\infty}\rangle$ and $|{\rho}^{(-)}_{\infty}\rangle$, respectively.

First we examine the parameter sets for Case 1 in Table 1(a). In accordance with the boson number symmetry argued earlier, here we find that at every value of $\Omega$ the local and non-local observables $n_{1}$ and $n_{1}n_{2}n_{3}$ are invariant in the NESS under the collective sign change, as shown in Fig. 1. More precisely, we see that the full statistical distribution of $n_{1}$ is the same.

Next we examine Case 2 where the sign of $J$ is kept fixed while the sign of $U$ and $\Delta$ are changed. This is not a number-conserving transformation and we do not expect the observables will remain the same after the transform. In Fig. 2, we plot the expectation value of the observable $n_{1}$ for both $|{\rho}^{(+)}_{\infty}\rangle$ and $|{\rho}^{(-)}_{\infty}\rangle$. We find that $\langle n_{1}\rangle$ differs between $|{\rho}^{(+)}_{\infty}\rangle$ and $|{\rho}^{(-)}_{\infty}\rangle$ at each value of $\Omega$, and therefore conclude that the number statistics is not invariant under only $U_{l}\rightarrow-U_{l}$ and $\Delta_{l}\rightarrow-\Delta_{l}$. This case is similar to the interaction sign change in the equilibrium BHM , where the hopping energy remains fixed, and the equilibrium phase changes.

.2 Disordered trimer

To further demonstrate that the invariance is very general, we now test the boson number symmetry in the presence of strong disorder. As before, $|{\rho}^{(+)}_{\infty}\rangle$ and $|{\rho}^{(-)}_{\infty}\rangle$ respectively denote the upper and lower sign choices of the parameters. We consider two specific cases analogous to those for the uniform trimer. In Case 1, the change in the NESS is examined when the hopping energy, interaction strength, and detuning all change sign; In Case 2, the sign of the hopping energy is kept fixed while the signs of the interaction strength and detuning are changed.

We first examine parameters for Case 1 as listed in Table 1(b). In this case the upper sign choice and lower sign choice of the parameters are related by the boson number symmetry transformation. Consequently, Fig.  3 reveals that the local and non-local observables $n_{1}$ and $n_{1}n_{2}n_{3}$ are the same between $|{\rho}^{(+)}_{\infty}\rangle$ and $|{\rho}^{(-)}_{\infty}\rangle$ at any given $\Omega$. We see in fact that the entire statistical distribution of $n_{1}$ is the same as in the uniform trimer case.

On the other hand, the parameter transformation in Case 2 is not of the type with boson number symmetry discussed earlier. Consequently, Fig. (4) shows that $\langle n_{1}\rangle$ is different between $|{\rho}^{(+)}_{\infty}\rangle$ and $|{\rho}^{(-)}_{\infty}\rangle$ at each value of $\Omega$.

Finally, we note that although we only consider the observables in the NESS, the invariance under the number-conserving transformation is at the level of EOM, and the dynamical observables should also remain invariant.

Discussion

We have given an analytical argument and provided numerical evidence for a boson number symmetry of the DDBHM. Specifically, the symmetry is that the boson number statistics of the system state are invariant to collective changes in the sign of the interaction energies, detunings, and hopping energies. In other words, simultaneously changing the sign of all of the parameters of the number-conserving terms of the system Hamiltonian does not observably change the state. On the other hand, we have also numerically shown that keeping the sign of the hopping energy fixed while changing the signs of the detunings and interaction energies does not leave the number statistics invariant.

We have therefore shown two contrasts to the case of equilibrium phases of the BHM: 1) the number statistics of the NESS of the DDBHM can exhibit a strong dependence on the sign of the hopping energy, and 2) it is possible for the number statistics of the NESS to be exactly the same for opposite signs of the interaction energy with the same magnitude.

These theoretical predictions are experimentally testable with existing superconducting circuit technology, and the symmetry is applicable beyond the DDBHM to any situation where the Lindbladian jump operators are real.

For self-interactions of bosons in superconducting circuits, it is strong attractive interactions that are experimentally accessible Bourassa et al. (2012); Hacohen-Gourgy et al. (2015) rather than strong repulsive interactions. Therefore, the equivalence between attractive and repulsive interactions that we have shown here for driven-dissipative bosonic phenomena indicates that superconducting circuits with strong attractive interactions are a viable platform for investigating predictions made for driven-disspative bosonic phenomena involving strong repulsive interactions, such as repulsively induced photon superbunching Grujic et al. (2013), fermionized photons Carusotto et al. (2009), polariton crystalization Hartmann (2010), photon transport resonances Biella et al. (2015), first-order dissipative quantum phase transitions Weimer (2015), and diffusive-insulator transport phase transitions Debnath et al. (2017).

Note added: After our paper was written, we noticed a preprint by Li and Koch Li and Koch (2017) reaching similar conclusions.

Methods

.3 Effective Equation of Motion

The open boundary Bose-Hubbard chain under homogeneous coherent driving is described by the following Hamiltonian ($\hbar=1$):

[TABLE]

$\omega_{l}$ denotes the bare frequency of the $l$th site, $J_{l,l+1}$ denotes the hopping amplitude between the $l$th and ($l+1$)th site, $U_{l}$ denotes the boson interaction energy on the $l$th site, $\Omega$ denotes the drive amplitude, and $\omega_{d}$ denotes the drive frequency. With on-site dissipation to a Markovian bath, the density matrix $\rho$ of the chain is governed by the following Lindblad master equation: $\frac{d}{dt}\rho=\mathcal{L}\rho=-i[H,\rho]+\gamma\sum_{l}\mathcal{D}[b_{l}]\rho$, where $\mathcal{D}[b]\rho=\frac{1}{2}(2b\rho b^{\dagger}-b^{\dagger}b\rho-\rho b^{\dagger}b)$ and $\gamma$ is the local dissipation rate. To eliminate the time dependence, the master equation is multiplied from the left by $U$ and from the right by $U^{\dagger}$, where $U=e^{i\omega_{d}t\sum_{l}b_{l}^{\dagger}b_{l}}$. The resulting effective EOM is (setting $\Omega$ real)

[TABLE]

where $\tilde{\rho}=U\rho U^{\dagger}$ and $\Delta_{l}=\omega_{l}-\omega_{d}$ is the site-dependent drive detuning, which plays the role of a chemical potential. For simplicity, in the main text we write $\tilde{\rho}$ as $\rho$ and $\tilde{H}$ as $H$.

.4 Numerical Simulations

The numerical simulation is performed by employing a matrix product density operator (MPDO) representation of ${\rho}$ Zwolak and Vidal (2004); Verstraete et al. (2004), which amounts to a quantum mechanical treatment characterized by a refinement parameter $\chi$ that designates the maximum size of the tensors that represent each site, and therefore the maximum amount of total correlations (classical plus quantum) between bipartitions of the chain that can be captured by the MPDO. Linking each site tensor with its neighbor in the MPDO is a diagonal matrix of $\chi$ “singular values” that represents these correlations.

In the MPDO picture the system density matrix ${\rho}$ becomes a vector, denoted $|{\rho}\rangle$, and the superoperator $\mathcal{{L}}$ becomes a regular operator $\mathcal{{L}}_{\sharp}$ such that $\langle{\rho}|\mathcal{{L}}_{\sharp}|{\rho}\rangle=0$ at the NESS. To obtain an approximation for ${\rho}_{\infty}$ under a given set of system parameters $U_{l}$, $J_{l}$, $\Delta_{l}$, $\Omega$, and $\gamma$, we first use the hybrid evolution method of Ref. Gangat et al. (2016) to evolve the MPDO representation of a random initial state ${\rho}$ under a desired choice of parameters until convergence in achieved. We then sweep the value of $\Omega$ in increments, converging the MPDO with real time evolution at each increment. Convergence is considered complete when $\langle\mathcal{{L}}_{\sharp}\rangle\lesssim 10^{-3}$ and the singular values between the first two sites of the MPDO are converged on a logarithmic scale. We find that $\chi=15$ and a timestep size of $10^{-1}$ is sufficient to achieve this for all of the cases that we consider. We verify uniqueness of the NESS by performing the sweep of $\Omega$ in both directions. We truncate the Hilbert space on each site at four quanta, and always choose $\gamma=1$.

Acknowledgements

This work is partially supported by Ministry of Science and Technology, Taiwan, under Grants No. MOST 104-2112-M-002 -022 -MY3, MOST 105-2112-M-002-023-MY3, MOST 106-2811-M-002-054 (A.A.G., Y.J.K.).

Author contribution statement

A.A.G proposed the project, performed the numerical simulations and wrote the manuscript; I.P.M. proposed the general proof; Y.J.K. supervised the project and wrote the manuscript. All the authors discussed the results and the manuscript.

Additional information

The authors declare that they have no competing financial interests. Correspondence and requests for materials should be addressed to Y.J.K. (email: yjkao@phys.ntu.edu.tw) .