Non-equilibrium quantum thermodynamics in Coulomb crystals

TL;DR

This paper investigates the non-equilibrium work statistics during sudden parameter changes in ion Coulomb crystals near a structural phase transition, revealing divergence and critical behavior linked to quantum Ising models.

Contribution

It provides an analytical and numerical analysis of irreversible work in Coulomb crystals near the linear-zigzag transition, connecting it to quantum critical phenomena.

Findings

Irreversible work diverges near the transition

Work fluctuations exhibit non-analytic behavior

Results are experimentally verifiable with current ion trap technology

Abstract

We present an in-depth study of the non-equilibrium statistics of the irreversible work produced during sudden quenches in proximity to the structural linear-zigzag transition of ion Coulomb crystals in 1+1 dimensions. By employing both an analytical approach based on a harmonic expansion and numerical simulations, we show the divergence of the average irreversible work in proximity to the transition. We show that the non-analytic behaviour of the work fluctuations can be characterized in terms of the critical exponents of the quantum Ising chain. Due to the technological advancements in trapped ion experiments, our results can be readily verified.

Equations

$$H_0=$$

$$\frac{Q^2}{8\pi {ϵ}_0}\sum _{i\ne j}{\left[(X_i-X_j{)}^2+(Y_i-Y_j{)}^2\right]}^{-1/2}$$

$$H(\omega )=\frac{1}{2}\sum _{j=1}^L\left[-g^2\frac{{\partial }^2}{\partial y_j^2}+({\omega }^2-h_1)y_j^2+h_2(y_j+y_{j+1}{)}^2+h_3{y}_j^4\right],$$

$$P_F(W)\equiv \sum _{n,\bar{m}}{p}_n^0{p}_{\bar{m}|n}^{t_f}\delta [W-({ϵ}_{\bar{m}}-{ϵ}_n)],$$

$${\omega }_k^2=g^2\left[{\omega }^2-h_1+4h_2{\cos }^2\frac{k}{2}\right]$$

$${\chi }_F(t)\equiv \int dW{e}^{iWt}{P}_F(W)=\mathrm{Tr}[e^{i{H}_{f}t}U(t_f,0{)}^{†}{e}^{-i{H}_{i}t}U(t_f,0){\rho }_i]$$

$$\displaystyle W_{\text{IRR}}=\sum_{k}\Bigg{[}\frac{1}{2}\left(\Omega_{k}\omega_{k}^{f}-\omega_{k}^{i}\right)\coth\frac{\beta\omega_{k}^{i}}{2}-\frac{1}{\beta}\ln\frac{\sinh{(\frac{\beta\omega_{k}^{f}}{2})}}{\sinh(\frac{\beta\omega_{k}^{i}}{2})}\Bigg{]},$$

$${\sigma }_W^2=\sum _k\frac{{{\omega }_k^f}^2\cosh (\beta {\omega }_k^i)\left({\Omega }_k^2-1\right)+{\left({\omega }_k^f{\Omega }_k^2-{\omega }_k^i\right)}^2}{4{\sinh }^2(\beta {\omega }_k^i/2)},$$

$$W_{\text{IRR}}^{\mathrm{soft}}$$

$${\sigma }_{\mathrm{soft}}^2$$

$$W_{\text{IRR}}(T=0)$$

$$⟨W⟩$$

$$1-e^{\frac{W_{\text{IRR}}-W_{\text{IRR}}^{(\max )}}{L}}=f\left[L^{\frac{1}{\nu }}({\omega }^2-{\omega }_{\star }^2)\right]$$

Full text

I Introduction

The recently renovated interest in non-equilibrium thermodynamics of quantum systems, spurred from tremendous advances in experimental techniques, has found a plethora of interesting developments and applications Campisi et al. (2011); Esposito et al. (2009); Goold et al. (2016); Vinjanampathy and Anders (2016); Millen and Xuereb (2016). From a theoretical perspective, several important achievements are already available in the literature. Prominent examples are the quantum generalisation of fluctuation relations such as the celebrated Jarzynski equality Jarzynski (1997); Tasaki (2000); Crooks (1999); Talkner and Hänggi (2007); Talkner et al. (2007) and the design of a single-atom thermal machine Abah et al. (2012), recently realised with a trapped ion Roßnagel et al. (2016).

In recent years, much interest has been devoted to the analysis of the quantum work extracted from, or absorbed by, a quantum system. While many definitions of work in a quantum setting have been proposed, the most popular one in the literature, based on two-time measurements Talkner et al. (2007); Campisi et al. (2011), fulfils the Jarzynski equality but not the first law of thermodynamics when the system exhibits energy coherences Perarnau-Llobet et al. (2017). Furthermore, several proposals have been put forward to estimate work in quantum systems without the need of realising energy projections. These include schemes based on Ramsey interferometry Mazzola et al. (2013); Dorner et al. (2013); Campisi et al. (2013) experimentally realised in an NMR setting Batalhão et al. (2014). Other methods employ the aid of auxiliary continuous quantum systems Roncaglia et al. (2014); De Chiara et al. (2015); Villa and De Chiara (2017). Some others are based on population imbalance and coherence in bosonic Josephson junctions Lena et al. (2016).

In the many-body scenario, fluctuations of the work and of its irreversible contribution have been calculated mainly for spin chains in proximity to a quantum phase transition Silva (2008); Dorner et al. (2012); Fusco et al. (2014); Bayat et al. (2016); Paganelli and Apollaro (2017). In this context it has been shown that all the moments of the work probability distribution are singular when the system is dynamically driven close to the phase transition Bayat et al. (2016); Paganelli and Apollaro (2017).

In this context, most studies so far have been limited to archetypal examples of strongly correlated systems in condensed matter Mascarenhas et al. (2014); Sindona et al. (2014); Sotiriadis et al. (2013); Shchadilova et al. (2014); Russomanno et al. (2015); Marino and Silva (2014); Bayocboc and Paraan (2015); Zhong and Tong (2015); Beau et al. (2016). However, testing such predictions on experimental platforms poses additional challenges. In fact, the full probability distribution seems to be unaccessible to observation because of the complexity of energy projection and because Ramsey schemes would involve, unrealistically, an ancilla coupled to the whole many-body system.

In this work, we analyse the out-of-equilibrium quantum thermodynamics of a model specifically tailored to an experimental setup. Precisely, we consider ion Coulomb crystals (ICC): Many-body quantum systems of cold atomic ions confined to highly anisotropic traps and mutually interacting via Coulomb repulsion Dubin and O’Neil (1999). Here we estimate the statistics of the irreversible work production during sudden quenches in proximity to the phase transitions in such ICC. Our results can be experimentally tested in current experiments by measuring, as we show in this work, the transverse displacement distribution of the ions positions.

At equilibrium, ICC exhibit different structural arrangements, generally depending on the spatial properties of the trapping potential Dubin (1993). The structural transitions occurring between these configurations are typically phase transitions of the first order Piacente et al. (2004). Here, we focus on quasi one-dimensional arrangements obtained for strongly anisotropic traps, which exhibit a transition from a linear to a zigzag configuration Fishman et al. (2008). In the limit of ultracold ions, quantum fluctuations become relevant and the linear-zigzag transition becomes a quantum phase transition of the Ising universality class in 1+1 dimensions, at the thermodynamical limit Shimshoni et al. (2011a, b); Cormick and Morigi (2012); Silvi et al. (2013); Fishman et al. (2015); Silvi et al. (2014); Podolsky et al. (2014). The production of defects during a quick change of the trap anisotropy has been studied theoretically del Campo et al. (2010); De Chiara et al. (2010); Nigmatullin et al. (2016); Silvi et al. (2016) and experimentally Ulm et al. (2013); Pyka et al. (2013).

With this setup in mind we analyse the fluctuations of the work performed upon the ICC system by changing suddenly the transverse confinement frequency near the linear-zigzag transition point. We compare analytical results from calculations based on a harmonic expansion with numerical calculations based on the density matrix renormalization group algorithm (DMRG) White (1992, 1993) in the matrix product state formalism Östlund and Rommer (1995); Schollwöck (2011). We show that when approaching the critical point such fluctuations display a singularity, and they exhibit a universal scaling compatible with the quantum Ising model.

II Linear-Zigzag model for ICC

Ions confined to a fully anisotropic 3D trap and interacting via repulsive Coulomb interaction undergo a structural phase transition in which their spatial geometry changes from a one-dimensional linear chain to a planar zigzag configuration Fishman et al. (2008). The control parameter of such a transition is the frequency $\omega$ of the transverse harmonic trapping. This transition is driven by a mechanical instability of the chain that is associated with a soft mode at the boundary of the Brillouin zone whose frequency vanishes at the critical transverse trapping frequency. In the following paragraph we briefly review the analytical approximations needed to recast the linear-zigzag transition into a simple short-range model. For convenience, we restrict the motion of the ions to the $XY-$plane in which the zigzag structure develops. Here $X$ is the direction parallel to the trap axis and $Y$ is the direction perpendicular to the trap axis that emerges as a result of spontaneous symmetry breaking or because of a small anisotropy in the transverse confinement.

The strong repulsion between the ions makes them practically distinguishable particles, which in turn allows us to write a Hamiltonian in a first quantization

[TABLE]

in which $Q$ is the ion charge, $M$ is the mass, $(X_{j},Y_{j})$ and $(P_{x,j},P_{y,j})$ are the position and momentum of the $j-$th ion, respectively, and $V_{L}(x)$ is the longitudinal component of the confining potential. The quantum nature of this model stems from the commutator $[X_{i},P_{x,j}]=[Y_{i},P_{y,j}]=i\hbar\delta_{i,j}$. As shown in Fishman et al. (2008); Shimshoni et al. (2011a), when the chain is sufficiently close to criticality, the longitudinal and transversal components of $H$ effectively decouple. One can therefore fix the average equilibrium positions of the ions along the longitudinal direction to $x_{j}=ja$, with $a$ the effective Wigner lattice spacing. Afterwards, the transverse dynamics Hamiltonian $H_{y}$ can be Taylor-expanded at fourth order in the displacements $y_{j}$. The resulting theory, still long-range, can then be recast into a short-range model through an expansion of the scattering matrix of the harmonic modes, at second order in $\delta k$ around the soft mode ($\delta k=k-\pi/a$) Silvi et al. (2014). This mapping effectively simplifies the Hamiltonian of Eq. (1) into

[TABLE]

in which $g=\sqrt{\hbar^{2}/Ma^{2}E_{0}}$ plays the role of an effective Planck constant measuring the impact of quantum fluctuations Silvi et al. (2013); Podolsky et al. (2014). Here all quantities have been expressed in dimensionless scales, according to: $H=H_{0}/E_{0}$ with $E_{0}=Q^{2}/(4\pi\epsilon_{0}a)$, $y_{j}=Y_{j}/a$, and $\omega=\omega_{0}/\sqrt{E_{0}/Ma^{2}}$. Finally, $h_{1}=7\zeta(3)/2,h_{2}=\ln 2$, and $h_{3}=93\zeta(5)/8$, with $\zeta$ being the Riemann function, are universal constants Silvi et al. (2013, 2014). The Hamiltonian of Eq. (2) captures accurately the dynamics of the linear-zigzag quantum phase transition, and the critical point, identified by transverse frequency $\omega=\omega_{C}(g)$, can be computed as a function of $g$, and it was estimated to scale as $\omega_{C}(g)\approx h_{1}-3h_{3}g|\ln g|/2\pi+\mathcal{O}(g)$ for small $g$ Podolsky et al. (2014). In what follows we are going to extensively study the non-equilibrium statistics of the irreversible work generated after sudden changes of the transverse frequency from $\omega_{i}=\omega$ to $\omega_{f}$ such that $|\omega_{i}^{2}-\omega_{f}^{2}|=\Delta\omega$. First, we are going to present analytical results obtained using an approximated harmonic version of Eq. (2) and then compare them with numerical simulations, based on the DMRG algorithm.

III Non-equilibrium quantum thermodynamics

In quantum mechanics work is not a quantum observable Talkner et al. (2007) but a generalised measurement Roncaglia et al. (2014). As such, it is strongly affected by quantum fluctuations arising in the measurement process. The key figure in this respect is the probability distribution of the work generated when the system is subject to a time-dependent Hamiltonian, but is otherwise isolated by sources of heat or dissipation.

In this paradigm the time-dependent Hamiltonian $H_{i}=H(\omega_{i})$ is controlled via the frequency $\omega$. In turn, $\omega$ is quenched in time according to a certain time-dependent protocol $\omega(t)$, within the time window $[t_{i},t_{f}]$ [and accordingly, $\omega_{i}=\omega(t_{i})$ and $\omega_{f}=\omega(t_{f})$]. At the beginning of the protocol the system is assumed at equilibrium in the Gibbs state $\rho_{i}=e^{-\beta H(\omega_{i})}/\mathcal{Z}_{i}$, where $\mathcal{Z}_{i}=\mathrm{Tr}[e^{-\beta H(\omega_{i})}]$ is the partition function, and the inverse temperature $\beta$ is also expressed in dimensionless units. For later convenience we also define the final equilibrium partition function $\mathcal{Z}_{f}=\mathrm{Tr}e^{-\beta H(\omega_{f})}$. Two sets of energy measurements are then performed, the first prior to the protocol and corresponding to the eigenstates of $H_{i}$, and the second one right after the protocol and corresponding to the eigenstates of $H_{f}$. One can define the work distribution performed during the $H_{i}\to H_{f}$ transformation as

[TABLE]

in which $\epsilon_{n}$ and $\epsilon_{\bar{m}}$ are the eigenvalues of the initial and final Hamiltonian respectively, $p^{0}_{n}=e^{-\beta\epsilon_{n}}/\mathcal{Z}$ is the initial probability distribution in the energy levels, and $p^{t_{f}}_{\bar{m}|n}=|\mathinner{\langle{\epsilon_{\bar{m}}}|}U\mathinner{|{\epsilon_{n}}\rangle}|^{2}$ is the transition probability for the system to evolve from the state $\mathinner{|{\epsilon_{n}}\rangle}$ to $\mathinner{|{\epsilon_{\bar{m}}}\rangle}$, after the time evolution $U=U(t_{i}\to t_{f})=\mathcal{T}\!\exp\int_{t_{i}}^{t_{f}}-iH(t)dt$.

For a sudden quench ($U=\mathbbm{1}$), the average work is simply $\langle W\rangle={\rm Tr}[\rho_{i}(H(\omega_{f})-H(\omega_{i}))]$ while the free energy difference is $\Delta F=-\beta^{-1}\ln(Z_{f}/Z_{i})$. Because of the relation $\langle W\rangle\geq\Delta F$, we define the irreversible work as the extra work needed to perform the transformation: $W_{\text{IRR}}\equiv\langle W\rangle-\Delta F$.

IV Harmonic approximation

For small quantum fluctuations, corresponding to small values of $g$, the dynamics of the Wigner crystal described by the Hamiltonian (2) can be expressed in terms of small quantum displacements, coupled harmonically, around the classical equilibrium ion positions. In the linear phase the classical equilibrium positions are $y_{j}=0$ while in the zigzag phase these are $y_{j}=(-1)^{j}b/2$ where the zigzag width $b$ is determined by $\omega$ Fishman et al. (2008). In this regime one can find the normal frequencies associated with a normal mode at momentum $k\in[-\pi,\pi]$ of the harmonic chain of oscillators. In the linear phase these read

[TABLE]

At zero temperature, this semiclassical model predicts a critical transverse frequency at $\omega_{C}=\sqrt{h_{1}}$, for which the frequency of the soft-mode at $k=\pi$ and the quadratic term in Eq. (2) vanish. For $\omega>\omega_{C}$ the chain spatial configuration is linear, while it is zig-zag in the opposite case.

Instead of calculating directly the work probability distribution of Eq. (3), we analytically compute its Fourier transform, namely the characteristic or moment-generating function:

[TABLE]

similarly to the methods reported in Deffner and Lutz (2008); Carlisle et al. (2014). We extract the average irreversible work $W_{\text{IRR}}$ and its statistical variance $\sigma^{2}_{W}=\langle W^{2}\rangle-\langle W\rangle^{2}$ by computing the first two moments of $\chi_{F}(t)$. For an instantaneous quench within the same phase one finds

[TABLE]

[TABLE]

in which $\omega_{k}^{i(f)}$ is the initial (final) frequency of the $k$ mode and $\Omega_{k}=(\omega_{k}^{i2}+\omega_{k}^{f2})/2\omega_{k}^{i}\omega_{k}^{f}$. The irreversible work and the work variance, at zero temperature, are shown in Fig. 1 for small quenches within the same phase and for chain lengths ranging from 60 to 144 ions. It is interesting to note the extensiveness of both the irreversible work and its variance. In agreement with previous results Silva (2008); Dorner et al. (2012); Fusco et al. (2014); Bayat et al. (2016) both $W_{\text{IRR}}$ and $\sigma^{2}_{W}$ diverge at the critical frequency as a consequence of the vanishing of the lowest eigenfrequency.

While the harmonic approximation works well far from the critical point, the vanishing soft mode frequency causes an unphysical divergence, even for a finite number of ions, as the chain approaches criticality, as evidenced in Fig. 1. On both the linear and the zigzag sides of the transition, the irreversible work and the statistical variance are monotonically increasing functions for $\omega\to\omega_{C}$ respectively. Due to the vanishing excitation gap at $k=\pi$ (soft mode), any quench close to the critical point, no matter how small, will always require an amount of work much larger than the mere energy difference between the two equilibrium configurations: $\langle W\rangle\gg\Delta F$.

In order to understand better this divergence, we isolate the contribution to the irreversible work and to the variance of the soft mode. Limiting the sums in Eqs. (6) and (7) to $k=\pi$ and expanding up to second order in $\Delta\omega^{2}$ we obtain

[TABLE]

Both expressions are valid in the linear phase by taking $\gamma_{W}=2$ and $\gamma_{\sigma}=8$ and in the zigzag phase by taking $\gamma_{W}=\sqrt{2}$ and $\gamma_{\sigma}=4$. The expression for $\sigma_{\rm soft}^{2}$ is exact in the linear phase. These results are shown in the inset of Fig. 1.

The scaling found for the soft mode is however modified when adding the other modes. In the thermodynamic limit and close to criticality we can expand Eq. (6) in the linear phase, at $T=0$, as

[TABLE]

in which the constant $A$ depends on a small momentum cut-off introduced when turning the summation in Eq. (6) into an integral. A similar expression holds in the zigzag phase De Chiara et al. (2008).

We remark that the results of this section have been obtained assuming the short range effective model of Eq. (2). We obtain similar expressions, within the harmonic approximation, for the long-range model of Eq. (1) finding the same scaling with renormalised parameters, e.g. the critical frequency is modified and has a weak finite-size correction $\sim 1/L^{2}$. So far, however, we have neglected the non-linear couplings between the normal modes. To overcome this, in the next section we solve numerically the full anharmonic problem.

V Full anharmonic model

In this section we present numerical results from the treatment of the full short-range Hamiltonian (2). These results rely on the assumption that the initial state is prepared at zero temperature $\beta\to\infty$, and thus can be found via variational methods. In fact, the quantum many-body ground states $|\Psi_{G}(\omega)\rangle$ of the Hamiltonian from Eq.  (2) are simulated with the DMRG algorithm, using a numerical technique for continuous-variables quantum systems analogous to Refs. Iblisdir et al. (2007); Silvi et al. (2013, 2014). For any given $\omega$, we evaluate the corresponding ground state energy $E_{G}(\omega)=\langle\Psi_{G}(\omega)|H(\omega)|\Psi_{G}(\omega)\rangle$ and the total fluctuation of the transverse displacement operators $\mathcal{Y}^{2}(\omega)=\sum_{j}\langle\Psi_{G}(\omega)|y^{2}_{j}|\Psi_{G}(\omega)\rangle$. Such data are sufficient to evaluate the average work $\langle W\rangle$ arising from a sudden quench where $\omega$ is instantaneously driven from $\omega_{i}$ to $\omega_{f}$. In fact, we can simplify

[TABLE]

which we can easily compute from the equilibrium data we acquired. Moreover, this expression for the average work shows how to measure it in an experiment by estimating the average quadratic transverse displacement of the ions.

Fig. 2 displays the irreversible work generated by a small quench of the Hamiltonian (2), at $T=0$ and for several values of $L$. The first feature we notice is the disappearance of the divergence at the critical point. This is a clear signature of the finite-size effects in the quantum many-body system: at finite size $L$ the energy gap remains finite for all $\omega$, thus actually smearing out the non-analyticity. On the other hand, the harmonic approach exhibits a critical behaviour even at finite size. Moreover, by increasing the number of ions two features appear: the peak in $W_{\text{IRR}}$ becomes increasingly sharper and its position slowly shifts towards larger frequencies, in contrast to the harmonic theory.

As for the value of the peak itself, we expect it to grow linearly in $L$, $W_{\text{IRR}}$ being an extensive quantity and as already derived from the harmonic theory, see Eq. (10). Fig. 2 fully confirms this prediction. In this case we can also draw a direct connection between our findings and recently results in the study of infinitesimal quenches and ground state fidelity susceptibility for the quantum Ising model Paganelli and Apollaro (2017). In this respect, one can investigate the finite-size effect on the maxima positions by plotting $\omega_{\textrm{max}}$ as a function of $1/L$, see Fig. 2, left inset.

To gain further insight into the behaviour of the irreversible work, we adopt the finite-size scaling ansatz with a model-dependent scaling function $f$

[TABLE]

which is typical of quantities that diverges logarithmically with the control parameter as in the case of the specific heat in the 2D classical Ising model Barber (1983).

In Fig. 3 we show the rescaled data and we obtain a collapse of the irreversible work as in the Ising model with the same critical exponent $\nu=1$. We remark that given the magnitude of the data shown in Fig. 3, a similar collapse plot (not shown) is obtained by expanding the exponential function: $1-\exp[(W_{\text{IRR}}-W_{\text{IRR}}^{\rm(max)})/L]\simeq(W_{\text{IRR}}^{\rm(max)}-W_{\text{IRR}})/L$.

VI Conclusions

In this work we have investigated the irreversible work production associated with infinitesimal quenches of the transverse frequency around a structural phase transition of Coulomb crystals. We have employed different approaches starting from a harmonic approximation that allows us to obtain analytical results but in turn generates an unphysical vanishing gap at $k=\pi$ at the critical point, which is known to drive the structural phase transition. We isolated the contribution of the soft mode to the irreversible work finding a power law scaling of the irreversible work. Finally, we have studied the full anharmonic model through DMRG $T=0$ simulations. With these different approaches we have observed the extensiveness of the irreversible work, and how they generate different scaling properties. Interestingly, the scaling laws recently found for the irreversible work of the Ising model, in terms critical exponents and collapse ansatz, are recovered in the full anaharmonic model when the $g$ parameter is big enough to appreciate shifts from the classical critical point due to pure quantum effects. Beyond the fundamental interest, it has been found that such critical behaviours could be used to design a quantum Otto engine where it has been showed that, for a working substance around criticality, the Carnot point can be reached Campisi and Fazio (2016). In this sense our setup is of particular experimental interest, since we have also shown how the calculated quantities can be related to the fluctuations in the transverse displacement of the ions.

Acknowledgements.

We acknowledge support from the Horizon 2020 EU collaborative projects QuProCS (Grant Agreement 641277) and TherMiQ (Grant Agreement 618074). PS gratefully acknowledges support from the EU via UQUAM and RYSQ, the DFG via SFB/TRR 21, and the Baden-W rttemberg Stiftung via Eliteprogramm for Postdocs.