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.


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): l l l l l l l where J l,l+1 denotes the hopping amplitude between the lth and (l + 1)th site, U l denotes the boson interaction energy on the lth site, γ is the local dissipation rate, Ω denotes the drive amplitude (assumed real), and Δ l = ω l − ω d , which plays the role of a chemical potential, is the site-dependent drive detuning when ω l is the bare frequency of the lth site and ω d 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 We therefore conclude that Ω → −Ω is unnecessary to preserve the boson number statistics in the NESS; the invariance only requires U l → −U l , J l,l+1 → −J l,l+1 , and Δ l → −Δ l . 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 ρ ρ 〈 | | 〉 = 0 #  at the NESS. To obtain an approximation for ρ ∞ under a given set of system parameters U l , J l , Δ l , Ω, 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  〈 〉 − 10 # 3  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  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 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 Δ 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 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 n 1 and n 1 n 2 n 3 are the same between ρ | 〉 ∞ + ( ) and ρ | 〉 ∞ − ( ) at any given Ω. 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 〈n 1 〉 is different between ρ | 〉 ∞ + ( ) and ρ | 〉 ∞ − ( ) at each value of Ω.
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 . ω ω + + +