Optimal Persistent Currents for Interacting Bosons on a Ring with a Gauge Field Marco Cominotti,1, 2 Davide Rossini,3 Matteo Rizzi,4 Frank Hekking,1, 2 and Anna Minguzzi1, 2 1 arXiv:1310.0382v2 [cond-mat.quant-gas] 8 Jul 2014 Universit´ Grenoble Alpes, LPMMC, F-38000 Grenoble, France e 2 CNRS, LPMMC, F-38000 Grenoble, France 3 NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, I-56126 Pisa, Italy 4 Institut f¨r Physik, Johannes Gutenberg-Universit¨t Mainz,Staudingerweg 7, D-55099 Mainz, Germany u a We study persistent currents for interacting one-dimensional bosons on a tight ring trap, subjected to a rotating barrier potential, which induces an artificial U (1) gauge field. We show that, at intermediate interactions, the persistent current response is maximal, due to a subtle interplay of effects due to the barrier, the interaction and quantum fluctuations. These results are relevant for ongoing experiments with ultracold atomic gases on mesoscopic rings. PACS numbers: 67.85.-d, 03.75.Lm, 71.10.Pm, 73.23.Ra A quantum fluid confined on a ring and subjected to a U (1) gauge potential displays a periodicity in the particle current as a function of the flux of the corresponding classical gauge field. This persistent current phenomenon is a manifestation of the Aharonov-Bohm effect, and reflects the macroscopic coherence of the manybody wave function along the ring. Such currents were observed more than 50 years ago in bulk superconductors [1] and, more recently, in normal metallic rings, overcoming the challenges of the decoherence induced by inelastic scattering [2]. The most recent developments in the manipulation of ultracold atoms on ring traps [3, 4] have disclosed a novel platform for the study of persistent currents, which can be induced by the application of a rotating localized barrier or, alternatively, by inducing suitable artificial gauge fields [5]. Tunable localized barriers in toroidal Bose-Einstein condensates have been realized, using well-focused, repulsively tuned laser beams [4]. Also, recently the engineering of an atomic superconducting quantum interference device was demonstrated [6]. The unprecedented variety of interaction and barrier strength regimes paves the way to applications such as high-precision measurements, atom interferometry and quantum information, e.g., by the construction of macroscopic superposition of current states and flux qubits [7–9]. The scenario becomes particularly intriguing if the transverse section of the ring is sufficiently thin to effectively confine the system in one dimension (1D): the rich interplay between interactions, quantum fluctuations, and statistics acquires a role of primary relevance. In absence of any obstacle along the ring, the persistent currents display an ideal sawtooth behavior as a function of the flux, i.e., perfect superfluidity for any interaction strength at zero temperature [10, 11]. Diamagnetic or paramagnetic response depending on the population parity is expected [12] for fermions but not for bosons. If a localized barrier is added, persistent currents are smeared—their shape taking a sinusoidal form in the case of large-barrier or small-tunneling limit—as obtained for thin superconducting rings from a Luttinger-liquid ap- proach [13]. Beyond these limiting regimes the physics of bosonic persistent current remains unexplored. The aim of the present work is to provide a complete characterization of persistent currents for 1D bosons, in all interaction and barrier strength regimes. By combining analytical as well as numerical techniques suited for the 1D problem, we show that the current amplitude is a nonmonotonic function of the interaction strength and displays a pronounced maximum in all regimes of barrier height. The presence of an optimal regime illustrates the highly nontrivial combination of correlations, quantum fluctuations and barrier effects. Our results demonstrate that, in a large range of interaction strengths, unwanted impurities or imperfections on the ring only weakly affect the system properties. For the application to quantum state manipulation, the regimes of choice should be either very weak or very strong interactions, where the response to a localized external probe is stronger. Our predictions are readily amenable to experimental testing with quasi1D ultracold atomic gases confined in mesoscopic uniform and lattice rings. Indeed all the interaction regimes, from weakly interacting quasicondensate to the impenetrableboson Tonks-Girardeau limit, have been experimentally demonstrated for 1D bosonic wires [14]. Persistent currents for bosons under a gauge field.— Let us consider a 1D ring of circumference L, containing N bosons of mass M interacting with each other via a contact potential v(x − x ) = g δ(x − x ). The ring contains a localized barrier modeled as Ub (x, t) = U0 δ(x − V t), moving along its circumference at constant velocity V . This induces an effective gauge field with Coriolis flux Ω = M V L/2π in the rotating frame. In the corotating frame, the Hamiltonian reads N H= j=1 2 2 2M −i ∂ 2π g − Ω +U0 δ(xj )+ ∂xj L 2 N δ(xl −xj ) . j,l=1 (1) This generalizes the Lieb-Liniger model [15] to the rotating case and is nonintegrable due to the presence of the barrier. The corresponding many-body energy spectrum is periodic in Ω with period 1, giving rise to the 2 0.5 N 1 N 0.25 weak barrier strong barrier 0 0.25 0.5 1 Α 0.5 Α a 0.25 0.5 0.75 1 0 weak barrier strong barrier b 0.25 0.5 0.75 1 FIG. 1. Persistent current (a) in the rotating frame and (b) in the nonrotating frame, in units of I0 = 2π /M L2 for N = 18, for zero barrier (black-dashed line), weak barrier (λ = 0.2, red-solid line) and strong barrier (λ = 10, blue dotted line) at fixed interaction strength (γ = 0.004). Aharanov-Bohm effect. We consider the stationary regime, reached after the barrier has been adiabatically switched on at early times, such that no high-energy excitations are created. At zero temperature, the spatially averaged particle current I(Ω), or persistent current, is obtained from the groundstate energy E(Ω) via the thermodynamic relation [16] I(Ω) = − 2 2 1 ∂E(Ω) . 2π ∂Ω (2) In absence of a barrier, for any interaction strength, the ground-state energy is given by a series of parabolas with well-defined angular momentum, shifted with respect to each other by Galilean translation [17] and intersecting at the frustration points Ωj = (2j + 1)/2. The corresponding persistent current (Fig. 1) is a perfect sawtooth in the rotating frame and a staircase in the nonrotating frame, corresponding to states with well-defined angular momentum (dashed lines). The addition of a barrier breaks the rotational invariance and induces coherent superpositions of states of different angular momentum. This gives rise to gap openings in the many-body energy spectrum at the frustration points, thus forming the working points for a qubit based on a superposition of current states. The level mixing is visible in the persistent current. In the rotating frame, I(Ω) is a smeared sawtooth for weak barriers (solid lines) and a sinusoid for large barriers (dotted lines). In order to characterize the qubit, we focus on the amplitude α of the persistent current in the rotating frame, determining its value for all regimes of dimensionless interaction strength γ = M g/ 2 n0 , with n0 = N/L, and dimensionless barrier height λ = M U0 L/π 2 . Noninteracting and impenetrable boson limits.—Both for zero and infinitely large repulsive interactions, it is possible to find an exact solution to the many-body Schr¨dinger equation HΨ(x1 , ... xN ) = EΨ(x1 , ... xN ). o For a noninteracting (NI) Bose gas, the many-body N wave function ΨNI (x1 , ... xN ) = i=1 ψ0 (xi ) is simply given by the product of N identical single-particle wave functions ψ0 (xi ), which are ground-state solutions of the corresponding one-body Schr¨dinger equao tion, 2M −i∂x − 2π Ω ψn + U0 δ(x)ψn = εn ψn , and L has energy ENI = N ε0 . In the infinitely interacting limit of impenetrable bosons, or Tonks-Girardeau (TG) gas, the solution is obtained by mapping the system onto a gas of non-interacting fermions subjected to the same external potential [18], ΨTG (x1 , ... xN ) = 1≤j< ≤N sgn(xj − x ) det[ψk (xi )]. The corresponding N −1 energy ETG = k=0 εk and density profile n(x) = N −1 |ψk |2 directly reflect the fermionization properties k=0 of the strongly interacting Bose gas. The Friedel oscillations of the density profile [Fig. 2(d)] are indeed a signature of the strongly correlated regime [19]. The persistent current amplitude in the zero and infinitely interacting limits, obtained from Eq. (2) (see the Supplemental Material [20]), is shown in Fig. 3. We find that the current amplitude depends on the interaction regime. In particular, for all values of barrier strength, α is always larger in the strongly interacting regime than in the noninteracting one. This behavior is explained by noticing that the barrier affects the lowest-lying energy levels more than the high-energy ones [8]. As a result, the current amplitude is smaller for noninteracting bosons, occupying only the lowest level, than for the TG gas, where the levels are filled up to the Fermi energy. Weak interactions.—In order to explore the role of interactions for quantum state manipulation, we start from the noninteracting result and consider the effect of weakly repulsive interactions. In this regime we neglect quantum fluctuations [21] and describe the fluid as a BoseEinstein condensate, making use of the mean-field GrossPitaevskii (GP) equation. In the corotating frame, this 2 takes the form 2M (−i∂x − 2π Ω)2 Φ + U0 δ(x)Φ + g|Φ|2 Φ = L µΦ, where Φ is the condensate wave function and µ the chemical potential. This equation admits stable darksoliton √ solutions, with a size given by the healing length ξ = / 2M gn0 . We have found an analytical solution for Φ in terms of Jacobi elliptic functions [20], thus extending Refs. [22] [see also Ref. [23]]. In particular, by comparing it with the numerical solution of the GP equation, we find that the barrier pins the soliton and turns it to the ground-state solution. The resulting density profile has a minimum at the barrier position. The depth and the healing length of the soliton depend on the interaction strength, the barrier height, and the rotation velocity [see Fig. 2(a) and [20]]. The solitary suppression of the density is accompanied by a winding of the superfluid phase, known as a phase-slip [inset of Fig. 2(a)]. We note that a finite-width barrier, as could be realized experimentally, yields very similar profiles [Fig 2(b)]. Its extension to multiple finite-width barriers is known to give rise to a wealth of nonlinear modes [24]. The soliton energy is obtained from the GP energy functional EGP [Φ] = dx Φ∗ 2 /2M (−i∂x − 2π Ω)2 Φ + L g|Φ|4 /2 + U0 δ(x)|Φ|2 , which encodes the dependence on Ω also in the condensate wave function. Computing the 3 Ψ ^2/N | |2 N 0.2 0.2 0.3 1 0.5 0 0.5 1 Π 0 0.1 a 0 0 Θ Π 0.2 Π Λ 95.5 Λ 1.9 Λ 0.5 Λ 0.1 Π ‰ ‰‰ ‰ ‰ ‰ ‰ ‰ ‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰ ‰‰ ‰ ‰‰ ‰‰ ‰‰ ‰ ‰‰ ‰‰‰‰‰‰‰ ‰ ‰‰ ‰ ‰‰ ‰ ‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰ ‰‰‰ ‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰ ‰ ‰ ‰‰ ‰‰ ‰‰ ‰ ‰‰‰ ‰ ‰‰ ‰‰ ‰ ‰ ‰ ‰ ‰‰ ‰ ‰‰ ‰ ‰ ‰ ‰‰ 0.1 ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰ ‰ ‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ ‰‰ ‰ ‰‰ ‰ ‰ l=1.9 l= -p ‰‰ ‰ ‰‰‰ l=19.1 c) l= bL ‰ L d) ‰ ‰ 0 -p 0 p -p p q q †y§^2 êN Φ 0.3 0.1 b) 0 -p 0 q d-barrier s=0.5 s=2.0 p FIG. 2. Density profiles vs angular coordinate θ along the ring. (a) GP analytical solution at weak interaction strength γ = 0.01 for a delta barrier for various values of the barrier strength λ; the inset shows the soliton phase vs θ for the same values of barrier strength. (b) GP numerical solution (solid lines) at γ = 0.01 for a Gaussian barrier of strength λ = 1.9 for various values of the barrier width σ (in units of n−1 ), and corresponding barrier potential (dashed lines, in units of 10 × 2π 2 2 /mL2 ). 0 (c) MPS solution at intermediate interaction γ = 3.3 and filling 0.15. (d) Analytical TG solution (lines) vs MPS (squares) in the hard-core limit. The other parameters used are N = 18 and Ω = 0.4 in all the curves. persistent current amplitude α we find that it increases monotonically with the interaction strength γ, as illustrated in Fig. 3. This is due to the fact that the healing length ξ decreases at increasing interaction strength, thus yielding a more effective screening of the barrier, thereby restoring superfluidity. Strong interactions.— We now turn to the effect of quantum fluctuations on the persistent current at strong interparticle interactions, up to the TG limit. This can be achieved with the help of Luttinger liquid (LL) theory [25], a low-energy, quantum hydrodynamics description of the bosonic fluid, in terms of the canonically conjugate fields θ and φ corresponding to density fluctuations and superfluid phase, respectively. In the rotating frame, the effective LL Hamiltonian for a uniform ring is H0 = vs 2π L dx K ∂x φ(x)− 0 2 1 2π Ω + (∂x θ(x))2 . (3) L K The microscopic interaction strength enters through the Luttinger parameter K and the sound velocity vs . In the case of repulsive contact interactions their dependence on the interaction strength is known (see, e.g., Ref. [26]): at vanishing interactions, K tends to infinity and vs vanishes, while, in the TG limit, K = 1 and vs corresponds to the Fermi velocity of the fermionized Bose gas. In the presence of a barrier, nonlinear terms appear in the Hamiltonian. We treat these perturbatively in the limits of weak and strong barrier strength. In the regime of a weak barrier, its contribution to the Hamiltonian Hb = dx U0 δ(x)ρ(x) is obtained by keeping only the lowest harmonics in the density field expansion +∞ ρ(x) = [n0 + ∂x θ(x)/π] l=−∞ e2ilθ(x)+2ilπn0 x , yielding as the most important term Hb ∼ 2U0 n0 cos[2θ(0)]. This backscattering term breaks angular momentum conservation, allowing for transitions between angular momentum states [27]. At the same time, quantum fluctuations renormalize the barrier strength down according to Ueff = U0 (d/L)K , where d is a suitably chosen shortdistance cutoff [20]. As a result, the persistent current close to the frustration point Ω = 1/2 is a smeared saw- tooth I(Ω)/N = −I0 δΩ 1 − 4δΩ2 + λ2 /π 2 eff −1 , (4) where δΩ = Ω − 1/2, λeff = M Ueff L/π 2 , and I0 = 2π /M L2 . The corresponding amplitude α is shown in Fig. 3. For decreasing interactions down from the TG limit, the quantum fluctuations of density increase, suppressing the barrier more strongly, hence, increasing α. In the opposite case of strong barrier strength, we model the transport across the barrier with a tunneling Hamiltonian Ht = −2tn0 cos[φ(L) − φ(0) + 2πΩ], where t is the tunneling amplitude. In this limit, quantum phase fluctuations lead to a renormalization down of the tunnel amplitude [28], i.e., I(Ω)/N = −(2t/L )(d/L)1/K sin(2πΩ) . (5) Note the dual nature of this result (power 1/K) in comparison to the one obtained above for a weak barrier (power K). Indeed, the quantum phase fluctuations increase with increasing interactions, thereby suppressing the tunneling amplitude more, thus, decreasing α. The duality between the two models can also be used to establish a link between the tunnel amplitude and the barrier height [7, 20]. The persistent current amplitude for large barrier is also shown in Fig. 3. Interestingly, the renormalization by quantum fluctuations at intermediate interactions is effective enough to turn a relatively large barrier into a weak one. This quantum healing phenomenon completely changes the physical scenario and illustrates the dramatic effect of interplay of interactions and quantum fluctuations. Both for weak and large barriers, the LL description breaks down for sufficiently weak interactions when the short-distance cutoff required in the theory increases, until it becomes comparable with the system size. Intermediate interactions and barrier strengths.— Away from the weakly and the strongly interacting regime and for arbitrary barrier strength it is difficult to tackle the many-body Schr¨dinger equation correo 4 1.0 Λ 0.1 0.8 Λ 1.9 0.6 Α Λ 9.5 0.4 Λ 19.1 0.2 Λ 38.2 0.0 0.001 0.01 0.1 1 Γ 10 Λ 95.5 100 1000 FIG. 3. Persistent current amplitude α, in units of I0 = 2π /M L2 , as a function of the interaction strength γ at varying barrier strength λ, for N = 18, from GP equation (black dotted lines), LL approach (blue solid lines), numerical MPS calculations (red squares, red thin lines are guides to the eye) and NI and TG exact solutions (green dashes). Orange diamonds and cyan circles are for a Gaussian barrier of width σ = 0.5, 2 (n−1 ) respectively. 0 sponding to the Hamiltonian of Eq. (1) with analytical approaches; therefore, we employ numerical simulations based on the density-matrix renormalization-group (DMRG) [30]. After discretizing the space on Ns sites, and using the standard Peierls construction, we map Eq. (1) onto a Bose-Hubbard model on a 1D ring lattice at low filling: Ns −tBH e− Hlat = 2πiΩ Ns b† bj+1 +H.c. + j j=1 + (β δj,1 nj − µ nj ) , UBH nj (nj −1) 2 (6) where b† (bj ) are bosonic creation (annihilation) operj ators at site j, tBH is the tunnel energy for bosons on adjacent sites, and UBH is the on-site repulsion energy. The phase twist Ω/Ns accounts for the effect of the Coriolis flux, β is the barrier strength, and µ is the chemical potential. An experimental implementation of the Hamiltonian of Eq. (6) with a localized barrier on one lattice site has been recently proposed [9]. We compute the ground-state wave function and energy using a matrix product state (MPS) [31], optimized for periodic boundary conditions [20]. The resulting density profiles are shown in Fig. 2 (c, d). As compared to Figs. 2 (a), the healing effect is much more pronounced and might be efficiently used to screen unwanted impurities. As anticipated above, for large interactions the density profiles display Friedel oscillations, in very good agreement with the analytical predictions in the TG limit. The oscillations are strongly damped at intermediate and weak interactions. We also compute the amplitude α of the persistent current for a large range of barrier heights λ and interaction strengths γ (Fig. 3) [32]. We find a nonmonotonic behavior connecting the weak- and strong- interaction results as a consequence of the subtle interplay between backscattering and interaction effects. This result allows us to confirm the expected regimes of validity of the analytical estimates. We note that the position of the maximum, i.e., the optimum persistent current, depends on the barrier strength. Experimental considerations.— Our results are also valid for finite-width barriers, used in realistic implementations. We have checked that the use of a Gaussian barrier does not lead to qualitative changes of our results in the mean-field regime [33] (see Fig. 3 for a choice of barrier strength; similar results hold for all the barrier strength regimes), provided that the barrier width remains comparable to the interparticle distance, a condition achievable e.g., with a microscope-focused stirring beam [34]. Thermal fluctuations give rise to additional smearing of persistent currents starting from a typical temperature kB T ∼ N E0 = 2π 2 2 n0 /M L [11, 12, 35]. Within typical experimental constraints on energy and time scales, persistent currents might be observed in the next-generation experiments of mesoscopic ring traps on a chip (i.e., with a ring diameter of ∼ 5 µm), corresponding to a typical energy N E0 ∼ 550 Hz for 87 Rb atoms taking N = 18, as in Fig. 3. Conclusions.—We disclosed the role of quantum fluctuations and correlations in the emergence of persistent currents of rotating bosons in reduced dimensionality, in the presence of a localized barrier. Our results, summarized in Fig. 3, evidence the presence of a maximum as a function of the interaction strength, which is due to the competition between correlations and fluctuations. While at increasing interactions a classical bosonic field screens the barrier more and more, going towards the strongly correlated Tonks-Girardeau regime quantum fluctuations screen the barrier less and less, due to the increasingly fast spatial decay of phase-phase correlations [36]. These predictions can be tested via timeof-flight measurements, similar to those used to probe circulation in atomic rings [3, 4]. In the present work we studied a thermodynamic quantity—the persistent current—but the interplay of barrier and quantum fluctuations will have a similar impact on out-of-equilibrium properties involving wave emission, ranging from collective modes to transport phenomena, including shock wave generation and Hawking radiation of superfluids [37]. Our results are also relevant to other mesoscopic bosonic quantum fluids, as thin superconducting rings, photons in nonlinear twistedpipe waveguides, and solid-state photonic or polaritonic nanocavities etched on a ring-necklace shape. Acknowledgments. We are indebted to L. Amico, R. Citro, R. Fazio, L. Glazman, N. Pavloff, I. Safi, P. Silvi and W. Zwerger for discussions. This work is supported by the ERC Handy-Q grant N.258608, Institut universitaire de France, Italian MIUR through FIRB Project RBFR12NLNA. MPS simulations were performed on the 5 TQO cluster of the Max-Planck-Institut f¨r Quantenopu tik (Garching) and on MOGON cluster in Mainz. [1] B. S. Deaver and W. M. Fairbank, Phys. Rev. Lett. 7, 43 (1961); N. Byers and C. N. Yang, Phys. Rev. Lett. 7, 46 (1961); L. Onsager, Phys. Rev. Lett. 7, 50 (1961). [2] L. P. L´vy, G. Dolan, J. Dunsmuir, and H. Bouchiat, e Phys. Rev. Lett. 64, 2074 (1990); D. Mailly, C. Chapelier, and A. Benoit, Phys. Rev. Lett. 70, 2020 (1993); H. Bluhm, N. C. Koshnick, J. A. Bert, Julie M. E. Huber, and K. A. Moler, Phys. Rev. Lett. 102, 136802 (2009); A. C. Bleszynski-Jayich, W. E. Shanks, B. Peaudecerf, E. Ginossar, F. von Oppen, L. Glazman, and J. G. E. Harris, Science 326, 272 (2009). [3] S. Gupta, K. W. Murch, K. L. Moore, T. P. Purdy, and D. M. Stamper-Kurn, Phys. Rev. Lett. 95, 143201 (2005); C. Ryu, M. F. Andersen, P. Clad´, V. Natarae jan, K. Helmerson, and W. D. Phillips, Phys. Rev. Lett. 99, 260401 (2007); W. H. Heathcote, E. Nugent, B. T. Sheard, and C. J. Foot, New J. Phys. 10, 043012 (2008); K. Henderson, C. Ryu, C. MacCormick, and M. G. Boshier, New J. Phys. 11, 043030 (2009); B. E. Sherlock, M. Gildemeister, E. Owen, E. Nugent, and C. J. Foot, Phys. Rev. A 83, 043408 (2011); S. Moulder, S. Beattie, R. P. Smith, N. Tammuz, and Z. Hadzibabic, Phys. Rev. A 86, 013629 (2012); S. Beattie, S. Moulder, R. J. Fletcher, and Z. Hadzibabic, Phys. Rev. Lett. 110, 025301 (2013). [4] A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 106, 130401 (2011); K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 110, 025302 (2013). ¨ [5] J. Dalibard, F. Gerbier, G. Juzeli¯nas, and P. Ohberg, u Rev. Mod. Phys. 83, 1523 (2011). [6] C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G. Boshier, Phys. Rev. Lett. 111, 205301 (2013). [7] D. W. Hallwood, T. Ernst, and J. Brand, Phys. Rev. A 82, 063623 (2010); D. Solenov and D. Mozyrsky, Phys. Rev. A 82, 061601 (2010); A. Nunnenkamp, A. M. Rey, and K. Burnett, Phys. Rev. A 84, 053604 (2011); D. Solenov and D. Mozyrsky, J. Comput. Theor. Nanosci. 8, 481 (2011). [8] C. Schenke, A. Minguzzi, and F. W. J. Hekking, Phys. Rev. A 84, 053636 (2011). [9] L. Amico, D. Aghamalyan, F. Auksztol,H. Crepaz, R. Dumke, and L.-C. Kwek, Sci. Rep. 4, 4298 (2014). [10] A. J. Leggett, in Granular Nanoelectronics, edited by D. K. Ferry et al., Nato ASI Series B, Vol. 251 (Plenum, 1991), p. 297; A. Mueller-Groeling, H. A. Weidenmueller, and C. H. Lewenkopf, Europhys. Lett. 22, 193 (1993). [11] D. Loss, Phys. Rev. Lett. 69, 343 (1992). [12] A. A. Zvyagin and I. V. Krive, Low Temp. Phys. 21, 533 (1995). [13] F. W. J. Hekking and L. I. Glazman, Phys. Rev. B 55, 6551 (1997); K. A. Matveev, A. I. Larkin, and L. I. Glazman, Phys. Rev. Lett. 89, 096802 (2002). [14] S. Dettmer, D. Hellweg, P. Ryytty, J. J. Arlt, W. Ertmer, K. Sengstock, D. S. Petrov, G. V. Shlyapnikov, H. [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] Kreutzmann, L. Santos, and M. Lewenstein, Phys. Rev. Lett. 87, 160406 (2001); F. Gerbier, J. H. Thywissen, S. Richard, M. Hugbart, P. Bouyer, and A. Aspect, Phys. Rev. A 67, 051602 (2003); B. Paredes, A. Widera, V. Murg, O. Mandel, S. F¨lling, I. Cirac, G. V. Shlyapnikov, o T. W. H¨nsch, and I. Bloch, Nature (London) 429, 277 a (2004); T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 5687 (2004). E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). F. Bloch, Phys. Rev. B 2, 109 (1970). F. Bloch, Phys. Rev. A 7, 2187 (1973). M. Girardeau, J. Math. Phys. 1, 516 (1960). The Friedel oscillations weakly depend on the rotation velocity, as in [38]. See Supplemental Material for details. D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Phys. Rev. Lett. 85, 3745 (2000). R. Kanamoto, L. D. Carr, and M. Ueda, Phys. Rev. Lett. 100, 060401 (2008); Phys. Rev. A 79, 063616 (2012). M. Schecter, D. M. Gangardt, and A. Kamenev, Ann. Phys. 327, 639 (2012); K. Anoshkin, Z. Wu, and E. Zaremba, Phys. Rev. A 88, 013609 (2013); G. Arwas, A. Vardi, and D. Cohen, Phys. Rev. A 89, 013601 (2014). Y. Li, W. Pang, and B. A. Malomed, Phys. Rev. A 86, 023832 (2012); in Localized Excitations in Nonlinear Complex Systems, edited by R. Carretero-Gonz´lez a et al. (Springer, 2014), p. 171. F. D. M. Haldane, Phys. Rev. Lett. 47, 1840 (1981). M. A. Cazalilla, J. Phys. B 37, S1 (2004). R. Citro, A. Minguzzi, and F. W. J. Hekking, Phys. Rev. B 79, 172505 (2009). Close to the TG limit, a renormalization group calculation yields corrections to the LL result, i.e., an additional reduction of the tunnel energy [39]. U. Weiss, Solid State Commun. 100, 281 (1996). U. Schollw¨ck, Ann. Phys. 326, 96 (2011). o I. Danshita and A. Polkovnikov, Phys. Rev. B 82, 094304 (2010); Phys. Rev. A 85, 023638 (2012). The MPS approach generally performs worse in the weakly interacting regime, where single-particle wavefunctions are delocalized across the lattice. Moreover the soft-core nature of bosons increases the inaccuracies due to the cutoff in the maximum on-site occupancy [20]. A quantitative study for the TG regime has been reported in C. Schenke, Ph.D. thesis, Universit´ de Grenoe ble. 2012, p. 51. R. Desbuquois, L. Chomaz, T. Yefsah, J. L´onard, J. e Beugnon, C. Weitenberg, and J. Dalibard, Nat. Phys. 8, 645 (2012). K. K. Das, M. D. Girardeau, and E. M. Wright, Phys. Rev. Lett. 89, 170404 (2002); P. Vignolo and A. Minguzzi, Phys. Rev. Lett. 110, 020403 (2013). A non monotonous behavior of tunnel current was also found in fermionic systems in the BCS-BEC crossover [40]. ´ P.-E. Larr´, A. Recati, I. Carusotto, and N. Pavloff, Phys. e Rev. A 85, 013621 (2012). W. Zwerger, L. B¨nig, and K. Sch¨nhammer, Phys. Rev. o o B 43, 6434 (1991). D. Yue, L. I. Glazman, and K. A. Matveev, Phys. Rev. B 49, 1966 (1994). A. Spuntarelli, P. Pieri, and G. C. Strinati, Phys. Rev. Lett. 99, 040401 (2007). 6 Supplemental Material for: “Optimal Persistent Currents for Interacting Bosons on a Ring with a Gauge Field” Details of the exact solutions for the limiting cases of ideal and Tonks-Girardeau gas We outline the solution for the single particle problem used to obtain the ideal-gas and the Tonks-Girardeau (TG) solution of the main text. The single particle eigenfunctions of the one-body Schr¨dinger equation take the form o ψn (x; Ω) = 1 Nn e−iΩπ eikn (x−L/2) + An,Ω e−ikn (x−L/2) x ∈ [0, L/2) 1 Nn eiΩπ eikn (x+L/2) + An,Ω e−ikn (x+L/2) x ∈ [−L/2, 0) (7) + − By imposing twisted boundary conditions, unity normalization and the cusp condition ∂x ψn (0+ ; Ω) − ∂x ψn (0− ; Ω) = λψn (0; Ω), we obtain kn = ±λ π sin(kn L) , L cos(2πΩ) cos(kn L) where the ± sign refers to a number of particles N odd or even respectively [1], An,Ω = cos(kn L/2+Ωπ) An,Ω = − cos(kn L/2−Ωπ) for N even, while Nn = n L 1 + A2 + 2An,Ω sin(kLL) n,Ω kn sin(kn L/2+Ωπ) sin(kn L/2−Ωπ) for N odd and is the normalization factor. Weak barrier limit — In the regime λ 1 we determined perturbatively the persistent current amplitude, choosing (±) for simplicity N odd. The single-particle energy levels εn are obtained by degenerate perturbation theory around (2π )2 2 the unperturbed parabolas n = 2M L2 (Ω − n) . Close to the frustration point Ω = 1/2 the pairs of degenerate levels are n and 1−n . We immediately get ε(±) = n (2π )2 2M L2 [δΩ2 + (n − 1/2)2 ] ± 2(n − 1/2) ˜ δΩ2 + λ2 n (8) (−) (N −1)/2 (+) ˜ where δΩ = Ω − 1/2 and λn = λ/(2n − 1). The TG ground-state energy is given by ETG = n=0 (εn + εn ). In this sum the ± terms in the r.h.s. of Eq. (8) compensate except the highest-energy one [2], yielding Eq. (4) of the main text with λef f λ/N . This allows us to choose the short-distance cutoff of the Luttinger theory (see below). Details of the soliton-like solution at weak interactions 2 We discuss the soliton-like solution of the GP equation 2M (−i∂x − 2π Ω)2 Φ + U0 δ(x)Φ + g|Φ|2 Φ = µΦ. Our L approach extends the one of Ref. [3, 4] to the case where a delta barrier is present. We first recast the GP equation ˜ in dimensionless form by introducing Φ(θ) = L/2πΦ(2πx/L), µ = M L2 µ/(2π 2 2 ) and g = gM N L/(π 2 ), and take ˜ ˜ ˜ θ ∈ [0, 2π]. A parametrization of the condensate wavefunction in density-phase representation Φ(θ) = f (θ)eiφ(θ) yields − f + f (φ )2 − 2Ωf φ + (Ω2 − µ + g f 2 )f = 0 , ˜ ˜ (9) −2f φ − f φ + 2Ωf = 0 . (10) The effect of the delta barrier is replaced by the cusp condition f (0+ ) − f (0− ) = λf (0), where λ = M U0 L/(π 2 ), which, assuming a symmetric cusp f (0+ ) = −f (0− ) and introducing the density s = f 2 reads s (0+ ) = λs(0). We first integrate Eq. (10) to obtain φ , φ = C + Ω, f2 (11) where C is an integration constant. Substituting this result into Eq. (9) we get −f + C 2 /f 3 + (˜f 2 − µ)f = 0 , which, g ˜ upon integration and change of variables, yields s 2 = −4C 2 + 2˜s3 − 4˜s2 + 4As , g µ A being an integration constant. Introducing the potential U (s) = 2C 2 − 2As + 2˜s2 − g s3 , we see that the problem µ ˜ is equivalent to the one of a classical particle of unitary mass with position s and velocity s having zero total energy. 7 Denoting U (s) = −˜(s − s1 )(s − s2 )(s − s3 ) = 0 with s1 ≤ s2 ≤ s3 , an allowed trajectory is possible in the interval g s1 < s < s2 if A > 0. In the presence of the barrier, in order to satisfy the the cusp condition, the soliton trajectory starts at initial position smin > s1 . The soliton is then found upon integration, s(θ) θ ds smin −2U (s ) = dθ = θ 0 Introducing the change of variable y 2 = (s − s1 )/(s2 − s1 ), the integral corresponds to the Jacobi elliptic function sn(u|m) with m = (s2 − s1 )/(s3 − s1 ). By imposing periodic boundary conditions for the condensate density (i.e., ˜ requiring that at half-period the soliton solution reaches its maximum value s = s2 ), we find π g (s3 − s1 )/2 + α = K and s(θ) = s1 + (s2 − s1 )sn2 [(K − α)θ/π + α]. (12) Here α = F [arcsin( (smin − s1 )/(s2 − s1 ))|m] with F [φ|m] being the incomplete elliptic integral of the first kind, and K is the corresponding complete elliptic integral. 2π Imposing the normalization condition 0 dθs(θ) = 1 we obtain g 2πs1 + 4(K − α)(K − α − E + α )/(π˜) = 1, where E[φ|m] and E are respectively the incomplete and complete elliptic integrals of the second kind, and α = E[arcsin( (smin − s1 )/(s2 − s1 ))|m] . Substituting Eq. (12) for s in equation s 2 = −2U (s) and equating the terms with the same power of sn, yields the following parameter identification: µ = 3˜s1 /2 + (1 + m)[(K − α)/π]2 ; ˜ g A= C2 = 2˜s1 − 3˜s2 /2 − 2m[(K µ g 1 As1 − µs2 + g s3 /2 . ˜ 1 ˜ 1 (13) 4 − α)/π] /˜ ; g (14) (15) Finally, the cusp condition yields an equation for smin , λsmin = −2U (smin ). (16) Equations (13 – 16) form a coupled set of equations which can be solved as a function of m and smin . Using Eqs. (11) and (12), we then obtain the solution for the phase of the soliton solution: φ(θ) = C(Π[(s1 − s2 )/s1 ; (K − α)θ/π + K|m] − Π[(s1 − s2 )/s1 ; α|m])/(K − α)s1 + Ω(θ + π) , (17) where Π[n; u|m] is the incomplete elliptic integral of the third kind. π Imposing 2π-periodicity of the condensate wave function, integration of Eq. (11) yields 2πn = −π 2π(n − Ω)/C = 2π(Π[(s1 − s2 )/s1 ; K|m] − Π[(s1 − s2 )/s1 ; α|m])/((K − α)s1 ) . dθ C + 2πΩ, hence s (18) Fig. 4 shows the soliton density and phase at fixed barrier strength at varying Coriolis flux Ω, both for the n = 0 case of no winding, and of one quantum of angular momentum, ie n = 1. In summary, we find the soliton solution according to the following strategy: for given values of the interaction constant g , the barrier strength λ and the Coriolis flux Ω, we express s1 , s2 and s3 (hence µ, A, and C) as a function ˜ ˜ of m and smin , then solve simultaneously Eq. (16) and the one obtained after equating Eqs. (18) and (15). This uniquely fixes m and smin , and hence the entire soliton solution. By performing imaginary-time numerical integration of the GP equation, we have checked that the analytical soliton solution given by Eqs. (12), (17) coincides with the numerical ground state, see again Fig. 4. 8 1.5 0.1 b Π 0 Π 0.5 0.4 0.2 0.0 Π D o t t e0 0.2 0 0 Π 0.2 Π 1.0 0.8 0.6 0.5 Π 0.1 0 0 Θ Π Θ Φ Π Φ Π 0.3 2Π 0 Ψ ^2 N Φ Ψ ^2/N | |2 N 0.2 0 0.3 1.5 Ψ ^2/N | |2 N 0.3 0.2 0 0.2 Π 0 Π 0.1 0 Π 0 Θ Π FIG. 4. Left and central panel, soliton solution for the density (main figure) and phase (inset) as a function of the spatial coordinate along the ring for various values of the reduced Coriolis flux Ω as indicated in each figure, at fixed dimensionless barrier height λ = 1.9, N = 18 and g = 1. Right panel, comparison of the analytical soliton solution of the GP equation (black ˜ solid line), with the numerical solution obtained via imaginary-time integration (red dashed line) taking N = 18, Ω = 0.4, ˜ num λ = 0.5, and g = 1. The numerical ground state energy in units of 2π 2 2 /M L2 is EGS = 0.2965(3) to be compared with the ˜ sol ˜GS = 0.2966. one of the soliton solution E Details of the Luttinger-liquid solution We outline the derivation yielding Eqs.(4) and (5) of the main text. In the weak barrier case, we start from the mode expansion for the density and phase fields for the uniform ring, θ(x) = θ0 + φ(x) = φ0 + 1 2 q=0 1/2 2πK qL [eiqx bq + e−iqx b† ] , q 2πx 1 (J − Ω) + L 2 q=0 2π qLK (19) 1/2 sgn(q)[eiqx bq + e−iqx b† ] , q (20) where q = 2πj/L with j integer, J is the angular momentum operator, K is here the Luttinger parameter and the following commutation relations hold: [bq , b† ] = δq,q ; [J, e−2iθ0 ] = e−2iθ0 . The latter property implies that the q zero-mode θ0 acts as a (normalized) raising operator for the states |J of given angular momentum: e−2iθ0 |J = |J +1 . The lowest-order relevant term in the barrier Hamiltonian induces transitions of one quantum of angular momentum [5] due to the zero-mode part in θ(x) ≡ θ0 + δθ(x), i.e. |J − 1 J|e2iδθ(0) + |J J + 1|e−2iδθ(0) . Hb ∼ 2U0 n0 cos(2θ(0)) = n0 U0 J The calculation is performed in two steps. First, using the mode expansion of the fields and averaging the total Hamiltonian over the nonzero modes, we obtain an effective Hamiltonian for the angular momentum operator [6] HJ = E0 (J − Ω)2 + n0 Ueff |J + 1 J| + H.c. . (21) J The effective barrier strength Ueff is obtained by averaging over the density fluctuations Ueff = U0 e±i2δθ(0) = U0 (d/L)K , d being a short-distance cutoff of the order of the interparticle distance. By choosing d = K/n0 , this expression coincides with the exact TG result at K = 1 (see above) and takes into account the shrinking of the linear region of the excitation spectrum once interactions are decreased away from the TG point. In a second step, we perform degenerate perturbation theory between states of given angular momentum, obtaining E(Ω)/N = E0 [δΩ2 − δΩ2 + λ2 /4π 2 ], with E0 = 2π 2 2 /M L2 , δΩ = Ω − J − 1/2, and λeff = M Ueff L/π 2 . Upon eff using the thermodynamic relation, Eq. (2) of the main text one readily obtains Eq. (4) of the main text. In the strong barrier or weak-tunnel limit the appropriate mode expansion is the one with open boundary conditions, θ(x) = θ0 + i q>0 φ(x) = φ0 + q>0 πK qL 1/2 π qLK 1/2 sin(qx)[bq − b† ] , q (22) cos(qx)[bq + b† ] , q (23) 9 where q = jπ/L with j integer. Note that no circulation is allowed in the limit of infinitely large barrier, and hence the angular momentum J does not enter here. As in the weak-barrier case, we average the barrier-contribution of the Hamiltonian over the fluctuation modes obtaining the Ω-dependent part of the energy, E(Ω)/N = −2(t/L) cos(φ(L)− φ(0) + 2πΩ) = −2(t/L)(d/L)1/K cos(2πΩ), where we have made the same choice for the short-distance cut-off length d as in the weak-barrier case above. This readily leads to Eq.(5) of the main text. The relation between tunneling amplitude t and dimensionless barrier strength λ = M U0 L/π 2 at given interparticle interactions is obtained using the result [7], t/L = Γ(1 + K) Γ(1 + 1/K)K ( ωc ) 1+K (U0 /L)−K , (24) where the model-dependent cut-off frequency ωc is determined using the exact TG result in the case K = 1. In essence, for large values of the barrier strength we perform exact calculations of the current I vs Coriolis flux Ω, then use Eq. (5) of the main text as fitting function to extract the value of the tunnel amplitude t. The resulting dependence of t on the barrier strength is shown in Fig. 5. For sufficiently large values of the barrier strength (i.e., λ 50) a very good agreement is found with the hyperbolic law (24), allowing therefore to extract ωc . tL 5 4 3 2 1 0 0 20 40 60 80 100 120 140 Λ FIG. 5. Tunneling amplitude t vs. barrier strength λ in the TG limit for N = 18. For the red dots, the value of λ is too weak for the tunneling approximation to hold, and a deviation from the hyperbolic behavior is observed. Details of the matrix-product-state numerical approach at intermediate interactions In order to tackle the ground-state many-body problem for intermediate interactions, we resorted to numerical simulations based on a discretized version of the Hamiltonian (1) in the main text. This can be written in terms of a Bose-Hubbard model for a 1D chain of Ns lattice sites with periodic boundary conditions (PBC) [see Hamiltonian (6) in the main text]. By taking the continuum limit of Eq. (6) in the main text, we link the parameters tBH , UBH and β of the BoseHubbard model to the parameters M, g and U0 of the continuum Hamiltonian (1) in the main text. Specifically we have 2 /2M = tBH a2 , g = UBH a and U0 = aβ, where a is the lattice spacing of the discrete model. Our simulations are based on the Density Matrix Renormalization Group (DMRG) method, (see, e.g., Ref. [8]). Specifically, we have used a Matrix Product State (MPS) representation of the many-body wavefunction and performed a variational minimization of the energy cost function, ruled by Hamiltonian (6) in the main text, site by site. We controlled the accuracy of the simulation by adjusting the size m of each matrix (bond link) composing the MPS. It is worth mentioning that special care has to be taken for PBC systems with a DMRG-based algorithm, since it performs a factor m2 worse than for open-ended systems [9, 10]. Here we adopted an improved method based on the controlled factorization procedure for long products of MPS transfer matrices, which takes into account only p m2 singular values, without compromising the accuracy and reducing the computational effort [11–13]. This is justified by the fact that, for large chains, the local physics of the system is weakly affected by the properties of the boundaries. In a similar fashion, we approximated the effective Hamiltonian on the MPS basis by expanding it via a singular value decomposition, keeping only the contributions associated to its s ∼ p largest eigenvalues [12]. Unfortunately the implementation of any symmetry in the MPS wavefunction, such as, in our case, the one accounting for the conservation of the total number of bosons, is less trivial with PBC rather than with OBC [14]. This forces us to work in the grand canonical ensemble and fix an average number of particles by suitably tuning the chemical potential. In our simulations we have chosen Ns = 120 sites and fixed the chemical potential µ = µ(tBH , UBH ) in such a way to have N ∼ 18 particles in the system. A fine-tuning of µ in order to have a well defined number of particles has been performed in absence of any barrier, in the case when the system is translationally invariant (the effect of a finite β 10 -5 10 5 Ω = 0. Ω = 0.2 Ω = 0.4 Ω = 0.46 Ω = 0.48 4 3 0.04 2 0.06 10 10 -4 |Em - E∞| E(Ω) - E(0) [× 10 ] 6 -6 -8 0.08 0.1 1/m m = 10 m = 15 m = 20 m = 25 1 0 0 0.1 0.2 Ω 0.3 0.4 0.5 FIG. 6. Ground-state energy per lattice site, in units of tBH , for the system described by the Hamiltonian (6) in the main text, as a function of the phase twist Ω and for different values of the bond-link dimension m (here we fixed p = s = 4m). The dashed line is obtained by extrapolating the values of the energies for m → ∞, according to a power-law fit of the numerical data (symbols): Em ∼ E∞ + A m−B . Here we simulated a Bose-Hubbard chain of length Ns = 120, with UBH = 10, µ = −1.81, β = 1 (parameters are expressed in units of tBH ). We also chose nmax = 2, checking that this is sufficient with this value of interaction. Inset: Scaling of the energies with m, for fixed values of Ω. is that of smearing such accuracy in N —see below). We thus kept an average low filling nj 0.15 (here nj = b† bj j counts the number of bosons on site j), mimicking the continuum limit of Eq.(2) of the main text and minimizing lattice effects. The simulations presented in Fig. (2) and Fig. (3) of the main text have been performed with m = 20, while we fixed p = s up to 120. We also cut the maximum allowed occupation per site to nmax bosons, going from nmax = 1 in the TG regime, to nmax = 4 for UBH = 0.1 (the lowest chosen value for UBH , in units of tBH ). We carefully checked that m 20 is already sufficient to get reliable results for the particle current I(Ω) ∝ −∂Ω E, where E(Ω) is the ground-state energy. An example of the typical performances of our MPS code as varying the bond-link dimension m is provided in Fig. 6, where we explicitly show the dependence of the ground-state energy E as a function of the phase twist Ω. For fixed Ω, we found a dependence of such energy on m consistent with a power-law behavior: Em ∼ E∞ + A m−B , from which we could extrapolate the asymptotic value E∞ (inset). We point out that the bond link values considered in this work are considerably smaller than those used in typical DMRG simulations for open-ended chains, where mobc ∼ 500 could be easily reached. Nonetheless they are sufficient to reach an accuracy of the order 10−7 in terms of absolute values of energies, thus leading to an error in the extrapolation of the particle current I(Ω) of the order of 1%, that is barely visible on the scale of Fig. 6. A further and most important source of inaccuracies comes from the determination of the average number of particles N , which in some cases (especially for small interactions) becomes less accurate when varying the barrier strength. As a combination of these two sources of errors, the bars in Fig. 3 of the main text have been computed from point to point. We found a very similar scenario for all the interaction regimes considered in this work. [1] This follows from the fact that the mapping function A = Π1≤j< ≤N sgn(xj − x ) used to build the TG wavefunction ΨTG (x1 , ... xN ) = A det[ψk (xi )] is periodic (antiperiodic) for odd (even) N respectively; correspondingly, the single-particle orbitals need to be periodic (antiperiodic) for odd (even) N . [2] This is typical of fermionic systems, where only the highest occupied level determines the transport properties. [3] R. Kanamoto, L. D. Carr, and M. Ueda, Phys. Rev. Lett. 100, 060401 (2008). [4] R. Kanamoto, L. D. Carr, and M. Ueda, Phys. Rev. A 79, 063616 (2009). [5] Higher harmonics in the expansion induce jumps of more than one unit of angular momentum, and this processes can be ignored in a weak barrier approximation. [6] A. O. Gogolin and N. V. Prokofev, Phys. Rev. B 50, 4921 (1994). [7] U. Weiss, Solid State Commun. 100, 281 (1996). [8] U. Schollw¨ck, Ann. Phys. 326, 96 (2011). o [9] F. Verstraete, D. Porras, and J. I. Cirac, Phys. Rev. Lett. 93, 227205 (2004). [10] F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57, 143 (2008). [11] P. Pippan, S. R. White, and H. G. Evertz, Phys. Rev. B 81, 081103(R) (2010). [12] D. Rossini, V. Giovannetti, and R. Fazio, J. Stat. Mech., (2011) P05021. [13] M. Weyrauch and M. V. Rakov, Ukr. Phys. J. 58, 657 (2013) [arXiv:1303.1333]. [14] This study is presently under investigation, and will be subject to a future work by us.