Emergent Time Scale in Entangled Quantum Dynamics of Ultracold Molecules in Optical Lattices

We derive a novel lattice Hamiltonian, the \emph{Molecular Hubbard Hamiltonian} (MHH), which describes the essential many body physics of closed-shell ultracold heteronuclear molecules in their absolute ground state in a quasi-one-dimensional optical lattice. The MHH is explicitly time-dependent, making a dynamic generalization of the concept of quantum phase transitions necessary. Using the Time-Evolving Block Decimation (TEBD) algorithm to study entangled dynamics, we demonstrate that, in the case of hard core bosonic molecules at half filling, the MHH exhibits an emergent time scale over which spatial entanglement grows, crystalline order appears, and oscillations between rotational states self-damp into an asymptotic superposition. We show that this time scale is a non-monotonic function of the physical parameters describing the lattice. We also point out that experimental mapping of the static phase boundaries of the MHH can be used to measure the molecular polarizability tensor.


Introduction
In recent years, ultracold atomic gases have provided near perfect realizations of condensed matter Hamiltonians, acting as quantum simulators [1,2] that allow the study of complex condensed matter phenomena in a clean and highly controllable environment. Ultracold polar molecular gases, which have recently been brought to the edge of quantum degeneracy in their absolute ground state [3,4], offer additional features over atomic gases, such as a large internal Hilbert space and a greater susceptibility to external fields via a permanent electric dipole. There have been a number of proposals on how to use ultracold molecular gases for mimicking well-known Hamiltonians such as spin-1 lattice models [5]. Ultracold molecules have also been suggested as a model system for the study of strongly correlated 2D quantum phases [6] or for quantum information processing schemes [7,8,9]. However, these proposals frequently involve complex and yet-to-be implemented experimental techniques. In this article, we instead focus on the completely new quantum many body physics which results naturally from the simplest quantum lattice experiments that can be performed in the immediate future with established techniques in ultracold molecular quantum gases.
Towards this end we derive a novel lattice Hamiltonian, which we refer to as the Molecular Hubbard Hamiltonian (MHH). The MHH describes the physics of an ultracold polar molecular gas in a 1D optical lattice that is oriented using a DC electric field, giving rise to a resonant dipole-dipole interaction, and is driven between rotational levels using a microwave AC field. In particular, new aspects of our derivation include explicit dependence of hopping energy on the molecular polarizability tensor. This in turn allows a determination of the tensor elements, an outstanding experimental issue, from the borders of the static phase diagram of the MHH, which are identical to those of the extended Bose-Hubbard Hamiltonian [10] when a single molecular rotational level is occupied.
Beyond the statics, the MHH naturally has a dynamical component due to the AC driving fields, as well as an internal structure in terms of rotational modes which is inherently different from spinor atomic systems [11,12]. We study this dynamical aspect with Time-Evolving Block Decimation (TEBD) [13,14], a newly developed entangled quantum dynamics algorithm which takes spatial entanglement (specifically, Schmidt number [15]) as a cut-off. We find an emergent time scale in the case of half-filling for hard core bosonic molecules. We emphasize that a quantum lattice model requires low filling (average number of particles per site), in contrast to a mean field lattice model, for which the filling would typically be quite high. Thus, although experiments can most easily access the mean field regime of hundreds of molecules per site with a single pair of counter-propagating laser beams, we look slightly ahead to the quantum regime, which will require two pairs of such beams in order to create an array of quasi-1D "tubes." A third pair is then used to create the lattice in each tube. This technique is already well established for ultracold atoms [16].
Dynamical aspects of quantum phase transitions are just beginning to be considered [17,18], and have so far been a limited area of study restricted to mean field considerations, due to lack of numerical tools. With the recent advent of entangled quantum dynamics algorithms, namely TEBD, dynamical properties of many-body systems are becoming amenable to numerical study. For example, TEBD has been used to address key questions such as the dynamics of a quantum quench [19,20] or the speed at which correlations propagate in a lattice [21]; these are not issues which can be studied with other dynamical methods such as dynamical mean field theory (DMFT) [22]. We give a brief review of TEBD in Sec. 3. The reader interested in computational details can find them in Ref. [23].
The first main contribution of this paper is to present a careful derivation of the Molecular Hubbard Hamiltonian. This is done in Sec. 2, with some previously known aspects of molecular physics relegated to Appendix A. The second main contribution is to present an emergent time scale for half filling; although we treat the case of hard core bosons, the MHH can also be applied to fermionic molecules. To this end, in Sec. 3 we first give a brief explanation of TEBD and the quantum measures we use. Then, in Sec. 4 we present and analyze our simulations, with an accompanying convergence study in Appendix B. Finally, in Sec. 5 we summarize.
In the following subsections and Appendix A we justify Eq. (1) with a careful derivation and present the energy scales of each term.

Derivation of the Molecular Hubbard Hamiltonian
The full molecular Hamiltonian in second quantization iŝ H = d 3 rψ † (r) Ĥ kin +Ĥ rot +Ĥ DC +Ĥ AC (t) +Ĥ opt (r) ψ (r) The terms on the first line correspond to single-molecule effects: kinetic energy, rotation, the DC electric field which orients the dipole, the AC microwave field which drives transitions between rotational levels, and the far off-resonant optical lattice potential, respectively. The second line is the two-molecule resonant dipolar energy. The field operatorsψ can be either bosonic or fermionic. We focus on the bosonic case for brevity. There are five key assumptions underlying our derivation, as follows. We consider all five assumptions to be reasonable for present and near-future experiments.
(i) We consider ultracold closed-shell polar heteronuclear diatomic molecules, characterized by permanent dipole moment d and rotational constant B. The most experimentally relevant bosonic species in this category are SrO, RbCs, and LiCs [6]. The individual molecules are assumed to be in their electronic and vibrational ground states, and it is assumed that none of these degrees of freedom can be excited at the large intermolecular separations and low temperatures/relative energies that we consider.
(ii) The molecule is assumed to have a 1 Σ ground state. The characteristic trapping potential length is chosen large enough compared to the internuclear axis to assume spherical symmetry, i.e. a locally constant potential.
(iv) We consider only the lowest three rotational levels. All AC fields will be sufficiently weak to allow this assumption.
(v) We work in the "hard-core" limit where at most one molecule is allowed per site. This is enforced by strong dipole-dipole interactions on-site. We consider the lattice spacing large enough to include only nearest-neighbor dipole-dipole interactions. Other short-range interactions such as exchange or chemical reactions or long range interactions such as dispersion and quadrupole-quadrupole interactions are not considered.
We proceed to follow the usual procedure [25] of expanding the field operators of our secondquantized Hamiltonian in a Wannier basis of single-molecule states centered at a particular discrete position r i :ψ where i is a site index and the sum is over all lattice sites. For our Wannier Basis we choose the singlemolecule basis that diagonalizes the rotational and DC electric field Hamiltonians, spanned by kets |E; JM . In this basis, which we refer to as the "dressed basis" (the DC field "dresses" the rotational basis) we have the field operator expansion We note that such a basis, while highly efficient for the hard core limit we consider, becomes progressively worse for higher filling factors, till in the mean field limit the single-molecule basis, whether dressed or not, is so poor that many bands must be considered. Here we do not include a band index for simplicity, although the generalization of Eq. (1) to include multiple bands is straightforward. This choice of Wannier basis associates the terms in Eq. (1) to the terms in Eq. (2) as follows: (8) where the operators H kin , H opt , etc., are taken to be in position space representation. For the derivation of the single-molecule terms (rotational, DC electric field, and AC electric field) and discussion of the properties of our Wannier basis we refer the reader to Appendix A. In the following sections we present the derivation of the tunneling (hopping) and dipole-dipole terms, which have new aspects not heretofore appearing in the literature [26].

Tunneling
The tunneling term represents the sum of the molecular kinetic energy with the potential energy of the lattice. After expanding in the Wannier basis of Eq. (4), we find the effective tunneling Hamiltonian where t J,J ′ ,M was defined in Eq. (5). To understand why this operator mixes states of different J, we note that the kinetic energy and (far off-resonant) optical lattice potential do not mix rotational eigenstates. Because our Wannier basis states are dressed and therefore superpositions of rotational eigenstates with different J, the tunneling operator in the dressed basis will mix J. Although the dressed basis makes the tunneling more complex to analyze, it simplifies other terms in the MHH, such as the DC term, and is in any case a more standard basis for analysis of the diatomic molecules we study here. Comparable basis changes are sometimes made in other quantum many body systems, where, for instance, particles and holes are mixed, or particles are paired. Note that, because we assume z-polarized fields, M is still a good quantum number. To discuss the actual form of the tunneling energies {t J,J ′ ,M } we must first examine the interaction of a diatomic molecule with the optical lattice.

Interaction with an Optical Lattice
The charge redistribution that occurs when a molecule is subjected to a static, spatially uniform electric field E is reflected in its dipole moment d via the polarizability series where the first, second, and third order coefficients α jk , β jkl , and Γ jklm are elements of the polarizability, hyperpolarizability, and second hyperpolarizability tensors, respectively. The polarizability tensor is a symmetric rank-two tensor with no more than six independent elements (less if molecular symmetry is greater), and characterizes the lowest order dipole moment induced by an applied electric field. From this tensor we can form the scalar invariants referred to as the polarizability and the polarizability anisotropy, respectively. Note that we use the tilde to clarify thatα with elements α jk is a tensor, not a scalar -we reserve the accent circumflex (the "hat" symbol) for quantum operators. In linear molecules, such as diatomic molecules, the presence of only two distinct moments of inertia allows for the classification ofα according to its components along and perpendicular to the internuclear axis, denoted α and α ⊥ , respectively. In the presence of AC electric fields with frequency ω we speak of the dynamic polarizability tensorα (ω), with the series of Eq. (10) being the zero frequency limit. The tensorα (ω) is, in general, complex, with the real part inducing a dipole moment and the imaginary part accounting for power absorption by the dipole and out-of-phase dipole oscillation. In the case of Σ diatomic molecules in their electronic and vibrational ground states [27] α where the e ′ q are molecule-fixed spherical basis vectors. The parallel and perpendicular dynamic polarizabilities are respectively. In these expressions d νΛ(v)−XΣ(0) is the transition dipole moment from the ground state to the νΛ (v) state (following the usual diatomic molecular notation, Λ ∈ {Σ, Π} ≡ {0, 1} is the quantum number associated with the projection of the total electronic orbital angular momentum along the internuclear axis, i.e., in the molecule-fixed basis) and the sum over ∓ accounts for the near-resonant and typically far off-resonant terms. Transformingα from the molecule-fixed basis to the space-fixed basis using the transformation discussed in Appendix A, we find where C (j) m is an unnormalized spherical harmonic, (. . .) denotes the Wigner 3-j coefficient [28], and the e p are space-fixed spherical basis vectors.
The interaction of the lattice with the molecule is represented by the Hamiltonian If the electric field has polarization p in the space-fixed spherical basis then we find

Entangled Quantum Dynamics
For light linearly polarized in thex-direction we obtain whereas for light linearly polarized in theŷ-direction we find Since C (0) 0 = 1, these terms give a state-independent energy shift. The C (2) q terms produce a tensor shift. Because the depth (in energy) of a typical optical lattice is much smaller than the energy of transitions between rotational levels (of order B, as defined in Appendix A), we can ignore far off-resonant Raman coupling between different J manifolds and use only the diagonal matrix elements. The C (2) 2 term and the C (2) −2 will both mix M in the J ≥ 2 manifolds, but do not affect the lowest two rotational levels, again, because we neglect Raman couplings. Thus x, y, and z polarizations all have the same Hamiltonian in this approximation. We can calculate the matrix elements of C (2) 0 in the field free basis using the Wigner-Eckart theorem to find In our effective Hamiltonian we choose right circular polarization for the z lattice, x polarization for the x lattice, and y polarization for the y lattice, where each "lattice" refers to a pair of counter-propagating laser beams used to create a standing wave. We consider the fields making up the optical lattice to have sinusoidal spatial profiles, resulting in sine-squared intensity profiles. In addition, we assume that the y and z lattices are tight, meaning that the molecules are strongly confined at the potential minimum (for a red-detuned trap). This tight confinement allows us to approximate them via a Taylor series, e.g., sin 2 (k z z) ≃ k 2 z z 2 in the vicinity of the molecule. Using the above results, the matrix elements of the Hamiltonian for the optical lattice can be written or, more compactly, as by defining We now define, as is customary, the "lattice heights" in the x, y, and z directions, respectively, as The tight confinement in the transverse (y and z) directions strongly suppresses tunneling in these directions, making the overall lattice effectively 1D along x. From Eqs. (24)- (25), it is apparent that different rotational levels experience different trapping frequencies and different tunneling energies. To make this clearer, we parse our full field-free tunneling matrix element as we proceed to compute each piece separately.
In the evaluation of the first integral, Eq. (29) we assume that the Bloch function of a molecule in the sinusoidal optical lattice is a Mathieu function along x. This may seem to contradict our assumption of spherical symmetry in the above derivation. However, the assumption of spherical symmetry (i.e. a locally constant potential) need only hold on the order of an internuclear axis (∼ 5Å) near the molecule. In contrast, on the order of the characteristic lattice length h/µω opt the rigid-rotor molecule is indistinguishable from a point particle (such as an alkali atom), and so spherical symmetry is not required. With this understanding, we recognize t (0) JM as the expression for the hopping energy for point particles in optical lattices [29] with the additional feature that the lattice height along the quasi-1D direction V 0 = V (JM ) x is dependent on J through the polarizability tensor. Thus, altering the expression from the theory of point particles in optical lattices, we obtain the result where A = 1.397, B = 1.051, C = 2.121, and is the recoil energy. For the second integral, Eq. (30), we approximate the Wannier functions with the ground state of a simple harmonic oscillator where the harmonic oscillator lengths are given by  Then where λ is the wavelength of the optical lattice. Because we consider tight traps such that the lattice height in the y and z directions is much greater than the lattice height in the JM , and so we neglect it. Thus, This is equivalent to the array of tubes we discussed in Sec. 1, where each tube is isolated from its neighbors.
Using tabulated values of the polarizabilities for LiCs [30] as given in Table 1, we find that, for a reasonable lattice height V (00) x /E R ≃ 10, the tunneling term for the |11 state is only about 20% of that in the |00 state, as shown in Fig. 1. For LiCs in a red-detuned optical lattice of wavelength λ = 985nm, We reiterate that the above matrix elements and tunneling energies {t JM } have been computed in the field-free basis for simplicity. To transform to the dressed basis, we use the unitary matrix with dressed eigenvectors as columns, recovering Eq. (9), where the tunneling matrix element is no longer diagonal in J.

Dipole-Dipole Interactions
The induced dipoles from the DC field give rise to a resonant dipole-dipole interaction. The Hamiltonian for this interaction in the two-site dressed basis spanned by |E; where we have defined and for notational simplicity we have suppressed the E subscripts. Note that because of our choice of polarizations of the optical lattice and AC and DC electric fields, The resonant dipole-dipole interaction between two permanent dipoles d 1 and d 2 whose respective centers of mass are separated by a vector R in the space-fixed frame iŝ where e R is a unit vector in the direction of R. Using standard angular momentum recoupling we recast this in spherical tensor notation as where (T ) (k) q denotes the component of the rank-k spherical tensor T that has projection q along R, C (j) m (R) is an unnormalized spherical harmonic in the polar coordinates defined with respect to R, and we have defined the tensor product of the vector operatorsd 1 andd 2 as In the last line, j 1 , m 1 , j 2 , m 2 |J, M is a Clebsch-Gordan coefficient. We now take matrix elements of Eq. (41) in the two dressed-molecule basis |E; J 1 M 1 , J 2 M 2 , where molecule 1 is on site i and molecule 2 is on site i + 1, yielding Because our DC field is polarized along z, only (d 1 ) 0 matrix elements are nonzero, enforcing µ = 0, m = 0. With this in mind, the interaction takes the particularly simple form The intermolecular axis plays a crucial role in the sign of the interaction. Two molecules oriented along the intermolecular axis attract if their dipoles are parallel and repel if their dipoles are antiparallel. Two molecules oriented perpendicular to the intermolecular axis, on the other hand, repel if their dipoles are parallel and attract if their dipoles are antiparallel. The DC field that orients the molecules in our setup is polarized along z, perpendicular to the intermolecular quasi-1D axis x. This gives rise to repulsive interactions for positive dipole matrix elements. With this geometry the dipole potential becomes yielding where λ is the wavelength of the optical lattice.

Energy Scales
We proceed to clarify the energy scales associated with each term in Eq. (1). Between previous discussion in Sec. 2 and that of Appendix A, all terms in Eq. (1) are now clearly defined. The energy scales of the dressed basis are B, the rotational constant, which is roughly 60h GHz, and dE DC , which is of order 1 − 10B. The DC term has no length scale associated with it because the field is uniform, and the length scale of the rotational term is the internuclear separation, on the order of angstroms. The relative contribution of the DC electric field and rotational terms in Eq. (1) are expressed through the dimensionless parameter the ratio of the DC field energy to the rotational level splitting. The energy scales of the AC term arehω, where ω is the angular frequency of the driving field, and dE AC . The scalehω is of order 2B for small β DC ≪ 1, and of order B √ β DC for large β DC ≫ 1. The AC field energy dE AC is of order 0.5hω. The single-molecule time scale associated with dE AC is the Rabi period, the time it takes for the population of a two-level system to cycle once, as seen in Figure 3(a). In real time, this is on the order of 10ps for the parameters in the preceding paragraph. The time scale associated with ω is the time scale on which the small oscillations in Figure 3(a) occur, of order 0.5ps. The length scale of the AC field is on the order of centimeters, and so we can neglect this in light of the micron length scale of the trap.
The tunneling term has several scales. The optical lattice near the point of confinement has a length scale given by the harmonic oscillator length l (00) ho,x ∼100nm and an energy scale of E R ≈1.4h kHz. The energy scales of the tunneling operator proper are given by the {t JJ ′ M } which are of order 10 −1 -10 −2 E R ∼100h Hz for the given recoil energy.
There are also many scales for the dipole term. For the B and d specified in the first paragraph of this section and β DC = 1.9, the characteristic length scale where the dipole-dipole energy becomes comparable to the rotational energy is approximately 348 Bohr radii (18.4nm). Outside this region the Born-Oppenheimer adiabatic approximation is easily fulfilled [6]. Since the length scale of our optical lattice is of order µm, we Term Length scale Energy scale are justified in working within the Born-Oppenheimer framework. For the same parameters, the length scale where the off-resonant van der Waals potential C 6 /r 6 ≈ −d 4 /(6Br 6 ) becomes comparable to the dipole-dipole interaction is This length is very small, on the order of tens to hundreds of Bohr radii. Outside of this region the resonant dipole potential dominates and the intermolecular force is repulsive. This repulsion enforces the hard-core limit. The energy scale of the dipole-dipole force is E; 00|d|E; 00 2 /λ 3 ∼1.2h kHz, with higher J being an order of magnitude or so lower for small β DC , and of the same order for large β DC (see Fig. 1(a)).
To summarize, the scales of the problem are shown in Table 2.

Novel Features of the Molecular Hubbard Hamiltonian
The MHH, Eq. (1), has a number of novel features which distinguish it from the Hamiltonians typically considered in the quantum lattice and condensed matter literature [32,25]. First, the tunneling energies {t J,J ′ M } not only depend on the rotational level J, M but even change rotational states from J to J ′ . This is due both to the polarizability tensor's dependence on rotational level, and to the dressed basis. This differs from other Hubbard models which consider spin degrees of freedom, as tunneling does not occur between spin states -hopping does not cause spin transitions. If we consider populating a single mode (e.g. J = 0, M = 0) in the Ω → 0 limit, then Eq. (1) becomes the extended Bose-Hubbard Hamiltonian, and the phase diagram is known [33,10]. This gives ideas of how to characterize the static phases of the MHH. However, because the tunneling energy depends on J, the borders of the phase diagram will depend on the rotational state of the system. We will discuss this property and provide an application in Sec. 4.
Second, the Hamiltonian is fundamentally time-dependent because it is a driven system. This allows for the study of dynamic quantum phases, requiring the concept of a quantum phase diagram to be generalized to an inherently time-dependent picture. In a case study for hard core bosonic molecules at half filling presented in Sec. 4, we show that the MHH has an emergent time scale.

Time-Evolving Block Decimation
The Time-evolving Block Decimation algorithm (TEBD) is a new method [34,35] designed to study the dynamics of entangled quantum systems. The essential idea of TEBD is to provide a moving "spotlight" in Hilbert space which tracks a dynamical system. The portion of the Hilbert space so illuminated is an exponentially small fraction of the full Hilbert space; this is justified by the fact that real, physical quantum many body systems, especially in real materials, typically explore only a small, lowly-entangled part of the total Hilbert space.
In fact, TEBD moves the full quantum many-body problem from the NP-complete complexity class to the P class through an exponential reduction in the number of parameters needed to represent the many body state. We can understand the possibility of this reduction through an analogy to image compression. Present digital cameras are capable of producing a roughly 3000 × 3000 array of pixels. Downloading the images from such a camera, one notices that there are far less than 10 Megapixels worth of data per image. Image compression algorithms such as JPEG produce images of remarkable quality with only a small fraction of the raw data. The reason that these algorithms are so effective is that a physical image, as opposed to a random 2D pixel array, is not the "most common" or most probable image; it contains a great deal of structure and regularity. In the same way, physical states in Hilbert space tend to be lowly entangled (by some entanglement measure), even though a general state in Hilbert space has a much larger probability of being highly entangled. There is no general proof of this fact, just as there is no guarantee that an image will come out perfectly crisp after JPEG compression; it is simply a trend observed in many-body quantum systems.
To be slightly more specific, TEBD performs a partial trace over a particular bipartite splitting of the lattice, and then keeps the χ largest eigenvalues of the resulting reduced density matrix. The cut-off parameter χ is based on the Schmidt measure [15], and so it also serves as a measure of the degree of spatial entanglement. This idea is not unique to TEBD. In fact, the density matrix renormalization group (DMRG) method first proposed by White [36] did something analogous years before. TEBD's innovation is that at each time step it re-optimizes the truncated basis (thus the "moving spotlight"). The Schmidt number is just the number of non-zero eigenvalues in the reduced density matrix, and so is an entanglement measure natural to quantum many body systems. The parameter χ is the number of non-zero eigenvalues in the reduced density matrix that TEBD retains. It is the principal convergence parameter of the algorithm, both in entanglement and in time. Although the time-propagation method we use is Trotter-Suzuki [37], it turns out that, due to a normalization drift, χ controls convergence at long times.
With χ interpreted as an entanglement measure, we can say that TEBD treats the system not as a wavefunction in a d L -dimensional Hilbert space (L is the number of lattice sites), but as a collection of wavefunctions in d 2 -dimensional two-site spaces that are weakly entangled with the environment created by the rest of the system. To facilitate this viewpoint, we replace the d L coefficients of the full many-body wavefunction with L sets of (dχ 2 + χ) coefficients corresponding to the wavefunctions of each bipartite splitting. The most computationally expensive portion of the TEBD algorithm is typically the diagonalization of these local coefficient matrices at a cost of O (d 3 χ 3 ). Looping over all L − 1 bipartite splittings and evolving the system for a total time t f in time steps of length δt, one obtains an asymptotic scaling of O L t f δt d 3 χ 3 . This scaling can be greatly improved by the presence of conserved quantities. When a conserved quantity exists in the system we are able to diagonalize reduced density matrices corresponding to distinct values of this conserved quantity independently, which can result in significantly smaller reduced density matrices to diagonalize. Implementing this idea, scalings of O (χ 2 ) have been reported for fixed d [38]. In addition, conserved quantities in the presence of selection rules can reduce the local dimension. For example, in the case of the MHH, z-polarized electric fields disallow transitions from a particular M to any other. If we begin with all molecules in a particular M state, this allows us to restrict our attention only to states with this M. In our numerics we conserve both the projection M, and the total number N. Furthermore, to match our hard core requirement, we allow only zero or one molecules per site, so that the local dimension is d ≤ R + 1, R being the magnitude of the greatest angular momentum that we consider (note that the local dimension d, mentioned only here in Sec. 3.1, bears no relation to the permanent electric dipole moment d used throughout the rest of our treatment).

Quantum Measures
We use a suite of quantum measures to characterize the reduced MHH, Eq. (55) below. The fewbody measures we use are n J i , the number in the J th rotational state on the i th site, E ≡ Ĥ , the expectation of the energy, and 1 L n J , the average number in the J th rotational state per site (L is the number of lattice sites). The latter is a J-dependent filling factor. The many body measures we use include the density-density correlation between rotational modes J 1 and J 2 evaluated at the middle site where ⌊q⌋ is the floor function, defined as the greatest integer less than or equal to q. As an entanglement measure we use the Meyer Q-measure [40,41,42] whereρ (k) is the single-site density matrix obtained by tracing over all but the k th lattice site, and the factor outside of the bracket is a normalization factor (d is the on-site dimension). This gives an average measure of the entanglement of a single site with the rest of the system. The Q-measure can also be interpreted as the average local impurity (recall that the Tr(ρ 2 ) = 1 if and only ifρ is a pure state).
To determine what measures we can use to ascertain the static phases of our model we reason by analogy with the extended Bose-Hubbard Hamiltonian where we know that the possible static phases are charge density wave, superfluid, supersolid, and Bose metal [10]. The charge density wave is an insulating phase appearing at half integer fillings which has a wavelength of two sites. Like the Mott insulating phase, it has an excitation gap and is incompressible. While the extended Bose-Hubbard Hamiltonian has only one charge density wave phase due to the presence of only one species, the MHH has the possibility of admitting several charge density wave phases due to the presence of multiple rotational states. As such, we define the structure factor where N is the total number of molecules. We recognize this object as the spatial Fourier transform of the equal-time density-density correlation function between rotational states J 1 and J 2 , evaluated at the edge of the Brillouin zone. This measure is of experimental interest because it is proportional to the intensity in many scattering experiments, e.g. neutron scattering [43]. Crystalline order between rotational states J 1 and J 2 is characterized by a nonzero structure factor S (J 1 J 2 ) π . The charge density wave is the phase with crystalline order but no off-diagonal long-range order as quantified by the superfluid stiffness of rotational state J (note that ρ s bears no relation to the density matrixρ). If both the structure factor and the superfluid stiffness are nonzero, the phase is called supersolid. If both the structure factor and the superfluid stiffness are zero, the phase is called Bose Metal. Finally, if the structure factor is zero and the superfluid stiffness is nonzero, the phase is superfluid. In one dimension the entire superfluid phase is critical, and so there is no order parameter [10].

Case Study: Hard Core Bosonic Molecules at Half Filling
In the following, we consider a particular case of Eq. (1) for dynamical study. We choose the hard core case, which can occur naturally due to strong on-site dipole-dipole interactions, and half filling, which is an interesting point in a number of models, including the repulsive Fermi-Hubbard Hamiltonian and the extended Bose-Hubbard Hamiltonian discussed in Sec. 3.2. For example, in the latter case, the charge-density-wave phase requires a minimum of half-filling [10].
If we assume that our system begins in its ground state (J = 0, M = 0) we need only include states which have a dipole coupling to this state. For z-polarized DC and AC fields, this means we only consider M = 0 states, yielding the reduced Hamiltonian This is the specific case of the MHH that we study using TEBD. A matter of practical concern, as apparent in Table 2, is the large disparity between the timescales of the first three (Rotational, DC, and AC) and the last three (kinetic, tunneling, and Dipole-Dipole) terms. The accumulation of error resulting from truncating the Hilbert space at each TEBD timestep causes the algorithm to eventually fail after a certain "runaway time," making studies over long times intractable [44]. This invites a multiscale approach in the future [45,46]. In our current numerics we artificially increase the recoil energy and dipole-dipole potential to be of the order of the rotational constant in order to study Eq. (55) using TEBD. In particular, we take t J = 10B η 1 + 2 ∆ᾱ α J (J + 1) (2J + 1) (2J + 3) where the dimensionless variable η becomes an ersatz "lattice height." To see the scaling more explicitly, we compare the above with the actual expressions for the MHH parameters If we now scale E R to be 10B/1.397 and set d such that 2mE R d 4/3 / h 2 π 2 3 2 = 10B for this E R , we recover Eqs. (56) and (57) provided we make the definition Since this dimensionless parameter plays the same role as the quasi-1D lattice height scaled to the recoil energy did in the actual MHH, we refer to it as the lattice height. For the polarizability tensor, we choose ∆α/ᾱ = 165.8/237, corresponding to LiCs [30]. This rescaling does not change the qualitative static and dynamical features of Eq. (55); it only makes Eq. (55) treatable directly by TEBD, without multiscale methods. First, we point out that if we consider populating a single rotational state (e.g. J = 0, M = 0) in the Ω → 0 limit, then Eq. (55) becomes the extended Bose-Hubbard Hamiltonian, and the phase diagram is known [33,10]. Because the tunneling energy is different for different rotational states (see Eq. (31)) and this difference depends only on the properties of the polarizability tensor, we can relate the borders of the phase diagram for different rotational states to properties of the polarizability tensor. The MHH thus gives a means to measure the polarizability tensor, a standing issue in experiments [47]. Our calculations in Sec. 2 can be used to compare directly to the phase diagram from the literature [33,10]. In fact, this aspect of our work, unlike the simulations below, is not restricted to 1D.
However, our main focus at present is on the dynamics of the MHH. In the following numerical study, we explore dynamics as a function of the physical characteristics of the lattice, namely, number of sites L and effective lattice height η. Specifically, we study L = 9, 10, and 21 lattice sites with N=4, 5, and 10 molecules, respectively, and η ranging from 1 to 10. We fix the dipole-dipole term as in Eq. (56), and fix the DC field parameter to be β DC = 1.9. While β DC = 1.9 may not correspond to a physically realizable situation, its exploration provides insight into the MHH.
The Rabi oscillations between the J = 0 and the J = 1 states damp out exponentially in the rotational time t r ≡ Bt/h as with some characteristic time scale τ , as seen in Fig. 2. We note that an exponential fit has a lower reduced chi-squared than a power-law, or algebraic fit. We also tried fit functions where the oscillations do not decay to zero, but rather persist with some asymptotic nonzero amplitude. We find that the fit functions Eqs. (63) and (64) above fit the data better as quantified by the convergence properties of the algorithms used, as discussed in Appendix B. (c) Site-averaged population vs. rotational time for 10 sites. Note that there is no significant difference between an odd and even number of sites. Note that there is no significant difference between this and the smaller system sizes.  (a) Structure factors vs. rotational time for 9 sites. Note the similar asymptotic behavior to the populations in Fig. 2(a).
(b) Squared modulus of Fourier transform of S (00) π vs. rotationally scaled frequency for L = 9 sites. Note the similarity with Fig. 2(b) above. vs. rotationally scaled frequency for L = 9 sites. Note the absence of the Rabi frequency.
(e) Structure factors vs. rotational time for 21 sites. Note the lack of significant difference with the smaller odd system size.
(f) Comparison of the S (01) π correlation structure factor for odd and even numbers of sites. Note that the even site (exactly half filling) structure factor grows faster and larger than the odd site (slightly less than half filling) structure factor. Figure 3. Dependence of structure factors within and between rotational states J on the number of lattice sites. We do not consider the off-resonant J = 2 and higher rotational states because they have a very small occupation; J = 2 is shown explicitly in Fig. 2. (a) Site-averaged population vs. rotational time for 21 sites with η = 5. Note that the J = 0 and J = 1 states now appear to converge to different fillings.
(b) Squared modulus of Fourier transform of n 00 vs. rotationally scaled frequency for L = 21 sites and η = 5. Note the presence of several new frequencies not observed in the η = 1 case ( Fig. 2(f)). In particular, Ω 00 , 2Ω 00 , and 3Ω 00 , are denoted by arrows.
(c) Site-averaged population vs. rotational time for 21 sites with η = 10. Note the similarity to the η = 1 case (Fig. 2(e)) and the difference from the η = 5 case (Fig. 4(a))-the asymptotic behavior is not a monotonic function of the lattice height.
(d) Squared modulus of Fourier transform of n 00 vs. rotationally scaled frequency for L = 21 sites and η = 10. Note that the frequencies that emerged during η = 5 have persisted. The time scale τ also describes the decay of physically measurable quantities, for example the structure factors as defined in Eq. (53) and illustrated in Fig. 3. We show the emergent time scale τ for various lattice heights and systems sizes in Table 3.
Examining Fig. 2, one observes that the driven system approaches a dynamical equilibrium that is a mixture of rotational levels. The time scale with which the system relaxes to this equilibrium, τ , cannot be determined from the single-molecule physics, and so we refer to τ as an emergent time scale. For the low lattice height η = 1, the populations of the first two rotational states appear to oscillate around and asymptotically converge to roughly quarter filling, with J = 1 being lower due to contributing to population of J = 2 via an off-resonant AC coupling (Fig. 2(a)). For η = 5, the asymptotic equilibrium is an uneven mixture of rotational states that favors occupation of the J = 0 to the η = 1 case (Fig. 3(e)). Note also that S (01) π is now nonzero, and is periodic with the Rabi frequency Ω 00 at short times and twice the Rabi frequency at long times (see also Figs. 5(d) and 5(b)).  state ( Fig. 4(a)), and the emergent time scale for reaching this equilibrium is shorter than it was for η = 1 by roughly a factor of four. As the lattice height is then increased to η = 10, the populations return to the trend of η = 1, again converging to quarter filling with a time scale comparable to that of η = 1 (Fig. 4(c)). This illustrates the fact that the emergent time scale τ is not, in general, a monotonic function of the parameters of the lattice.
While the dynamics of the site-averaged rotational state populations are superficially similar for η = 1 and η = 10, the underlying physics is not identical, as can be seen by comparing Figs. 2(f), 4(b), and 4(d). These figures display the squared modulus of the Fourier transform of the site-averaged number in the J = 0 state. The only significant frequency observed for η = 1 is the Rabi frequency Ω ∼ 0.064B/h. In contrast, the η = 5 case has numerous other characteristic frequencies. As we raise the lattice height to η = 10, the frequencies that arose for η = 5 remain, even though the overall visual trend of the site-averaged number reflects that of the single-frequency η = 1 behavior. While we do not explicitly see the new frequencies in the site-averaged number, we do see them in the structure factors.  Table 3. Emergent time scales τ and τ Q and their fit asymptotic standard errors for various lattice heights and system sizes.
An example is Fig. 5(b), which clearly displays the 2Ω frequency behavior of the correlation structure factor S (01) π for η = 10. This frequency, which we easily pick out in the site-averaged number's Fourier transform, can also be seen in the Fourier transform of S (01) π , see Fig. 5(d). We find that the emergent time scale τ does not depend strongly on the size of the system L, even though the distribution of molecules on the lattice is, in general, quite different for different numbers of sites, as can be seen by comparing Figs. 2(a) and 2(e). Examining Fig. 2(c) and Table 3, the L = 10 case has a smaller τ than either of the odd L cases. We think this has to do with the filling being exactly 1/2 and not, strictly speaking, with the number of lattice sites, as the L = 9 and L = 21 cases have fillings less than 1/2. We see this clearly by comparing Fig. 4 with Figs. 2(a), 2(c), and 2(e). Fig. 4 displays n 0 0 /N, a quantity which is independent of filling but dependent, in general, on the number of lattice sites. There is a weak dependence on the number of lattice sites. On the other hand, Figs. 2(a), 2(c), and 2(e) display n 00 /L, a quantity which is independent of the number of lattice sites but dependent, in general, on the filling. There is a marked difference between L = 10, which has filling of 5/10 = 1/2 and the others, which have fillings< 1/2, but there is not a significant difference between L = 9 and L = 21, which have fillings of 4/9 and 10/21, respectively.
The dependence of τ on the filling is also evidenced by the correlation structure factor S (01) π in Fig. 3(f), which shows that there is a stronger correlation between the J = 0 and J = 1 states for exactly half filling than for fillings less than half, regardless of the system size. Half filling is known to be important in the extended Bose Hubbard model, where it marks the introduction of the charge density wave phase. We thus interpret this greater correlation structure factor as the appearance of a dynamic charge density wave phase between rotational states at half filling. This is in contrast to the usual behavior, where the structure factors S (00) π and S (11) π are nonzero whenever there is nonzero occupation of the particular rotational state and the structure factor S (01) π is much smaller-essentially zero, see Figs. 3(a) and 3(e). These results for the structure factors means that the J = 0 and J = 1 states tend to lie on top of one another, and not to "checkerboard" with a different rotational state occupying alternating sites. This is due to the fact that the Rabi flopping time scale is much shorter than the dipole-dipole time scale, meaning that the population cycles before there is sufficient time for the molecules to rearrange to a configuration which is energetically favorable with respect to the dipole-dipole term. However, because the population in each rotational level asymptotically reaches some nonzero value, we do see a small amount of rearrangement after many Rabi periods for any filling, corresponding to a nonzero S (01) π . Note that this rearrangement does not affect the site-averaged numbers, but rather the distribution of rotational states among the lattice sites. This asymptotic distribution emerges on time scales longer than we have considered, and is more prone to finite size effects than the site-averaged quantities, so we do not make a conjecture about it here.
We find that the Q-measure saturates as with a different time scale τ Q , see Fig. 4 and Table 3. We also find that the saturation time scale of the Q-measure is not, in general, a monotonic function of the lattice height η, as shown in Fig. 4.
(a) Dependence of the population damping time scale τ on the number of lattice sites. When we remove the dependence on the filling by dividing through by the total number, we see that there is little difference in the time scales with which systems of different size approach dynamic equilibrium. Contrast Figs. 2(a), 2(c), and 2(e), which display a profound dependence on filling when the dependence on lattice sites has been removed.
(b) Dependence of spatial entanglement on number of lattice sites. We see that systems of different size have different spatial entanglement in their static ground state. The time scale of the Q-measure saturation, τ Q , is shorter for L = 10 than it is for the odd L cases. This follows the general trend of τ and τ Q responding correspondingly to changes in the Hamiltonian parameters, and so we associate this shorter time scale partially with the filling, not entirely with the system size. This time scale is different from the time scale τ at which the populations approach an asymptotic equilibrium, though both time scales respond similarly to changes in the Hamiltonian parameter, see Table 3. For example, if τ Q gets larger as a parameter is changed then τ also gets larger, as illustrated in Figs. 4 and 7(b). The time scale τ Q displays a stronger dependence on the number of lattice sites L than τ , as can be seen in Figs. 4 and 4. This is because τ describes a quantity that has been averaged over sites, while τ Q does not.

Conclusions
We have presented and derived a novel lattice Hamiltonian, the Molecular Hubbard Hamiltonian (MHH). The MHH is a natural Hamiltonian for connecting theoretical studies of the dynamics of quantum phase transitions to near-term experimental setups using ultracold molecular gases. We presented a case study of this new Hamiltonian for hard core bosonic molecules at half filling. Starting from an initial condition of half filling in the J = 0, M = 0 state, we found that initial large (a) Dependence of spatial entanglement on lattice height. Note that the spatial entanglement and its associated time scale are not monotonic functions of the lattice height. Note also that the entanglement of the static ground state appears to be largely insensitive to the lattice height.
(b) Dependence of the site-averaged number on the lattice height. Note that the emergent time scale τ is not a monotonic function of the lattice height. Note also that τ responds in the same way that τ Q does to changes in the lattice height. oscillations in the system self-damp to an asymptotic equilibrium which consists of a lattice height and filling-dependent spatially entangled superposition of dressed states. This occurs on an emergent time scale τ which can not be predicted from the single molecule theory. We showed that τ depends non-monotonically on lattice height, weakly on lattice size, and strongly on filling (as apparent in simulations with odd and even numbers of sites). We also discovered a separate emergent time scale τ Q which describes how quickly the many body spatial entanglement saturates. We demonstrated that τ Q and τ respond similarly to changes in the Hamiltonian parameters and that τ Q depends on the filling, the lattice size, and, non-monotonically, on the lattice height. In addition to these emergent time scales, we studied the time-dependent structure factors and their frequency-domain Fourier transforms.
In future studies we will consider different filling factors, DC field strength to rotation ratios β DC , and initial conditions, as well as polarized and unpolarized spin-1/2 fermionic molecules. In addition, we will use multiscale methods to study how the emergent time scale demonstrated above compares to experimental time scales for physical systems, and thereby make quantitative predictions for experiments.
We acknowledge useful discussions with Deborah Jin, Heather Lewandowski, and Jun Ye. This work was supported by the National Science Foundation under Grant PHY-0547845 as part of the NSF CAREER program.

Relationship between operators in space-fixed and molecule-fixed coordinate systems
It is well known that the representation of the angular momentum operators in a molecule-fixed coordinate frame lead to the anomalous commutation relations [J i , J k ] = −ihǫ ijk J k [48]. The simplest way to avoid this trouble is to transform all expressions into the space-fixed frame where the angular momentum operators satisfy the normal commutation relations [J i , J k ] = ihǫ ijk J k [49]. If the moleculefixed axes are obtained by rotation of the space-fixed axes through the Euler angles {φ, θ, χ} [28] (which we collectively abbreviate as (R)), then the component of a k th -rank spherical tensor T that has projection p along the space-fixed z axis, denoted (T ) (k) p , can be expressed in terms of the molecule fixed components as where D (k) pq (R) ⋆ is the complex conjugate of the pq element of the k th -rank rotation matrix (Wigner D-matrix). To avoid confusion, we will label all space-fixed components with the letter p and all molecule-fixed components with q. From the orthogonality of the rotation matrices we have the inverse relationship

Rotational Hamiltonian
In the rigid rotor approximation the rotational Hamiltonian is simplŷ where we have defined the rotational constant B ≡ 1/2µr 2 e , with µ the molecule's reduced mass and r e its equilibrium internuclear separation. Typical values of B are ∼ 60h GHz [50]. This Hamiltonian has eigenvalues BJ (J + 1) and eigenstates |JM , with J the total angular momentum and M its projection along the internuclear axis.

DC Field Term
The dipole moment of a polar molecule in a rotational eigenstate is zero in an average sense due to the spherical symmetry of the rotational Hamiltonian. We break this symmetry by introducing a DC electric field along the space-fixed z axis, with Hamiltonian where E DC is the electric field amplitude. The field defines the spherical space-fixed axis p = 0, and the molecule-fixed internuclear axis defines q = 0. We transform between them using a first-rank rotation matrix as outlined above: The matrix elements of the DC Hamiltonian in the basis which diagonalize the rotational Hamiltonian Eq. (A.5) are where we use the notation (. . .) for the Wigner 3-j symbol [28]. Note that the symbol d refers to the permanent dipole moment of a molecule, and is not to be confused with the dipole operator denoted bŷ d. We refer to the basis which simultaneously diagonalizes the Rotational and DC Hamiltonians as the "dressed basis," and we denote the kets that span this basis by |E; JM , where the labels J and M are the zero field values of the corresponding quantum number and the symbol E is a reminder that these kets are superpositions of field free rotational states and DC field.
The effects of the DC field can be clearly seen by considering the dressed state wavefunctions, energies, and dipole moments to lowest order in perturbation theory in the dimensionless parameter β DC ≡ dE DC /B, the ratio of the field energy to the rotational level splitting: , (A.9) where ∆E (2) JM is the lowest non-zero shift in the energy. The DC field mixes states of different J, breaking the (2J + 1)-fold degeneracy of the rotational Hamiltonian, and so J is no longer a good quantum number. In the case of a z-polarized field, M remains a good quantum number, and a degeneracy persists for all states with the same |M|. This mixing aligns the molecule with the field, inducing a nonzero dipole moment. This means of orienting polar molecules, known as "brute force" orientation, works well for molecules that both have a large dipole moment and can be efficiently rotationally cooled [51]. While more effective means of orienting molecules using intense laser fields are known [52], they complicate the theoretical discussion and the experimental setup, and so we do not consider them here.
In larger fields the rotational levels become deeply mixed, which allows states that are weak-field seeking in low fields to become high-field seeking in high fields [53]. The actual mixing of rotational levels vs. β DC is depicted in Fig. A2 for the lowest three dressed levels. We note that there always exists a field E R such that the lowest R dressed states' dipole moments are all positive, as this is important to ensure the stability of a collection of dipoles. The universal curve of the induced dipole moments (in units of d) vs. β DC of the first two dressed rotational manifolds are shown in Figure 1(a). The universal curve of the dressed state energies energies (in units of B) vs. β DC is shown in Figure 1(b). For reference, β DC = 1 corresponds to a field of roughly 1.93 kV cm for B ∼ 60h GHz and d ∼ 9 D. Expanding the field operators in Eq. (1) in a Wannier basis of dressed states centered at a particular discrete position r i as described in Eq. (4), we find where E JM is the energy of the |E; J, M dressed state (see Fig. 2(a)) andn E,JM is the number operator associated with this same state. If the DC field were aligned at a small angle θ a to the z field of the trap (say, in the xz plane), then small dipole moments mixing M ′ = M ± 1 states would arise and the M ′ = M dipoles would decrease slightly (we can view them as being in an effective field of E eff = cos θ a E DC ). Treating the new contribution perturbatively in the small parameter sin θ a β DC , we find the lowest order couplings to the ground state  Figure A1. Dressed state dipole moments and energies. Note that the J = 1, M = 0 resonant dipole moment changes from weak-field seeking to high-field seeking at β DC ≈ 5. All rotational states have a field where this transition occurs, and the dipole tends monotonically towards unity after this field. The 10|d|00 dipole moment (and all transition dipole moments, generically) tends towards zero monotonically as β DC increases. Note also that the energetic differences between rotational levels are smallest at zero field and grow monotonically thereafter.
(a) Composition of 1 st dressed state.
(b) Composition of 2 nd dressed state.
(c) Composition of 3 rd dressed state. Figure A2. Compositions of dressed states vs. scaled rotational energy. The states become deeply mixed in large fields, and that the dressed state |E; JM whose zero field value is |JM does not always have the greatest overlap with |JM for all β DC . The field strength where the first dressed state changes from weak-field to high-field seeking, β DC = 5, is also roughly the place where its overlap with the |00 field-free level is greater than the overlaps with all other field-free levels.
and associated timescale τ θa for occupation of M = 0 states from the ground state,

AC Field Term
An AC microwave field of frequency ω resonantly drives transitions between two DC dressed states |E; J ′ M ′ and |E; JM with energy difference (E J ′ M ′ − E JM ) /h ≈ ω provided the induced dipole moment E; J ′ M ′ |d|E; JM is nonzero. Two states separated by an energy difference ∆E that is offresonant from the driving field (i.e. ∆E ≫ ω) will also be coupled, albeit much more weakly. In our system we resonantly couple the lowest two dressed rotational levels, |E; 10 and |E; 00 . We consider the case of z polarization, in which the effective Hamiltonian in the dressed Wannier basis iŝ H AC (t) = −π sin (ωt) JM Ω JM â † E;J,Mâ E;J+1,M + h.c , (A.14) where is the Rabi frequency. This is the frequency with which the populations of a two-level system cycles.
In experiments, the AC field has spatial curvature on the order of cm which is negligible on the µm system size scale.
In the absence of couplings between sites, the physics of the system is determined by the on-site, single-molecule physics. The percentage population of each component in both the |E; J, M dressed and |JM field-free bases are shown below for one Rabi period. In these plots only the |E;   Numerically, we must have a finite upper bound to the sum in Eq. (B.1), which we call J cut . This does not cause difficulty in practice, as the overlap of a dressed state |E; JM with a field-free state |J ′ M diminishes rapidly as J ′ differs more greatly from J. We find the coefficients in Eq. (B.1), as well as the dressed state energies and dipole moments by simultaneously diagonalizing the rotational and DC field Hamiltonians in a basis consisting of the first J cut rotational levels. Because TEBD scales poorly with the on-site dimension, we form as small an on-site basis as possible by keeping the eigenvectors corresponding to the R lowest dressed levels. To form a proper basis, we must renormalize these eigenvectors (which, for z-polarized field, does not change their orthogonality). We now demonstrate the convergence of these two procedures To show convergence of the first procedure, we plot the difference between the energy of the J th rotational state calculated for a particular value of J cut = i and one higher value, ∆E J (i) as a function of i. The results for various field strengths are shown in Figures 1(a)-1(b). We see very fast convergence for the low fields (e.g. β DC = 1.9) of interest. In our numerics we use J cut = 25, which ensures convergence for any of the β DC considered. Figure B1. Convergence with respect to DC dressing rotational state cutoff. As few as 7 field-free levels are needed for the weak field β DC = 1.9 to have the dressed state energies of interest converge to machine precision (left panel), and even a large DC field β DC = 20 requires only 12 field-free levels for the energy to converge (right panel).
To determine convergence with respect to the second procedure, examine Figs. 2(a)-2(b), which show the amount of the total dressed wave function norm | E; J0|E; J0 | 2 that lies outside of the first R fieldfree rotational levels for R = 3 and R = 4, respectively. For R = 4 the renormalization of the first three rotational levels is a very small effect for the β DC we consider, and the fourth level is not populated to any appreciable extent during time evolution for any β DC (see Fig.3(a)), so we expect that keeping the R = 4 lowest levels will give sufficient accuracy. By direct simulation, we find six digit accuracy in the suite of quantum measures defined in Sec. 3.2; specifically, we compare R = 3 to R = 4.

Many Body Considerations
There are also convergence issues that are inherent to the TEBD algorithm. The first, called the Schmidt error, is the error that arises from truncating the Hilbert space at each time step. We can Figure B2. Convergence with respect to local dimension cutoff. Dressed states with greater J lose more of their norm in truncation, as mixing occurs most strongly with adjacent J. Also, as the field is increased, the states become more deeply mixed, and so all states lose more of their norm. Truncating the local basis at the J = 3 dressed level incurs at most a 1% loss of norm for any of the states that are appreciably populated during time evolution (right panel).
parameterize the error per step in terms of the entanglement cutoff parameter χ as where λ [l] is a vector containing the eigenvalues of the reduced density matrix obtained by tracing over all sites but l, and α l is the local index that entangles the site l with the rest of the system, with smaller α l states having greater weight. We find that, among the measures we use, the one that is the most sensitive to χ is the Q-measure, which we plot for four values of χ in Fig. B3. Increasing χ improves the accuracy over longer times, but there is always a time after which the measure begins to deviate. This is the normalization drift alluded to in Sec. 3.1. The χ-dependent time after which the Schmidt error dominates is referred to as the runaway time [44]. In the case study of Sec. 4, we used χ = 50 for all simulations, which gives the Q-measure accurately to within four decimal places over the time scales considered. The second intrinsic source of error in TEBD is due to the Trotter-Suzuki expansion of the propagator [37]. We parameterize this error in terms of δt, the time step. When we halve the time step from that used in the simulations above (= 2π/(133ω)), we find no change in the measures to the ninth digit. It is clear that the Schmidt error discussed above is the chief source of error in our simulations.
To extract the emergent time scales defined in Eqs. (63) and (65), we used two different methods. The first is the nonlinear curve fitting routine "fit" in gnuplot. The second is the "NonlinearRegression" package in Mathematica 6.0. Both methods use nonlinear regression, which fits the data to a specified nonlinear function of the model parameters. The goodness of the fit is quantified by the asymptotic standard errors of the model parameters, which gives the standard deviation of each parameter. A low percent asymptotic error means that the model parameters cannot be adjusted very far without noticeably changing the goodness-of-fit. Both gnuplot and Mathematica returned the same values for the emergent time scales to within the stated asymptotic standard error. Figure B3. Convergence with respect to entanglement cutoff parameter. The left figure shows the spatial entanglement measure Q for various values of the TEBD entanglement cutoff parameter χ. As χ is increased, Q remains close to its true value for longer. In the right figure we plot the log of the absolute difference in Q for two values of χ divided by its arithmetic mean. We see at least four-digit accuracy for the largest values of χ we consider. Note also that even small values of χ are accurate for short times.