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

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.

Introduction

The non-equilibrium behaviour of Bose-Hubbard systems has received considerable theoretical attention recently[1-8]. 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[9-11], strong interactions are more accessible with attractive interaction energies than with repulsive interaction energies[12,13]. 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 non-equilibrium steady state (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[13-15].

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 (ħ = 1):where J denotes the hopping amplitude between the lth and (l + 1)th site, U denotes the boson interaction energy on the lth site, γ is the local dissipation rate, Ω denotes the drive amplitude (assumed real), and Δ = ω − ω, which plays the role of a chemical potential, is the site-dependent drive detuning when ω is the bare frequency of the lth site and ω is the drive frequency.

The NESS of the DDBHM, denoted ρ∞, is defined as the fixed point of the evolution given by equation (1), . We observe that the EOM for ρ given by equation (1) is the same as the EOM for ρ* if the Hamiltonian is negated (H → −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 ρ → ρ* 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. This is discussed and demonstrated further in[16], an independent work whose findings overlap with the ones presented here.

For the DDBHM there is a further simplification of the symmetry. The transformation H → −H entails Ω → −Ω, which is equivalent to b → −b. However, b → −b itself does not change the boson number statistics. To see this, note that if b → −b, then

We therefore conclude that Ω → −Ω is unnecessary to preserve the boson number statistics in the NESS; the invariance only requires U → −U, J → −J, and Δ → −Δ. We note that a dimer system is a special case, where the sign of J also need not change (see ref.[16]).

Numerical Simulation

The numerical simulation is performed by employing a matrix product density operator (MPDO) representation of ρ[17,18], which amounts to a quantum mechanical treatment characterized by a refinement parameter χ 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 χ “singular values” that represents these correlations.

In the MPDO picture the system density matrix ρ becomes a vector, denoted |ρ〉, and the superoperator becomes a regular operator such that at the NESS. To obtain an approximation for ρ∞ under a given set of system parameters U, J, Δ, Ω, and γ, we first use the hybrid evolution method of Ref.[19] to evolve the MPDO representation of a random initial state ρ under a desired choice of parameters until convergence in achieved. We then sweep the value of Ω in increments, converging the MPDO with real time evolution at each increment. Convergence is considered complete when and the singular values between the first two sites of the MPDO are converged on a logarithmic scale. We find that χ = 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 Ω in both directions. We truncate the Hilbert space on each site at four quanta, and always choose γ = 1.

Results

Uniform trimer

To test the arguments set forth above, we perform numerical investigations on a DDBHM trimer system (see ref.[16] for results with a uniform dimer 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 Δ 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 Δ are changed. The simulation parameters are summarized in Table 1. Both cases are examined at several different values of the drive strength Ω. At each value of Ω we denote the NESS for upper and lower sign choices by and , respectively.

Table 1

Simulation parameters for the DDBHM trimer.

Case	(a) Uniform			(b) Disordered
	J	U	Δ	J1,2, J2,3	U1, U2, U3	Δ1, Δ2, Δ3
1	±1	±10	±1	±1, ∓3	±8,0, ∓10	0, ∓10, ±1
2	1	±10	±1	1, −3	±8,0, ∓10	0, ∓10, ±1

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) is kept fixed while the sign of the interaction strength (U) and detuning (Δ) are changed. γ = 1 for all cases. (a) Uniform trimer, (b) Disordered trimer with non-uniform parameters.

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 Ω the local and non-local observables n1 and n1n2n3 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 n1 is the same.

Figure 1

Boson number statistics for site 1 and three-site correlators under a number-conserving transformation in a uniform trimer. Gray bars correspond to the lower sign choice in the parameters listed in Table 1(a), Case 1; and orange bars correspond to the upper sign choice. (a) Boson number statistics as a function of drive strength Ω. (b) Three-site correlator as a function of drive strength. Although the NESS changes between the two different sign choices at each drive strength, the number statistics and the correlator remain the same.

Next we examine Case 2 where the sign of J is kept fixed while the sign of U and Δ 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 n1 for both and . We find that 〈n1〉 differs between and at each value of Ω, and therefore conclude that the number statistics is not invariant under only U → −U and Δ → −Δ. This case is similar to the interaction sign change in the equilibrium BHM, where the hopping energy remains fixed, and the equilibrium phase changes.

Figure 2

Boson number expectation value on site 1 as function of varying drive strength (Ω) for a uniform trimer. Gray bars correspond to the lower sign choice in the parameters listed in Table 1(a), Case 2; and orange bars correspond to the upper sign choice. The transformation from the lower sign choice to the upper sign choice is not number-conserving and the boson number expectation value is not invariant.

Disordered trimer

To further demonstrate that the invariance is very general, we now test the boson number symmetry in the presence of strong disorder. This is important for experimental tests, where at least some disorder is inevitable, and complements the numerical tests in ref.[16], which considers only uniform systems. As before, and 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 n1 and n1n2n3 are the same between and at any given Ω. We see in fact that the entire statistical distribution of n1 is the same as in the uniform trimer case.

Figure 3

Boson number statistics for site 1 and three-site correlators under a number-conserving transformation in a disordered trimer. Gray bars correspond to the lower sign choice in the parameters listed in Table 1(b), Case 1; and orange bars correspond to the upper sign choice. (a) Boson number statistics as a function of drive strength Ω. (b) Three-site correlator as a function of drive strength. Although the NESS changes between the two different sign choices at each drive strength, the number statistics and the correlator remain the same even in the presence of strong disorder.

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 〈n1〉 is different between and at each value of Ω.

Figure 4

Boson number expectation value on site 1 as function of varying drive strength (Ω) for a disordered trimer. Gray bars correspond to the lower sign choice in the parameters listed in Table 1(b), Case 2; and orange bars correspond to the upper sign choice. The transformation from the lower sign choice to the upper sign choice is not number-conserving and the boson number expectation value is not invariant.

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 the EOM, and the dynamical observables should also remain invariant. This is numerically demonstrated in the independent work of ref.[16] for the case of a DDBHM dimer as well as a spin system.

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, strong attractive interactions are more readily accessible[12,13] 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[3], fermionized photons[1], polariton crystalization[2], photon transport resonances[5], first-order dissipative quantum phase transitions[6], and diffusive-insulator transport phase transitions[8].

Our findings confirm and complement the results in the recent independent work by Li and Koch[16].

Methods

Effective Equation of Motion

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

ω denotes the bare frequency of the lth site, J denotes the hopping amplitude between the lth and (l + 1)th site, U denotes the boson interaction energy on the lth site, Ω denotes the drive amplitude, and ω denotes the drive frequency. With on-site dissipation to a Markovian bath, the density matrix ρ of the chain is governed by the following Lindblad master equation: , where and γ 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†, where . The resulting effective EOM is (setting Ω real)where and Δ = ω − ω is the site-dependent drive detuning, which plays the role of a chemical potential. For simplicity, in the main text we write as ρ and as H.

Numerical Simulations

The numerical simulation is performed by employing a matrix product density operator (MPDO) representation of ρ[17,18], which amounts to a quantum mechanical treatment characterized by a refinement parameter χ 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 χ “singular values” that represents these correlations.

In the MPDO picture the system density matrix ρ becomes a vector, denoted |ρ〉, and the superoperator becomes a regular operator such that at the NESS. To obtain an approximation for ρ∞ under a given set of system parameters U, J, Δ, Ω, and γ, we first use the hybrid evolution method of ref.[19] to evolve the MPDO representation of a random initial state ρ under a desired choice of parameters until convergence in achieved. We then sweep the value of Ω in increments, converging the MPDO with real time evolution at each increment. Convergence is considered complete when and the singular values between the first two sites of the MPDO are converged on a logarithmic scale. We find that χ = 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 Ω in both directions. We truncate the Hilbert space on each site at four quanta, and always choose γ = 1.