Free Fermions with a Localized Source

TL;DR

This paper investigates the dynamics of free fermions on an infinite lattice coupled to a localized source, revealing how the particle growth rate depends on lattice geometry, dimensionality, and input rate, with quantum Zeno effects influencing the behavior.

Contribution

It provides explicit formulas for particle growth rates in free fermion systems with localized sources, highlighting non-monotonic dependence and quantum Zeno effects.

Findings

Growth rate depends on lattice geometry and dimensionality.

Non-monotonic relationship between input rate and growth rate.

Large input rate leads to inverse proportionality due to quantum Zeno effect.

Abstract

We study an open quantum system of free fermions on an infinite lattice coupled to a localized particle source. In the long time limit, the total number of fermions in the system increases linearly with growth rate dependent on the lattice geometry and dimensionality. We express the growth rate in terms of lattice Green functions and derive explicit formulae in one dimension and for the square lattice. The interplay between the dynamics and the coupling to the environment leads, in contrast to classical systems, to a non-monotonic dependence of the particle growth rate on the input rate. We show that for all lattices the particle growth rate is inversely proportional to the input rate when the latter becomes large. This is a manifestation of the quantum Zeno effect.

Equations

$$N(t)\simeq \{\begin{array}{ll}4\sqrt{\frac{t}{\pi }}& d=1\\ \frac{2\pi t}{\ln t}& d=2\\ \frac{\Gamma }{1+\Gamma W_d}t& d>2\end{array}$$

$$W_d=\frac{d}{2}{\int }_0^{2\pi }\dots {\int }_0^{2\pi }\prod _{i=1}^d\frac{d{q}_i}{2\pi }\frac{1}{{\sum }_{i=1}^d\left(1-\cos q_i\right)}$$

$${\partial }_t\rho =-\mathrm{i}[H,\rho ]+2L\rho L^{†}-\{L^{†}L,\rho \}$$

$$H=\sum _{n=-\infty }^{\infty }\left(c_n^{†}{c}_{n+1}+c_{n+1}^{†}{c}_n\right)$$

$$N(t)\simeq C_d(\Gamma )t$$

$$C_1(\Gamma )=2\Gamma -\frac{2{\Gamma }^2}{\pi }\left[{\gamma }^{-2}+\left({\gamma }^{-1}-{\gamma }^{-3}\right){\tan }^{-1}\gamma \right]$$

$$\gamma =\sqrt{(\Gamma /2{)}^2-1}$$

$$C_1(\Gamma )=2\Gamma +\frac{4{\Gamma }^2}{\pi }\left[\ln \frac{\Gamma }{4}+\frac{1}{2}\right]+O({\Gamma }^4\ln \Gamma )$$

$$C_1(\Gamma )=4{\Gamma }^{-1}-\frac{128}{3\pi }{\Gamma }^{-2}+36{\Gamma }^{-3}-\frac{4096}{15\pi }{\Gamma }^{-4}+\dots$$

$$C_{2}(\Gamma)=2\Gamma-\tfrac{16}{\pi}\,\Gamma^{2}\big{[}I_{1}(\Gamma)+I_{2}(\Gamma)\big{]}$$

$$I_1(\Gamma )$$

$$I_2(\Gamma )$$

$$K(x)={\int }_0^{\pi /2}\frac{d\theta }{\sqrt{1-x^2{\sin }^2\theta }}$$

$$C_2(\Gamma )=2\Gamma -A_2{\Gamma }^2+A_3{\Gamma }^3+\dots$$

$$\begin{array}{rl}{A}_2& =\frac{4}{{\pi }^3}{\int }_0^{1}dx\left[2K^2(x)+K^2(x^{\prime })\right]=1.8054571566\dots \\ A_3& =\frac{4}{{\pi }^3}{\int }_0^{1}dx\left[K^2(x)K(x^{\prime })+K^3(x^{\prime })\right]=1.29401145979\dots \end{array}$$

$$C_2(\Gamma )\simeq \frac{8}{\Gamma }$$

$$\underset{\Gamma \to \infty }{\lim }\Gamma C_d(\Gamma )=B_d=4d$$

$${\sigma }_{i,j}(t)=⟨c_i^{†}{c}_j{⟩}_t=\mathrm{Tr}\left(\rho (t)c_i^{†}{c}_j\right)$$

$$\frac{d{\sigma }_{i,j}}{dt}$$

$${\sigma }_{i,j}(t)=S_i(t)S_j^{\ast }(t)$$

$$\frac{d{S}_n}{dt}=\mathrm{i}[S_{n+1}+S_{n-1}]-\Gamma {\delta }_{n,0}{S}_n$$

$$S_n(0)={\delta }_{n,0}$$

$${\hat{S}}_n(s)={\int }_0^{\infty }dt{e}^{-st}{S}_n(t),$$

$$s\widehat{S}_{n}-\delta_{n,0}={\rm i}\big{[}\widehat{S}_{n+1}+\widehat{S}_{n-1}\big{]}-\Gamma\delta_{n,0}\widehat{S}_{n}$$

$$S(s,q)=\sum _{n=-\infty }^{\infty }{\hat{S}}_n(s)e^{-\mathrm{i}qn},$$

$$sS-1=2\mathrm{i}S\cos q-\Gamma {\hat{S}}_0$$

$$S=\left[1-\Gamma {\hat{S}}_0(s)\right]\Psi (s,q),\Psi (s,q)=\frac{1}{s-2\mathrm{i}\cos q}$$

$${\hat{S}}_0(s)=\frac{1}{2\pi }{\int }_0^{2\pi }dqS(s,q)$$

$${\hat{S}}_0(s)=\frac{1}{\Gamma +\sqrt{s^2+4}}$$

$$S(s,q)=\frac{\sqrt{s^2+4}}{\Gamma +\sqrt{s^2+4}}\Psi (s,q)=\frac{1}{\Gamma +\sqrt{s^2+4}}\frac{\sqrt{s^2+4}}{s-2\mathrm{i}\cos q}$$

Full text

1 Introduction

Quantum phenomena are fragile and susceptible to decoherence by the interplay of an open quantum system with its environment through interactions and measurements. An environment that is ‘constantly watching’ [1] generates a competition between Hamiltonian dynamics and dissipation [2] that often leads to interesting non-equilibrium processes. Hence a suitably engineered coupling to a reservoir may be a resource to generate a non-equilibrium state that can be investigated experimentally and analytically [3, 4].

Classical interacting particle processes are an ideal laboratory to explore the differences between thermodynamics and non-equilibrium phenomena. In this respect, simple models retaining non-trivial physical features and amenable to analytical treatment play an important role in understanding of time-reversal breaking, rare events, large deviations and dynamical phase transitions [5, 6, 7]. Coupling to reservoirs represented by Markov processes have found numerous applications to non-equilibrium dynamics ranging from biophysical models to urban traffic flows [8]. One such paradigmatic model is an exclusion process—a lattice gas of random walkers subject to a hard-core constraint that prevents two walkers from being simultaneously at the same site [9]. This classical exclusion process has been instrumental to exploring non-equilibrium states generated by various means such as unbalanced reservoirs, internal driving field or relaxation from a specific initial condition. For example, the study of the exclusion process on an infinite lattice has led to a better understanding of current statistics, height distribution in the Kardar-Parisi-Zhang equation and their relations with Tracy-Widom distributions from random matrix theory (see [10, 11, 12] and references therein). Another example is the symmetric exclusion process driven by a localized source [13, 14]. The statistics of the total number $\mathcal{N}(t)$ of particles entering an initially empty system displays a strong dependence on dimension [13]. For instance, the average number of particles $N(t)=\langle\mathcal{N}(t)\rangle$ grows according to444We tacitly assume that the underlying lattice is hyper-cubic if not stated otherwise, so the Watson integral (2) corresponds to the lattice ${\mathbb{Z}}^{d}$ with $d\geq 3$. Equation (1) accounts for simple exclusion processes with total hopping rate equal to 2 for all $d$.

[TABLE]

where $\Gamma$ is the injection rate and $W_{d}$ the Watson integral defined via

[TABLE]

Interestingly, in low dimensions, viz. for $d=1$ and $d=2$, the leading behavior in (1) is independent of the injection rate.

In the present work we study a similar situation in the quantum setting, namely we analyze an open quantum system of fermions driven by a source that injects fermions at the origin. We consider the simplest non-interacting (spinless) fermions; the Pauli exclusion makes this quantum system somewhat analogous to the exclusion process.555This analogy is intuitively appealing, but somewhat misleading: the one-dimensional symmetric exclusion process can be mapped exactly to the Heisenberg XXX spin-chain in imaginary time and not to free fermions (see also [15]). Our goal is to study how the number $N$ of fermions in the system increases with time, as a function of the injection rate $\Gamma$ and the dimensionality $d$ of the lattice.666A closely related quantum process defined in continuous space was analyzed in [16]. The growth of the number of fermions differs significantly from the corresponding behavior of the classical exclusion process (1). Some qualitative behaviors are simpler and more intuitive in the quantum case. In particular, the average number of fermions increases linearly with time, $N\simeq C_{d}(\Gamma)t$, regardless of the spatial dimension $d$. This is easy to appreciate—quantum walks are ballistic in any dimension [17], and so the effective jamming of the lattice in the vicinity of the origin by classical particles that leads to the sub-linear growth when $d\leq 2$, see (1), is avoided. The injection rate affects the asymptotic behavior in all $d$ and the particle growth rate has an expected asymptotic $C_{d}(\Gamma)\sim\Gamma$ when $\Gamma\to 0$. A counter-intuitive behavior occurs in the opposite limit, $\Gamma\to\infty$, where $C_{d}(\Gamma)\sim 1/\Gamma$. This is a signature of the quantum Zeno effect.

The outline of this work is as follows. In the next section  2, we define the model using the Lindblad equation formalism and present a few chief results. In section 3, we calculate the growth rate and the density profile for the one-dimensional lattice. In section 4, a general approach valid for an arbitrary dimension is presented and explicit calculations are performed for the square lattice. Various extensions are discussed in section 5. Details of the calculations and an overview of the lattice Green functions are presented in the appendices.

2 The model and some exact results

We study an open quantum system on an infinite lattice driven by a localized source of fermions. The system is initially empty; at time $t=0$, the source is turned on and fermions are injected at the origin at a constant rate.

The evolution of open quantum systems is described by the Lindblad equation [18, 19, 20, 21, 22]

[TABLE]

for the density matrix $\rho(t)$. The set of operators $L$ describe the interactions with reservoirs. In our situation when the fermions are injected at the single site which we set to be the origin, we have a single operator $L=\sqrt{\Gamma}\,c_{\bf 0}^{\dagger}$ modeling the insertion with intensity $\Gamma$. The curly brackets in equation (3) denote the anti-commutator, ${\rm i}=\sqrt{-1}$ and we set $\hbar=1$. We consider identical non-interacting spinless lattice fermions described by the Hamiltonian $H=\sum_{\langle{\bf i},{\bf j}\rangle}(c_{\bf i}^{\dagger}c_{\bf j}+c_{\bf j}^{\dagger}c_{\bf i})$, where the sum is taken over pairs of neighboring lattice sites ${\bf i}$ and ${\bf j}$. We consider uniform systems, so the hopping rates are equal and set to unity. Fermions satisfy the Pauli exclusion principle: $\{c_{{\bf i}}^{\dagger},c_{{\bf j}}\}=\delta_{ij}$ and all the other anticommutators vanish. In one dimension, the Hamiltonian reads

[TABLE]

A basic characteristic of this open quantum system is the total number of particles $\mathcal{N}(t)$ at time $t$. The full statistics of this evolving random quantity is described by the probability distribution $P(N,t)=\text{Prob}[\mathcal{N}(t)=N]$. We focus on the average $N(t)=\langle\mathcal{N}(t)\rangle$ and show that it exhibits a linear growth

[TABLE]

when $t\gg 1$. The particle growth rate $C_{d}(\Gamma)$ depends on the strength $\Gamma$ of the source and on the spatial dimension $d$ (and also on the geometry of the underlying lattice).

It is possible to obtain a formal expression for the growth rate in any dimension in terms of lattice Green functions (see Appendix B), and we succeeded in deriving explicit tractable formulae for $d=1$ and $d=2$ (which are experimentally most relevant). We now state some of the results derived in the following sections.

In one dimension (section 3), the fermion growth rate $C_{1}(\Gamma)$ is given by

[TABLE]

with

[TABLE]

Note that equation (6) is formally valid when $\Gamma>2$ and $\gamma$ is real and positive. For $\Gamma<2$, the result remains valid after analytical continuation.

From equation (6), one can derive extremal behaviors. When $\Gamma\ll 1$, we have

[TABLE]

The particle growth rate $C_{1}(\Gamma)$ increases with $\Gamma$ for sufficiently small insertion rates and then decreases asymptotically approaching to zero when $\Gamma\to\infty$, see Fig. 1 (a). The vanishing of the particle growth rate when $\Gamma\to\infty$ is the manifestation of the quantum Zeno effect (QZE), one of the simplest and most paradigmatic differences between classical and quantum open systems. The QZE refers to the counter-intuitive phenomenon that a frequent measurement of a quantum system does not allow the state to evolve unitarily [23, 24, 25]. For weak measurements, the system is only slightly perturbed by the measurement back-action. When the measurements are strong, the system gets projected to an eigenstate of the measurement operator (the wavefunction collapses) and if one keeps on measuring, the system is unable to evolve. Fundamentally, there is little difference between continuous measurements and a dissipative process as the former entails connecting the system to a large reservoir (the measurement device). The only real difference is that measurements are typically recorded. It is thus quite natural that the dynamics of simple lattice problems with local injection or removal of particles exhibits the QZE.

The QZE has been demonstrated in various experiments with ultra-cold atoms, see e.g. [26, 27, 28]. Some recent experiments [29, 30, 31] are particularly close to our model. The experimental setup [29, 30, 31] involves a local loss process induced by shining a focused electron beam onto an atomic Bose-Einstein condensate residing on the one-dimensional optical lattice. The QZE manifests itself in the number of atoms lost from the condensate: linear growth was found for small $\Gamma$ and an inverse $1/\Gamma$ scaling was observed when $\Gamma$ is large. A recent theoretical analysis [32] of this one-dimensional system with local losses combines the Lindblad framework with continuum modeling based on the Luttinger liquid description and emphasizes the similarity with the Kane-Fisher barrier problem [33, 34].

In higher dimensions, exact results for the Laplace transforms are given in Section 4. For the square lattice

[TABLE]

where $I_{1}$ and $I_{2}$ are the integrals

[TABLE]

and $K(x)$ is the complete elliptic integral of the first kind of modulus $x$:

[TABLE]

In (10a) we have used the standard notation $x^{\prime}=\sqrt{1-x^{2}}$ to denote the complementary modulus.

The dependence of the particle growth rate $C_{2}(\Gamma)$ on the insertion rate $\Gamma$ is qualitatively the same as in one dimension, see Fig. 1 (b). When $\Gamma\ll 1$, the particle growth rate can be written as

[TABLE]

with

[TABLE]

The asymptotic expansion (12) is obtained by expanding (10a)–(10b) when $\Gamma\ll 1$ and it is in excellent agreement with the exact prediction for sufficiently small $\Gamma$. When $\Gamma\to\infty$, the leading behavior is given by

[TABLE]

The universal decay, $C_{d}\simeq B_{d}\Gamma^{-1}$ as $\Gamma\gg 1$, occurs in all spatial dimensions and for hyper-cubic lattices the amplitude has a simple value

[TABLE]

The universal $\Gamma^{-1}$ decay in all dimensions is a signature of the QZE, which is a local dynamical effect, due to the measurement, and is not related to the global geometry of the system under study. Note that the dimension $d=2$ which separates transient and recurrent behavior for classical random walks plays no special role in the quantum case.

3 Free fermions on $\mathbb{Z}$ with a source

The basic quantities that we wish to determine are the average total number of particles as a function of time and the particle density profile. Both observables are encoded in the two-point correlation functions, defined as

[TABLE]

Using the Lindblad equation (3) with Hamiltonian (4) and $L=\sqrt{\Gamma}\,c_{0}^{\dagger}$ we deduce that the two-point correlations (16) satisfy a closed system777A similar property holds for the symmetric exclusion process [12] and for some of its quantum generalizations [15]. of differential equations

[TABLE]

These equations are linear and exactly solvable.

It is more convenient to deal with equations

[TABLE]

which are identical to (17) apart from the source term, the last term on the right-hand side of (17), which is absent in (18). Interestingly, (18) describe the dual problem in which particles are removed from the origin at rate $\Gamma$; this process is modeled by the Lindblad equation (3) with the same Hamiltonian (4) and $L=\sqrt{\Gamma}\,c_{0}$. In order to calculate the growth rate we shall analyze this dual problem starting with a single particle at the origin.

A crucial observation simplifying the analysis of (18) is that the solution of (18) factorizes, that is the ansatz

[TABLE]

is consistent with (18) if the functions $S_{n}(t)$ satisfy

[TABLE]

with initial condition

[TABLE]

The initial-value problem (20)–(21) represents a continuous-time quantum walk, namely a quantum particle propagating on a chain, in the presence of a trap of strength $\Gamma$ at the origin. The survival of a quantum particle subject to a complex optical potential [35, 36, 37, 38] was analyzed in [39]. The factorization (19) shows that, in the present case, this non-Hermitian potential is equivalent to the more rigorous framework provided by the Lindblad equation that also includes the measurement process (see [40, 41, 42, 43, 44, 45] for related approaches in the context of first detection time of a quantum walk).

3.1 Solution of the one-body problem

Equation (20) has been analyzed in the past, see [39] and references therein. The approach is based on the Laplace-Fourier transform. Performing the Laplace transform with respect to time,

[TABLE]

we recast (20) into

[TABLE]

Performing the Fourier transform in space,

[TABLE]

we rewrite equation (23) as

[TABLE]

yielding

[TABLE]

The definition (24) implies

[TABLE]

which in conjunction with (26) fixes

[TABLE]

and

[TABLE]

Inverting the Fourier transform, we deduce

[TABLE]

The inverse Laplace transform of this expression determines the fermion density $S_{n}(t)$ at a given time and at position $n$ (see Appendix A).

3.2 Calculation of the growth rate

We now calculate the growth rate $C_{1}(\Gamma)$. The average total number of particles can be expressed in terms of the correlation functions

[TABLE]

where the $\sigma_{i,j}(t)$’s satisfy the differential equations (17) with initial condition $\sigma_{i,j}(0)=0$ for all $i,j$. Note that the linear system (17) can formally be rewritten as

[TABLE]

Here $\sigma=\sigma(t)$ is an infinite vector with components $\sigma_{i,j}(t)$; $M$ is the matrix that encodes the homogeneous terms888The same $M$ arises in the dual problem of removing a particle from the origin.; and $V$, the source term in (17), is a vector with a unique non-zero component, $2\Gamma\delta_{i,0}\delta_{j,0}$. Equation (31) is solved to yield

[TABLE]

The vector ${\rm e}^{(t-\tau)M}V$ represents two-point correlations in the dual problem (18) with a fermion starting at the origin. Hence, using the ansatz (19), we deduce

[TABLE]

where the functions $S_{i}$ satisfy the dynamics (20). Combining (30) and (33) we obtain

[TABLE]

where $\Pi(t)$ represents the survival probability of a single fermion starting at the origin, that can be removed at the origin with rate $\Gamma$. This survival probability

[TABLE]

approaches, when $t\to\infty$, to

[TABLE]

Using the Parseval-Plancherel identity

[TABLE]

for the Fourier transforms and setting $f(t)=S_{0}(t)$ and $g(t)=S_{0}^{*}(t)$, we express the integral in (36) through the Laplace transform $\widehat{S}_{0}(s)$. This approach yields [39]

[TABLE]

When $\Gamma>2$, we compute the integrals and arrive at

[TABLE]

with $\gamma\equiv\sqrt{(\Gamma/2)^{2}-1}$; for $\Gamma<2$, the survival probability can then be obtained by analytical continuation of (39).

Thus, we conclude that the average number of fermions grows linearly in time and

[TABLE]

leading to the announced result (6) for the particle growth rate.

3.3 The density profile

We wish to determine the density profile of the particles at time $t$ starting with an initially empty system. It is equivalent to studying the dual problem of a system initially full of holes, where holes are removed at the origin (and replaced by particles). In this formulation, the governing equations are given by (18) and the initial conditions become

[TABLE]

The answer to this problem can be written as a linear combination

[TABLE]

where $S_{n}^{(k)}(t)$ solves the one-body non-Hermitian dynamics (20) with initial condition

[TABLE]

This one-body problem is solved by the Laplace transform and we obtain

[TABLE]

Performing the inverse Laplace transform yields

[TABLE]

where $J_{|n-k|}(2t)$ is a Bessel function and $\Gamma{{\rm i}}^{|n|+|k|}G_{|n|+|k|}(t)$ denotes the inverse Laplace transform of the expression in the second line on the right-hand side of Eq. (44).

We conclude that the density of particles at site $n$ is given by

[TABLE]

with $S_{n}^{(k)}(t)$ given by (45). To extract the long time behavior, we re-parametrize $k$ and $n$:

[TABLE]

Due to symmetry, it suffices to consider the $n\geq 0$ range, so $v\geq 0$. In the $t\to\infty$ limit, the density profile is obtained by applying the saddle-point method for the functions $J_{n}$ and $G_{m}$. The asymptotic formulae read

[TABLE]

Equation (48) is standard, see [46]; equation (49) is derived in a similar fashion. The asymptotic formulae (48)–(49) are valid in the region $-1\leq u-v\leq 1$ and $-1\leq|u|+v\leq 1$. Recalling that $v\geq 0$ one gets $|u|\leq 1-v$.

In the long time limit, the profile is obtained by substituting (48) and (49) into (45) and (46). After converting the discrete sum over the integer $k$ into integrals with respect to the continuous variable $v$ we find the density at the position $n$ to be

[TABLE]

The profile ${\mathcal{N}}(v)$ describes the density at the bulk where $v=\frac{n}{2t}>0$ and it has the form

[TABLE]

Near the origin there is a spike, namely an inner layer of width $n=O(1)$ where the density profile deviates from the continuous density profile ${\mathcal{N}}(v)$ that varies on the scale $n=O(t)$. The deviation is given by

[TABLE]

This correction decays rapidly when $|n|\gg 1$. Using integration by parts one extracts the asymptotic approach of the spike near the origin to the bulk density:

[TABLE]

Figure 2 shows a comparison of the asymptotic bulk density profile $\mathcal{N}(v)$ to exact numerical results at finite time. Note the excellent agreement in the bulk. The spike near the origin described by $N_{n}(t)\simeq{\mathcal{N}}(0)+I_{n}$ with $I_{n}$ given by (52) is also in good agreement with numerical results for finite time. There are also oscillations near the causal horizon caused by the dispersion of the particles. These oscillations are not accounted by our asymptotic formulas, but they can be deduced from the exact results (44)–(46) using a more accurate treatment in the region $|n|-2t=O(t^{1/3})$. This $t^{1/3}$ scaling of the width of the front region and its staircase like internal structure arises in many problems involving free fermions, see e.g. [47, 48].

Setting $n=0$ in equations (50)–(52), we find the density at the origin:

[TABLE]

Computing the integral gives

[TABLE]

with $\gamma=\sqrt{(\Gamma/2)^{2}-1}$ when $\Gamma>2$; when $\Gamma<2$, the result follows from (55) by analytical continuation. As expected, $1-N_{0}=\Pi_{\infty}$ with $\Pi_{\infty}$ given in (39). As another consistency check we compute the total average number of fermions

[TABLE]

and recover the expected answer.

4 Free fermions on ${\mathbb{Z}}^{d}$ with a source

An open quantum system of free fermions performing independent continuous time quantum walks and coupled to a localized source can be studied on an arbitrary graph. Explicit calculations are possible for lattices admitting sufficiently simple lattice Green functions. Thus, after describing the general formalism, we focus on the square lattice.

4.1 General Formalism

In $d$ dimensions, we again solve equation (20) by performing Laplace and Fourier transforms. We first define

[TABLE]

where ${\bf n}=(n_{1},\ldots,n_{d})$, ${\bf q}=(q_{1},\ldots,q_{d})$ and ${\bf q}\cdot{\bf n}=q_{1}n_{1}+\ldots+q_{d}n_{d}$. To avoid cluttered notation, we write

[TABLE]

We also define

[TABLE]

The function $\Phi_{d}(s)$ is given through a $d-$fold integral which can be reduced to a single integral involving a Bessel function. Indeed, writing $\Psi_{d}(s)$ as

[TABLE]

and plugging (78) into the definition of $\Phi_{d}(s)$ we get

[TABLE]

Recalling an integral representation of the Bessel function

[TABLE]

we express $\Phi_{d}(s)$ through a single integral:

[TABLE]

Applying the Laplace-Fourier transform to the dual problem we arrive at

[TABLE]

From the definition (56b), we obtain

[TABLE]

Using (61) and (62) we fix the value of $\widehat{S}_{\bf 0}(s)$,

[TABLE]

with $\Phi_{d}(s)$ defined in (57b). Therefore (61) becomes

[TABLE]

The general relations (36) and (40) are valid in arbitrary dimension; thus the growth rate is formally given by

[TABLE]

with

[TABLE]

This integral can be expressed through the Laplace transform $\widehat{S}_{\bf 0}(s)$, given in (63), via the Parseval-Plancherel identity (37). Hence, to have more manageable results, one needs to compute $\Phi_{d}(s)$ explicitly to obtain an integral representation comparable to (38). The ‘euclidean’ versions of integrals (57a)–(57b) arise in numerous problems and are known as lattice Green functions (LGF). The LGF are reviewed in Appendix B where we present explicit results for $\mathbb{Z}^{d}$ with $d=2$ and $d=3$, and briefly discuss other lattices. The LGF is rather complicated even for the cubic lattice; therefore, in the following we focus on the square lattice.

4.2 Growth rate for the process on square lattice

For the square lattice, the integral in (60) can be expressed through the complete elliptic integral of the first kind. One gets (see [49] and Appendix B)

[TABLE]

leading to

[TABLE]

Using again the Parseval-Plancherel identity (37), we can write

[TABLE]

The last expression can be viewed as a complex integral along vertical contour, infinitely close to the imaginary axis. The function $K(z)$ is an even function, analytic in the complex plane with branch cuts along the real axis from $-\infty$ to $-1$ and from $1$ to $+\infty$. In the vicinity of the cuts, for $z=s\pm{\rm i}\epsilon$ with $\Re(s)\geq 1$ and $\epsilon\to 0^{+}$, we have the following relation [46, 50]

[TABLE]

(see e.g., formula 162.02, page 39 in [50]). We cut the integral appearing in (69) into four parts: $y\in]-\infty,-4],[-4,0],[0,4]$ and $[4,\infty[$. Inside the middle ranges, the function is analytic. In the boundary integrals, we use (70)999The correct sign is selected depending on which side of the cut the integrand lies.. This leads to

[TABLE]

This integral representation leads to the announced results (9), (10a) and (10b).

To derive the behavior for $\Gamma\to\infty$, we start from the asymptotic

[TABLE]

Keeping the leading term in (72) we recast $I_{2}(\Gamma)$ into

[TABLE]

Subtracting (73) from (10b) yields

[TABLE]

Combining this with (73) and keeping terms up to $O(\Gamma^{-3})$, we obtain

[TABLE]

Using (10a) we obtain a similar expansion

[TABLE]

Combining (74) and (75) we arrive at

[TABLE]

with

[TABLE]

A numerical estimate shows

[TABLE]

with very high precision. This is probably an identity which would be interesting to prove. If (78) is true, we get $B_{2}=8$, leading to the announced result (14). For completeness, we show the $2D$ density profile in Fig. 3. Like in one dimension it has a pronounced peak near the injection site. Moreover, in contrast to the analogous classical problem, it is very anisotropic. It clearly displays interference of particles along high symmetry lines in the lattice. In principle the 1D analyses from section 3.3 can be extended to 2D but analytic extraction of the asymptotics is more challenging and we leave it as an open problem.

5 Discussion

Theoretical studies of quantum gases with sinks, sources, and other dissipative mechanisms have been mostly carried out using one-dimensional models. Attractive features of one-dimensional systems include integrability of some toy models, and realizability in experiments with cold atoms. In higher dimensions, approximations and simulations are required. For the simple quantum process that we studied, an open system of free fermions with injection from a reservoir, we showed how to calculate the particle growth rate for simple lattices in any dimension. Analytical formulae have been derived for the one-dimensional lattice and for the square lattice. With more effort, explicit results for the particle growth rate could be derived for the triangular lattice, the cubic lattice, and perhaps even for the body-centered hyper-cubic lattices as for such lattices the lattice Green function can be expressed through hypergeometric functions. More detailed characteristics can also be studied analytically, e.g., we have calculated the density profile in one dimension. While the latter depend on the exact dispersion relations, linear growth of the average particle number is to be expected on any lattice for free fermions. The linear growth is simply induced by the fact that particles spread ballistically on any local translational invariant lattice Hamiltonian.

In addition to the average total number of injected particles, one could explore higher cumulants and study the entire distribution function $P(N,t)$, or the full counting statistics. The large deviation (LD) formalism developed for classical systems [51, 6, 52, 53] has been recently applied to simple quantum systems such as single and two-qubit systems [54, 55, 56, 57]. Large open systems of non-interacting fermions are also tractable in stationary one-dimensional settings, for example non-equilibrium steady states of boundary driven systems have been established [58, 59], and the large-deviation function has been subsequently calculated [60]. The stationarity is crucial in the procedure [61] employed in the calculation of the large-deviation function [60]. For other applications of the LD formalism to boundary-driven non-interacting quantum chains, see Refs. [62, 63, 64, 65, 66, 67, 68].

Our system is manifestly non-stationary as the number of particles increases with time. In some evolving free fermion systems, fluctuations have been probed analytically; for instance, in a system of lattice fermions starting from the domain wall initial condition, the number of fermions in the initially empty half-line grows linearly in time, while the variance increases logarithmically [69]. Recently, the full statistics of the current in this problem has been calculated using powerful techniques inspired by previous studies of classical stochastic systems [70]. We believe that system could be also explored using integrable probability methods.

We totally ignored interactions. A natural generalization of the present model is the XXX chain with all spins originally pointing down in which the spin at the origin can be flipped upwards with a certain rate. The equations for the correlators in this model are hierarchical rather than recurrent and the factorization property that played a key role in the present study also does not hold. Various matrix representation techniques have been applied to open interacting finite spin chains coupled to reservoirs [71, 72, 73] and it would be interesting to extend these results to the infinite chain case (see [74] for a review and references).

There are many other open quantum systems with local sinks or source which can be potentially treated. The most obvious extension is to bosons, see [16, 75, 76] for work in this direction. For example, in contrast to fermions, a phase transition was reported in [75] in the case of weakly interation bosons with a local sink.

Acknowledgments. We are grateful to Jean-Marc Luck, Michel Bauer and Tomaž Prosen for stimulating discussions, and to S. Mallick for a careful reading of the manuscript. PLK thanks the Institut de Physique Théorique for hospitality. DS acknowledges support from the FWO as post-doctoral fellow of the Research Foundation – Flanders.

Appendix A Single particle case

In this appendix we determine the density profile in the dual problem starting with a single particle at the origin; we consider the most symmetric setting when the sink is also at the origin. The inverse Laplace transform of (27) admits an integral representation [77]

[TABLE]

In the long time limit, the integral term dominates. Computing its asymptotic behavior we find

[TABLE]

The density at the origin

[TABLE]

oscillates indefinitely reflecting quantum interference. Further, the density at the origin vanishes infinitely often: $N_{0}(t_{k})=0$ at times $0<t_{1}<t_{2}<t_{3}<\ldots$ significantly depending on $\Gamma$ for $k=O(1)$ and behaving universally, $t_{k}\simeq\tfrac{\pi}{8}+\tfrac{\pi}{2}k$, for $k\gg 1$.

After averaging (81) over oscillations we arrive at

[TABLE]

When $n>0$, the inverse Laplace transform of (29) admits an integral representation

[TABLE]

In the scaling limit

[TABLE]

we have

[TABLE]

and therefore (83) yields

[TABLE]

Using the asymptotic of $J_{n}(2t)$ in the scaling limit (84) we find that $N_{n}=|S_{n}|^{2}=0$ when $|v|>1$ and

[TABLE]

when $|v|<1$. Equation (86b) has been obtained by averaging over oscillations. The scaled density profiles (86b) are shown in Fig. 4 for three values of the strength $\Gamma$ of the trap. As a useful consistency check, we recover the value (39) of the asymptotic survival probability

[TABLE]

Appendix B Lattice Green functions

One can relate the Euclidean version of the integrals (57a)–(57b) to formulae arising in the context of lattice random walks [49]. Let $P_{k}$ be the probability that a random walker starting at the origin returns back after $k$ steps. The generating function $P(z)=\sum_{k\geq 0}P_{k}z^{k}$, known as the lattice Green function (LGF), admits an integral representation

[TABLE]

with

[TABLE]

for the hyper-cubic lattice. The general formula (87) is valid for many other lattices, the specificity of the lattice is reflected in the function $\Lambda({\bf q})$. For instance,

[TABLE]

for the triangular lattice, the body-centered cubic (bcc) lattice and the face-centered cubic (fcc) lattice. The LGFs for many low-dimensional lattices are derived in a book by Hughes [49]; for more recent results, see [78, 79, 80, 81] and references therein. Here we outline some remarkable special cases.

In two dimensions, for the square lattice, the LGF admits a remarkably simple expression through the complete elliptic integral (or as an hypergeometric function)

[TABLE]

We have used this expression to derive (9) and (10a)–(10b). For a two-dimensional triangular lattice, the LGF has a more complicated form

[TABLE]

with $k=\sqrt{\frac{2(\zeta_{+}-\zeta_{-})}{c}}$ and $\zeta_{\pm}=\frac{3}{z}+1\pm\sqrt{3+\frac{6}{z}}\,.$

In three dimensions, for the body-centered cubic lattice, a neat expression for the LGF has been found long ago [82]

[TABLE]

For the three-dimensional cubic lattice, the LGF reads [79]

[TABLE]

with $k_{1}^{2}=\frac{16\xi^{3}}{(1-\xi)^{3}(1+3\xi)}$ and $\xi^{2}=\frac{1-\sqrt{1-z^{2}/9}}{1+\sqrt{1-z^{2}}}\,.$

In four and higher dimensions, the LGF of the body-centered hyper-cubic lattices can be expressed [81] through the hypergeometric function in arbitrary dimension

[TABLE]

In contrast, the calculations of the LGF of the hyper-cubic lattices require rather abstract mathematics such as Calabi-Yau ordinary differential equations [78, 80, 81].