Optimal stochastic unraveling of disordered open quantum systems: application to driven-dissipative photonic lattices

TL;DR

This paper introduces an efficient stochastic method for averaging observables in disordered open quantum systems, simplifying the sampling process by pairing each trajectory with a single disorder configuration, and applies it to study disorder effects in driven-dissipative photonic lattices.

Contribution

The paper presents a novel optimal sampling technique for disordered open quantum systems, reducing computational complexity by linking each trajectory to one disorder configuration.

Findings

Disorder influences dissipative physics of fermionized bosons in 1D chains.

On-site inhomogeneities affect phase transitions in 2D driven-dissipative lattices.

The method efficiently computes configuration averages in complex quantum systems.

Abstract

We propose an efficient numerical method to compute configuration averages of observables in disordered open quantum systems whose dynamics can be unraveled via stochastic trajectories. We prove that the optimal sampling of trajectories and disorder configurations is simply achieved by considering one random disorder configuration for each individual trajectory. As a first application, we exploit the present method to the study the role of disorder on the physics of the driven-dissipative Bose-Hubbard model in two different regimes: (i) for strong interactions, we explore the dissipative physics of fermionized bosons in disordered one-dimensional chains; (ii) for weak interactions, we investigate the role of on-site inhomogeneities on a first-order dissipative phase transition in a two-dimensional square lattice.

Figures

Equations

$${\partial }_t\hat{\rho }(t)=\mathcal{L}\hat{\rho }=-i\text{commutator}\hat{H}\hat{\rho }+\sum _j\frac{{\gamma }_j}{2}\mathcal{D}[\widehat{L_j}]\hat{\rho }(t),$$

$$\mathcal{D}[\widehat{L_j}]\hat{\rho }(t)=2{\hat{L}}_j\hat{\rho }{\hat{L}}_j^{†}-\text{anticommutator}{\hat{L}}_j^{†}{\hat{L}}_j\hat{\rho },$$

$$\mathcal{L}(\mathbf{w})={\mathcal{L}}_{\mathrm{clean}}+{\mathcal{L}}_D(\mathbf{w}).$$

$$\hat{\rho }(t,\mathbf{w})={\mathrm{e}}^{\mathcal{L}(\mathbf{w})t}{\hat{\rho }}_0,$$

$$⟨\hat{O}⟩(t,\mathbf{w})=\text{Tr}\left[\hat{O}\rho (t,\mathbf{w})\right].$$

$$\overline{⟨\hat{O}⟩}$$

$$=\underset{R\to \infty }{\lim }\frac{1}{R}\sum _{i=1}^R⟨\hat{O}⟩({\mathbf{w}}_i),$$

$$\overline{\hat{\rho }}(t)=\left[\int d\mathbf{w}{p}_D(\mathbf{w})e^{\mathcal{L}(\mathbf{w})t}\right]{\hat{\rho }}_0,$$

$$\exp [{\mathcal{L}}_{\mathrm{ave}}t]=\int d\mathbf{w}{p}_D(\mathbf{w})e^{\mathcal{L}(\mathbf{w})t},$$

$$\hat{\rho }={\int }_{\mathcal{H}}d\psi p(\psi )|\psi ⟩⟨\psi |,$$

$${ϵ}_{\mathrm{clean}}=\alpha \sqrt{\frac{{\mathrm{Var}}_{\psi }[\hat{O}]}{T}},$$

$$\overline{⟨\hat{O}⟩}=\int d\mathbf{w}{p}_D(\mathbf{w})\int d\psi p(\psi |\mathbf{w})o(\psi ),$$

$${ϵ}^2$$

$$={\alpha }^2\frac{\overline{⟨{\hat{O}}^2⟩-⟨{\hat{O}⟩}^2}+\left[\overline{⟨\hat{O}{⟩}^2}-{\left(\overline{⟨\hat{O}⟩}\right)}^2\right]}{R}.$$

$$ϵ=\alpha \sqrt{\frac{\overline{{\mathrm{Var}}_{\psi }(\hat{O})}+{\mathrm{Var}}_D(⟨\hat{O}⟩)}{R}},$$

$$ϵ(R,T)=\alpha \sqrt{\frac{\overline{{\mathrm{Var}}_{\psi }(\hat{O})}+T{\mathrm{Var}}_D(⟨\hat{O}⟩)}{RT}}.$$

$${\hat{H}}_{\mathrm{BH}}=\sum _{j=1}^L\left[-{\Delta }_j{\hat{a}}_j^{†}{\hat{a}}_j+\frac{U}{2}{\hat{a}}_j^{†2}{\hat{a}}_j^2+F\left({\hat{a}}_j^{†}+{\hat{a}}_j\right)\right]-J\sum _{⟨i,j⟩}{\hat{a}}_i^{†}{\hat{a}}_j,$$

$${\partial }_t\hat{\rho }(t)=-i\text{commutator}{\hat{H}}_{\mathrm{BH}}\hat{\rho }+\frac{\gamma }{2}\sum _{j=1}^L\left[2{\hat{a}}_j\hat{\rho }{\hat{a}}_j^{†}-\text{anticommutator}{\hat{a}}_j^{†}{\hat{a}}_j\hat{\rho }\right].$$

$$\begin{gathered} \frac{d|\psi (t)⟩}{dt}=-i\left(\hat{H}(\mathbf{w})-i\gamma \sum _j{\hat{a}}_j^{†}{\hat{a}}_j\right)|\psi (t)⟩+ \\ +\sum _j{N}_j(t){\hat{a}}_j|\psi (t)⟩, \end{gathered}$$

$${\Psi }_B(i_1,i_2,\dots ,i_N)=\prod _{k<j}^N\mathrm{sgn}(i_k-i_j){\Psi }_F(i_1,i_2,\dots i_N),$$

$$\mathcal{E}(k_1,\dots ,k_N)=-N{\Delta }_0-2J\sum _{i=1}^N\cos (k_i).$$

$$\begin{gathered} \frac{d{\alpha }_j}{dt}=\left[-i({\Delta }_j-U(|{\alpha }_j{|}^2-1)-\gamma /2)\right]{\alpha }_j- \\ -iJ\sum _{⟨j^{\prime },j⟩}{\alpha }_{j^{\prime }}+iF+\sqrt{\gamma /2}\chi (t), \end{gathered}$$

$${\mathcal{S}=\frac{\partial n_{ss}(F/\gamma )}{\partial (F/\gamma )}|}_{F=F^{\ast }},$$

$$\mathcal{S}\sim W^{-\beta },$$

$$g^{(1)}(r=|\mathbf{l}-\mathbf{m}|)=⟨{\hat{a}}_{\mathbf{l}}^{†}{\hat{a}}_{\mathbf{m}}⟩-⟨{\hat{a}}_{\mathbf{l}}^{†}⟩⟨{\hat{a}}_{\mathbf{m}}⟩,$$

$$g^{(1)}(r)\sim r^{-\alpha }.$$

$$g^{(1)}(r)\sim {\mathrm{e}}^{-r/\lambda (W)},$$

Full text

I Introduction

The simulation of open quantum many-body systems is a formidable task. This is due to the exponential growth of the Hilbert space with the system size combined with the need of representing mixed states Breuer and Petruccione (2007). While systems in contact with the surrounding environment are expected to relax towards a (generalized) Gibbs ensemble Ilievski et al. (2015); Langen et al. (2015), in the most general scenario a system may not thermalize and evolve towards a non-equilibrium state. This situation can be reached by driving the system away from its thermodynamic equilibrium Lee et al. (2013a); Le Boité et al. (2013); Carusotto and Ciuti (2013) or by properly engineering the system-bath couplings in order to let the system evolve toward some target attractor of the dynamics Kapit (2017); Müller et al. (2012); Diehl et al. (2008); Wolf et al. (2011).

In the last decade, impressive experimental advances made possible to study the open nonequilibrium dynamics of large ensambles of quantum particles. There exist several experimental platforms that allow to explore this physics, such as photonic cavities Jacqmin et al. (2014); Rodriguez et al. (2017); Carusotto and Ciuti (2005), superconducting resonators Houck et al. (2012); Fitzpatrick et al. (2017); Labouvie et al. (2016), optomechanical resonators Aspelmeyer et al. (2014) and Rydberg atoms Lee et al. (2013b); Viteau et al. (2012); Saffman et al. (2010). Their high versatility allows to simulate the physics of several interesting Hamiltonians Georgescu et al. (2014), where most of the dynamical, spectral and correlation properties would be otherwise difficult to access. Indeed, their lossy nature has been repeatedly exploited in order to directly probe their state by observing the emission in the environment. For this reason, the stabilization and control of non-trivial, nonequilibrium many-body states in open systems has attracted an increasing number of theoretical works (see e.g. Jin et al. (2018); Rota et al. (2018a); Ma et al. (2018); Biella et al. (2017); Lebreuilly et al. (2017); Schiró et al. (2016); Jin et al. (2016); Wilson et al. (2016); Bardyn et al. (2013); Roncaglia et al. (2010) and references therein). Upon the variation of an Hamiltonian parameter and/or of the coupling rate with the surrounding environment, these systems can be driven across a (dissipative) phase transition in a properly defined thermodynamic limit Minganti et al. (2018). The peculiar behavior at criticality and the mechanism triggering the transition do not always follow the paradigms of equilibrium statistical mechanics Rota et al. (2017); Marino and Diehl (2016); Scarlatella et al. (2018); Schütz et al. (2016).

This scenario can be further enriched if one considers the effects of disorder. For closed systems, the interplay of interaction and disorder in a quantum many-body system can lead to remarkable effects such as the breakdown of ergodicity due to quantum effects, a phenomenon known as many-body localization Basko et al. (2006); Oganesyan and Huse (2007); Nandkishore and Huse (2015). Similar studies have recently been performed in the context of open systems Xu et al. (2018); Vakulchyk et al. (2018); Máximo et al. (2015); Celardo et al. (2017); Biella et al. (2013); Celardo et al. (2012); See et al. (2018) in order to unveil the competition between dissipation and localization phenomena. However, the theoretical investigations of the nature of far-from-equilibrium correlated phases in the presence of disorder have only scratched the surface Maghrebi (2017). Moreover, while the non-equilibrium extension of stochastic mean-field theory has been developed Kulaitis et al. (2013), the scarcity of effective numerical and analytical tools to tackle the problem beyond the mean-field paradigm has limited the study of the effects of disorder on extended systems. Indeed, dealing with disordered systems introduces an additional layer of complexity with respect to the clean problem because it requires to average over a number of disorder realisations to extract the properties of the system. This motivated us to develop an efficient numerical approach to simulate the dynamics of generic open many-body quantum systems in presence of disorder.

For clean homogenous systems, a few theoretical and computational methods have been put forward in the last years to simulate the physics of open systems under generic nonequilibrium conditions, such as renormalisation group calculations exploiting the Keldysh formalism Sieberer et al. (2016); Maghrebi and Gorshkov (2016), corner-space renormalization Finazzi et al. (2015), permutation-invariant solvers Shammah et al. (2018), full configuration-interaction Monte Carlo Nagy and Savona (2018), tensor-network techniques Mascarenhas et al. (2015); Cui et al. (2015); Jaschke et al. (2019); Werner et al. (2016); Kshetrimayum et al. (2017) and cluster methods Biella et al. (2018); Jin et al. (2016).

A convenient way to model open quantum systems is via stochastic quantum trajectories Carmichael (1998); Plenio and Knight (1998); Gambetta et al. (2008); Daley (2014), which can be used to compute expectation values and to reveal important dynamical features. According to the specific stochastic protocol, pure quantum states evolve according to stochastic jump-diffusion processes. The density matrix of the system can be retrieved by a proper average over an ensemble of trajectories. Quantum trajectories can also be combined with many-body methods, such as those based on matrix product states Grujic et al. (2012), renormalization group approaches Finazzi et al. (2015), Gutzwiller mean-field Casteels et al. (2018), semiclassical approaches Jäger et al. (2017), or with an analytical ansatz on the nature of the states Verstraelen and Wouters (2018). Critical phenomena in extended bosonic Vicentini et al. (2018) and spin Mascarenhas (2017); Rota et al. (2018b) lattices have been studied thanks to the Monte Carlo sampling of those trajectories. The main computational advantage of this approach is that one only needs to evolve pure states instead of the full density matrix of the system, reducing the memory cost.

In this paper we propose an optimal stochastic trajectory approach to efficiently deal with disordered systems. We show how the additional computational cost needed to evaluate disorder configurations can be strongly reduced by sampling the disorder distribution with stochastic trajectories, evolving according to different disorder configurations. We prove, under general assumptions, that our approach provides the best estimate of the exact disorder-averaged expectation values at a fixed computational cost.

Our method can be applied to generic open quantum many-body systems. In the present work, we focus on the paradigmatic driven-dissipative Bose-Hubbard model Vicentini et al. (2018); Foss-Feig et al. (2017); Le Boité et al. (2013) on a square lattice. We test our method against exact simulations and we combine it with several trajectory protocols to explore both the strongly and the weakly interacting regime.

In the strongly interacting limit, we studied the impact of disorder on the spectral properties of a driven-dissipative fermionized gas of bosons in one dimension Girardeau (1960). Tuning the frequency of the pump, one can excite different many-body resonances which, in the clean system, respect a selection rule based on the conservation of the quasi-momentum Carusotto et al. (2009). We investigate how the addition of disorder implies the violation of this rule and allows to access new resonances with spectroscopic experiments.

The driven-dissipative Bose-Hubbard model in two dimensions is known to also undergo a first-order dissipative phase transition from a low- to a high-density phase as the driving strength is increased Vicentini et al. (2018); Foss-Feig et al. (2017). At the critical point, the system is bistable and the whole lattice site population jumps simultaneously between two metastable states Minganti et al. (2018). We study the fate of this transition in the presence of on-site disorder in the weakly interacting regime. For closed systems, the effect of disorder on a classical Proctor et al. (2014); Imry and Ma (1975); Imry and Wortis (1979) and quantum Hrahsheh et al. (2012) first-order phase transition is not as well known as the second-order. In the case of a dissipative criticality this issue, to the best of our knowledge, has never been addressed so far.

Beyond the fundamental theoretical interest, the method and the findings of our work are relevant for ongoing experiments in extended driven-dissipative photonic lattices Hartmann (2016); Noh and Angelakis (2017). Despite the impressive degree of controllability at single and few-site level, the presence of inhomogeneities when assembling several fundamental units is, at present, unavoidable. In the case of photonic lattices realised via semiconductor cavities and superconducting circuits the main source of disorder is local, i.e. the bare cavity frequency varies from site to site. So that, in order to predict the outcome of experiments (often analyzed by means of quantum and homodyne trajectories) one need to access the role of inhomogeneities and its relation to collective phenomena.

The paper is organized as follows: In Sec. II we introduce the general framework to study disordered open quantum systems while in Sec. III we describe our method. The driven-dissipative Bose-Hubbard model is introduced in Sec. IV, together with numerical verification of the optimality of our strategy. We then exploit the optimal stochastic unraveling method to investigate the strongly interacting limit of this model in Sec. IV.1, and the first-order dissipative transition in the weakly interacting regime in Sec. IV.2. Finally, in Sec. V we draw our conclusions and discuss future perspectives.

II Disordered open quantum systems

In this paper we will consider quantum systems coupled to memoryless environments which leads to a Markovian dynamics for the system density matrix $\hat{\rho}$. In this case, the dynamics is governed by a master equation in the Lindblad form ($\hbar=1$) Breuer and Petruccione (2007)

[TABLE]

where $\mathcal{L}$ is the Liouvillian superoperator. The commutator on the right hand side of Eq.(1) accounts for the unitary part of the dynamics ruled by the Hamiltonian $\hat{H}$ while the incoherent dissipative processes are given by

[TABLE]

where $\hat{L}_{j}$ is the set of the jump operators representing different dissipative channels acting at a rate $\gamma_{j}$. At long times, the system will approach the steady-state $\hat{\rho}^{\rm SS}=\lim_{t\to\infty}\rho(t)$ (and thus $\partial_{t}\hat{\rho}^{\rm SS}=0$) regardless of the initial conditions. In what follows we will consider situations where $\hat{\rho}^{\rm SS}$ is guaranteed to be unique Nigro (2018); Albert and Jiang (2014).

In a system with static disorder, we take into account random variations of the parameters describing the system, which we encode in the vector of random variables $\mathbf{w}$. Those inhomogeneities satisfy a certain probability distribution $p_{D}(\mathbf{w})$, which will determine how likely a certain configuration $\mathbf{w}$ is to occur. The open system dynamics of the disordered system is obtained by adding to the Lindbladian of the clean system $\mathcal{L}_{\rm clean}$ a second Lindblad term $\mathcal{L}_{D}(\mathbf{w})$, which accounts for both coherent and incoherent disorder-induced processes

[TABLE]

Let us note that $\mathcal{L}(\mathbf{w})$ is still a well-defined Lindbladian according to Eq.(1) and can be formally integrated to give the time-evolution of the density matrix for the system in the presence of static disorder,

[TABLE]

where $\operatorname{\hat{\rho}}_{0}$ is the state at $t=0$. The expectation value of an observable $\hat{O}$ at a given time $t$ for a specific disorder realization $\mathbf{w}$ reads as

[TABLE]

In order to access the statistical properties of the system we need to take the average of observables over the possible disordered configurations $\mathbf{w}$, weighted by their probability $p_{D}(\mathbf{w})$. We will refer to this average as the configuration average, and will represent it with an overline $\overline{\bullet}$. Dropping the time-dependance we get:

[TABLE]

where, in the second line, we introduced the sum over the set of $R$ disorder configurations $\{\mathbf{w}_{i}\}_{i=1,\dots,R}$ sampled according to the distribution $p_{D}(\mathbf{w})$.

Note that it is always possible to define a disorder-averaged density matrix as

[TABLE]

such that $\Tr[\overline{\operatorname{\hat{\rho}}}(t)\hat{O}]=\overline{\langle\hat{O}\rangle}$ for any given observable $\hat{O}$. Through Eq.(8), it is also possible to define the averaged (or quenched) propagator

[TABLE]

which governs the averaged dynamics of the system. This would allow us to obtain directly the averaged dynamics of the system in the particular cases when the formal expression (9) can be exactly integrated.

III Optimal stochastic unraveling of disordered open quantum systems

In this section we present our approach, sketched in Fig.1, to the problem of disordered open quantum systems. In Sec. III.1 we briefly review trajectory-based methods for clean systems. In Sec. III.2 we introduce our extension to efficiently compute disorder-averaged expectation values. A numerical demonstration of the optimality of our method will be given in Sec. IV.

III.1 Stochastic trajectories

Let us start by briefly reviewing the main ideas of the stochastic trajectory approach to clean open quantum systems ($\mathcal{L}_{D}=0$). To compute numerically the state of the system at a given time, one needs to evolve the density matrix of the system according to Eq. (1). This method is limited by very high memory requirements and it becomes unfeasible as the number of sites is large enough. Indeed, one need to store and evolve the full density matrix $\hat{\rho}$ which is a Hermitian matrix of dimension $d^{N}\times d^{N}$ for a system made of $N$ sites with local Hilbert space having dimension $d$.

As already mentioned in the introduction, a possible approach with a lower memory cost involves the stochastic unravelling of the density matrix by means of trajectories evolving according to a given protocol Daley (2014). The unravelling can be formally written as an integral over the space of trajectories

[TABLE]

where we omitted the time-dependance. Such procedure allows us to interpret the density matrix at each time $t$ as the mixture of several pure states $\ket{\psi}$, and to recast the master equation (1) into a jump-diffusion stochastic differential equation evolving a set of pure states Carmichael (1998). In practice, the distribution $p(\psi)$ is sampled by independently evolving a finite number of trajectories. As this is a stochastic sampling of the distribution, the error of the estimated expectation value $\hat{O}$ scales as 111This assumes that all trajectories are independent. Because the initial state is fixed, for small times below the correlation time the trajectories are not independent variables, and a covariance correction to the variance should be included through the substitution $\operatorname{\mathbb{V}ar}_{\psi}\rightarrow\operatorname{\mathbb{V}ar}_{\psi}+\operatorname{\mathbb{C}ov}_{\psi}$. As this is related to the correlation time, it must decay exponentially in time. Note that this only affects the efficiency of the sampling, but not the expectation value, which is unaffected.

[TABLE]

where $T$ is the number of trajectories and $\operatorname{\mathbb{V}ar}_{\psi}[\hat{O}]$ is the statistical variance of the observable across the trajectories which depends on the specific evolution protocol and on the observable $\hat{O}$. The prefactor $\alpha=\epsilon/\sigma$ arises from the definition of the absolute deviation $\epsilon=\langle|\hat{O}-\langle\hat{O}\rangle|\rangle$, as opposed to that of the standard deviation $\sigma=\sqrt{\langle\hat{O}^{2}-\langle\hat{O}\rangle^{2}\rangle}$. As it is obvious from the definition, both are measures of the deviation (spread) of the distribution $\mathcal{O}$, but while the first has a more direct interpretation when sampling, the latter has much better mathematical properties. It is possible to show that $\alpha\in[0,1]$, and we will simply treat it as a trivial proportionality factor with a value usually around that of the normal distribution, namely $\alpha=\sqrt{2/\pi}$ (Kenney and Keeping, 1947, pgs. 81-85 and 113).

III.2 Stochastic trajectories for disordered systems

We can now introduce our method to perform an efficient estimate of disorder-averaged expectation values. Our aim is to efficiently compute the expectation value defined by Eq. (6) when the density matrix is unravelled stochastically. By inserting Eq.(10) into Eq.(6) we get:

[TABLE]

where $p(\psi|\mathbf{w})=p(\psi,\mathbf{w})/p_{D}(\mathbf{w})$ is the conditional probability of obtaining the state $\ket{\psi}$ given the disorder configuration $\mathbf{w}$ and $o(\psi)=\bra{\psi(\mathbf{w})}\hat{O}\ket{\psi(\mathbf{w})}$. Sampling the disorder-averaged expectation value can therefore be interpreted as performing the sampling of multivariate distribution $p(\psi,\mathbf{w})$ which depends on $\psi$ and $\mathbf{w}$.

There exist several ways to perform such a sampling. Indeed, we have freedom to choose how many trajectories we want to evolve according to a given disorder configuration. In the general case one evolves $T$ trajectories for each of the $R$ disorder realizations for a total computational cost $\mathcal{C}=R\times T$. We will show that, when computing disorder-averaged expectation values, the most computationally efficient sampling strategy is realized for $T=1$, i.e. one needs to evolve only one trajectory for each disorder configuration. This strategy is sketched in Fig. 1. In this case, the error committed when computing the expectation value is

[TABLE]

The first term $\overline{\langle\hat{O}^{2}\rangle-\langle{\hat{O}\rangle}^{2}}=\overline{\operatorname{\mathbb{V}ar}_{\psi}[{\hat{O}}]}$ is the disorder-averaged statistical variance of the observable [see Eq.(11)] and accounts for the uncertainty due to the trajectory protocol. This is the same contribution found in the clean case described by Eq.(11). The second term $\overline{\langle\hat{O}\rangle^{2}}-\left(\ \overline{\langle\hat{O}\rangle}\ \right)^{2}=\operatorname{\mathbb{V}ar}_{D}[\langle\hat{O}\rangle]$, instead, is the variance of the expectation value of the observable $\langle\hat{O}\rangle$ due to the presence of disorder. Exploiting these identifications, we can rewrite Eq.(14) in a more compact form as

[TABLE]

which is valid when sampling one trajectory for each disorder realization. The generalization of Eq.(15) to the $T>1$ case is obtained by performing the substitution $\operatorname{\mathbb{V}ar}_{\psi}[\hat{O}]\rightarrow\operatorname{\mathbb{V}ar}_{\psi}[\hat{O}]/T$ which reflects the reduction in the statistical error when evolving multiple trajectories. After some algebra, we obtain the general formula for the error as a function of $R$ and $T$

[TABLE]

We point out that in the absence of disorder, $\operatorname{\mathbb{V}ar}_{D}=0$ and the only term playing a role is the uncertainty associated with the trajectories. In this case, the distinction between different disorder configurations becomes meaningless, and Eq.(16) reduces to Eq. (11) where the total number of trajectories is $T\rightarrow\mathcal{C}$. Equivalently, when integrating the master equation exactly for $R$ different configurations of $\mathbf{w}$ one obtains a similar formula where the contribution of the trajectories is absent ($\operatorname{\mathbb{V}ar}_{\psi}=0$) and only $\operatorname{\mathbb{V}ar}_{D}$ contributes to the error.

From Eq. (16) it is easy to deduce that the most efficient sampling strategy, i.e. the choice of $T$ which minimizes the error for a fixed computational cost $\mathcal{C}=R\times T$, is realised for $T=1$. We remark that the only assumption we made to obtain Eq.(16) is to consider that the dynamics of the $T$ trajectories we evolve for each disorder realization is not correlated 222If the evolution is started from a pure state, then the trajectories will be correlated for times $t<\tau$ below the correlation time. Nevertheless, our method still yields correct results, but the error will scale slightly worse, as there are cross-correlations that are harder to average out., as required for an effective stochastic unraveling Gardiner and Zoller (2004). Thus Eq.(16) holds for arbitrary evolution protocols and disorder distributions.

IV Numerical applications to the driven-dissipative Bose-Hubbard model

As a first application of the optimal stochastic unraveling of disordered open quantum systems, we study the driven-dissipative Bose-Hubbard model with on-site static disorder. The coherent part of the evolution (in the frame rotating at the driving frequency and setting $\hbar=1$) is ruled by the following Hamiltonian

[TABLE]

where $\Delta_{j}=\omega_{p}-\omega_{j}$ is the local frequency detuning ($\omega_{p}$ and $\omega_{j}$ being the drive and the local bare frequency, respectively), $U$ is the on-site interaction, $F$ the strength of the coherent drive (the phase is set in such a way that $F\in\mathbb{R}$) and $J$ sets the hopping rate between nearest-neighbouring sites. The local detunings $\{\Delta_{j}\}$ are assumed to be random variables obeying a gaussian distribution $\mathcal{N}(\Delta_{0},W)$ with mean $\Delta_{0}$ and variance $W^{2}$. Dissipative processes, i.e. local boson leakage, are modelled with a set of local jump operators $\hat{L}_{j}=\hat{a}_{j}$ ($j$ being the site index) with uniform rate $\gamma$ leading to the following expression for the master equation

[TABLE]

Let us start by presenting a numerical verification of the optimality of our sampling method by considering a driven-dissipative Bose-Hubbard chain consisting of $L=5$ sites with periodic boundary conditions and interaction strength $U=\gamma$. First, we have computed the observables via a brute force integration of the master equation averaged over $2000$ disorder configurations, sufficient to reduce the statistical error below $5\cdot 10^{-3}$. We take the value of the average boson density $\hat{n}=\sum_{j=1}^{L}\hat{a}^{\dagger}_{j}\hat{a}_{j}/L$ computed in this way as the exact value $n_{exact}$ against which we will compare the sampling performed with trajectories.

We then unravel the dynamics of the master equation (Eq.18) in terms of quantum trajectories Carmichael (2008) within the so-called wavefunction Monte Carlo technique Dalibard et al. (1992); Mølmer et al. (1993). In this framework, a pure quantum state $\ket{\psi(t)}$ is evolved according to a non-Hermitian Hamiltonian in the presence of stochastic quantum jumps. This jump-differential equation reads

[TABLE]

where $N_{j}(t)$ is a Poissonian counter which is equal to 1 with probability $p_{j}(t)=\gamma\expectationvalue{\hat{a}_{j}^{\dagger}\hat{a}_{j}}{\psi(t)}/\innerproduct{\psi(t)}{\psi(t)}$ and 0 otherwise, encoding the stochastic nature of quantum jump trajectories. Note that the above equation does not preserve the norm of $\ket{\psi(t)}$, which should be re-normalized when computing observables. The disorder configuration is encoded in the vector of on-site energies $\mathbf{w}=\{\Delta_{i}\}$. For every configuration $\mathbf{w}$ we evolve $T=1,\ 10,\ 100,\ 1000$ trajectories. Considering $R$ different disorder configurations, we extract the steady-state average boson density $n_{i}$ for every $i=1,\dots,\mathcal{C}=RT$ trajectory.

In Fig.2 we show the mean absolute error $\epsilon(R,T)=\frac{1}{\mathcal{C}}\sum_{i=1}^{\mathcal{C}}|n_{i}-n_{exact}|$ computed both numerically and analytically (Eq.(16) with $\alpha=\sqrt{2/\pi}$), finding a perfect agreement. The fact that at a fixed computational cost $\mathcal{C}$ the lowest sampling error $\epsilon$ is achieved by considering $T=1$ further confirms the optimality of our method.

In what follows we will adopt this method to compute disorder averages and study the effect of on-site disorder on a first-order dissipative phase transition in the driven-dissipative Bose-Hubbard lattice.

IV.1 Strongly interacting (fermionized) bosons in 1D chains

In this section we will focus on the strongly interacting limit $U/J\gg 1$, where the boson-boson interaction is strong enough to forbid the double occupation of each site. In this hard-core limit, the many-body Hamiltonian (17) can be exactly diagonalized in one spatial dimension with periodic boundary conditions via fermionization of the bosonic wavefunction. Indeed, any $N$-body bosonic wavefunction $\Psi_{B}(i_{1},i_{2},\dots,i_{N})$ ($i_{1},i_{2},\dots,i_{N}$ being the positions of the $N$ particles) can be one-to-one mapped onto a fermionic one Girardeau (1960). For $N$ bosons, the wavefunction in real-space representation is given by

[TABLE]

where $\prod_{k<j}^{N}{\rm sgn}(i_{k}-i_{j})$ ensures that $\Psi_{B}$ respects the bosonic statistics. In the limit of vanishing pumping (i.e., $F=0$), the eigenstates of (17) can be simply labelled by the occupation number of single-particle eigenstates of a free-fermion system. The notation $\ket{k_{1},\dots,k_{N}}$ indicates the fermionic eigenstates obtained via the symmetrization of the wavefunction with one particle in each $k_{1},\dots,k_{N}$ orbital Carusotto et al. (2009). The energy of this state with $N$ particles, according to Eq. (17) is

[TABLE]

In a spectroscopic experiment, as the pump frequency is varied, one do expect an $N$-body resonance to appear whenever $\Delta_{0}=-2J\sum_{i=1}^{N}\cos(k_{i})/N$. In particular, when the system is homogeneously driven, only many-body states with total momentum $\mathcal{K}=\sum_{i}^{N}k_{i}=0$ can be accessed. We also remark that, because the pump term $F(\hat{a}+\hat{a}^{\dagger})$ only couples the vacuum with 1-body states $\ket{k_{n}}$ at first order in $F$, resonances with $N\geq 2$ will only arise when the pump strength $F$ is non-perturbative with respect to $\gamma$.

While the spectroscopy of the clean system has been studied in Refs.Biella et al. (2015); Carusotto et al. (2009), here we want to characterize the effect of local disorder on the resonances of the fermionized system. To this aim, we studied the average population in the steady-state $\mathcal{N}_{ss}=\sum_{j}\braket{\hat{a}_{j}^{\dagger}\hat{a}_{j}}$, a witness of the presence of those $N-$particle resonances, as a function of the detuning $\Delta_{0}$. In Fig. 3 (a) and (b) we show $\mathcal{N}_{ss}$ for a chain of five and nine resonators, respectively. The position of one- and two-body resonances which respect the selection rule $\mathcal{K}=0$ is marked by black and orange (light gray) dotted lines, respectively. The position of the others one- and two-body resonances with $\mathcal{K}\neq 0$ is denoted by dashed lines. In the panel (a) the continuous line has been obtained by averaging the exact solution obtained via exact diagonalization of the Liouvillian over $R=300$ disorder realizations for each value of $\Delta_{0}$. To show the advantage of our method introduced in Sec.III, we computed the same observables through a wavefunction Montecarlo algorithm Mølmer et al. (1993); Carmichael (1993); Daley (2014) where every quantum trajectory evolves according a different realization of disorder. The data obtained with our method (filled dots) are in very good agreement with the exact simualtions, and can be obtained (for a grid of equal spacing) with a computational cost that is roughly $2^{5}=32$ times lower.

Let us now focus on the position of the resonances. If $W=0$, only the resonances for which $\mathcal{K}=0$ are excited by the homogeneous drive. For small disorder (i.e. $W/\gamma<1$) the intensity of the $\mathcal{K}=0$ resonances is only affected perturbatively. However, one-body resonances with $\mathcal{K}\neq 0$ start to be visible in the spectrum. Indeed, the disorder perturbatively changes the shape of one-body states allowing to populate them even with an homogeneous pump profile.

Among the several two-body $\mathcal{K}\neq 0$ resonances (most of which are not shown), the one at $\Delta/\gamma=10\gamma$ is particularly interesting because it is the only one significantly populated. An intuitive explanation is obtained by considering that the drive couples one- and two-body states with the same total momentum $\mathcal{K}$. Because the coupling is inversely proportional to their energetic mismatch, only $\ket{\mp\frac{1}{5}\pi,\ \pm\frac{3}{5}\pi}$ which has $\mathcal{K}=\pm\frac{2}{5}\pi$ couple to the nearby $\ket{\pm\frac{2}{5}}$ one-body state since the pump conserves the total momentum $\mathcal{K}$. For large disorder ($W/\gamma>1$) this one- and two-body resonances start to overlap, eventually merging into a unique broad peak for larger values of $W$.

This disorder-enhanced emergence of resonances is expected to become more relevant for longer chains, where more $\mathcal{K}\neq 0$ resonances are now accessible. We exploited our technique to investigate a 9 site chain, which would otherwise be computationally much more demanding. The results are presented in the bottom panel of Fig.3. In this case, all the five one-body resonances becomes visible (black vertical lines) and the nearby two-body states with $\mathcal{K}\neq 0$ with the same $\mathcal{K}$ emerge.

We can conclude that, for the quite wide range of disordered considered, the resonances of the clean system are still observables when a spectroscopy of the system is performed. Furthermore, some of the many-body resonances unaccessible because of the symmetry of the pump, emerge in presence of disorder.

IV.2 Role of disorder on a first-order dissipative phase transition in 2D lattices

In this subsection we will study the fate of a first-order criticality in a weakly interacting driven-dissipative Bose-Hubbard lattice with disorder.

The master equation Eq.(18) can be recasted into a third-order partial differential equation for the Wigner quasiprobability distribution Carmichael (1998). Considering the regime of weak interactions, we can perform the truncated Wigner approximation by neglecting the third order terms, obtaining a well-defined Fokker-Planck equation. The time-evolution of the Wigner distribution can be efficiently sampled by solving the associated stochastic differential equation for the Wigner trajectories Vogel and Risken (1989)

[TABLE]

where $\expectationvalue{\chi(t)\chi^{\star}(t^{\prime})}=\delta(t-t^{\prime})$ and $\expectationvalue{\chi(t)\chi(t^{\prime})}=0$ is a delta-correlated noise. Expectation values are obtained by averaging over a large number of such trajectories. This approach has been shown to be very effective in the clean case for a range of nonlinearities up to $U\lesssim 0.5\gamma$ Vicentini et al. (2018), allowing for the simulation of large lattices up to $22\times 22$ sites. We solved Eq.(22) with a combination of recently developed stiff solvers Rackauckas and Nie (2017a, 2018) implemented in Julia Bezanson et al. (2017) within the DifferentialEquations.jl package Rackauckas and Nie (2017b).

In the homogeneous case ($W=0$), for some values of $\Delta_{0}$, the system is known to undergo a first-order phase transition in the steady-state from a low- ($F<F_{crit}$) to a high-density ($F>F_{crit}$) phase as the driving strength is increased Foss-Feig et al. (2017); Vicentini et al. (2018). At criticality the density matrix is bimodal: the steady-state density matrix is the symmetric mixture of the two phases Minganti et al. (2018). As a consequence, when $F\approx F_{crit}$, the system displays bistability: the whole lattice simultaneously jump between two metastable states (the low- and high-density phases). Such mechanism is the fingerprint of a dissipative first-order transition. The transition rate depends on the asymptotic decay rate towards the steady state $\tau^{-1}$ Vogel and Risken (1988), which vanishes in the thermodynamic limit leading to the critical slowing down in the system dynamics at the critical point Vicentini et al. (2018). With this in mind, we exploited the stochastic unraveling method presented in this article to efficiently simulate disordered lattices and investigate the interplay of disorder and criticality.

IV.2.1 Boson density

Let us start our analysis by studying the bosonic population across the transition. In Fig. 4 we computed the probability distribution of the average steady-state density $n_{ss}$ for several values of the disorder $W$ on a $N=L\times L$ lattice with $L=14$. In the insets we show typical snapshots of the site-resolved density along a single trajectory, where the formation of domains can be seen for stronger disorders. We note that the abrupt jump present in the clean case [see Fig. 4, panel (a)] becomes progressively smoother when influenced by the local disorder. Interestingly, even for disorder distributions well within the single-site linewidth ($W\ll\gamma$) the transition is rounded and replaced by a smooth crossover.

The underlaying mechanism responsible for the suppression of the criticality is the depletion of the homogeneous $\mathbf{k}=0$ mode of the lattice which is macroscopically occupied across the transition for $W=0$. In Fig. 4, we show the homogeneous fraction $f_{0}=n_{\mathbf{k}=0}/n=\sum_{i,j=1}^{N}\braket{\hat{a}^{\dagger}_{i}\hat{a}_{j}}/\sum_{i=1}^{N}\braket{\hat{a}^{\dagger}_{i}\hat{a}_{i}}$ in the steady state for different disorder strengths $W$. Local inhomogeneities in the energy landscape force the system to populate $\textbf{k}\neq 0$ modes, competing against the transition. We point out that the strongest depletion takes place for $F<F_{crit}$. Indeed, in the linear regime the term that dominates the local energy is proportional to the random local detuning. When $F>F_{crit}$ the leading term due to the nonlinearity $U$ grows as $n_{ss}^{2}$, making this phase less sensitive to the presence of disorder.

To better characterise the effect of disorder on the criticality, we studied the behavior of the slope at the transition

[TABLE]

where $F^{*}$ is the value of $F$ for which $\mathcal{S}$ is largest. In Fig. 6 we show $\mathcal{S}^{-1}$ as a function of $W^{-1}$ for different system sizes. As the size of the array is increased, the data progressively display a power-law behavior

[TABLE]

with $\beta=0.98\pm 0.05$. This result proves that the discontinuous population jump occurring in the thermodynamic limit is smoothened by disorder. Thus only a perfectly homogeneous system ($W=0$) would display a true criticality $\mathcal{S}\to\infty$ when $L\to\infty$ if one drives the $\mathbf{k}=0$ mode. This does not exclude that the criticality could be restored by engineering a space-dependent driving field which couples more effectively to the modes of the disordered system.

The value of the exponent $\beta$ has been extrapolated by fitting the data for $L=22$. Interestingly, any finite-size system made of $L\times L$ sites, follows the clean-system behavior in Eq. (24) up to a certain value $1/W^{*}$ which increases as $L$ is increased. This is due to finite-size effects and, in particular, the relation between $W^{*}$ and $L$ is connected with the typical correlation length of the systems. To get further insight about this point, in the following section we study the behavior of the correlation functions at the transition in clean and disordered systems.

IV.2.2 Correlation functions

The rounding of the transition due to disorder is also witnessed by the connected part of the one-body correlation functions

[TABLE]

where $\mathbf{l}$ and $\mathbf{m}$ are the $2D$ coordinates of lattice sites. In a coherently-driven system it is necessary to consider the connected part of the correlators in order not to account for the coherence imprinted by the external driving field 333Given that $\braket{\hat{a}_{\mathbf{i}}}\neq 0,\ \forall\mathbf{i}$ because of the local driving term one would simply get $\lim_{|\mathbf{l}-\mathbf{m}|\to\infty}\braket{\hat{a}^{\dagger}_{\mathbf{l}}\hat{a}_{\mathbf{m}}}\sim\braket{\hat{a}^{\dagger}_{\mathbf{l}}}\braket{\hat{a}_{\mathbf{m}}}\to{\rm const}$..

Results are shown in Fig. 7 444Each point is obtained considering the mean value of the correlator Eq.(25) averaged over all the pairs of sites at a distance $r$. In the clean case ($W=0$), the system displays long-range order and the correlation function decays as a power-law

[TABLE]

The exponent $\alpha=0.16\pm 0.01$ is obtained by fitting the data for $L=22$. This result is somehow unexpected since in proximity of a first-order transition one would not expect a divergent correlation length Sachdev (2001). However, this is a peculiar feature of a perfectly clean system and disorder suppresses long-range order.

Indeed, in a disordered system ($W\neq 0$), the correlation function (25) decays exponentially as

[TABLE]

where $\lambda(W)\sim W^{-\eta}$ (with $\eta=1.1\pm 0.2$) is the correlation length which decreases as $W$ is increased. This is a consequence of the formation of coherent domains with typical size $\lambda(W)$. This means that as long as one considers arrays with $L\lesssim\lambda(W)$, the system will behave cooperatively. Indeed, for a system of size $L$, we do expect departure from the critical thermodynamic behavior when $L>\lambda(W)\sim W^{-\eta}$. This is consistent to what has been shown in Fig. 6, since we find that $1/W^{*}\sim L^{1/\eta}$. Let us conclude this section remarking that all the exponent obtained here are not expected to be universal since the transition is of first order Vicentini et al. (2018).

V Conclusions and outlook

In this work we have presented a general sampling method to efficiently compute configuration averages of observables in disordered open systems, whose dynamics has been unraveled with stochastic trajectories. We have proven that the optimal strategy is realized when each trajectory evolves according to a different disorder configuration. Our approach allows us to drastically reduce the computational cost needed to perform disorder-averaged quantities, because the scaling of the statistical error with the number of trajectories and disorder configurations is independent of the specific stochastic evolution protocol and of the particular disorder distribution Remarkably, the present method can be applied to any Markovian non-equilibrium quantum system numerically solved via stochastic differential equations. Our method can represent an invaluable tool to address the dynamical and steady-state properties of a wide class of driven-dissipative systems.

As a first application, we have been able to study the role of disorder in the driven-dissipative Bose-Hubbard model both in the strongly and weakly interacting regime. In the case of strong interactions we have analyzed the emergence of spectral resonances due to fermionization of bosons in a disordered 1D chain. In the weakly interacting case we have explored the role of disorder on a dissipative first-order phase transition occurring in 2D lattices, showing how the criticality is suppressed. Thanks to the optimal unraveling method, we have been able to show that the correlation functions of the disordered system decay exponentially in space while for the clean case they decay algebraically at criticality. We have also shown that the derivative of the boson density scales critically (as in the clean system) up to sizes comparable with the typical size of the domains. The emerging picture is directly relevant for the ongoing experiments in state-of-art photonic quantum simulators based on superconducting circuits and semiconductor microcavities. In these systems, the presence of on-site disorder in the local frequency of resonators is unavoidable. Experimentally, by tuning the size (and/or the disorder) of a sample, one can explore different regimes where collective critical dynamics or density domains dominates the physics.

Many interesting perspective are open for the future. The case of spatially correlated noise can be treated within our formalism. The combination between our stochastic trajectory approach with cluster techniques Biella et al. (2018); Jin et al. (2016) and tensor-network ansatz Orús (2014) represents an intriguing direction. Furthermore, developing exact methods to calculate exactly the average Liouvillian dynamics can be useful to shade light on general properties on open many-body disordered systems. It would be also interesting to study the impact of disorder on a second-order dissipative criticality and to study models where long-range interactions competes with local disorder Santos et al. (2016).

Acknowledgements.

We thank N. Bartolo, N. Carlon-Zambon, I. Carusotto for stimulating discussions. We acknowledge support from ERC (via Consolidator Grant CORPHO No. 616233). This work was granted access to the HPC resources of TGCC under the allocations 2018-AP010510493 and 2018-A0050510601 attributed by GENCI (”Grand Equipement National de Calcul Intensif”).