{'fraction_non_alphanumeric': 0.03611757985851751, 'fraction_numerical': 0.008436223032646425, 'mean_word_length': 4.597638957206099, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Engineering quasi-steady-state correlations in uncorrelated thermal states using stochastic driving

Nonequilibrium quantum dynamics can give rise to the emergence of novel steady states. We propose a scheme for driving an initially uncorrelated thermal state to generate customized correlation functions by determining and reverse engineering the steady-state two-point functions for a class of Markov processes. We also extend the formalism to the calculation of four-point functions. We then apply our method to generating power-law correlated fermionic Green's functions. Furthermore, we find that the power-law patterns emerge at much shorter times than the convergence to the steady state, at which point the disorder in the two-point correlations disappears. On the other hand, the density-density correlations exhibit steady-state disorder while following a power-law trendline. These ideal steady states appear as intermediate-time quasi-steady states in the presence of perturbations.


I. INTRODUCTION
There has been significant progress in driving nonequilibrium quantum matter to prepare desired quantum states, particularly in the context of quantum simulations [1,2]. While the objective is often creating states with equilibrium counterparts, such as the ground state of interacting Hamiltonians [3,4], it has been realized that quantum driving may give rise to the emergence of states without equilibrium counterparts. Floquet engineering, i.e., controlling the evolution by periodic driving, has been a fertile ground for generating such states [5][6][7]. Furthermore, dissipation engineering has been similarly effective for open systems [8][9][10][11][12][13].
In this paper, we focus on stochastic drives [14][15][16][17][18] corresponding to pure dephasing for noninteracting fermions [19], with formal similarities to open systems [20][21][22][23]. We propose a scheme for creating almost arbitrary two-point correlation functions starting from uncorrelated thermal states. Our method utilizes a single temporally uncorrelated noise source coupled to a fixed drive Hamiltonian, resulting in spatially correlated noise. The Lindblad master equation governing the dynamics only contains one double commutator, analogous to a fully dissipative open system.
We show that, with this scheme, we can engineer a broad class of custom-ordered equal-time two-point Green's functions (also known as the single-particle density matrix [24]). Long-range couplings [25][26][27] are necessary for creating long-range power-law correlations in the shortest possible time for a drive Hamiltonian with a fixed norm (See Ref. [28] for general results on limitations of creating critical correlations in open systems). As two-point correlations do not fully characterize a quantum system's many-body wave function, we also calculate the long-time limit of the four-point densitydensity correlation function to further investigate the nature of the steady states. We find that while the steady-state two-point function decays as a clean power law, the density-density correlator remains highly disordered but follows a power-law trendline with a different exponent.
Under ideal conditions, our pure dephasing dynamics allow the system to escape the generic fate of a noise-driven quantum state; namely, a featureless infinite-temperature fixed point originating from a noise-induced heating rate that constantly deposits energy into the system [17,29,30]. This scheme relies on the predominance of a single noise source, engineered using high-frequency pulse control. To achieve pure dephasing, we need to eliminate other noise sources or nonstochastic contributions to the Hamiltonian. In practice, however, even if we engineer our spatially correlated noise to have the largest energy scale, perturbations will ultimately force the system to flow to the infinite-temperature fixed point. However, for infinitesimally small perturbations (compared with the energy scale of the engineered noise, which can be made as large as possible within the experimental constraints), the timescales associated with the engineered noise and the perturbations separate. Thus, the steady states predicted in this paper for the ideal case are expected to emerge as slightly perturbed quasi-steady states at intermediate timescales before the system heats up to infinite temperature at much longer evolution times.
A noise term dominating the dynamics can only be generated in a controllable synthetic platform, with a strong random signal driving the system and naturally occurring noise sources reduced as much as possible. Two experimental platforms are potential candidates for our proposed correlation engineering method: cold fermionic atoms in optical lattices and digitized quantum processors with pulse level control. For a short-range drive Hamiltonian, fermionic atoms are the most convenient but creating direct long-range couplings is challenging with cold atoms. However, there has been recent progress in simulating long-range hopping with temporal driving [31][32][33]. On the other hand, general-purpose digitized devices can incorporate nonlocal couplings (either directly designed into the hardware or by using swap gates). Although qubits in these devices are not fermionic, the fermionic exchange statistics can be implemented through Jordan-Wigner or Braviy-Kitaev mappings [34][35][36][37][38].
The outline of this paper is as follows. In Sec. II, we extend the results of Ref. [30] on the steady-state two-point functions to steady-state four-point functions under stochastic dephasing dynamics. In Sec. III, we present our method for reverse engineering the drive to provide the desired matrix of fermionic equal-time Green's functions. In Sec. IV, we apply our method to the case of power-law correlated Green's functions. In Sec. V, we examine the steady-state densitydensity correlation functions for the drive used in engineering the Green's function. Finally, we present our conclusions in Sec. VI.

II. STEADY-STATE TWO-AND FOUR-POINT FUNCTIONS FOR FULLY STOCHASTIC QUADRATIC DRIVING
Here we present the theoretical framework for stochastic dynamics capable of generating customized noise-averaged correlations in fermionic systems, initially in an uncorrelated thermal state. For a system of N spinless fermions, where the creation/annihilation operators satisfy the anticommutation relation {c † i , c j } = δ i j , consider a quadratic Hamiltonian that stochastically fluctuates around zero. We write H(t) as Here, J is an N × N Hermitian matrix, and ξ(t) is Gaussian white noise with ξ(t) = 0 and ξ(t)ξ(t ) = w 2 δ(t − t ). The overline symbol represents averaging over the realizations of noise. The noise strength w sets the equilibration timescale associated with this noise term, but does not affect the nature of the steady state in the absence of perturbations. Although the couplings J i j are general, fixing the geometry is necessary for engineering correlations decaying as a power law of the distance between sites. We assume the sites as arranged on a one-dimensional line. The evolution by H(t) describes an isolated system with a vanishing average Hamiltonian driven by temporally uncorrelated and spatially correlated noise, with J encoding the spatial correlations. The system will reach a steady-state density matrix in the long-time limit t → ∞, given by where ρ(0) is the initial density matrix. Since the Hamiltonians at different times commute, time ordering is not necessary. We can view the dynamics as a Markov process described by a stochastic Schrödinger equation. The long-time limit of the dynamics can be alternatively viewed as a time average over random evolution times after a quench to Hamiltonian V, as the white noise is simply creating a random walk for the total evolution time with V. The equation of motion for the density matrix is then a Lindblad master equation . Similar to Ref. [30], we write a Heisenberg-picture analog of the master equation for a noise-averaged Heisenberg operator O(t) as Our analysis is focused on the unperturbed evolution above. However, in general we can have both stochastic and nonstochastic perturbations to H(t), which give rise to where H 1 is nonstochastic and j are the strengths of other noise terms. These perturbations introduce a longer timescale than the timescale associated with w, in which the system generically flows to the infinite-temperature fixed point with a density matrix proportional to the identity. In the remainder of this paper, we assume our timescales are shorter than the time needed to reach infinite temperature, neglect the perturbations, and analyze the problem using the ideal master equation (3). We also set w to 1, which amounts to a simple rescaling of evolution time.
An important observation is that under evolution with Eq. (3), generic quadratic (quartic) operators remain quadratic (quartic). In contrast, the evolution can spread operators into nontrivial linear combinations of other quadratic (quartic) operators. The above observation follows from the fact that in the commutator of two-quadratic forms, the quartic terms cancel out and we are left with another quadratic form . We now consider a general quadratic operator where column i of the unitary matrix U is the eigenvector of J with eigenvalue D i . In other words, We note that although all derivations in the paper are for white noise, the nature of the long-time limit also applies to a wider class of noise spectra with finite zero-frequency spectral power, e.g., the Ornstein-Uhlenbeck process effectively only rescales the evolution time compared to white noise [17]. Assuming J has a nondegenerate spectrum, we can then write the long-time limit of the noise-averaged Heisenberg operator for t → ∞ as We now extend this formalism to quartic operators. Of particular interest are density-density correlation functions. To find the steady-state correlation functions, again we use the Heisenberg-picture approach. Due to the cancellation of the sextic terms in the commutator of a quartic form and a quadratic from V = i j J i j c † i c j , we find where (9) Therefore, all quartic operators of the form above retain this structure when evolving with the master equation (3). By using Eq. (9) twice, we can write the equation of motion for the tensor R as where We now use the diagonalization of J [Eq. (6)] to find an exact solution to the equation above. We write the Kronecker deltas as elements of U † U to obtain In the absence of any accidental vanishing of the exponent, we can calculate the long-time limit of the R tensor by noting that we must have either b = a and d = c or b = c and d = a (which double counts the case with a = b = c = d) [40,41]. Thus we find For the density-density correlator However, we need to work with the more general tensors above as they appear in the Heisenberg-picture operator.

III. ENGINEERING STEADY-STATE GREEN'S FUNCTIONS IN UNCORRELATED THERMAL STATES
An important class of correlation functions, i.e., the equaltime single-particle Green's functions, are given by the expectation values where ρ(t) is the density matrix of the system at time t. We can collect these correlators into an N × N Hermitian matrix G, with elements G i j = G * ji . Let us consider a system initially prepared in a trivial thermal state at inverse temperature β with the density matrix with the partition function Z = tr[exp (−βH 0 )]. As H 0 is diagonal in the occupation-number basis, the different sites are uncorrelated in the thermal state above and we can write the initial Green's function as where f represents the Fermi-Dirac function. Suppose we drive this system by a stochastic quadratic Hamiltonian as in Eq. (1). Can we engineer custom-ordered Green's functions in the steady state of the system that emerges in the long-time limit? More precisely, suppose our goal is choosing µ i and J i j such that for given (perhaps even arbitrary) Hermitian correlation matrix G and initial inverse temperature β.
For the c † i c j operator, we have O ηγ (0) = δ ηi δ γ j , which leads to the following matrix elements when inserted in Eq. (5): Since we are engineering the matrix J, we can require it to have a nondegenerate spectrum. Thus, analogously to Eq. (7), the Heisenberg-picture matrix corresponding to the long-time limit of the c † i c j operator is given by The above expressions lead to the relationships below between the time-dependent and steady-state correlation functions and the initial correlation functions through the eigen-vectors of J: The above equation implies that Two conditions follow immediately: (i) G T and J must have the same eigenvectors since the same unitary operator U diagonalizes both of them, and (ii) the eigenvalues g σ of G T are the diagonal elements of A and related to the initial correlation functions as follows: where the matrix W of transition probabilities is defined through We have inserted Eq. (16) into Eq. (22) to obtain the condition (23). The matrix W is a doubly stochastic (it has nonnegative elements and the sum of every row and every column is equal to unity), which implies the sum rule α g α = α f (µ α ). Given a matrix G of desired steady-state correlators, we can find g σ and U by diagonalizing G. The matrix W can then be explicitly constructed. For a given initial inverse temperature β, we can then write the chemical potentials in the initial Hamiltonian as The columns of the matrix U = (|g 1 , |g 2 , . . . ) are eigenvectors |g α of G = α g α |g α g α |. The matrix J can then be any linear combination as long as the eigenvalues {d α } are nondegenerate. It may appear that the Hermitian matrix G can be completely arbitrary, and we can create any correlation functions we choose. Hermiticity is not the only requirement of a physical matrix of Green's functions, however. For example, all diagonal elements are the occupation number of a site and must be between 0 and 1, so creating states with nonphysical occupation numbers should be impossible. Furthermore, the method fails for many physically allowed Green's functions, as discussed below.
First, the matrix W must be invertible, which is not guaranteed. Indeed we have found that W tends to be singular for uniform systems, and some disorder in the occupation numbers is necessary for its invertibility. Second, the Fermi-Dirac functions are bounded as 0 f (µ α ) 1, which is equivalent to the condition that the chemical potentials written in Eq. (25) are real numbers. Therefore, the eigenvalues and eigenvectors of G must be consistent with the following constraint: for all j.
The constraint of Eq. (27) is not too restrictive if we are only interested in a pattern of correlations between different sites (off-diagonal elements of G) without simultaneously specifying on-site occupation numbers (diagonal elements of G). This observation is because we can change the eigenvalues of G by adding a multiple of the identity matrix to G, which only changes the occupation numbers, or multiplying G by a constant, which only rescales the Green's functions, without changing their overall pattern. More specifically, suppose for some matrix G, a = min σ W −1 jσ g σ < 0 and b = max σ W −1 jσ g σ > 1. Transforming G → 1 b−a (G − aI) by a shift to the diagonal elements and a rescaling guarantees the constraint (27) is satisfied. Here, I is the identity matrix.
While the initial chemical potentials µ j are uniquely determined by G and β as shown in Eq. (25), we have many choices for J i j . Essentially, the only constraint on J is through its eigenvectors, and we can choose the eigenvalues at will as long as they are nondegenerate. The absence of degeneracy is essential for our method. For example, in the extreme case where all eigenvalues d α are degenerate, the matrix J becomes proportional to the identity matrix, which will not change the initial Green's functions. The hopping matrix J is also entirely local in the limit of all degenerate eigenvalues. There is indeed competition between locality and equilibration time. To speed up the emergence of the steady state, we spread out these eigenvalues and make the gaps d α − d σ as large as possible.
The timescale for the equilibration goes as according to Eq. (5). Of course, we can arbitrarily shorten the timescale if we allow the energy scale of the Hamiltonian to run to infinity. To find an optimal J with fixed energy scale, we set the spectral norm of J, i.e., the square root of the largest eigenvalue of the matrix J † J to unity, ||J|| = 1, which implies −1 d α 1. To minimize t eq , we can make all the gaps identical by choosing While using all permutations of the above values of d α give the same equilibration time, we have found that the smoothest site dependence in matrix J corresponds to the case where eigenvectors of G are ordered according to their eigenvalues, i.e., g α g α+1 .

IV. ENGINEERING POWER-LAW TWO-POINT FUNCTIONS
We now consider an example of creating off-diagonal Green's functions for an arbitrary exponent ν. An example of a target G i j for i j is shown in Fig. 1(a) for ν = 1.5. We found that requiring a uniform average density gives a singular W matrix. Thus, to make the matrix W invertible, we choose random values drawn from a uniform distribution for the diagonal elements G ii (note that G ii is the expectation value of the occupation number of fermionic mode and must be between 0 and 1 for any physical system). Our results on the structure of the Hamiltonian J, the emergence of the steady-state correlations and the behavior of the four-point functions are not sensitive to the realization of diagonal disorder.
To apply our method, we find the eigenvalues and eigenvector of G, which in turn yield the initial chemical potentials according to Eq. (25). As long as the initial Fermi-Dirac functions f (µ i ), which only depend on the eigenvalues and eigenvectors of G, are between 0 and 1, we can find real chemical potential for any β. As an example, we show in Fig. 1(b) the f (µ i ) that yield the Green's functions shown in Fig. 1, together with the corresponding G ii . The initial occupation numbers f (µ i ) are close to the target occupation numbers G ii but not identical. The smallest difference for i = 47 is 0.13%,  while the largest difference for i = 28 is 53%. We note that while the initial occupation numbers are redistributed in the final steady state, the total charge is conserved, and we have In Fig. 2, we show the elements of the J matrix of Eq. (26) obtained by eigenvalues (29). We find that the elements of J exhibit similar critical correlations as the correlations we seek to create. We examined many different permutations of the eigenvalues. While the J can become highly disordered, it appears that the correlation structure does not change from critical as long as we keep the eigenvalues separated for optimal equilibration.
It is illuminating to examine the time dependence of the two-point function as it approaches the engineered critical correlations. As shown in Fig. 3, they start from zero and quickly acquire a noisy pattern aligned with the overall power-law dependence. The evolution continuously shifts the correlators up in a log-log plot while reducing the noise until it converges to the target disorder-free power law. The exponent of the critical correlation emerges much faster than the quantitative emergence of those correlation functions.

V. DENSITY-DENSITY CORRELATIONS IN THE STEADY STATE
To understand the nature of the steady state we have generated, we compute the density-density correlation functions, as the two-point function does not fully characterize the steady state. A widely studied class of correlation functions is the connected density-density correlator We treat the two terms in the above expression separately. To find n in j , we use the formalism of Sec. II to obtain the longtime limit of the rank-4 tensor R for the density-density operator. We can then obtain the time dependence as well as the long-time limit of the n in j correlator in terms of the initial four-point-function expectation values c † α c β c † γ c δ (0) at t = 0, which using the Wick's theorem is given by The long-time limit for t → ∞ follows from Eq. (13) as shown below (16). Inserting Eq. (32) into the above expression, upon some matrix algebra, leads to where the W matrices are defined in Eq. (24). The subtraction of the background n i n j is subtle. It might appear that each n i in the long-time limit is simply a diagonal element of G. However, that is only correct if we first take the noise average of each expectation value of density and multiply the noise-averaged results. It is more natural to take the quantum averages of density for various realizations of noise first, multiply them, and then take the noise average, keeping the noise averaging consistently in the last step. The long-time limit can be obtained by time averaging the product of expectation values upon a quantum quench to Hamiltonian V. In this approach, the background term yields two terms that are identical to terms in n in j (∞) . We find where the matrixF = U † diag(F 1 , . . . , F N )U . Using these expressions we calculated the long-time limit of R i j for the same drive that produces power-law two-point functions of Eq. (30). As shown in Fig. 4, while the long-time limit of G has no disorder in the off-diagonal terms, the density-density correlation functions are highly disordered. The minima in |R i, j | decay as | j − i| −3 for ν = 1.5. It seems that the steady states emerging upon engineering clean x −ν two-point functions have highly disordered density-density four-point functions with a trendline scaling as x −2ν .

VI. CONCLUSIONS
In this paper, we extended the Heisenberg-picture formalism of stochastic driving to the case of fermionic four-point functions. We then proposed a scheme of reverse-engineering the drive to generate customized fermionic two-point functions from a completely uncorrelated initial thermal state. Finally, we examined the case of generating critical correlations and found that the steady state exhibits clean engineered twopoint functions and highly disordered density-density correlation functions. Such steady states with custom-ordered correlations provide novel examples of drive engineering.

VII. ACKNOWLEDGMENTS
I thank Michael Kolodrubetz and Zohar Nussinov for helpful discussions. This work was supported by NSF Award No. DMR-1945395. I thank the Kavli Institute for Theoretical Physics, which is supported by the NSF Grant No. PHY-1748958, for hospitality. I also thank the Aspen Center for Physics, which is supported by NSF Grant No. PHY-1607611, for hospitality.