Open quantum systems with local and collective incoherent processes: Efficient numerical simulations using permutational invariance

The permutational invariance of identical two-level systems allows for an exponential reduction in the computational resources required to study the Lindblad dynamics of coupled spin-boson ensembles evolving under the effect of both local and collective noise. Here we take advantage of this speedup to study several important physical phenomena in the presence of local incoherent processes, in which each degree of freedom couples to its own reservoir. Assessing the robustness of collective effects against local dissipation is paramount to predict their presence in different physical implementations. We have developed an open-source library in python, the Permutational-Invariant Quantum Solver (PIQS), which we use to study a variety of phenomena in driven-dissipative open quantum systems. We consider both local and collective incoherent processes in the weak-, strong-, and ultrastrong-coupling regimes. Using PIQS, we reproduce a series of known physical results concerning collective quantum effects and extend their study to the local driven-dissipative scenario. Our work addresses the robustness of various collective phenomena, e.g., spin squeezing, superradiance, and quantum phase transitions, against local dissipation processes.

Figure 1

Open quantum system dynamics. An ensemble of identical two-level systems is collectively coupled to a bosonic cavity through a coherent dynamics. The action of dissipative processes on the two-level system dynamics is quantified by different rates for homogeneous local (dashed box) and collective processes (dot-dashed box), given by ${\mathcal{L}}_{\mathrm{loc}}\left[\rho \right]$ and ${\mathcal{L}}_{\mathrm{col}}\left[\rho \right]$, respectively: local and collective emission, set by ${\gamma }_{\downarrow }$ and ${\gamma }_{\Downarrow }$, local and collective dephasing, set by ${\gamma }_{\varphi }$ and ${\gamma }_{\mathrm{\Phi }}$, and local and collective pumping, set by ${\gamma }_{\uparrow }$ and ${\gamma }_{\Uparrow }$, respectively.

Figure 2

Dicke space structure. (a) Block-diagonal structure of the density matrix in the Dicke basis, $|j,m\rangle \langle j,m^{\prime }|$, for $N=6$. Each block has a different value of $j$, decreasing from the top-left one for which $j=\frac{N}{2}$. The size of each block is set by the multiplicity of $m$. Red equals 1, white equals 0. The blue region outside of the diagonal blocks (darker shading) is shown to highlight sectors that cannot be populated. (b) The density matrix of the Dicke state $\left|\frac{N}{2},\frac{N}{2}\right.〉$, which is the fully excited state (top left matrix element). (c) The Dicke state $\left|\frac{N}{2},-\frac{N}{2}\right.〉$ is the ground state if the Hamiltonian $H=\hslash {\omega }_0{J}_z$ (bottom right matrix element in the first block). (d) Qualitative representation (saturated palette) of the steady state for $H=\hslash {\omega }_0{J}_z$, and ${\gamma }_{\downarrow }=0.3{\gamma }_{\uparrow }$, ${\gamma }_{\uparrow }={\omega }_0$. Due to the choice of the Hamiltonian, only the matrix elements $p_{jmm}$ on the main diagonal can be populated, and the choice of the relative strength of the dissipation determines the fact that the matrix elements $p_{jjj}$ are the most populated ones. (e) The density matrix of the maximally symmetric Dicke state $\left|\frac{N}{2},0\right.〉$ (central matrix element in the first block). (f) The density matrix of the subradiant Dicke state $\left|0,0\right.〉$. All these states are diagonal in the Dicke space, $p_{jm{m}^{\prime }}=0$ for $m\ne m^{\prime }$ (bottom right matrix element).

Figure 3

Nondiagonal density matrices. Representation of quantum states in the Dicke basis $|j,m\rangle \langle j,m^{\prime }|$. We show for $N=10$ only the block with $j=\frac{N}{2}$, which corresponds to the largest one in Fig. 2, since it is the only block with nonzero matrix elements. (a) The GHZ state (four red matrix elements at the corners of the first block), (b) the symmetric coherent spin state (CSS) with $a=b=\frac{1}{\sqrt{2}}$ in Eq. (13), and (c) the antisymmetric CSS, with $a=-b=\frac{1}{\sqrt{2}}$ in Eq. (13).

Figure 4

Dicke space for $N=8$ showing the Dicke states $|j,m\rangle$ arranged in Dicke ladders, each with a given degeneracy $d_j^N$. Sketch of the dynamics experienced by a Dicke state, expressed in terms of the rates ${\mathrm{\Gamma }}^{\left(i\right)}$, where we omit the subscripts for clarity ${\mathrm{\Gamma }}^{\left(i\right)}={\mathrm{\Gamma }}_{j,m,m^{\prime }}^{\left(i\right)}$. The effect of different driven-dissipative collective and incoherent processes is shown with colored arrows. We show the Dicke space [56] and the action of each rate onto a given Dicke state in the center of the Dicke space. All processes contribute to the rate ${\mathrm{\Gamma }}_{j,m,m^{\prime }}^{\left(1\right)}$. Local pure dephasing, local emission, and local pumping connect neighboring Dicke ladders. Collective emission and collective pumping connect vertically subsequent rungs in the same Dicke ladder, while collective pure dephasing destroys correlations among different Dicke states. In the dashed rectangle, the contribution to the time evolution of the state probability is expressed in terms of incoming probability from the other eight neighboring states according to the rates ${\mathrm{\Gamma }}^{\left(i\right)}$, with $i>1$, and population depletion at a rate ${\mathrm{\Gamma }}^{\left(1\right)}$ (dashed black arrow).

Figure 5

PIQS performance, small $N$. Running time as a function of the number $N$ of TLSs in the system, up to $N=10$. (a) Semilogarithmic plot giving the time (in seconds) required to construct the Liouvillian superoperator, with $H={\omega }_0{J}_z$ in the Dicke basis (blue circles joined by solid lines) and in the uncoupled basis (orange circles joined by dashed lines). We set ${\omega }_0={\omega }_x=1$, the local processes ${\gamma }_{\downarrow }={\gamma }_{\varphi }={\gamma }_{\uparrow }=0.1$, and the collective processes ${\gamma }_{\Downarrow }={\gamma }_{\mathrm{\Phi }}={\gamma }_{\Uparrow }=0.01$. (b) Semilogarithmic plot giving the time (in seconds) required to solve the dynamics given by the Liouvillian in the Dicke basis (blue circles joined by solid lines) and in the uncoupled basis (orange circles joined by dashed lines). In both cases we use QuTiP's master equation mesolve for $H={\omega }_0{J}_z$, where ${\omega }_0=1$, with local processes ${\gamma }_{\downarrow }={\gamma }_{\varphi }={\gamma }_{\uparrow }=0.1$, and collective processes ${\gamma }_{\Downarrow }={\gamma }_{\mathrm{\Phi }}={\gamma }_{\Uparrow }=0.01$. We use 1000 integration time steps from $t_0=0$, when the system is initially fully excited, to $t_{\text{max}}=4ln\left(N\right)/N{\gamma }_{\Downarrow }$.

Figure 6

PIQS performance, large $N$. Semilogarithmic plots of the running time as a function of the number of TLSs $N$ in the system, up to $N=100$. We compare the time (in seconds) required to build the matrix governing the dynamics and to solve it. We compare two methods that both leverage PIQS. The first one is the most flexible one, which builds the Liouvillian (blue circles joined by dashed lines) and solves the dynamics using QuTiP's mesolve (blue circles joined by solid lines). The second one, applicable in the subset of problems diagonal in the Dicke basis, uses an ad hoc permutational-invariant solver, $\mathtt{pisolve}$, which further decreases the system's size, building the matrix $M$ (red stars joined by dashed lines) and solving the relative rate equations (red stars joined by solid lines). The dynamics is set by Eq. (1) with $H=\hslash {\omega }_0{J}_z$, where $\hslash {\omega }_{\text{0}}=1$, the local processes are set by ${\gamma }_{\downarrow }={\gamma }_{\varphi }={\gamma }_{\uparrow }=0.1$, and the collective processes are set by ${\gamma }_{\Downarrow }={\gamma }_{\mathrm{\Phi }}={\gamma }_{\Uparrow }=0.01$. We use 1000 integration time steps from $t_0=0$, when the system is initially fully excited, to $t_{\text{max}}=4t_{\mathrm{D}}=4ln\left(N\right)/N{\gamma }_{\Downarrow }$.

Figure 7

Superradiant light emission. Superradiant decay of the total inversion, $\langle J_z\rangle \left(t\right)$, and emitted light, $\langle J_{+}{J}_{-}\rangle \left(t\right)$, for $N=20$ TLSs, for different initial states, with $H={\omega }_0{J}_z$, ${\omega }_0\gg {\gamma }_{\Downarrow }$, ${\gamma }_{\Downarrow }/{\gamma }_{\varphi }=1$. The time is expressed in units of the superradiant delay time, $t_{\text{D}}=lnN/N{\gamma }_{\Downarrow }$.

Figure 8

Steady-state superradiance. (a) Value of the steady-state light emission ${\langle J_{+}{J}_{-}\rangle }_{\text{ss}}/N$ as a function of the local pumping strength ${\gamma }_{\uparrow }$, which is in units of $N{\gamma }_{\Downarrow }$. We set ${\gamma }_{\Downarrow }$ and tune ${\gamma }_{\uparrow }$ for $N=10$ (thinnest curves), 20, 30, 40 (thickest curves). Two different cases are considered: in the first one, only collective emission is present, ${\gamma }_{\downarrow }=0$ (solid black curves), see also Fig. 2 of Ref. [10]; in the second one (dashed red curves), the local emission rate ensures the detailed balance condition with respect to local pumping, ${\gamma }_{\downarrow }={\gamma }_0(1+n_{\text{T}})$ and ${\gamma }_{\uparrow }={\gamma }_0{n}_{\text{T}}$; we consider the high-temperature limit, $n_{\text{T}}\gg 1$, thus setting ${\gamma }_{\downarrow }={\gamma }_{\uparrow }$. (b) We study the steady-state light emission for fixed ${\gamma }_{\Downarrow }={\omega }_0$ and $N=40$, varying both ${\gamma }_0$ and $n_{\text{T}}$.

Figure 9

Spin squeezing. (a) Time evolution of the spin squeezing parameter ${\xi }^2$ for two different dynamics and two different initial states: the solid curve corresponds to ${\gamma }_{\downarrow }=\mathrm{\Lambda }/5$ and ${\gamma }_{\Downarrow }=0$, the dashed curve to ${\gamma }_{\downarrow }=0$ and ${\gamma }_{\Downarrow }=\mathrm{\Lambda }/5$. Here we consider $N=50$ TLSs. The thin curves qualitatively reproduce Fig. 3 of Ref. [61] and correspond to an initially fully excited state, $\left|\frac{N}{2},\frac{N}{2}\right.〉$. The thick curves instead correspond to the nonsymmetric Dicke state with longest time under spin-squeezed evolution, ${\xi }^2<1$, the $\left|j,j\right.〉$ Dicke state with $j=14$. (b) Study of the minimum spin-squeezing parameter ${\xi }^2$, for any initial Dicke state $|j,m\rangle$: We consider ${\gamma }_{\downarrow }=\mathrm{\Lambda }/5$ and ${\gamma }_{\Downarrow }=0$ and $N=20$; the size of the circles and the shading of the filling give the strength of the maximum spin squeezing obtained. (c) Trade-off between the spin-squeezing time $\tau$ and the minimum spin-squeezing parameter ${\xi }^2$ for different initial states, fixing $N=20$ for two different dynamical conditions, either in the presence of collective dissipation only (circles joined by the solid orange segmented line) or local dissipation only (dashed blue segmented line). (d) The minimum value of the spin-squeezing parameter ${\xi }^2$ (maximum squeezing) is explored for different values of ${\gamma }_{\downarrow }$ and ${\gamma }_{\Downarrow }$, with initial state the fully excited state for $N=20$. (e), (f) Color plots for the minimum spin-squeezing parameter value ${\xi }^2$ and the time at which it is reached $\tau$ as a function of ${\gamma }_{\downarrow }$, fixing ${\gamma }_{\Downarrow }=\mathrm{\Lambda }/5$. The dashed lines help identify qualitatively the region with spin squeezing (parameter space to the right of the line).

Figure 10

Open Dicke model. Here we consider the Wigner function of the photonic part of the steady state; $x$ and $p$ are the conjugate operators of the photonic mode operators. We fix the cutoff of the photonic Hilbert space at $n_{\text{ph}}=20$ and study the influence of local and collective processes. For $N=10$ TLSs, $H$ is given by Eq. (25), with ${\omega }_{\text{cav}}={\omega }_0$, $g=2{\omega }_0/\sqrt{N}$, $\kappa ={\omega }_0$, and ${\gamma }_{\varphi }=0.01{\omega }_0$. (a) Adding ${\gamma }_{\uparrow }=0.1{\omega }_0$ does not restore the superradiant phase transition. (b) Adding ${\gamma }_{\downarrow }=0.1{\omega }_0$ restores the superradiant phase transition, and reproduces one panel of Fig. 1 of Ref. [78]. (c) Adding both ${\gamma }_{\downarrow }=0.1{\omega }_0$ and the collective pumping ${\gamma }_{\Uparrow }=0.1{\omega }_0$ is detrimental to the superradiant phase. (d) Adding to ${\gamma }_{\downarrow }=0.1{\omega }_0$, also the collective emission ${\gamma }_{\Downarrow }=0.1{\omega }_0$ still allows us to resolve the superradiant phase.

Figure 11

Boundary time crystal. The dynamics of Eq. (28) is shown in the case of weak dissipation (${\omega }_x\gg \frac{N}{2}{\gamma }_{\Downarrow }$) for $N=30$, with the system initialized in the excited state $\left|\frac{N}{2},\frac{N}{2}\right.〉$, for different dephasing rates ${\gamma }_{\varphi }/{\omega }_x=0,0.01,0.1,1$ (thicker to thinner curves). (a) The total inversion, $\langle J_z\rangle \left(t\right)$, normalized by $N/2$, is plotted in time. These results extend those obtained for ${\gamma }_{\varphi }=0$ in Ref. [85]. (b) The total spin length, $\langle J^2\rangle \left(t\right)$, normalized by $j_2=\frac{N}{2}(\frac{N}{2}+1)$, is shown. In panels (c) and (d), the same quantities are plotted for different collective dephasing rates, ${\gamma }_{\mathrm{\Phi }}/{\omega }_x=0,0.01,0.1,1$ (thicker to thinner curves).

Figure 12

Spin-excitation exchange in multiple TLS ensembles. The robustness of the spin-excitation exchange given by the dynamics of Eq. (29) is studied for two dissipatively coupled ensembles of TLSs. Black and red curves display the time evolution of the total spin inversion of the first and second ensembles, respectively, normalized by $N_1/2$ and $N_2/2$, respectively. Here we set $N_1=5$, $N_2=15$, the system is initially in the state $\left|\frac{N_1}{2},-\frac{N_1}{2}\right.〉\otimes \left|\frac{N_2}{2},\frac{N_2}{2}\right.〉$, and $n_{\text{T}}\ll 1$, thus neglecting collective pumping, ${\gamma }_{\Uparrow }$, in Eq. (29). (a) For ${\gamma }_{\Downarrow }={\omega }_0$ (solid curves), the first ensemble, initially unexcited, starts at $\langle J_z^{\left(1\right)}\left(0\right)\rangle /(N_1/2)=-1$ (black thin solid curve) and ends up in an excited steady state, while the initially excited second ensemble, starting at $\langle J_z^{\left(2\right)}\left(0\right)\rangle /(N_2/2)=1$ (red thick solid curve), ends up deexcited; these solid curves qualitatively reproduce the results of Ref. [86]. We investigate the effect of local incoherent emission, ${\gamma }_{\downarrow }$, adding it to both ensembles, showing that the excitation exchange becomes only transient, as the excitation lost by the second ensemble (red thick dashed curve) is partly acquired by the first ensemble (black thin dashed curve) at short times; eventually also the second ensemble relaxes to the ground state. (b) We show the effect of collective emission from each of the two ensembles, ${\gamma }_{\Downarrow ,i}={\gamma }_{\Downarrow }$ (dot-dashed curves), as well as that of local pure dephasing, for ${\gamma }_{\varphi }={\gamma }_{\Downarrow }$ (dotted curves). As in panel (a), the system is initialized in an antisymmetric state, with the first ensemble not excited at $t=0$ (black thin curves) and the second ensemble fully excited (red thick curves). The delay time is $t_{\text{D}}=ln\left(N_2\right)/N_2{\gamma }_{\Downarrow }$.

Figure 13

Ultrastrong-coupling regime. We compare the solution of the master equation valid below the ultrastrong-coupling (USC) regime and for effective models (black dashed curve) with the correct USC master equation (blue solid curve) for $N=10$ TLSs in a lossy cavity. The parameters are ${\omega }_{\text{cav}}={\omega }_0$, $g=0.1{\omega }_0$, ${\gamma }_{\downarrow }=\kappa =0.01{\omega }_0$, and the photonic Hilbert subspace has a cutoff of $n_{\text{ph}}=3$. (a) The intracavity photon population, $\langle a^{†}a\rangle \left(t\right)$, is calculated with two different Liouvillians, relative to the USC regime (solid blue curve) and to the non-USC regime (dashed black curve). We subtract to these quantities the photon population in the steady state of the USC Liouvillian ${\langle a^{†}a\rangle }_{\text{GS}}$. The initial state is the tensor state $\left|\frac{N}{2},-\frac{N}{2}\right.〉\otimes \left|0_{\text{ph}}\right.〉$. The different dynamics prompt different time evolutions tending to the USC steady state (dashed yellow line) and non-USC steady state (solid green line). (b) The extracavity photon emission spectrum $S\left(\omega \right)$ is shown, according to the two calculations, here initializing the system in the ground state of the light-matter Hamiltonian. The dotted vertical red lines correspond to the polariton frequencies $\omega ={\omega }_0\pm \sqrt{N}\frac{g}{2}$.