Quantum quench dynamics of the attractive one-dimensional Bose gas via the coordinate Bethe ansatz

We appreciate the time that the referee has invested in reading our manuscript and writing the report. We would liked to have performed the calculations for larger numbers of particles, but even going to $N=5$ is prohibitive, due to the numerical difficulties as described in appendix B. Our further responses to the requested changes are below.

Requested changes:

Minor questions/comments

• I find the notation for the eigenstates of the Hamiltonian slightly confusing, in particular in Section 4.1. Maybe it would be better to use a notation which makes more transparent the presence of bound rapidities in the state.

Response: We appreciate the referee's point, but it would not be so straightforward to introduce a new notation to address this, particularly as we are using the same quantum number convention to label states of the attractive and repulsive system. We would have to somehow retain the current quantum numbers and introduce an additional notation indicating, e.g., bound states also, and we think that on balance doubling up the notation like this would probably make the presentation less, rather than more, clear. Where the structure of the state (i.e., as a two- or three-body bound state, or a gas-like state) is important we do indicate this in the text alongside the quantum number pair labelling the state. It is perhaps not as elegant as it could be, but we believe all the important information is communicated unambiguously.

• Why do the authors not consider the time evolution of the two point function g(1)(x,t) after the quench?

Response: We do not show the time evolution of g(1)(x,t) mainly because of the computational expense of calculating this quantity. As compared to g2(x,t), the first-order correlation function requires an additional integration over space, and thus an additional combinatorial increase in the number of terms involved. It would therefore be very expensive to obtain a high-resolution coordinate-space representation of g(1)(x,t) analogous to the one we present for g(2)(x,t). The most important features of the first-order correlations are shown in our plots of the momentum distribution n(k,t), which is perhaps a more useful physical picture anyway. The calculation of this quantity is more tractable as we only plot the first six momentum modes, which can be evaluated exactly from a comparatively sparse spatial representation of g1(x,t).

We appreciate the time that the referee has invested in reading our manuscript and writing the report, and the positive comments Our responses to the requested changes are below.

Requested changes: 1- The authors should make a communication effort, putting in value in the text the most important results (and maybe put some benchmarking results in appendices?). In this way the paper would become accessible to a broader audience, and directly target experimentalists.

Response: We respond to this first point in the requested changes with reference to the issue raised by the referee in "Weaknesses": that it is not easy to grasp the main physical conclusion of the research. We would like to point out that our purpose in this paper is, as is often the case in theoretical/computation papers of this nature, primarily to report the results of our calculations. We submit, respectfully, that it is not essential that such papers have a "main physical conclusion". We believe the contribution of the paper should be clear enough: we have applied an essentially exact methodology to investigate the nonequilibrium correlation functions resulting from this particular quench scenario, in the limiting case of small particle numbers to which our approach is applicable. We have, where possible, contrasted the behaviour of this system to the more well-established behaviour of the Lieb-Linger gas following a quench to repulsive interactions, and compared our results to those of other approaches to this problem where such exist -- namely, the results for the post-quench values of g2 found by Piroli et al. (Refs. [35,36]).

We note that all three referees have indicated the clarity with which our results are reported, and referees one and two have both remarked positively on the level of detail with which our results are described. While we may be missing an opportunity to promote our work more vigorously to, e.g., experimentalists in this area, we believe that the quality of communication of results in the paper is of a standard sufficient for publication, consistent with many of the remarks made by the three referees.

Having said that, we do agree that it would be worthwhile to include more text to summarise and highlight the key results of Sections 4 and 5 at the ends of these sections. We have expanded on the summary at the end of Section 4 to write:

"In summary, quenches from the ideal-gas ground state to attractive values of \gamma result in the occupation of energy eigenstates containing bound states in addition to the gas-like scattering states of the attractively interacting model, which are analogous to the eigenstates of the repulsively interacting Lieb–Liniger gas. As the magnitude |\gamma| of the final interaction strength is increased, the postquench occupations of the gas-like excited states approach those of their counterparts following a quench to the corresponding repulsive interaction strength, and the occupations of bound states eventually decrease. However, these bound states significantly influence the dynamics of postquench correlation functions for all final interaction strengths we have considered, causing large oscillations in local correlations and in the occupation of the zero-momentum mode. For large attractive values of \gamma, bound states are highly localized and thus influence the second-order correlation function only at small separations, whereas at larger separations this function exhibits postquench dynamics similar to those observed following quenches to repulsive interactions [32].

We have added also the following summary paragraph at the end of Section 5:

"In summary, our results for the time-averaged local second-order correlation function g(2)DE(0) are consistent with an enhancement of this quantity over the initial ideal-gas value by a factor of \simeq 4 in the limit of strong final interaction strengths, and thus with the predictions of Refs. [35, 36] in this limit. Our calculations also reveal an enhancement of the local third-order correlation function g(3)DE(0) over the ideal-gas value by a factor of \simeq 20 for strong interactions, suggesting that the postquench state would be susceptible to large three-body recombination losses in practice. Results for time-averaged correlation functions at interparticle separations larger than the characteristic extent of bound states are comparable to those obtained previously [32] for quenches to repulsive interactions."

We believe that these revisions to the text, together with the additions made in response to the referee's second point, described below, make the central results of our work clearer.

2- In the introduction, I suggest to make a more explicit connection between the current work and the previous work by the same authors Ref 33, as well as to explain more in detail their contribution wrt to previous works in Refs [35-36].

Response: Thank you for this suggestion. We have revised and extended the introduction. It now reads in part:

"In Refs. [32, 33] we developed a methodology for the calculation of equilibrium and nonequilibrium correlation functions of the repulsively interacting Lieb-Liniger gas based on the semi-analytical evaluation of matrix elements between the eigenstates of the Lieb-Liniger Hamiltonian given by the coordinate Bethe ansatz. Here we extend this approach to the attractively interacting gas, for which the Bethe rapidities that characterize the eigenstates are in general complex-valued, indicating the presence of multiparticle bound states. We apply our method to calculate results for the time evolution of correlation functions following a quench to attractive interactions from the ideal-gas ground state for a system of four particles. As in our previous studies of quenches to repulsive interactions [32, 33], we find that finite-size effects are significant for quenches to weak final interaction strengths. For strong final interaction strengths our results for the time-averaged local second-order correlation function are consistent with the stationary values in the thermodynamic limit calculated in Refs. [71, 72]. In contrast to that work, however, our approach allows us to also calculate the time-averaged value of the postquench third-order correlation function, which we find to be dramatically enhanced over the ideal-gas value, implying that three-body recombination losses would be significant in experimental realizations of the quench. Our approach also allows us to calculate the dynamical evolution of correlation functions following the quench, and for a quench to strong attractive interactions we observe behaviour similar to that following a quench to repulsive interactions of the same magnitude, superposed with characteristic contributions of bound states at small interparticle separations."

We hope the referee agrees that this both makes more explicit the connection between the current work and our previous Refs. [32,33], and explains the contribution of this work with respect to the results of Piroli et al. [35,36].

3- Page 9, the authors mention that the lowest nonzero momentum modes start to deviate from the $k^{-4}$ scaling. Why is this surprising? The scaling holds only at large k. The authors should explain what they meant.

Response: We simply meant to point out that for the weakest interaction strengths we consider, this scaling effectively persists for all nonzero k in the finite-size geometry we consider, and that as interactions are increased less trivial structure of n(k) is observed. We have revised the text to state: "For the case of $\gamma = −0.1$ (red empty circles), interactions are sufficiently weak that no visible deviation from this scaling is visible at the smallest nonzero momenta $k_j$ resolvable in our finite geometry. By contrast, for $\gamma = −0.21$ (green triangles), less trivial behaviour of the momentum distribution can be seen, with the lowest nonzero momentum modes deviating visibly from the $\propto k^{-4}$ scaling."

4- Page 11, the authors state that at stronger interactions the agreement with mean-field theory are in closer agreement with their results. Can the authors state why? Is it coincidental? One would expect that mean-field theory to hold only at weak interactions.

Response: We do not have a precise argument for this point, so we have not added anything to the manuscript text. We do not think it is simply coincidental, however. In the case of a repulsive Bose gas mean-field theory degrades as the interaction strength is increased, due to particles being scattered out of the mean-field wavefunction (i.e., a uniform condensate). For the correlated cluster-like states of the attractive Bose gas we consider here, this picture no longer applies. The mean-field wavefunction in this case corresponds to all particles residing in the same single-particle state -- as in the usual Gross-Pitaevskii/Hartree approximation -- except that this state is now localized in space (in a symmetry-broken picture). As the interaction strength is increased the particles in the exact N-particle wavefunction become more tightly bound together, and this wavefunction apparently comes to more closely resemble the Hartree product of localized states, as the attractive interactions cause the particles to effectively "condense" together in real space. (In the finite system the wavefunction is of course in both cases a symmetry-restored superposition of these localized states, but apparently the exact symmetric wavefunction also agrees more closely with the symmetrized mean-field wavefunction as the interaction strength is increased.)

5- Page 15 The authors mention a monochromatic oscillatory behaviour for a quench to $\gamma=-0.5$. This is not easy to check since the figure is in log time scale. Can the authors provide an inset at short times in linear scale?

Response: With respect, we do not think it is worthwhile to add an inset showing the short-time behaviour here, given that the figure panels are already quite dense and that this is a relatively unimportant point. We hope the referee would agree that there would be little to be gained by attempting to mislead the reader on a point of such little consequence (indeed, all it really serves to show is that the dynamics of our finite-size system are uninteresting for such weak interactions). Moreover, we were careful to say that the behaviour is "nearly" monochromatic, and that this is the case should be sufficiently clear from the log-scale plot as it is. We have, however, revised the text later in the paragraph from "determines the frequency" to "determines the dominant frequency" for consistency.

6- Page 16 The authors notice that $g^3$ is less noisy than $g^2$. Is there any explanation for this?

Response: We have rewritten the text to say: "For larger attractive values of the postquench interaction strength, on the other hand, the shape of g(3)(0,t) is more regular compared to g(2)(x=0,t), reflecting the fact that only one three-body bound state contributes to the postquench wavefunction, whereas multiple states containing bound pairs are present."

7- Fig 7: the long-times behaviour look very noisy. Is it an artefact of the log scale? Or what is the reason for that? Should we trust the curves at such times?

Response: We believe that this appearance results, as the referee suggests, simply from the log scale, and probably the particular linestyles used. It is perhaps not pretty but we believe it communicates the important points clearly enough. We do expect that the curves can be trusted at these late times, as only weakly occupied (and high energy) eigenstates are neglected in calculating these curves. The effects of these omitted contributions would therefore be very small, high frequency perturbations to the plotted curves, and it is unlikely these corrections would even be visible on the scales of these plots.

We appreciate the time that the referee has invested in reading our manuscript and writing the report. Our responses to the requested changes are below.

Requested changes: 1. I think the presentation would benefit from dividing subfigures of figures 1,2 and 11 into separate figures and placing them closer to where they are referred to in the text. It's just a suggestion.

Response: While we appreciate the suggestion, all the subplots in each of these figures are closely related to one another, and therefore we believe they are best presented together. Currently they are all within a page of where they are discussed in the text - breaking them up wouldn’t help this a lot. Additionally, there are already quite a few figures in the manuscript, and we believe that creating 8 more would potentially make the manuscript more difficult to read.

1. On page 9 the authors write: “In Fig. 1(a) we compare our exact results to the mean-field solution just on the localized side of the crossover at $\gamma = −0.21$ (green crosses), and find that the exact many-body solution (green dotted line) is slightly more localized.” Is there an explanation why the exact solution is more localized. Naively one could expect the quantum fluctuations to delocalize the exact solution as compared with the mean-field result.

Response: The referee is correct that we would expect the mean-field solution to be more localised in general. The opposite behaviour observed here for $\gamma=-0.21$ appears to be due to finite-size broadening of the crossover, as can be seen from close inspection of the entropies plotted in the inset to Fig. 1(d). We have added the following text on this point: "We note that this behaviour is consistent with that of the entanglement entropy [inset to Fig. 1(d)], which is smaller for the exact solution than for the mean-field approximation for $|\gamma| \gtrsim 0.23$. By contrast, at weaker interaction strengths finite-size rounding of the crossover yields an entropy for the exact system larger than the mean-field value."

1. On page 10 the authors write “We note that the momentum distributions for the most strongly interacting systems considered here are much broader than those of ground states in the strongly repulsive regime (cf., e.g., Ref. [33]).” Could authors provide a bit more details here instead of just referring to the previous work?

Response: Thank you for pointing this out. We have amended the text to state: "We note that the momentum distributions for the most strongly interacting systems considered here are much broader than the "hump" that forms in the ground-state momentum distribution of the repulsive gas in the strongly interacting Tonks limit, which extends to $\simeq 2k_F$ [33, 87, 88]."

1. On page 24 the authors write “For comparison, we also plot the constant ideal-gas form g(2)(x,t=0)=1−1/N of this function in the initial state”. Wouldn’t it be worth to compare the post-quench results with the g(2) of the ground state with final value of the interactions instead?

Response: The purpose of plotting this ideal-gas value is largely so the reader can see where g2(x) has been enhanced or suppressed relative to the initial state. As the initial value is 1-1/N=0.75 rather than simply 1, we think it is worthwhile to have it indicated here. Comparison to the ground state at the final $\gamma$ would perhaps not be so useful -- recall from Fig. 2 that at $\gamma=-40$ we have g2(0)=100! It might be more enlightening to compare to thermal values of the correlation function corresponding to the conserved energy, as we have done in Refs. [32,33]. However, calculating these quantities was computationally prohibitive in the attractive case. We have amended the text to read: "In Fig. 11(b), we plot the second-order correlation function gDE(2)(x) in the diagonal ensemble for several postquench interaction strengths $\gamma$ and compare these to the initial-state form g(2)(x,t=0)=1−1/N of this function (horizontal line)."