The dynamics and prethermalization of one dimensional quantum systems probed through the full distributions of quantum noise

Quantum noise correlations have been employed in several areas in physics including condensed matter, quantum optics and ultracold atom to reveal non-classical states of the systems. So far, such analysis mostly focused on systems in equilibrium. In this paper, we show that quantum noise is also a useful tool to characterize and study the non-equilibrium dynamics of one dimensional system. We consider the Ramsey sequence of one dimensional, two-component bosons, and obtain simple, analytical expressions of time evolutions of the full distribution functions for this strongly-correlated, many-body system. The analysis can also be directly applied to the evolution of interference patterns between two one dimensional quasi-condensates created from a single condensate through splitting. Using the tools developed in this paper, we demonstrate that one dimensional dynamics in these systems exhibits the phenomenon known as"prethermalization", where the observables of {\it non-equilibrium}, long-time transient states become indistinguishable from those of thermal {\it equilibrium} states.


I. INTRODUCTION
Probabilistic character of Schrödinger wavefunctions manifests itself most directly in quantum noise. In manybody systems, shot-to-shot variations of experimental observables contain rich information about underlying quantum states. Measurements of quantum noise played crucial role in establishing nonclassical states of photons in quantum optics [1], demonstrating quantum correlations and entanglement in electron interferometers [2], and verifying fractional charge of quasi-particles in quantum Hall systems [3][4][5]. In atomic physics so far, noise experiments focused on systems in equilibrium. Recent work includes analysis of counting statistics in atom lasers [6], establishing Hanbury-Brown-Twiss effect for both bosons and fermions [7], analysis of quantum states in optical lattices [8][9][10][11][12], observation of momentum correlations in Fermi gases with pairing [13] and investigation of thermal and quantum fluctuations in one and two dimensional condensates [14][15][16][17][18][19][20][21].
In this paper we demonstrate that analysis of quantum noise should also be a powerful tool for analyzing non-equilibrium dynamics of strongly correlated systems. Here we study the two equivalent dynamical phenomena; one given by the interaction induced decoherence dynamics in Ramsey type interferometer sequences for two component Bose mixtures in one dimension [22], and another given by the evolution of interference patterns of two one dimensional condensate created through the splitting of a single condensate [23,24]. We obtain a complete time evolution of the full distribution function of the amplitude of Ramsey fringes or interference patterns. In the case of Ramsey fringes, the average amplitude of Ramsey fringes measures only the average value of the trans-verse spin component. On the other hand, full distribution functions are determined by higher order correlation functions of the spins. Hence full distribution functions contain considerably more information about the time evolution of the system [24][25][26][27] and provide a powerful probe for the nature of the quantum dynamics under study. In particular, we use the simple expressions of full distribution functions to demonstrate the phenomena of "prethermalization" in these one dimensional systems, where observables in non-equilibrium long-time transient states become indistinguishable from those in thermal equilibrium states.
One dimensional systems with continuous symmetries, including superfluids and magnetic systems, have a special place in the family of strongly correlated systems. Quantum and thermal fluctuations are so extreme that long range order is not possible in equilibrium. Such systems can not be analyzed using standard mean-field approaches, yet they can be studied through the application of methods specific to one dimension such as exact Bethe ansatz solutions [28][29][30][31][32][33][34][35][36][37][38][39], effective description using Tomonaga-Luttinger and sine-Gordon models [40][41][42][43][44][45], and numerical analysis using density-matrix renormalization group(DMRG) and matrix product state (MPS) methods [46]. Such systems are often considered as general paradigms for understanding strongly correlated systems. One dimensional systems also give rich examples of integrable systems, where due to the existence of infinite number of conserved quantities, equilibration does not take place [47][48][49]. Hence the problem we consider in this paper is important for understanding fundamental issues such as the quantum dynamics of strongly correlated systems and equilibration/non-equilibration of many-body systems, as well as for possible applications of spinor arXiv:1104.5631v2 [cond-mat.quant-gas] 15 Jun 2011 condensates in spectroscopy, interferometry, and quantum information processing [50][51][52].
Our work is motivated by recent experiments of Widera et al. [22] who used two hyperfine states of 87 Rb atoms confined in 2D arrays of one dimensional tubes to perform Ramsey type interferometer sequences. They observed rapid decoherence of Ramsey fringes and the near absence of spin echo. Their results could not be explained within the single mode approximation which assumes a macroscopic Bose condensation into a single orbital state, but could be understood in terms of the multi-mode Tomonaga-Luttinger type model. Yet the enhanced decoherence rate and suppression of spin echo do not provide unambiguous evidences for what the origin of decoherence is. In this paper, we suggest that the crucial evidence of the multi-mode dynamics as a source of decoherence should come from the time evolution of the full distribution functions of the Ramsey fringe amplitude. Such distribution functions should be accessible in experiments on Atom Chips [14,[53][54][55][56] from the analysis of shot-to-shot fluctuations.
This paper is organized as follows. In Section II, we describe two physically distinct, but yet mathematically equivalent, dynamics in one dimensional systems, namely, the dynamics of spins in a Ramsey sequence and the dynamics of phase and contrast in interference patterns between two split condensates created from a single condensate. We start with illustrating the basic physics governing the dynamics studied in this paper, and give a summary of the central results including the prediction of prethermalization phenomena. The formal descriptions of the details of the theory for spin dynamics in Ramsey sequence start from Section III, where we give the Hamiltonian of the one dimensional system based on Tomonaga-Luttinger approach. In Section IV, we derive the analytical expression for the time evolution of the full distribution function for a simple case in which charge and spin degrees of freedom decouple. This decoupling limit gives a good approximation to the experimental situation of Widera et al [22]. A short summary of the result in this decoupling limit has been already reported in Ref [52]. More general case in which spin and charge degrees of freedom mix is studied in Section V. Such mixing introduces the dependence of spin distribution functions on the initial temperature of the system. All the results obtained in Sections IV and V can be extended to the study of the dynamics of interference patterns, using the mapping described in Section II. The details of the dynamics of phases and contrasts in interference patterns between split condensates is studied in Section VI. We demonstrate the prethermalization phenomena and show that the interference contrasts of split condensates in a steady state have indistinguishable distributions from those of thermal condensates at some effective temperature T eff . We conclude in Section VII with a discussion of possible extensions of this work.
(1) z x π/2 (2)  (2) π/2 pulse is applied to rotate each atom into the x direction; (3) spins freely evolve for time t. In actual experiments, final π/2 pulse is applied to measure the x component of spin. The imaging step (4) is omitted in the illustration. In this paper, spin operators refer to the ones before the final π/2 pulse.

A. Ramsey Dynamics
In this paper, we study the dynamics of one dimensional, interacting two-component Bose mixtures in the Ramsey-type sequence. In analogy with spin-1/2 particles, we refer to one component to be spin-up and the other component to be spin-down. In the experiment of cold atoms in Ref [22], two hyperfine states are used for these two components. In the following, we consider a generic situation where there is no symmetry that relates spin-up and spin-down. In particular, unlike fermions with spin-1/2, there is no SU(2) symmetry. In a typical experimental setup with cold atoms, there is a harmonic confinement potential along the longitudinal direction of condensates, but here we assume the absence of such a harmonic trap potential. Our consideration gives a good approximation for the central region of cold atom experiments in the presence of such potentials.
Ramsey-type sequence is described as follows( Figure  1): 1. All atoms are prepared in spin up state at low temperature 2. π/2 pulse is applied to rotate the spin of each atom into the x direction 3. Spins evolve for time t 4. Spins in the transverse direction (x − y plane) are measured In a typical experimental situations [22], the last measurement step is done by applying a π/2 pulse to map the transverse spin component into z direction, which then can be measured. In the following discussions, we describe the dynamics in the rotating frame of Larmor frequency in which the chemical potentials of spin-up and spin-down are the same in the absence of interactions. In this frame, the evolution of spins in the third step is dictated by the diffusion dynamics coming from interactions. Unlike the conventional use of the Ramsey sequence in the context of precision measurements, here we employ the Ramsey sequence as a probe of correlation functions in one dimensional system. The description of the spin dynamics starts from the highly excited state prepared after the π/2 pulse of step 2. The subsequent dynamics during step 3 crucially depends on the nature of excitations in the system. In particular, the dynamics of two-component Bose mixture in one dimension is quite different from those in three dimension. In three dimensions, bosons form a Bose-Einstein condensate(BEC) at low temperature, and particles occupy a macroscopic number of k = 0 mode. Then, the spin diffusion of three dimensional BEC is dominated by the spatially homogeneous dynamics coming from the single k = 0 mode at sufficiently low temperatures. On the other hand, bosonic systems in one dimension do not have the macroscopic occupancy of the k = 0 mode, and their physics is dominated by the strong fluctuations, to the extent that the system cannot retain the long range phase coherence even at zero temperature [41]. Thus, the spin dynamics of one dimensional bosonic system necessarily involves a large number of modes with different momenta and the spin becomes spatially inhomogeneous during the step 3 above.
Such dynamics unique to one dimension can be probed through the observation of transverse spin components in the fourth step. Since we aim to capture the multi-mode nature of the dynamics in one dimension, we consider the observation of spins at length scale l, given bŷ whereŜ a (r, t) with a = x, y are the transverse components of spin operators after time evolution of step 3 of duration t. We assume that l is much larger than the spin healing length ξ s , and much smaller than the system size L to avoid finite size effects. Furthermore, we assume that the number of particles within the length l, N l , is large, so that the simultaneous measurements ofŜ x l (t) andŜ y l (t) are in principle possible. For large N l , the non-commutativity ofŜ x l andŜ y l gives corrections of the order of 1/ √ N l compared to the average values. In this situation, it is also possible to measure the magnitude of transverse spin components,Ŝ ⊥ l = Ŝx l 2 + Ŝ y l 2 , which we will extensively study in the later sections. Due to quantum and thermal fluctuations, the measurements ofŜ a l (t) give different values from shot-to-shot. After the π/2 pulse of step 2, the spins are prepared in x direction, so the average value yields Ŝ x l (t = 0) ≈ N l /2 and Ŝ y l (t = 0) ≈ 0. In the rotating frame of Larmor frequency, the subsequent evolution does not change the expectation value of the y component so that Ŝ y l (t) ∼ 0 throughout. The decay of the average Ŝ x l (t) during the evolution in step 3 tells us the strength of spin diffusion in the system. The behaviors of Ŝ a l (t) due to spin diffusion are similar for one and three dimensions, and the difference is quantitative, rather than qualitative. On the other hand, a richer information about the dynamics of one dimensional system is contained in the noise ofŜ a l (t). Such noise inherent to quantum systems is captured by higher moments (Ŝ a l (t) ) n . In this paper, we obtain the expression for the full distribution function P a l (α, t) which can produce any moments ofŜ a l (t) through the relation where P a l (α, t)dα represents the probability that the measurement ofŜ a l (t) gives the value between α and α + dα. We will see in Section IV that it is also possible to obtain the joint distributions P x,y l (α, β, t) ofŜ x l (t) and S y l (t) as well as the distribution P ⊥ l (α, t) of the squared transverse magnitude Ŝ ⊥ l 2 . Now we summarize the main results of this paper, and give a qualitative description of spin dynamics in the Ramsey sequence. Elementary excitations of spin modes in the system are described in terms of linearly dispersing spin waves with momenta k and excitation energies c s |k|, where c s is the spin wave velocity. When certain symmetry conditions are satisfied(see discussion in Section IV), spin and charge degrees of freedom decouple, and these spin waves are free and they do not interact among themselves in the low energy descriptions within so-called Tomonaga-Luttinger theory [41,42]. Here, we describe the result in this decoupling limit, but the qualitative picture does not change even after the coupling between spin and charge is introduced, as we will see in Sec. V.
The initial state prepared after π/2 pulse in step 2 in which all spins point in the x direction is far from the equilibrium state of the system because interactions of spins are not symmetric in terms of spin rotations. Thus, the initial state contains many excitations and the subsequent dynamics of spins is determined by time evolution of the spin waves. A spin wave excitation with momentum k rotates spins with length scale ∼ 2π/k and time scale ∼ 1/(c s |k|). The amplitude of fluctuations coming from the spin wave with momentum k is determined by the initial state as well as the nature of spin wave excitations. We find that the energy stored in each mode is approximately the same(see discussion in Sec. IV B 5), thus the amplitude of fluctuations for wave vector k scales as 1/k 2 . Therefore, the fluctuation of spins is weak at short wave lengths and short times, and strong at long wave lengths and long times. In Fig.2, we illustrate such dynamics of spins due to fluctuations of spin wave excitations. It leads to the distributions presented in Fig.3 and Fig.4. Here, we have plotted distribution function of the squared transverse magnitude of spins Ŝ ⊥ K s is the spin Luttinger parameter, which measures the strength of correlations in 1D system (see Eq. (8) below), and ξ s is a spin healing length which gives a characteristic length scale in the low energy theory of spin physics.
The multi-mode nature of one dimensional system, in which spin correlations at different length scales are destroyed in qualitatively different fashion during the dynamics, can be revealed most clearly in the squared transverse magnitude of spins Ŝ ⊥ Here ξs is the spin healing length, and the x axis is scaled such that the maximum value of α is 1. Time is measured in units of ξs/cs where cs is the spin sound wave velocity. The evolution of the distribution crucially depends on the integration length. The steady state of the distribution of the squared transverse magnitude has a peak at a finite value for short integration length l/ξs = 20, whereas the peak is at 0 for long integration length l/ξs = 40.
spins strongly depends on the wavelength of the excitations. Spin excitations with momenta much smaller than ∼ 2π/l do not affect Ŝ ⊥ l 2 since these spin waves rotate the spins within l as a whole, while spin excitations with higher momenta lead to the decay of the magnitude. This is in stark contrast with the x component of the spin S x l , which receives contributions from spin waves of all wavelengths.
As a result of different contributions of spin wave excitations with different wavelengths to the integrated spin magnitude, there are two distinct behaviors of the distributions of Ŝ ⊥ l 2 ; one for short integration length l, which we call "spin diffusion regime" and another for long integration length l, which we call "spin decay regime." For short integration length l, the distribution func- is always peaked near its maximum value (ρl) 2 during the dynamics because the strengths of fluctuations coming from spin waves with high momenta are suppressed by 1/k 2 (Fig.3, l/ξ s = 20). While the magnitude of spins does not decay in this regime, fluctuations still lead to a diffusion ofŜ x l , thus, we call this regime the "spin diffusion regime." On the other hand, for long integration length, spin waves lead to fluctuations of the spins within the integration region, and the spins are randomized after a long time. This randomization of spins leads to the development of a Gaussian-like peak near S ⊥ Here axes are scaled such that the maximum value of α and β are 1. For for short integration length l/ξs = 20, the dynamics leads to the distribution with the "ring"-like structure, showing that the magnitude of spins does not decay much (spin diffusion regime). On the other hand, for longer integration lengths, the magnitude of spins decays quickly and the distribution forms a "disk"-like structure(spin decay regime).
0 and the maximum value (ρl) 2 are present, and one can observe the double peak structure. Because of the strong decaying behavior of the magnitude of spin, we call this regime the "spin decay regime." More complete behaviors of distribution functions can be captured by looking at the joint distribution functions from which we can read off the distributions of both Ŝ ⊥ l 2 andŜ x l , see Fig.4[57]. In the "spin diffusion regime" with short l, the joint distributions form a "ring" during the time evolution, whereas in the "spin decay regime" with long l, they form a "disk"-like structure in the long time limit. As we will see later, a dimensionless parameter given by l 0 ∼ π 2 l 4Ksξs determines whether the dynamics belongs to the "spin diffusion regime"(l 0 ≤ 1) or the "spin decay regime"(l 0 1). We emphasize that in three dimensions, spin waves are dominated by k = 0 mode and therefore, there is almost no decay in the magnitude of spins throughout the dynamics. Therefore, the existence of two qualitatively different behaviors of distribution functions of Ŝ ⊥ l 2 unambiguously distinguishes the dynamics in one and three dimensions.

B. Dynamics of interference between split condensates
The dynamics of Ramsey sequence considered above can be directly mapped to the dynamics of interference pattern of split one dimensional quasi-condensate [58]. More specifically we consider the following sequence of operations(see Fig.5). First, we prepare one-component 1D quasi-condensate in equilibrium. At time t = 0, the quasi-condensate is quickly split along the axial direction, This interference pattern contains the information about the local phase difference between the two quasi-condensate at the time of release. and the resulting two quasi-condensates are completely separated. The two quasi-condensates freely evolve for a hold time of t, and they are released from transverse traps to observe the interference pattern between these two quasi-condensates. Such dynamics of the interference patterns as a function of hold time t has been observed in the experiment by Hofferberth et.al. [23], where the average of interference patterns is analyzed in details. Here we study the unique one dimensional dynamics by looking at the full distributions of interference patterns (for an experimental study see: Gring et al [24]. ).
The dynamics of split condensates can be mapped to the dynamics of the Ramsey interferometer studied in this paper. The splitting of a quasi-condensate corresponds to the initial π/2 pulse in the Ramsey sequence. If we call one of the quasi-condensates to be L for left and another to be R for right, L(R)-condensate corresponds to spin-up (spin-down) component. Thus, the density difference between the two condensates corresponds to z component of the spin, namelyŜ z (r, is the creation operator of particles in i = L, R(R) condensate. Moreover, the local phase difference between the two condensates corresponds to the local spin direction in x − y plane. To see this, we first note that ψ † L (r, t)ψ R (r, t) corresponds to the spin raising operatorŜ + (r, t). This operator is expressed in terms of the phase difference between the condensatesφ s (r, t) asŜ + (r, t) ≡ ψ † L (r, t)ψ R (r, t) ∼ ρe iφs(r,t) where ρ is the average density of each condensate. Thus, for example, x and y spin operators are given byŜ x (r, t) ≡ ρ cos(φ s (r, t)) andŜ y (r, t) ≡ ρ sin(φ s (r, t)). Immediately after the splitting, the phases of the two quasi-condensates at the same coordinate along the axial direction are the same. Therefore, the splitting prepares spins in x direction in the language of the Ramsey sequence, and thus the splitting effectively amounts to the π/2 pulse.
The interference of two quasi-condensates measures the local phase difference at time t. If the phases of L and R condensates are the same, the interference pattern has a constructive peak at the center of two condensates, which we call x = 0. Thus, a shift in the interference pattern measures the local phase difference between the two condensates, which yields the information about S x (r) = ρ cos(φ s (r)) as well asŜ y (r) = ρ sin(φ s (r)). We note that the integrated interference contrast that can be extracted from experiments is given by the expres-sionĈ 2 = ρ l e iφs(r) 2 [26] and related to the transverse magnitude of spins as Ŝ ⊥ All the results obtained for Ramsey dynamics are directly applicable to the dynamics of interference patterns between split condensates. When the splitting process prepares two condensates with equal average number of particles, "spin" and "charge" degrees of freedom decouple, see Sec.VI for details. In particular, depending on the integration length of the interference patterns along the axial direction, there exists two regimes corresponding to "phase diffusion regime" and "contrast decay regime," analogous to "spin diffusion regime" and "spin decay regime" described in the previous section, respectively.

C. Prethermalization of one-dimensional condensates
The equilibration and relaxation dynamics of generic many-body systems are fundamental open problems. Among possible processes, it has been suggested that the time evolution of some systems prepared in nonequilibrium states results in the reaching of quasi-steady states within much shorter time than equilibration time. This quasi-steady state is often not a true equilibrium state, but rather it is a dephased state, and true equilibration takes place at much longer time scale. Yet in some cases, the physical observables in the quasi-steady states take the value corresponding to the one in thermal equilibrium at some effective temperature T eff , displaying disguised "thermalized" states. Such surprising nonequilibrium phenomena, called prethermalization, have been predicted to occur in quantum as well as classical many-body systems [59][60][61] and observed in Gring et al [24].
In particular, integrable one dimensional systems are known not to thermalize and indeed, experiments in Ref. [47] have observed an exceedingly long equilibration time. Yet even in this extreme case, we suggest in this section that many-body one dimensional systems can reach disguised "thermalized" states through prethermalization phenomena within a short time.
Ramsey dynamics and dynamics of interference patterns between split condensates described in the previous sections are particular examples of dynamics in which slow equilibration is expected because the system essentially consists of uncoupled harmonic oscillators in the low energy description (see Sec.IV). In the following, we give a heuristic argument that the distribution of the interference contrast amplitudes of the two non-equilibrium quasi-condensates are given by that of two equilibrium quasi-condensates at some effective temperature T eff . We give more details in Sec. VI C.
Long time after the splitting, the position of the interference peak becomes completely random, and yields no information. Therefore, we focus on the squared transverse magnitude of the spin Ŝ ⊥ l 2 , or equivalently, the interference contrastsĈ 2 . In the following, we describe the physics for the dynamics of split condensates, but the same argument can be applied to Ramsey dynamics. The interference contrast of the split condensates after a long time is determined by the "average" phase fluctuations present in the system. In the dephased limit, such fluctuations are determined by the total energy present in each mode labeled by momenta k. Now for sufficiently fast splitting, the energy E split contained in each mode is independent of momenta because the density difference of quasi-condensates along the axial direction is uncorrelated in the initial state beyond the spin healing length ξ s (see discussion below Eq. (10)). On the other hand, the interference contrast of thermal condensates at temperature T is determined by the thermal phase fluctua-tions caused by excitations whose energy is distributed according to equipartition theorem; each mode in the thermal condensates contains the equal energy of k B T . Thus from this argument, we find that the interference contrast of split condensates after a long time becomes indistinguishable from the one resulting from thermal condensates at temperature k B T eff = E split . We will show in Sec.VI C that in fact, the full distribution function of interference contrast becomes indistinguishable for these two states. In the case of splitting two dimensional condensates, equipartition of energy and existence of "nonequilibrium temperature" was pointed out in Ref. [61].
Here we propose the occurrence of such prethermalization within Tomonaga-Luttinger theory. We emphasize that within Tomonaga-Luttinger theory of low energy excitations, different modes are decoupled and therefore no true thermalization can take place. In realistic experimental situations, such integrability can be broken and relaxation and thermalization process are expected to occur after a long time dynamics. The requirement to observe the prethermalization phenomena predicted in this theory in experiments depends on the time-scale of other possible thermalization processes we did not consider in our model such as the effective threebody collisions [62,63], relaxation of high energy quasiparticles [64], or interactions among the collective modes through anharmonic terms we neglected in Tomonaga-Luttinger theory [65]. When all these processes occur at much slower time scale than the prethermalization timescale, which is roughly given by the integration length divided by the spin sound velocity ∼ l/c s for decoupling case, the observation of prethermalization should be possible. In one dimension, the dynamics is strongly constrained due to the conservation of energy and momentum, and therefore it is likely that the dynamics is dominated by the modes described by Tomonaga-Luttinger theory for long time for quasi-condensates with low initial temperatures.
Such prethermalizations are expected to occur even in higher dimensional systems [19,25,61,66]. We note that the conditions for the experimental observations of the phenomena might be more stringent because true thermalization processes are expected to take place much more quickly in two and three dimensional systems.

III. TWO COMPONENT BOSE MIXTURES IN ONE DIMENSION: HAMILTONIAN
In this paper, we study the dynamics of two-component Bose mixtures in one dimension through Tomonaga-Luttinger formalism [41,42]. As we have stated before, we assume the rotating frame, in which spin-up and spindown particles have the same chemical potential in the absence of interactions. The Hamiltonian of two compo-nent Bose mixtures in one dimension is given by Here ψ i with i =↑, ↓ describe two atomic species with masses m i and g ij are the interaction strengths given by [67] where ν ⊥ is the frequency of transverse confinement potential and a ij are the scattering lengths between spin i and j. System size is taken to be L, and we take the periodic boundary condition throughout the paper. In addition, we use the units in which = 1.
In the low energy description, the Hamiltonians for weakly interacting bosons after the initial π/2 rotation can be written in quadratic form, and given by where ρ is the average density of each species andn σ are variables representing the phase and density fluctuation for the particle with spin σ. These variables obey a canonical commutation relation [n σ (r),φ σ (r )] = −iδ(r − r ). In the Hamiltonian above, we included the kinetic interaction term g φ ↑↓ , which is zero for weakly interacting bosons, but allowed by inversion symmetry and non-zero for generic Tomonaga-Luttinger Hamiltonians.
We note that in the weakly interacting case, one can obtain the parameters of the Hamiltonian in Eq. (4) such as g ij and m i through hydrodynamic linearization of the microscopic Hamiltonian. In this case, we assume the small fluctuations of the densitiesn σ and phases ∇φ σ and expand the Hamiltonian in Eq. (3) to the second order in these variables through the expression ψ † σ ∼ √ ρ +n σ e iφσ , resulting in the form of the Hamiltonian in Eq. (4). Due to this assumptions of small spatial variations of the phase,n σ andφ σ represent the "coarsegrained" variables where collective modes have linear dispersions. The Hamiltonian of Tomonaga-Luttinger theory in Eq. (4) can also describe effective low-energy physics of strongly interacting systems, but in this case there is no simple relation between microscopic parameters and the parameters of the Hamiltonian in Eq. (4).
In order to describe the spin dynamics, we define spin and charge operators as the difference and the sum of spin up and down operators, i.e.φ s =φ ↑ −φ ↓ ,φ c =φ ↑ +φ ↓ , n s = 1 2 (n ↑ −n ↓ ),n c = 1 2 (n ↑ +n ↓ ). In this representation, the Hamiltonian in Eq. (4) becomes where interaction strengths are given by . The masses are given by the relations The spin variablesφ s andn s are "coarse-grained" in the sense that they represent the operators in the long wavelength beyond the spin healing length ξ s . ξ s is determined from microscopic physics and gives the length below which the kinetic energy of spins wins over the interaction energy, see Eq. (5) above. For weakly interacting bosons with m ↑ = m ↓ = m, it is given by ξ s = π/ √ mρg s .
In the following, we assume that the number of particles within the spin healing length is large, i.e. ξ s ρ 1. This condition is always satisfied for weakly interacting bosons.
In the next section, we consider the case g mix = 0 and g φ mix = 0, in which spin and charge degree of freedom decouple. Then the dynamics of spins is completely described by the spin Hamiltonian in Eq. (5). The general case in which g mix = 0 and g φ mix = 0 will be treated in Sec V.

A. Hamiltonian and initial state
The experiment of Widera et al. [22] used F = 1, m F = +1 and F = 2, m F = −1 states of 87 Rb for spin-up and spin-down particles, respectively. These hyperfine states have the scattering lengths a σσ such that a ↑↑ ≈ a ↓↓ . Consequently, the mixing Hamiltonian in Eq.(7) approximately vanishes for weak interactions. Motivated by this experiment, here we consider the decoupling of spin and charge degrees of freedom [52]. Spin dynamics in this case is completely determined by the spin Hamiltonian where K s is the spin Luttinger parameter representing the strength of interactions, and c s is spin sound velocity.
Other spin variables can be similarly defined in terms of coarse grained spin variablesn s andφ s . In the following, we consider the general transverse spins pointing in the direction (x, y, z) = (cos θ, sin θ, 0) integrated over l given bŷ where σ a with a = x, y are Pauli matrices. HereŜ θ l with θ = 0 corresponds to spin x operator and θ = π/2 corresponds to spin y operator. In order to explore the one dimensional dynamics resulting from Hamiltonian in Eq. (8), we analytically compute the mth moment of the spin operatorŜ θ l , Ŝ θ l m , after time t of the π/2 pulse of the Ramsey sequence. Then, the full distribution functions ofŜ x l andŜ y l , as well as the joint distribution of these will be obtained from Ŝ θ l m . In order to study the dynamics of Ramsey interferometer in terms of low energy variablesn s andφ s , we need to write down an appropriate state after the π/2 pulse in terms ofn s andφ s . If pulse is sufficiently strong, each spin is independently rotated into x direction after the π/2 pulse. Naively, this prepares the initial state in the eigenstate ofŜ x (r) = ρ cosφ s (r) with eigenvaluê φ s (r) = 0. However, due to the commutation relation betweenn s andφ s , such an initial state has an infinite fluctuation inn s and therefore, the state has an infinite energy according to Eq.(8). This unphysical consequence comes about because the low energy theory in Eq. (8) should not be applied to the physics of short time scale given by 1/E c where E c is the high energy cutoff of Tomonaga-Luttinger theory. During this short time dynamics, the initial state establishes the correlation at the length scale of spin healing length ξ s . The state after this short time dynamics can now be described in terms of the coarse-grained variablesn s (r) andφ s (r). The variableŝ n s (r) andφ s (r) are defined on the length scale larger than the spin healing length ξ s . Since the z component of spins are still uncorrelated beyond ξ s after the initial short time dynamics, the appropriate initial condition of the state is written as where the delta function δ(r − r ) should be understood as a smeared delta function over the scale of ξ s . Because the state after the short time dynamics is still close to the eigenstate of theŜ x (r) operator, spins are equal superpositions of spin-up and spin-down. Then the distribution ofŜ z l = l 0Ŝ z (r)dr is determined through random picking of the values ±1/2 for 2ρl particles. Due to the central limit theorem, the distribution ofŜ z l = l 0Ŝ = ρlη/2, which determines the magnitude of the fluctuation forŜ z (r) in Eq.(10). In Eq.(10), we also introduced the phenomenological parameter η which accounts for the decrease and increase of fluctuations coming from, for example, imperfections of π/2 pulse. The ideal, fast application of π/2 pulse corresponds to η = 1.
In the experimental realization of Ref. [22], η was determined to be between 0.8 and 1.3 through the fitting of the experiment with Tomonaga-Luttinger theory for the time evolution of the average x component of the spin, Ŝ x l . Through engineering of the initial state such as the application of a weak π/2 pulse, η can also be made intentionally smaller than 1.
A convenient basis to describe the initial state of the dynamics above is the basis that diagonalizes the spin Hamiltonian of Eq. (8). The phase and density of the spinsφ s (r) andn s (r) can be written in terms of the creation b † s,k and annihilation b s,k operators of elementary excitations for the spin Hamiltonian in Eq.(8) aŝ where we definedφ s,k andn s,k to be the Fourier transform of operatorsφ s (r) andn s (r). b † s,k creates a collective mode with momentum k and follows a canonical commutation relation [b s,k , b † s,k ] = 1. Note that k = 0 mode has no kinetic energy, and it naturally has different evolution from k = 0 modes. The Gaussian state determined by Eq.(10) takes the form of a squeezed state of operators b s,k , and it is given where 2W k = 1−α k 1+α k , α k = |k|Ks πρη . Here the state |n s,0 is the normalized eigenstate of the operatorn s,0 with eigenvalue n s,0 . The summation of k in the exponent has a ultraviolet cutoff around k c = 2π/ξ s . N is the overall normalization of the state. It is easy to check that ψ 0 |n s,kns,k |ψ 0 = ρη 2 δ k,−k , which corresponds to Eq. (10).

B. Moments and full distribution functions of spins
After free evolution of the initial state |ψ 0 for time t, the state becomes |ψ(t) = e −iHst |ψ 0 . We characterize the state at time t by the mth moments of spin operators, (Ŝ θ l ) m . As we will see below, the full distribution function can be constructed from the expression of (Ŝ θ l ) m [18]. We consider the evaluation of moments (Ŝ θ l ) m at time t, |ψ(t) . Each momentum k component of the initial state |ψ 0 independently evolves in time. Since k = 0 mode has a distinct evolution from other k = 0 modes, we separately consider k = 0 and k = 0 modes.

k = 0 mode
The Hamiltonian of k = 0 mode is given by H s,k=0 = πcs 2Ksn 2 s,0 in Eq. (11). Therefore, in the basis of n s,0 , k = 0 part of the state |ψ(t) is given by where N k=0 is the normalization of the state. The initial Gaussian state ofn s,0 stays Gaussian at all times, and any analytic operator of φ s,0 and n s,0 can be exactly evaluated through Wick's theorem. For example, k = 0 part contributes to the decay of the average of the x component of spin Ŝ This diffusion of the spin from k = 0 contribution is generally present in any dimensional systems, and not particular to one dimension. Physical origin of this diffusion is the interaction dependent on the total spin ,Ŝ 2 z . The eigenstate ofŜ x with eigenvalue ρl is the superposition of different eigenstates ofŜ z with eigenvalues m z , and they accumulate different phases e −itm 2 z in time. This leads to the decay of Ŝ x . In the thermodynamic limit L → ∞, the uncertainty ofŜ z becomes diminishingly small, and therefore, the decay of Ŝ x coming from k = 0 goes to zero. More interesting physics peculiar to one dimensional systems comes from k = 0 modes. In the case of three dimensional systems, macroscopic occupancy of a single particle state is absent in one dimension, so k = 0 momentum excitations have much more significant effect in one dimensional dynamics.

k = 0 contribution
The exact evaluation of spin moments (Ŝ θ l ) m for k = 0 is possible through the following trick. Consider the annihilation operator γ s,k (t) for the state |ψ(t) such that γ s,k (t) |ψ(t) = 0. If we write the operatorsφ s (r) in terms of γ s,k (t) and γ † s,k (t), then k = 0 part of the mth moment schematically takes the form (Ŝ θ l ) m ∼ exp(i k =0 C s,k γ s,k + C * s,k γ † s,k ) (Here and in the following, we drop the time dependence of γ s,k (t) from the notation). Using the property e γ s,k |ψ(t) = (1 + γ s,k + · · · ) |ψ(t) = |ψ(t) and the identity e A+B = e A e B e − k =0 |C s,k | 2 ). It is straightforward to check that γ s,k operator is given by the linear combination of b s,k and b † s,−k as follows, (15) γ s,k and γ † s,k obey a canonical commutation relation [γ s,k , γ † s,k ] = 1. In terms of these γ s,k , the expression ofφ s,k (t) becomes C s,k (t) measures the fluctuation, or variance, of phase in the kth mode at time t, given by |φ s,k (t)| 2 = φ s,k (t)φ s,−k (t) . Indeed, since γ s,k is the annihilation operator of our state at time t, we immediately conclude that |φ s,k (t)| 2 /L = |C s,k (t)| 2 .
Using the technique described above, mth moment of S θ l becomes (we include both k = 0 and k = 0 contributions in the expression below) where ξ . . + s m e ikrm ). s i takes either the value 1 or −1, and {si} sums over all possible set of values. Note that L is the total system size and l is the integration range.

Full Distribution Function
Calculation of the full distribution functions from moments in Eq. (17) is studied by the techniques introduced in Ref. [18] through mapping to the statistics of random surfaces. In this subsection, we provide the details of the calculation.
Eq. (17) is simplified if the integrations for each r i can be independently carried out. This is not possible in Eq. (17)  . To unentangle this, we introduce Hubbard-Stratonovich transforma- We apply a similar transformation for Imξ s,k . This removes the cross term between e ikri and e ikrj for i = j and allows us to independently integrate over r i 's. Associated with each transformation, we introduce auxiliary variables λ 1sk for Re(ξ s,k ), λ 2sk for Im(ξ s,k ). Then, mth moment becomes where we introduced ξ {ri} s,k = |φ s,k (t)| 2 L e ikri . Summation over {s i = ±1} can now be carried out. Furthermore, we introduce a new variables λ rsk and λ θsk , and replace λ 1sk and λ 2sk through the relation λ rsk = λ 2 1sk + λ 2 2sk and cos(λ θsk ) = λ 2sk / λ 2 1sk + λ 2 2sk . These operations result in the simplified expression, where with α k = |k|Ks πρη . The integration over λ rsk and λ θsk in Eq.(18) extends from 0 to ∞ and from −π to π, respectively.
Comparing the expression in Eq. (18) and the implicit definition of a distribution function in Eq.(2), it is easy to identify the distribution function as This function can be numerically evaluated through Monte Carlo method with weight λ rsk e −λ 2 rsk /2 for λ rsk and equal unity weight for λ θsk .
While we have assumed that the chemical potentials of spin-up and spin-down atoms are the same in the absence of interactions by going to the rotating frame, it is easy to obtain the expression for distribution functions in the lab frame. The energy difference E between spin-up and spin-down atoms results in the rotation of the spin in the x − y plane at a constant angular velocity E. Therefore, the distribution in the lab frame is obtained by replacing θ → θ + Et in Eq. (21).
In this section, we have focused on the distribution function of spins in x − y plane, but it is also possible to obtain the distribution function of z component of the spin, and we present the result in the Appendix A2.

Joint Distribution Function
From the expression for the spin operators in Eq. (9), we observe that the spin operators for the x and y directions commute in the low energy description. This is because spin operators in Tomonaga-Luttinger theory are coarse-grained over ∼ ρξ s particles, and since for weak interactions ρξ s 1, the uncertainty of measurements coming from non-commutativity ofŜ x l and S y l becomes suppressed. The possibility of simultaneous measurements of spin x and y operators implies the existence of joint distribution functions P x,y l (α, β), where P x,y l (α, β)dαdβ is the probability that the simultaneous measurements ofŜ x l and S y l give the values between α and α + dα, and β and β + dβ, respectively. Here we provide the expression for P x,y l (α, β) and proves that this is indeed the unique solution.
The joint distribution function P x,y l (α, β) is given by the following expression where the expression for χ(r, {λ jsk }) is given in Eq. (19).
To prove it, we first show that Eq.(22) reproduces the distribution function P θ l (α) in Eq.(21) for all θ. Then, we show that a function with this property is unique, and therefore the expression in Eq.(21) is necessarily the joint distribution function.
Given a joint distribution function P x,y l (α, β), we can determine the distribution function P θ l (γ) of a spin pointing in the direction (cos θ, sin θ, 0). Consider the spin S in the x−y plane with S = (α, β, 0) whose probability distribution is given by P x,y l (α, β). The projection of the spin S onto the axis pointing in the direction (cos θ, sin θ, 0) is given by |S| cos(φ − θ) where |S| = α 2 + β 2 is the magnitude of spin and φ is the angle Arg(α + iβ). After a simple algebra, we find |S| cos(φ − θ) = α cos θ + β sin θ. Then given a spin S = (α, β, 0), if one measures the spin along the direction (cos θ, sin θ, 0), the measurement result gives γ if and only if γ = α cos θ + β sin θ. From this consideration, the probability distribution that the measurement along the direction (cos θ, sin θ, 0) gives the value γ is given by Now, if we plug in the expression of Eq. (22) in Eq. (23), we see that P θ l (γ) agrees with Eq. (21) for all θ. Now we prove the uniqueness of a function with the above property, i.e. a function that reproduces Eq. (21) through the relation Eq.(23). Suppose you have another distributionP x,y l (α, β) that satisfies Eq.(23) for all θ. We define Q(α, β) = P x,y l (α, β) −P x,y l (α, β). Our goal is to show that Q(α, β) must be equal to zero. By definition, we have the equality for all θ and γ. If we take the Fourier transform of both sides of Eq. (24) in terms of γ, we obtain In the last line, we defined w = w(cos θ, sin θ) and α = (α, β). Notice that this equation holds for any w. Then this last expression is just like (two-dimensional) Fourier transform of Q. By taking the inverse fourier transform of the last expression in terms of w, we find thereby proving the uniqueness of the joint distribution P x,y l (α, β). From the joint distribution function in Eq. (22), one can also obtain other distributions, such as the distribution P ⊥ l (γ) of the square of the transverse spin magnitude, S ⊥ l (t) 2 , which is given by

Initial state
FIG . 6: The illustration of the dynamics for each harmonic oscillator mode, described by the Hamiltonian Eq.(8). The initial state contains a large fluctuation of density differencê n s,k given by n s,kns,−k = ηρ/2(see Eq. (10)), and its conjugate variable, the phase differenceφ s,k , has a small fluctuations. In the subsequent dynamics, such squeezed state evolves and energy oscillates between the fluctuations of the density difference and phase difference.

Interpretation of the distribution dynamics
The form of the distribution function in Eq. (22) encapsulates the interpretation in terms of dynamics originating from spin waves explained in Sec II. Here e iχ(r,{t jsk }) represents the spin direction at coordinate r, where the x − y plane of the spin component is taken to be a complex plane. Then Eq. (22) suggests that for a given instance of the set {λ jsk }, (S x l + iS y l ) is simply the sum of the local spin directions e iχ(r,{λ jsk }) over the integration length l. The local spin direction at position r are determined by the phase χ(r, {t jsk }), which receives contributions from each spin wave of momentum k with strength A k (t) = λ rsk |φ s,k (t)| 2 . Spin waves with momenta k rotate the spins as sin(kr + λ θsk ) (see the expression of χ(r, {λ jsk }) in Eq. (19)). The rotation strength A k (t) ∝ λ rsk has the distribution λ rsk e −λ 2 rsk /2 , which represents the quantum fluctuation of the spins. On the other hand, λ θsk is distributed uniformly between −π and π.
The dynamics of phase fluctuations |φ s,k (t)| 2 can in fact be easily understood by considering the Hamiltonian given by Eq.(8) as a harmonic oscillator for each k (Fig. 6). We first note that the initial state has a large fluctuation of densityn s,k because the initial π/2 pulse prepares the state in the (almost) eigenstate of S x = ρ cos(φ s,k ) with a small fluctuation ofφ s,k , the conjugate variable ofn s,k . The fluctuation ofn s,k is given by n s,kns,−k = ηρ/2 (see Eq. (10)). Because of this large fluctuation in the density, almost all the energy of the initial state is stored in the interaction term |n s,k | 2 in Eq. (8). Therefore, the total energy of each harmonic oscillator can be estimated as πcsρη 4Ks . During the dynamics dictated by the harmonic oscillator Hamiltonian, this energy oscillates between the density fluctuations and phase fluctuations in a sinusoidal fashion, see Fig. 6. In the dephased limit of the dynamics, approximately equal energy of the system is distributed to the phase and density fluctuations, and from the conservation of energy, we conclude that the characteristic magnitude of phase fluctuation is given by |φ s,k (t)| 2 ∼ π 2 ρη 4K 2 s k 2 . Such 1/k 2 dependence of |φ s,k (t)| 2 agrees with the more rigorous result in Eq. (20). Therefore the spin fluctuations dominantly come from spin waves with long wavelengths, as we have stated in Section II. Moreover, the weak dependence of spin dynamics on high momenta contributions justifies the use of Tomonaga-Luttinger theory for describing the dynamics. We will more carefully analyze the dependence of distributions on the high momentum cutoff in Sec IV D.
From the simple argument above, it is also clear that the spin fluctuations coming from spin waves with momenta k have the time scales associated with the harmonic oscillators given by 1 |k|cs . Again, this rough argument agrees with the more rigorous result presented in Eq. (20). Therefore, the fast dynamics is dominated by spin waves with high momenta and slow dynamics is dominated by low momenta. These considerations lead to the illustrative picture of Fig.2. Furthermore, this implies that the dynamics of the magnitude of spin Ŝ ⊥ l 2 reaches a steady state around the time l 4cs since spin waves with wavelength longer than l do not affect the magnitude. This should be contrasted with the evolution of the x component of spin which, in principle, keeps evolving until the time scale of ∼ L 4cs (see Fig.4). The strength of interactions and correlations are associated with Luttinger parameter, K s . K s influences the spin fluctuations |φ s,k (t)| 2 at all wavelength, and |φ s,k (t)| 2 depends on K s as 1/K 2 s for a fixed density. As is expected, in the limit of the weak interaction corresponding to large K s , the amplitude of spin fluctuation decreases. In Fig. 7 In order to illustrate the dynamics of the Ramsey sequence further, it is helpful to study the dynamics of the expectation value of the squared transverse magnitude, In Fig.8, we plot the evolution of Ŝ⊥ l (t) 2 with K s = 20, L/ξ s = 200 and l/ξ s = 20, 30, 40. We also plot- with the same parameters. It is easy to verify that S x l (t) is independent of integration length l [68]. As we have discussed in the previous section, Ŝ⊥ l (t)  It is interesting to ask if the long time limit of Ŝ⊥ l (t) 2 for sufficiently large integration length l attains the value which corresponds to the one expected from the randomization of spin patches of size ξ s . At low energies, spins within the length ∼ ξ s are aligned in the same direction, but spin waves can randomize the direction of the spin for each of l/ξ s patches. Since the magnitude of spin within ξ s is ξ s ρ, if the patches are completely randomized, the result of the random walk predicts that Ŝ⊥ l 2 ∼ (ξ s ρ) 2 (l/ξ s ). We will see below that, due to the properties of correlations in one dimension, the integrated magnitude of spin Ŝ⊥ l 2 never attains this form, albeit a similar expression is obtained (see Eq. (26)). Moreover, we identify the integration length l which separates the "spin diffusion regime" and the "spin decay regime" by finding the decaying length scale for Ŝ⊥ The results for the long time limit of Ŝ ⊥ l (t) 2 can be analytically computed. Following similar steps leading to Eq.(17), we find Here ξ {ri} s,k = |C s,k |(e ikr1 − e ikr2 ). We introduce dimensionless variables r i = r i /l, k = kl and the integration over k in the exponent can be carried out as In the second line, we approximated cos 2 (|k|c s t) ≈ sin 2 (|k|c s t) ≈ 1/2, which is appropriate for long time.
In the last line, we extended the upper limit of the integration for the second term to ∞ and the lower limit to 0. The former is justified because we know that high momentum contribution is suppressed by 1/k 2 , and the latter is justfied because we also know low momenta excitations with wavelengths larger than l do not affectŜ ⊥ l . Since ∞ 0 dy sin 2 (y) y 2 = π/2, we find, in the long time limit, where we expressed the result as a ratio of the asymptotic value and the value at shortest time scale of the theory given by t ∼ 1/µ. l 0 = π 2 ηρl . Therefore, l 0 ≈ 1 separates the "spin diffusion regime" and the "spin decay regime." An intuition behind the expression for l 0 can be explained through the following heuristic argument. The system enters the spin decay regime when the spins within the integration length l rotates by 2π across l. The angle difference between the spins at r = 0 and r = l in the long time limit is roughly given by ∆χ = is the characteristic magnitude of |φ s,k (t)| 2 in Eq. (20), which is given by the half of the maximum magnitude of |φ s,k (t)| 2 . Now the expectation of magnitude (∆χ) 2 over the quantum fluctuations represented by λ rsk can be computed, and it yields (∆χ) 2 ≈ π 2 ηρl 4K 2 s . When (∆χ) 2 becomes of the order of 1, the system enters the spin decay regime. This estimate gives the boundary between the two regimes l 0 = π 2 ηρl 8K 2 s ≈ 1 apart from an unimportant numerical factor.
It is notable that the Eq. (26) approaches the random walk behavior ∝ (ξ s ρ) 2 (l/ξ s ) very slowly, i.e. in an algebraic fashion. Therefore, even in the steady state, the system retains a strong correlation among spins. Moreover, Eq. (26) in the limit of l 0 → ∞ is not just the random walk value, but is proportional to K s , which measures the strength of fluctuations.
The calculation above shows that the spin diffusion regime and the spin decay regime are separated at the integration length scale ofl ≈ 8K 2 s π 2 ηρ . This length scale is nothing but the correlation length of spins in the long time limit. The calculation of the spin correlation length, for example, between S x (r) andŜ x (r ) can be done similarly to the calculation of Ŝ ⊥ l 2 . The result in the long time limit is where C = e −kc/(4πρη) is a small reduction of the spins due to the contributions from high energy sector. Thus, one expects qualitatively different behaviors of distribution functions for integration lengths l <l and l >l.

D. Momentum cut-off dependence
The description of dynamics presented above uses the low energy effective theory. In order to confirm the selfconsistency of our approach, we check that the distributions of spins are not strongly affected by high energy physics, i.e. they weakly depend on high momentum cutoff. We have seen an indication that this is indeed the case through the weak fluctuations of phases for large k, |φ s,k | ∝ 1/k 2 , in Sec IV B 5.
First of all, we analyze the high momentum cut-off k c ∼ 2π/ξ s dependence of the average value ofŜ x l . From the discussion in Sec IV B, it is straightforward to obtain that(here we ignore k = 0 contribution) where in the last line, we took the long time limit t ξ s /c s [68]. In this limit, only the first term in Eq.(28) depends on the cutoff k c , and moreover, the cutoff dependence is independent of time. The effect is to reduce the value of Ŝ x l through the multiplication of a number close to one in the weakly interacting limit. For example, if we take k c = 2π/ξ c , then the cut-off dependent term reduces the value by multiplying exp − kc 4πρη ≈ e −1/(4Ks) ≈ 1. In a similar fashion, higher moments of spin operators can be shown to have a weak dependence on the cutoff momentum k c , as long as the integration length is much larger than the healing length, l/ξ s 1. In this limit, m moments of, for example,Ŝ x l is reduced by exp −m kc 4πρη . Therefore, the full distribution function is simply reduced by the multiplication of a number close to one exp − kc 4πρη ≈ e −1/(4Ks) ≈ 1 in the weakly interacting regime. This gives the self-consistency check of our results in Sec IV B

V. DYNAMICS OF FULL DISTRIBUTION FUNCTION IN THE PRESENCE OF MIXING BETWEEN SPIN AND CHARGE DEGREES OF FREEDOM
In this section, we extend the analysis in Sec IV to a more general case, in which spin and charge degrees of freedom mix. We will see that the distribution functions even for this more general case have essentially the same structure as in Eq. (21), and are described by spin waves with fluctuations whose amplitude is determined by the fluctuations of phase |φ s,k | 2 . One important difference from the decoupling case is the dependence of spin distributions on the initial temperature of the system. The thermal excitations are present in the charge degrees of freedom in the initial state, and such thermal fluctuations increase the value of |φ s,k | 2 through the coupling between spin and charge during the evolution.

A. Hamiltonian and initial state
In a generic system of two component bosons in one dimension, spin and charge degrees of freedom couple through the mixing Hamiltonian in Eq. (7). Yet, Hamiltonian in Eq.(4) is still quadratic and it can be diagonalized. We define new operatorsφ 1 ,φ 2 ,n 1 ,n 2 by Mixing angle κ and scaling parameter s c are chosen so that the Hamiltonian is written in the following diagonal form Explicitly, κ and s c are given by where ± in the expression of tan κ is + when κ 0 > 0 and − when κ 0 < 0. We defined κ such that κ = 0 corresponds to decoupling of charge and spin, i.e. to g mix = 0 and g φ mix = 0 in Eq.(7). Parameters g 1 , g 2 , ρ 2m1 and ρ 2m2 are given by In the weakly interacting systems which we study in this paper, Luttinger parameters K i and sound velocities c i are determined for each Hamiltonian H i , i =↑, ↓, c, s, 1, 2 through At finite temperature, the state before the first π/2 pulse contains excitations, and these excitations are carried over to the charge degrees of freedom after the pulse.
Pulse only acts on the spin degrees of freedom, and the local sum density of spin-up and down is left untouched as long as the pulse is applied in a short time compared to the inverse of typical excitation energies, β = 1/(k B T ). In other words, the local density fluctuation of spinup,n ↑ (r), before π/2 pulse is converted to the sum of the local density fluctuation of spin-up and spin-down, n ↑ (r) +n ↓ (r) after π/2 pulse. In this strong pulse limit, then, the distribution ofn ↑ (r) before π/2 pulse is the same as the distribution ofn ↑ (r) +n ↓ (r) after π/2 pulse.
The distribution of the local density for spin-up atoms before π/2 pulse is determined by the density matrix for spin-up given by e −βH ↑ , where in the weak interaction regime we have (see Eq. (4)) Then, the density matrix which produces the distribution ofn ↑ (r) +n ↓ (r) required above is given by e −βH c↑ where where K c↑ = π 4 ρ m ↑ g ↑↑ and c c↑ = 2ρg ↑↑ m ↑ . The initial state for spins is determined by the π/2 pulse, and we obtained the state in Eq. (12). Then, the complete initial density matrix after the first π/2 pulse is given bŷ This density matrix evolves in time asρ(t) = e −itHρ 0 e itH . Since we assume that the preparation of the initial state is done through a strong, short pulse, the spin and charge degrees of freedom are unentangled in the initial state.

B. Time evolutions of operators
In order to calculate the distribution function ofŜ θ l , we again start from the calculation of mth moments, Tr ρ(t) Ŝ θ l m . Evaluation of moments can be done through a similar technique used in Sec IV B.
In the following, we describe convenient, timedependent operators γ s,k (t) and γ c,k (t) used to evaluate spin operators such as e iφ s,k . The first operator resides in the spin sector and it is again the annihilation operator of the initial spin state such that Trγ s,k (0)ρ 0 = 0. This operator is given in Eq. (15), which is with 2W k = 1−α k 1+α k and α k = |k|Ks πρη as before. The second operator is the operator of charge degrees of freedom, and it is given by where b c↑,k is an annihilation operator for the elementary excitations in H c↑ . Since γ s,k (t) and γ c,k (t) commute at t = 0, they commute at any time t. We will drop the time dependence of γ a,k (t) in the notation from now on.
From the expression of initial density matrixρ 0 in Eq. (34), it is easy to check that the density matrix at time t given byρ(t) = e −itHρ 0 e itH can be written as the tensor product of the density matrix of operators γ s,k (t) and that of γ c,k (t). This is becauseρ 0 is a tensor product of the density matrices of γ s,k (t = 0) and that of γ c,k (t = 0). This structure of the density matrices at time t allows the independent evaluation of γ s,k (t) and γ c,k (t) operators, and it is advantageous to express spin operators in terms of these operators. As we show in the Appendix B, we can writeφ s,k in terms of γ c,k (t) and γ s,k (t) as follows.
where explicit expression of C s,k and C c,k are given by Using Eq.(35), we find an expression for (Ŝ θ l ) m in terms of γ a,k with a = s, c as follows where ξ a,k = ( 2. k = 0, spin sector This calculation is analogous to Eq. (17) and the result can be directly read off from Eq. (17), and it is 3. k = 0, charge sector We first rewrite the density matrix at time t aŝ Then the trace of Ŝ θ l m for k = 0 spin sector is We can evaluate this by taking the trace in the basis of normalized coherent states |α k such that γ c,k |α k = α k |α k . The use of the identity 1 = 1 π d 2 α k |α k α k | as well as of an important equality e va † a =: e (e v −1)a † a : [18], where : O : is a normal ordering of O, leads to

Full distribution function
We can summarize the results above as Here, M c,k = 1+e −β|k|c c↑ 1−e −β|k|c c↑ . As before, we introduce the auxiliary variables to separate spatial integrations over r i . We can combine ξ * s,k ξ s,k and M c,k ξ * c,k ξ c,k so that we only need to introduce three sets of variables, λ 1,s,k , λ 2,s,k , λ 0 for Hubbard-Stratonovich transformation. Summing over {s i } simplifies the result, leading to the following expression for the full distribution function The last line can be confirmed by directly computing |φ s,k | 2 , using the expression in Eq. (35). The expression for |φ s,0 | 2 is given by Eq. (37). As before, the joint distributions as well as the distributions of squared transverse magnitude can be obtained through the same procedure as in Sec IV B 4.
The spin distribution in the presence of mixing between spin and charge degrees of freedom resembles the one in the absence of such mixing, and the only change is the additional contributions to phase fluctuations coming from the thermal excitations, represented by 1+e −β|k|c c↑ 1−e −β|k|c c↑ |C c,k | 2 in Eq. (41). |C c,k | is proportional to sin 2 κ as one can see from Eq. (36). Thus, for weak coupling of κ ∼ 0, the contribution is diminished by a factor of κ 2 .
In the experiment by Widera et al., they used Rb 87 in the presence of Feshbach resonance. They employed the theory which assumes the absence of mixing between spin and charge degrees of freedom to analyze the decay of the Ramsey fringes. The ratio of interaction strengths in their experiment can be roughly estimated as g c : g s : g mix ≈ 3.66 : 0.34 : 0.06 which leads to the value of κ ≈ 2 × 10 −2 . Therefore, the thermal contributions are diminished by about four order of magnitude and thus, their assumption of decoupling between spin and charge is justified.
In Fig.9, we have plotted the evolution of the joint distribution functions for different strength of the coupling g mix at a relatively large initial temperature k B T = 0.4 × 2πc c↑ /ξ s where 2πc c↑ /ξ s is approximately the high energy cut-off of Tomonaga-Luttinger theory. Here we took the system size L/ξ s = 400, the Luttinger parameter K s = 20, integration length l = 20ξ s . For Fig.9 a), the ratio of interaction is taken to be g c : g s : g mix = 1 : 1 : 0.1, and for For Fig.9 b), g c : g s : g mix = 1 : 1 : 0.3. One can see that with increasing strength of mixing, the large initial temperature affects the spin dynamics at earlier time more strongly. For comparison, also see Fig. 4.

VI. INTERFERENCE OF TWO ONE-DIMENSIONAL CONDENSATES
A. Dynamics of interference pattern As we have described in Sec.II B, the full distribution of interference patterns can be studied in exactly the same way as we have studied the full distribution of spins in previous sections. In the following, we more formally describe the dynamics of split condensates.
The low energy effective Hamiltonian of two quasi-condensates after splitting is given by (∇φ L (r)) 2 + g 2 (n L (r)) 2 , where we assumed weakly interacting bosons with a possible density difference ρ R − ρ L = 0 between the two condensates. Here and in the following, we consider the rotating frame and ignore the chemical potential difference g/2(ρ 2 R − ρ 2 L ) between left and right condensates arising from interactions.
The interference pattern measures the phase differencê φ L −φ R . We describe the system in terms of the "spin" variables that are the difference of left and right condensates and "charge" variables that are the sum of the two. Using the variablesφ s =φ R −φ L ,φ c =φ R +φ L , n s = (n R −n L )/2,n c = (n R +n L )/2 , we find the Hamiltonian of the system to be Therefore, when the splitting makes two identical quasicondensates with equal density, "spin" and "charge" degrees of freedom decouple and we can use a simpler theory derived in Sec IV B. On the other hand, when the splitting makes two condensates with unequal densities, more general theory of Sec V needs to be employed. In any case, the full time evolution of the distributions of interference patterns can be obtained, which in principle can be compared with experiments. It is notable that the mixing of the "spin" and "charge" degrees of freedom for small density difference ρ R − ρ L is not "small," in the sense that the mixing angle κ defined in section V takes the maximum value π/4. The spin decoupling in the limit of ρ R − ρ L → 0 is recovered not by taking κ → 0, but rather, by taking the time at which the effect of the coupling takes place to infinity. This is most explicitly shown in Eq. (36) where the charge contributions of fluctuations go to zero as c 1 → c 2 which is attained in the limit ρ R − ρ L = 0.

B. Interference patterns in equilibrium
The techniques to calculate the full distribution functions presented in previous sections are directly applicable to also obtaining a simple form of the full distribution functions of the interference patterns between two independent, thermal quasi-condensates. This problem has been previously analyzed in theory [69,70] as well as in experiments [14,21].
We consider the preparation of two independent one dimensional quasi-condensates. If they are prepared by cooling two independent quasi-condensates, the temperature of the left quasi-condensate T L and that of right quasi-condensate T R are generically different. The density matrix of the initial state is described byρ 0 = e −(β L H L +β R H R ) where β a = 1/(k B T a ) with a = L, R. It is important to note that the constant shift of phase φ a → φ a + θ ac does not change the energy of the system, so that for the average over thermal ensemble one has to integrate over θ ac . Physically, this simply means that the phases of independent condensates are random. Then the only interesting distribution here is the distribution of the interference contrast [14,18,26,69] given by,Ĉ . The analysis of the evaluation of distributions in the density matrix of thermal equilibrium state in Sec V can be directly extended to this case, and we obtain the distribution where c a and K a , a = L, R are the sound velocity and Luttinger parameters of left and right quasi-condensate.

C. Prethermalization of interference patterns
In Sec.II C, we gave a heuristic argument for prethremalization phenomena, where the distribution of the interference contrast amplitudes of the two non-equilibrium quasi-condensates are given by that of two equilibrium quasi-condensates at some effective temperature T eff . We identified the effective temperature to be the energy stored in each momentum mode. In Sec. IV B 5, we found this energy to be πcsρη 4Ks , thus we conclude k B T eff = πcsρη 4Ks . In the following, we formally derive the result above, using the expressions of full distributions of interference patterns. The distribution of the interference contrast is determined by |φ s,k | 2 given in Eq. (20). In the long time limit, we can take sin 2 (c s |k|t) ∼ cos 2 (c s |k|t) ∼ 1/2. Moreover, since the interference contrast is most affected by the excitations with small wave vectors k with α k = |k|Ks πρη < 1, we can approximate the expression as On the other hand, for two quasi-condensates in thermal equilibrium, the position of the interference peaks is again random. The interference contrast is determined by |φ s,k | 2 given in Eq. (49). Since the main contribution to the fluctuation comes from low momenta, we approximate e −β|k|c ≈ 1−β|k|c. It is easy to check that the sound velocity and Luttinger parameters for each condensate is related to those of the difference mode (see Eqs. (44)(45)(46)) as c L = c R = c s and K L = K R = 2K s . Thus we obtain Now the crucial observation is that our closed form expressions for distributions of interference contrasts of both split quasi-condensates and thermal quasicondensates are determined solely by |φ s,k | 2 , and they take precisely the same form in terms of |φ s,k | 2 . Moreover, the expressions given by Eq. (50) and Eq.(51) have the same dependence on wave vectors |k|. Therefore, the full distribution of interference contrast of split condensates become indistinguishable from that of thermal condensates with temperature where the second equality holds for weakly interacting bosons and the chemical potential of one quasi-condesate is given by µ = gρ. Thus, split one dimensional quasicondensates indeed display the prethermalization phenomenon. In Fig. 10, we plot the interference contrast P ⊥ l (γ) (see Eq.(25)) of split condensates in a steady state at time t = 60ξ s /c s for Luttinger parameter K s = 20, system size L = 400ξ s and two different integration lengths l/ξ s = 30, 40. Also we plot the interference contrast of the thermal quasi-condensates (see Eq.(48)) at temperature πcs 2ξs for the same integration length. This temperature corresponds to the effective temperature obtained in Eq. (52). Indeed we see only a small difference between the distributions of steady states and thermal states for The distributions of interference contrast for steady states of split quasi-condensates and two thermal quasicondensates. Here x axis is scaled such that the maximum value of interference contrast is 1. For the split condensates, we plot the distribution at time t = 60ξs/cs for Luttinger parameter Ks = 20, system size L = 400ξs and two different integration length l/ξs = 30, 40. The thermal quasicondensates are for temperature πcs 2ξs for the same integration length corresponding to the effective temperature obtained in Eq. (52).
both integration lengths. The small difference comes from the approximations made in obtaining the expressions given by Eq. (50) and Eq. (51).
In the previous paragraphs, we assumed that the splitting prepares quasi-condensates with identical average densities. Here we briefly consider the case in which the splitting process prepares two quasi-condensates with slightly different densities. In this case, the temperature of the initial quasi-condensates affects the interference contrast around the time scale of ξs (c L −c R )π ≈ µ( √ ρ L − √ ρ R ) , whereas the prethermalized, long-time transient state is reached around µ l ξs , where l is the integration length.
These analytic arguments can be confirmed through numerical simulations. In Fig. 11, we have plotted the evolution of the interference contrastĈ 2 for system size L = 500ξ s , integration length l = 40ξ s , and Luttinger parameter K s = 20 with initial temperature corresponding to the chemical potential µ. Here we consider a situation where the density of left quasi-condensate is different from that of right quasi-condensate by 20% such that ρ R = 1.2ρ and ρ L = 0.8ρ where ρ R(L) is the average density of the right (left) condensate, and 2ρ is the average density of the initial condensate before splitting. Also for comparison we have plotted the magnitude of interference contrast of two quasi-condensates in thermal states at temperature given by µ. From the plot, one can observe the existence of quasi-steady state plateau after short time. Notice that the magnitude ofĈ 2 in the steady state is larger than the value expected from thermalized states at the initial temperature. The subsequent slow decrease of the interference contrast is due to the effect of temperature in the initial state coming from the small difference of the two quasi-condensates. Such development of the plateau at larger value of interference  11: The evolution of the interference contrastĈ 2 for system size L = 500ξs, integration length l = 40ξs, and the effective spin Luttinger parameter Ks = 20 with initial temperature corresponding to the chemical potential µ. Time is measured in units of ξs/cs. Here we took the density ρR = 1.2ρ and ρL = 0.8ρ. The magnitude ofĈ 2 for two thermal quasicondensates at temperature kBT = µ is plotted as red dotted line for comparison.
contrast than the one for equilibrium state of the initial temperature indicates the phenomenon of prethermalization.

VII. CONCLUSION
In this work, we have shown how noise captured by full distribution functions can be used to study the dynamics of many-body system in one dimension. The analytical results of joint distribution functions obtained in Sec IV B allow not only the simple understanding of distribution functions from spin-wave picture, but also an intuitive visualization of the correlation in one dimensional system. Using this picture, we have also shown that the phenomena of prethermalization occur in one dimensional dynamics. The thermal-like behaviors of the prethermalized state is revealed through the full distribution functions that contain information about the correlation functions of the arbitrary order. For the experimental demonstration of such prethermalizations, see Gring et al [24].
The approach developed in this paper can be extended to other types of dynamics. While we focused on Ramsey type dynamics or dynamics of interference patterns for a split quasi-condensate, we can also change different physical parameters to induce the dynamics. For example, it is straightforward to apply our study to the sudden change (quench) of interaction strength [40,71].
In this paper, we focused on distribution functions obtained from Tomonaga-Luttinger Hamiltonian (4). It is of interest to extend our analysis to higher spins [72], and analyze them, for example, in the presence of magnetic field [73]. Since there are more degrees of freedom in these systems, distributions might capture the tendency towards various phases such as ferromagnetic ordering.
Then, we can apply the trick introduced in Section IV B to obtain e λŜ z l = e λ 2 l/2 −l/2 dr1dr2 k =0 |d s,k | 2 e ik(r 1 −r 2 ) + where G ij (k) are the matrix elements of G(k).
Appendix C: k = 0 contribution in the presence of mixing between charge and spin In this appendix, we evaluate Tr e i( i si)φs,0/ √ L ρ(t) .
We first obtain the operator e i( i si)φs,0 after time evo-lution as e iHt e i( i si)φs,0 e −iHt = exp i i s i √ L (Aφ s,0 + A n s,0 + B √ s c φ c,0 + B √ s c n c,0 ) .