Spectroscopic evidence for engineered hadron formation in repulsive fermionic SU ( N ) Hubbard Models

,


Introduction
Non-equilibrium dynamics already provided us with a plethora of interesting phenomena, ranging from negative temperature states [1,2] through universal scaling across quantum phase transitions [3][4][5][6] to thermalization dynamics [7][8][9].In this context, recently there has been a large surge of interest in understanding and analyzing confinement in spin chains [10][11][12][13][14].Not only is it related to many body localization and slow entanglement dynamics [15][16][17], but confinement is also responsible for creating bound states of several interacting particles [18][19][20].Such objects play a prominent role in quantum chromodynamics and can shed light on particle formation in the early universe [21].
In addition to confinement, short range attractive interactions also give rise to many-particle bound states [22,23].These are also responsible for Cooper pairing [24] and can be engineered in ultracold atomic setting, which provide a fantastic platform to simulate and investigate strongly interacting forms of matter in a laboratory framework [25][26][27][28].The hyperfine spin of fermionic atoms, in particular, can play the role of quark colors or spins in condensed matter, and thereby enable experimentalists to emulate phenomena appearing in quantum chromodynamics [29][30][31] in the context of hadronic matter as well as simulating SU(N ) generalizations of condensed matter systems [32][33][34][35][36][37][38].Attractive fermions such as 6 Li, e.g., have been proposed to display color superfluidity and baryon formation at sufficiently low temperatures [32,34].
Fermion gases with attractive interaction are, however, quite rare, and rather unstable against three-body losses, especially in the regime of large scattering length [39].Moreover, cooling down a fermion gas is a notoriously hard problem [25,40].It is mostly for these reasons that the state proposed in Ref. [34] has not been observed so far.
Due to these complications, these states are not only extremely difficult to create in a controlled fashion, but also their detection is elusive.As we discuss here, the above problems can be circumvented and hadron formation naturally appears for fermions with strongly repulsive interactions, too, for a set of very simple initial states.We also demonstrate that by using quench spectroscopy, we can detect these newly formed particles unambiguously.This observation opens the door towards simulating and capturing hadron formation with appropriately engineered initial state in a much more stable cold atomic environment.Indeed, cold atomic systems such as 137 Yb or 87 Sr are today almost routinely used to realize stable repulsive SU(N ) gases [41][42][43][44].In an op- e Sketch of the entropy density as a function of energy density.The white dashed line indicates the J = 0 limit.The entropy density has a maximum at energy density ∞, where the temperature diverges and changes sign.For states with < ∞ the temperature is positive, while for states with > ∞ it becomes negative.
tical lattice, they realize the repulsive SU(N ) Hubbard model [45,46], describing fermions of N different colors, ψ † rα , moving around a lattice, and interacting locally via a colorindependent interaction, U , with n r = α ψ † rα ψ rα the number of fermions at site r.
The initial state which we propose to engineer is a product of N -fermion states [47,48], where Ξ denotes a subset of sites.As we show, this state, depicted for N = 3 in Fig. 1, necessarily evolves to a negative temperature state for any U > 0, but its dynamical properties change dramatically upon increasing the ratio, U/J, and for U/J 2 a strongly interacting quantum gas of baryons and mesons emerges.Similar initial states for N = 2 have been investigated in Ref. [49].

Results
Quench spectroscopy and resonances.The formation of composite particles can be most easily detected by quench spectroscopy [10], i.e., by analyzing the time evolution of the state (2) and its Fourier spectrum.We focused on a one-dimensional chain of SU (3) fermions, and performed non-Abelian Time Evolving Block Decimation (TEBD) simulations on it [50][51][52].On the left panel in Fig. 2, we show the time evolution of the probabilities p n of having n fermions at a site for an initial state, where three fermions have been placed at every third site.
For U/J 1, the probabilities p n (t) display damped oscillations, and the initial occupations p = 0 relax to some asymptotic values, p ∞ n .The Fourier spectra of p n (t) consist of a broad band, consistent with a band of width 4J +U of fermionic weakly interacting excitations.A careful analysis reveals an exponential relaxation towards a stationary state with a rate ∼ U 2 /J/ [52].
This picture changes radically for U/J 2, where the probabilities p n (t) are pinned roughly to their initial values, indicating that the probability that three particles stay together remains close to 1, p 3 (t) ≈ p (0) 3 .The time averaged probabilities of having n = 2 or n = 1 fermions at a site are suppressed, but they are finite, and almost equal.Even more strikingly, small, almost undamped oscillations decorate the p n (t) curves.The Fourier transform of these signals reveals high energy spectral features around ω ≈ 2U (see Fig. 2b).As we demonstrate, the observed oscillations of p n (t) can be understood as a result of quantum oscillations between baryonic trion states which transform coherently into mesonic doublon states and single fermion states, corresponding to quarks.
Composite particles.In the large U limit, we can treat the hopping J as a perturbation [53,54].Isolated n-particle states have an energy J, these composite particles behave as quite stable entities: although they collide with each other, they can decay and transform into each other only under the condition that energy, charge, and SU (N ) spins are all conserved.
In the N = 3 case, e.g., the collision of a baryon (trion) of energy E 3 ≈ E (0) 3 = 3U with three quarks (free fermions) of energy E 1 ∼ J ≈ 0 and the subsequent decay into three mesons (doublons) of energy E 2 ≈ E (0) 2 = U provides the lowest order decay channel for baryons [55] (see Fig. 1d).Since quarks have a very small concentration for large U , this event is very unlikely, and baryons behave as very stable composite particles.They can, however, virtually transform back and forth into a meson and a quark, and oscillate between these states at a frequency ω ≈ E 1 ) = 2U .This process gives rise to the oscillations observed in Fig. 2b.
Composite particles move on the lattice with a suppressed hopping.Simple perturbation theory can be used to determine the effective hopping J n of an n-particle composite, yielding ) in all dimensions (see Methods).Composite particles are therefore extremely heavy in the large U limit.In case of N = 3, the meson band (n = 2), is quite narrow compared to the 'quark' band of ψ α particles, but the baryon band (n = 3) is even narrower.The spectral peaks in Fig. 2b at ω ≈ 2U ± 2J can thus be simply understood as a result of van Hove singularities associated with the edges of the quark band (see Fig. 2c).
Effective theory.The emergent composite particles interact very strongly.In the large U limit, heavy N -ions dominate, and the density of other composite particles is suppressed.The energy of N -ions is increased compared to its J = 0 value by an amount of δE N ≈ zN J 2 (N −1)U due to quantum fluctuations, where a fermion of color α jumps to one of z neighboring lattice sites.Placing two N -ions next to each other suppresses these quantum fluctuations, and gives rise to an attractive interaction, Introducing the operator Φ † r , creating a dressed N -ion at site r, we arrive at the following effective Hamiltonian in the large U limit with n Φ r = Φ † r Φ r the N -ion number operator.Notice that for N odd, the Φ particles are fermions, while for N even they are hard core bosons.Also notice that |V /J N | ∼ (U/J) N −2 , implying that N -ion -N -ion interactions become strong for any N > 2 for U J. Therefore, in the U/J → ∞ limit, a gas of spinless, strongly interacting Φ particles is recovered.For N = 3, in particular, these particles are dynamically bound fermions, analogous to baryons in QCD.
Although the concentration of other particles is suppressed in the large U limit, their presence is still essential.In particular, the concentration p N −1 of (N −1)-ions is non-negligible, ∼ J 2 /U 2 .These composite particles, which transform according to the conjugate representation, are created by the operators Θ † ᾱ.Although somewhat lighter than N -ions, they are also very heavy, and their interaction with the N -ions is even stronger than the N -ion -N -ion interaction V itself: neighboring Θ α and Φ particles can exchange an α fermion in a leading order process, yielding the effective Hamiltonian It is this interaction that is ultimately responsible for the motion and transport of (N − 1)-ions, which -as we demonstrate later -dominates mass diffusion on the background of almost immobile N -ions.
Negative temperature state.We now prove that the initial state (2) must thermalize to a negative temperature state, T < 0. For that we only need to show that the energy density of the state |Ξ , 0 ≡ Ξ|H|Ξ /N sites is larger than that of the infinite temperature state, ∞ ≡ H T =∞ /N sites .For simplicity, we assume a particle-hole symmetrical band, but the proof carries over to any lattice and any dimension.For any product state of the the form (2), then one has Ξ|H J |Ξ = 0, and one obtains 0 = ν U N (N − 1)/2, with ν the filling factor, i.e. the ratio of sites occupied by N -ions.Here H J denotes the hopping part (first term) in Eq. ( 1).
In the infinite temperature limit, the hopping part of the Hamiltonian also averages to zero in case of particle-hole symmetrical Hamiltonians.The interaction part can be averaged by observing that each α fermion state at a given site is occupied with probability ν.Thus the interaction averages to U α<β n α n β T →∞ = U ν 2 N (N − 1)/2, yielding ∞ = ν 0 < 0 .The latter inequality immediately implies that if the state |Ξ thermalizes, it must thermalize to a negative temperature state [1,56,57] (see Fig. 1e).
To verify the thermalization of the system [58], we have performed exact diagonalization on small chains of linear sizes L = 9 and L = 6, and extracted the effective temperature using the condition that the energy density of the thermal state, H L T /L, be equal to the energy density of the initial state, 0 .We then used the extracted negative temperature to evaluate p n T , and compared the predicted values with the asymptotic values determined from our simulations.
Fig. 3a shows the extracted inverse negative temperatures for both system sizes.For small interactions, a perturbative calculation in U yields k B T eff ∝ −J 2 /U , while for large U we obtain k B T eff ∝ −U/ ln(U/J), in excellent agreement with the finite size numerical results.The probabilities p n T eff agree very well with the dynamically determined values, thereby evidencing the relaxation to a thermal, negative temperature state.The extracted negative temperatures correspond to very "hot" states of the system (|k B T | J.In such large magnitude negative temperatures, eigenstates with energies much below the upper edge of the spectrum are excited.Low energy properties [59,60] that characterize the spectrum near the edge are therefore insufficient to describe the rich dynamics that simulations uncover.
Entropy growth and correlations.The emergence of slow, composite particles is clearly visible in the time evolution of von Neumann entropy.The initial state at time t = 0 is a product state, and has vanishing entanglement entropy, wherever we cut the system in two.Entropy is generated by particles traveling from one part of the system to another [61].As shown in Fig. 3c, in our one-dimensional simulations, the von Neumann entropy increases linearly with time, and the entropy growth is barely influenced by the interactions as long as U J. For U J, however, the entropy growth is rapidly suppressed, indicating that particles carrying the entanglement move very slowly.The linear in time entropy growth is expected to be a general features even in non-integrable models, which thermalize in the long time limit [62] and do not exhibit many-body localization.
The dynamics of composite particles is more directly captured through time dependent charge oscillations in Fig. 4. For small interactions, U J, the charge at the origin x = 0, n x=0 (t) exhibits weakly damped coherent oscillations with a frequency ∼ J in Fig. 4a.This picture changes entirely once we enter the regime U J; there charge oscillations slow down, and universal oscillations with a frequency ω ∼ J 3 /U 2 appear.Interestingly, these oscillations are different from simple composite fermion oscillations.Rather, our direct simulations with the Hamiltonian (3) show that for large U/J, n 0 (t) approaches a universal curve described by Eq. ( 3) with infinitely strong interaction, |V | → ∞ (see inset of Fig. 4b).
Equal time density-density correlations [63], C nn (x, t) ≡ n y (t)n x+y (t) − n y (t) n x+y (t) , exhibit an even more interesting picture.In the quark-dominated U J regime, a ballistic front propagation and a light-cone structure is observed [64] in Fig. 4c.In contrast, for U J, the ballistic front is suppressed, and C nn (x, t) consists of two distinct features (see Fig. 4d) .A large peak associated with heavy and hardly moving baryons is observed at x ≈ 0. In addition, a diffusive correlation profile appears for x ≈ 0, very well described by the expression in one dimension.The diffusion constant can be directly extracted from the numerical data, and it scales as This scaling is clearly related to the meson diffusion.
Mesons have a small and non-negligible concentration, and move much faster than baryons.They collide with the background trions after a collision time τ ∼ U/J 2 , and propagate diffusively due to these collisions with a diffusion constant D ∝ J 2 /U/ν 2 (see inset of Fig. 4d).

Discussion
In the previous sections we have shown that hadronic states can be engineered very easily in optical lattices: one just has to prepare an initial state, where SU(N ) fermions are placed in N -fermion groups to a subset of lattice sites.A repulsive interaction U J stabilizes composite particles in this case, and bound states and resonances of hadronic nature emerge dynamically as a result of strong interactions.These heavy composite particles are simply stabilized by conservation laws, especially by that of energy and charge conservation, and form a 'hadron' gas of strongly interacting bosons and fermions [47,48,65].
In the particular case of N = 3, the heaviest hadrons are baryons (trions), while composites of two particles behave as somewhat lighter mesons (doublons), in analogy with QCD.The original, bare ψ α fermions play the role of quarks in this case.Similarly rich picture emerges for larger values of N , where n-ions with n = 1, . . ., N exist, and become long-lived particles for U J.
These particles behave somewhat similar to resonances in particle physics: they can transform into each other via many-particle collisions, and their densities equilibrate with time to form a negative temperature gas of hadrons.Notice, however, that since the emergent composite particles are very heavy and the decay channels require many-particle collisions, thermalization takes place at time scales much longer than the characteristic time t ∼ /J of the bare ψ α fermions' propagation.
Interestingly, for bipartite lattices, one can also show that for the states |Ξ , the time evolution of the density of the repulsive gas is identical to that of the attractive gas.To prove this, we notice that on a bipartite lattice, the gauge transformation K : ψ rα → (−1) ||r|| ψ rα , flipping the sign of the fermion operators on one sublattice of the bipartite lattice, changes the sign of J, but leaves the initial state as well as the density operator invariant, K : H J,U → H −J,U and |Ξ → |Ξ .Combing the time reversal symmetry T with K then yields, T K e −i t H J,U / |Ξ = e −i t H J,−U / |Ξ .Since both the time reversal and the gauge transformation K leave the density invariant, this implies that, starting from the state, |Ξ , the time evolution of density correlations is identical for U and for −U .Therefore, the quantum quench protocol suggested here can also be viewed as effectively realizing a gas of strongly attractive ψ fermions [56].
Initial states for N = 2 with two fermions at every second site have been realized a long time ago [47][48][49].The preparation of trion states could be performed by following similar protocols, but for N > 2, three particle losses may become an important factor [39,66].To perform a cold atom experiment, one should therefore use atoms with a relatively short scattering length a s compared to the wavelength of the optical lattice, λ.This condition is satisfied by 173 Yb, having a scattering length a s ≈ 10 nm [38], much shorter than the wavelength of the confining laser.In fact, a simple calculation yields that the ratio of the interaction U and the three-body loss rate of a trion, γ scales as U/(γh) ∼ (λ/a s ) 3 .For the laser used in Ref. [38] with λ = 759 nm, e.g., and for a barrier height corresponding to U/J = 1, we obtain the estimate U/h = J/h ≈ 300 Hz, while the three-body loss rate remains γ ≈ 0.16 Hz J/h.The time scale of three-body losses is thus more than 10 3 times larger than that of the dynamical scale of the strongly interacting gas in the regime where baryons and mesons form, and should be experimentally accessible [38,67].
Our results demonstrate that not only is it possible to create easily a plethora of interesting particles in the repulsive SU(N > 2) Hubbard model in all dimensions, but their detection and distinction are also straightforward using quench spectroscopy.This could be useful for many other incarnations of particle production in condensed matter and cold atomic ensembles.Non-Abelian MPS simulations.MPS simulations have been performed for the one dimensional Hubbard model with N = 3 colors.The real time dynamics, generated by the Hamiltonian (1) has been simulated using the infinite chain non-Abelian Time Evolving Block Decimation (NA-TEBD) algorithm [52], while exploiting the full SU(3) × U(1) symmetry of the model.We have kept M mult = 2500 multiplets in the NA-MPS, which corresponds to a usual bond dimension M ≈ 15000.The time step in the second order Suzuki-Trotter approximation was set to J∆t = 0.01.

Methods
We have also performed TEBD simulations using the effective trion Hamiltonian (3) in the |V | → ∞ limit, by using the residual U(1) symmetry of the effective model.Since the infinitely strong interaction forbids trions to hop next to each other, we could map the strongly interacting limit to free fermions by using a method of Cheong and Henley [68].
Exact Diagonalization at finite T. Finite temperature simulations for the SU(3) Hubbard chain have been carried out for short chains of length L = 6 and L = 9 with periodic boundary conditions.Corresponding to ν = 1/3 filling, we have set the total number of fermions to n tot = L, and also restricted the calculations to the color symmetrical subspace with n 1 = n 2 = n 3 .We have determined the canonical density matrix 1 Z e −βH by using a Taylor expansion.The inverse temperature β has then been extracted by enforcing Ξ|H|Ξ ≡ H T .
Derivation of effective N -ion Hamiltonian.The effective Hamiltonians describing the motion of composite particles can be constructed by performing perturbation theory in J.The effective hopping for an n-ion, e.g., can be estimated in leading order [53,69,70] as where r and r are two neighboring sites, H U denotes the on site energy, and T r,r = −J N α=1 ψ † rα ψ r α the hopping operator from r and r.In the intermediate state, where l particles are on r and (n − l) at site r the denominator can be evaluated as E The product of the denominators can therefore be evaluated as U ((n − 1)!) 2 , while the order in which particles are transferred between the two sites amounts in a combinatorial factor, n!, yielding the result in the main text.
Estimation of effective temperature.Small U limit.To estimate the effective temperature for U J, we first notice that β in this limit is a very small negative number.We can therefore perform a perturbative calculation in this limit by keeping only the leading order terms in β and U , and by taking the expectation value H J + H U T with the unperturbed density operator, and equating that with the energy of the initial state.This procedure yields the equation where ρ( ) denotes the single particle density of states of the lattice.The prefactor (1 − ν)ν in the first integral is associated with the (diverging) chemical potential, e −βµ ≈ (1 − ν)/ν.On a D dimensional cubic lattice Eq. ( 7) yields In D = 1 dimension and for N = 3 and ν = 1/3, we obtain a prefactor 1/2, in perfect agreement with Fig. 3a.Large U limit.In the U J limit, we can consider each site as a grand canonical subsystem, weakly coupled to the rest of the system, and having a density matrix ρ ∝ e −β(U n(n−1)/2−µn) .Since N -ions are abundant, the chemical potential of the subsystem must be approximately µ ≈ U (N − 1)/2.The probability of having a single fermion is then given by p n=1 ≈ (1 − ν)N e βU (N −1)/2 , where the prefactor approximates the inverse partition function.
The value of p n=1 can, however, be also estimated quite simply by noticing that the state, where we have N particles at one site, say at the origin, |N , is not a true N -ion state at finite J/U .Quantum fluctuations allow each fermion forming the N -ion to move to neighboring sites with probability amplitudes ≈ J/((N − 1)U ).This implies that a state |N decomposes to a fermion and an (N − 1) -ion with a probability P N → (N −1)+1 ≈ zN J 2 /[(N − 1) 2 U 2 ], yielding p 1 ≈ p N −1 ≈ z ν N J 2 /[(N − 1) 2 U 2 ].Comparison with the thermal average then yields with z = 2D on a D-dimensional cubic lattice.
Data availability.The data that support the findings of this study are available from the corresponding author upon reasonable request.

FIG. 1 .
FIG. 1. a The initial state |Ξ for N = 3. Sites are either empty or occupied by N fermions of different colors.b Schematic motion for times t > 0: n-particle composites move between neighboring sites with hopping Jn, and have an on-site interaction energy n(n − 1)/2U .c Decay of N -ions (N = 3 in the figure) to smaller composites is forbidden by energy conservation.d Lowest order decay channel for N = 3.A a trion (baryon) and three fermions (quarks) transform into three doublons (mesons).eSketch of the entropy density as a function of energy density.The white dashed line indicates the J = 0 limit.The entropy density has a maximum at energy density ∞, where the temperature diverges and changes sign.For states with < ∞ the temperature is positive, while for states with > ∞ it becomes negative.

FIG. 2 .
FIG.2. a Evolution of the probabilities pn for U/J = 0.5.The number of empty and triply-occupied sites decreases rapidly, and a large fraction of sites becomes occupied by one or two fermions.The Fourier spectrum of p0(t) (right panel) consists of a spectrum of a broad frequency range up to ω/J 10. b Evolution of the probabilities pn for large interaction strength, U = 6J.The number of empty and triply-occupied sites remains almost constant apart from small oscillations.The Fourier spectrum of p0(t), is mostly restricted to a window 9.5 ω 14.5.The shading highlights the predicted spectral window, ω ≈ 2U ± (2J + 8J 2 /U ).Frequency peaks at the edges are due to van Hove singularities in the single fermion (quark) band.c Sketch of quasiparticle energies in the large U limit.Oscillation frequencies in panel b can be interpreted as a result of coherent oscillations between a baryon (trion) and a state decomposed into a meson (doublon) and a quark (fermion).

FIG. 3 .
FIG.3.a Extracted dimensionless inverse negative temperature −J/T /kB of the stationary state as a function of the interaction strength, U/J, as obtained for small SU(3) chains of lengths L = 6 and L = 9 and filling ν = 1/3.The fitted value of −J/T /kB has a clear maximum around U/J ≈ 3. The extracted temperature is negative for any repulsive interaction U .For small U we obtain kBT ≈ −2J 2 /U , while for large U , a fall-off kBT ∝ −U/ ln(U/J) is analytically predicted.b Stationary (longtime) charge distribution pn as a function of the interaction strength.Filled symbols represent data obtained from NA-TEBD simulations by calculating the time average of pn(t) for long times.Thermal predictions for a short chain of length L = 9 are shown by empty symbols.c Evolution of the von Neumann entropy of a half chain for different interaction strengths.We observe linear growth of entropy for all values of U but the growth rate is strongly reduced for large U values.

FIG. 4 .
FIG.4. a Average fermion number at a triply occupied site as a function of time for weak interactions, U/J ≤ 1.In the non-interacting limit, we observe algebraically decaying oscillations, which become weakly but exponentially damped for weak interactions.b Strong interactions, U/J 1, slow down the dynamics.Inset: The rescaled curves, nx=0(t J 3 /U 2 / ) , collapse to a single universal curve for U → ∞, well captured by the |V | → ∞ limit of the effective Hamiltonian (3) (continuous line).c Amplitude of equal time density-density correlations |Cnn(x, t)| for U = 0.5J.Correlations develop within a light-cone of a slope ∝ 4J/ , i.e., twice the maximal velocity of free fermions.d Amplitude of equal time density-density correlations for strong interaction, U = 6J.Apart from the central trion peak at x ≈ 0, correlations spread diffusively according to(5).Inset: The fitted diffusion constant D scales as D ∼ J 2 /U , evidencing meson (doublon) diffusion.