QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems. Part II: bosons, fermions and higher spins

We present a major update to QuSpin, SciPostPhys.2.1.003 -- an open-source Python package for exact diagonalization and quantum dynamics of arbitrary boson, fermion and spin many-body systems, supporting the use of various (user-defined) symmetries in one and higher dimension and (imaginary) time evolution following a user-specified driving protocol. We explain how to use the new features of QuSpin using seven detailed examples of various complexity: (i) the transverse-field Ising model and the Jordan-Wigner transformation, (ii) free particle systems: the SSH model, (iii) the many-body localized Fermi-Hubbard model, (iv) the Bose-Hubbard model in a ladder geometry, (v) nonlinear (imaginary) time evolution and the Gross-Pitaevskii equation, (vi) integrability breaking and thermalizing dynamics in the translationally-invariant 2D transverse-field Ising model, and (vii) out-of-equilibrium Bose-Fermi mixtures. This easily accessible and user-friendly package can serve various purposes, including educational and cutting-edge experimental and theoretical research. The complete package documentation is available under http://weinbe58.github.io/QuSpin/index.html.


What Problems can I Study with QuSpin?
Understanding the physics of many-body quantum condensed matter systems often involves a great deal of numerical simulations, be it to gain intuition about the complicated problem of interest, or because they do not admit an analytical solution which can be expressed in a closed form. This motivates the development of open-source packages [1,2,3,4,5,6,7,8,9,10,11,12,13], the purpose of which is to facilitate the study of condensed matter systems, without the need to implement the inner workings of complicated numerical methods which required years to understand and fully develop. Here, we report on a major upgrade to QuSpin [14] a Python library for exact diagonalisation (ED) and simulation of the dynamics of arbitrary quantum many-body systems.
Although ED methods are vastly outperformed by more sophisticated numerical techniques in the study of equilibrium problems, such as quantum Monte Carlo [15,16,17], matrix product states based density matrix renormalisation group [18,19,20], and dynamical mean-field theory [21,22,23], as of present date ED remains essential for certain dynamical non-equilibrium problems. The reason for this often times relies on the fact that the underlying physics of these problems cannot be explained without taking into consideration the contribution from high-energy states excited during the nonequilibrium process. Some prominent examples of such problems include the study of the many-body localisation (MBL) transition [24,25], the Eigenstate Thermalisation hypothesis [26], quantum quench dynamics [27], periodicallydriven systems [28,29,30], adiabatic and counter-diabatic state preparation [31], applications of Machine Learning to non-equilibrium physics [32,33,34], optimal control [35], and many more.
It is, thus, arguably useful to have a toolbox available at one's disposal which allows one to quickly simulate and study these and related nonequilibrium problems. As such, QuSpin offers easy access to performing numerical simulations, which can facilitate the development and inspiration of new ideas and the discovery of novel phenomena, eliminating the cost of spending time to develop a reliable code. Besides theorists, the new version of QuSpin will hopefully even prove valuable to experimentalists working on problems containing dynamical setups, as it can help students and researchers focus on perfecting the experiment, rather than worry about writing the supporting simulation code. Last but not least, with the computational processing power growing higher than ever before, the role played by simulations for theoretical research becomes increasingly more important too. It can, therefore, be expected that in the near future quantum simulations become an integral part of the standard physics university curriculum, and having easily accessible toolboxes, such as QuSpin, is one of the required prerequisites for this anticipated change.
2 How do I Use the New Features of QuSpin?
New in QuSpin 2.0, we have added the following features and toolboxes: (i) support for fermion, boson and higher-spin Hamiltonians with basis.ent_entropy and basis.partial_trace routine to calculate entanglement entropy, reduced density matrices, and entanglement spectrum (ii) general basis constructor classes for user-defined symmetries.
(iii) blocks class to automatically handle splitting the evolution over various symmetry sectors of the Hilbert space.
(iv) user-customisable evolve() routine to handle user-specified linear and non-linear equations of motion.
(vi) quantum_LinearOperator class which applies the operator "on the fly" which saves vast amounts of memory at the cost of computational time.
Before we carry on, we refer the interested reader to Examples 0-3 from the original QuSpin paper [14]. The examples below focus predominantly on the newly introduced features, and are thus to be considered complementary. We emphasize that, since they serve the purpose of explaining how to use QuSpin, for the sake of brevity we shall not discuss the interesting physics related to the interpretation of the results.
Installing QuSpin is quick and efficient; just follow the steps outlined in App. A.

The Spectrum of the Transverse Field Ising Model and the Jordan-Wigner Transformation
This example shows how to • construct fermionic hopping, p-wave pairing and on-site potential terms, and spin−1/2 interactions and transverse fields, • implement periodic and anti-periodic boundary conditions, • use particle conservation modulo 2, spin inversion, reflection, and translation symmetries, • handle the default built-in particle conservation and symmetry checks, • obtain the spectrum of a QuSpin Hamiltonian.
Physics Setup-The transverse field Ising (TFI) chain is paradigmatic in our understanding of quantum phase transitions, since it represents an exactly solvable model [36]. The Hamiltonian is given by where the nearest-neighbour (nn) spin interaction is J, h denotes the transverse field, and σ α j are the Pauli spin-1/2 matrices. We use periodic boundary conditions and label the L lattice sites 0, . . . , L − 1 to conform with Python's convention. This model has gapped, fermionic elementary excitations, and exhibits a phase transition from an antiferromagnet to a paramagnet at (h/J) c = 1. The Hamiltonian possesses the following symmetries: parity (reflection w.r.t. the centre of the chain), spin inversion, and (many-body) momentum conservation. In one dimension, the TFI Hamiltonian can be mapped to spinless p-wave superconducting fermions via the Jordan-Wigner (JW) transformation [36,37,38]: where the fermionic operators satisfy {c i , c † j } = δ ij . The Hamiltonian is readily shown to take the form In the fermionic representation, the spin zz-interaction maps to nn hopping and a p-wave pairing term with coupling constant J, while the transverse field translates to an on-site potential shift of magnitude h. In view of the implementation of the model using QuSpin, we have ordered the terms such that the site index is growing to the right, which comes at the cost of a few negative signs due to the fermion statistics. We emphasize that this ordering is not required by QuSpin, but it is merely our choice to use it here for the sake of consistency. The fermion Hamiltonian posses the symmetries: particle conservation modulo 2, parity and (many-body) "momentum" conservation.
Here, we are interested in studying the spectrum of the TFI model in both the spin and fermion representations. However, if one naïvely carries out the JW transformation, and computes the spectra of Eqs. (1) and (3), one might be surprised that they do not match exactly. The reason lies in the form of the boundary condition required to make the JW mapping exact -a subtle issue often left aside in favour of discussing the interesting physics of the TFI model itself.
We recall that the starting point is the periodic boundary condition imposed on the spin Hamiltonian in Eq. (1). Due to the symmetries of the spin Hamiltonian (1), we can define the JW transformation on every symmetry sector separately. To make the JW mapping exact, we supplement Eq. (2) with the following boundary conditions: (i) the negative spin-inversion symmetry sector maps to the fermion Hamiltonian (3)   maps to the fermion Hamiltonian (3) with anti-periodic boundary conditions (APBC) and even total number of fermions. Anti-periodic boundary conditions differ from PBC by a negative sign attached to all coupling constants that cross a single, fixed lattice bond (the bond itself is arbitrary as all bonds are equal for PBC). APBC and PBC are special cases of the more general twisted boundary conditions where, instead of a negative sign, one attaches an arbitrary phase factor. In the following, we show how to compute the spectra of the Hamiltonians in Eqs. (1) and (3) with the correct boundary conditions using QuSpin. Figure 1 shows that they match exactly in both the PBC and APBC cases discussed above, as predicted by theory. Code Analysis-We begin by loading the QuSpin operator and basis constructors, as well as some standard Python libraries.
1 from quspin.operators import hamiltonian # Hamiltonians and operators 2 from quspin.basis import spin_basis_1d, spinless_fermion_basis_1d # Hilbert space spin basis 3 import numpy as np # generic math functions 4 import matplotlib.pyplot as plt # plotting library First, we define the models parameters. 6 ##### define model parameters ##### 7 L=8 # system size 8 J=1. 0 # spin zz interaction 9 h=np.sqrt(2) # z magnetic field strength We have to consider two cases when computing the spectrum, as discussed in the theory section above. In one case, the fermionic system has PBC and the particle number sector is odd, while the spins are constrained to the negative spin inversion symmetry sector, while in the second -the fermion model has APBC with even particle number sector, and the spin model is considered in the positive spin inversion sector. To this end, we introduce the variables zblock ∈ {±1} and PBC ∈ {±1}, where PBC = −1 denotes APBC. Note that the only meaningful combinations are (zblock, PBC) = (−1, 1), (1, −1), which we loop over: 11  Within this loop, the code is divided in two independent parts: first, we compute the spectrum of the TFI system, and then -that of the equivalent fermionic model. Let us discuss the spins.
In QuSpin, operators are stored as sparse lists. These lists contain two parts: (i) the lattice sites on which the operator acts together with the coupling strength, which we call a site-coupling list, and (ii) the types of the operators involved, i.e. the operator-string. Notice the way we defined the periodic boundary condition for the spin-spin interaction using the modulo operator %, which effectively puts a coupling between sites L − 1 and 0. We mention in passing that the above procedure generalises so one can define any multi-body local or nonlocal operator using QuSpin.
In order to specify the types of the on-site single-particle operators, we use operator strings. For instance, the transverse field operator O = g L−1 j=0 σ x j becomes ['x',h_field], while the two-body zz-interaction is ['zz',J_zz]. It is important to notice that the order of the letters in the operator string corresponds to the order the operators are listed in the site-coupling lists. Putting everything into one final list yields the static list for the spin model: In QuSpin, the user can define both static and dynamic operators. Since this example does not contain any time-dependent terms, we postpone the explanation of how to use dynamic lists to Sec. 2.5, and use an empty list instead. The last step before we can construct the Hamiltonian is to build the basis for it. This is done using the basis constructors. For spin systems, we use spin_basis_1d which supports the operator strings 'z','+','-','I', and for spin-1/2 additionally 'x','y'. The first and required argument is the number of sites L. Optional arguments are used to parse symmetry sectors. For instance, if we want to construct an operator in the spin-inversion block with quantum number +1, we can conveniently do this using the flag zblock=1. Having specified the static and dynamic lists, as well as the basis, building up the Hamiltonian is a one-liner, using the hamiltonian constructor. The required arguments in order of appearance are the static and dynamic lists, respectively. Optional arguments include the basis, and the precision or data type dtype. If no basis is passed, the constructor uses spin_basis_1d by default. The default data type is np.complex128. The Hamiltonian is stored as a Scipy sparse matrix for efficiency. It can be cast to a full array for a more convenient inspection using the attribute H.toarray(). To calculate its spectrum, we use the attribute H.eigenvalsh(), which returns all eigenvalues. Other attributes for diagonalisation were discussed here, c.f. Ref. [14]. 25 # calculate spin energy levels 26 E_spin=H_spin.eigvalsh() Let us now move to the second part of the loop which defines the fermionic p-wave superconductor. We start by defining the site-coupling list for the local potential Let us focus on the case of periodic boundary conditions PBC=1 first. 31 if PBC==1: # periodic BC: odd particle number subspace only In the fermion model, we have two types of two-body terms: hopping terms c † i c i+1 − c i c † i+1 , and pairing terms c † i c † i+1 − c i c i+1 . While QuSpin allows any ordering of the operators in the static and dynamic lists, for the sake of consistency we set a convention: the site indices grow to the right. To take into account the opposite signs resulting from the fermion statistics, we have to code the site-coupling lists for all four terms separately. These two-body terms are analogous to the spin-spin interaction above: To construct a fermionic operator, we make use of the basis constructor spinless_fermion_basis_1d. Once again, we pass the number of sites L. As we explained in the analysis above, we need to consider all odd particle number sectors in the case of PBC. This is done by specifying the particle number sector Nf. 37 # construct fermion basis in the odd particle number subsector The construction of the basis is the same, except that this time we need all even particle number sectors: 50 # construct fermion basis in the even particle number subsector 51 basis_fermion = spinless_fermion_basis_1d(L=L,Nf=range( 0,L+1,2)) As before, we need to specify the type of operators that go in the fermion Hamiltonian using operator string lists. The spinless_fermion_basis_1d class accepts the following strings '+','-','n','I', and additionally the particle-hole symmetrised density operator 'z' = n − 1/2. The static and dynamic lists read as Constructing and diagonalising the fermion Hamiltonian is the same as for the spin-1/2 system. Note that one can disable the automatic built-in checks for particle conservation check_pcon=False and all other symmetries check_symm=False if one wishes to suppress the checks.

Free Particle Systems: the Fermionic SSH Chain
This example shows how to • construct free-particle Hamiltonians in real space, • implement translation invariance with a two-site unit cell and construct the singleparticle momentum-space block-diagonal Hamiltonian using the block_diag_hamiltonian tool, • compute non-equal time correlation functions, • time-evolve multiple quantum states simultaneously. Physics Setup-The Su-Schrieffer-Heeger (SSH) model of a free-particle on a dimerised chain was first introduced in the seventies to study polyacetylene. Today, this model is paradigmatically used in one spatial dimension to introduce the concepts of Berry phase, topology, edge states, etc. The Hamiltonian is given by where {c i , c † j } = δ ij obey fermionic commutation relations. The uniform part of the hopping matrix element is J, the bond dimerisation is defined by δJ, and ∆ is the staggered potential. We work with periodic boundary conditions. Below, we show how one can use QuSpin to study the physics of free fermions in the SSH chain. One way of doing this would be to work in the many-body (Fock space) basis, see Sec. 2.3. However, whenever the particles are non-interacting, the exponential scaling of the Hilbert space dimension with the number of lattice sites imposes an artificial limitation on the system sizes one can do. Luckily, with no interactions present, the many-body wave functions factorise in a product of single-particle states. Hence, it is possible to study the behaviour of many free bosons and fermions by simulating the physics of a single particle, and populating the states according to bosonic or fermionic statistics, respectively.
Making use of translation invariance, a straightforward Fourier transformation to momentum space, a k = 2/L L−1 j:even e −ikj c j and b k = 2/L L−1 j:odd e −ikj c j , casts the SSH Hamiltonian in the following form where the reduced Brillouin zone is defined as BZ = [−π/2, π/2). We thus see that the Hamiltonian reduces further to a set of independent 2 × 2 matrices. The spectrum of the SSH model is gapped and, thus, has two bands, see Fig. 2a.
Since we are dealing with free fermions, the ground state is the Fermi sea, |FS , defined by filling up the lowest band completely. We are interested in measuring the real-space non-equal time density autocorrelation function To evaluate the correlator numerically, we shall use the right-hand side of this equation.
As we are studying free particles, it is enough to work with the single-particle states. For instance, the Fermi sea can be obtained as |FS = k≤k F c † k |0 . Denoting the on-site density operator by n i , one can cast the correlator in momentum space in the following form: If we want to consider finite-temperature β −1 , the above formula generalises to where n FD (k, β) = 1/(exp(βE k ) + 1) is the Fermi-Dirac distribution at temperature β −1 , with E k the SSH dispersion. Figure 2b) shows the time evolution of C ij (t, β) for two sites, separated by the maximal distance on the ring: L/2. Code Analysis.-Let us explain how one can do all this quickly and efficiently using QuSpin. As always, we start by loading the required packages and libraries.
1 from quspin.operators import hamiltonian,exp_op # Hamiltonians and operators 2 from quspin.basis import spinless_fermion_basis_1d # Hilbert space fermion basis 3 from quspin.tools.block_tools import block_diag_hamiltonian # block diagonalisation 4 import numpy as np # generic math functions 5 import matplotlib.pyplot as plt # plotting library 6 try: # import python 3 zip function in python 2 and pass if already using python 3 7 import itertools.izip as zip 8 except ImportError: 9 pass After that, we define the model parameters 1 ##### define model parameters ##### 2 L=1 0 0 # system size 3 J=1. 0 # uniform hopping 4 deltaJ= 0.1 # bond dimerisation 5 Delta= 0.5 # staggered potential 6 beta=1 0 0. 0 # inverse temperature for Fermi-Dirac distribution In the following, we construct the fermionic SSH Hamiltonian first in real space. We then show how one can also construct it in momentum space where, provided we use periodic boundary conditions, it appears bock-diagonal. Let us define the fermionic site-coupling lists. Once again, we emphasise that fermion systems require special care in defining the hopping terms: Eq. (4) is conveniently cast in the form where all site indices on the operators grow to the right, and all signs due to the fermion statistics are explicitly spelt out. To define the static list we assign the corresponding operator strings to he site-coupling lists. Since our problem does not feature any explicit time dependence, we leave the dynamic list empty. Setting up the fermion basis with the help of the constructor spinless_fermion_basis_1d proceeds as smoothly as in Sec. 2.1. Notice a cheap trick: by specifying a total of Nf=1 fermion in the lattice, QuSpin actually allows to define single-particle models, as a special case of the more general many-body Hamiltonians. Compared to many body models, however, due to the exponentially reduced Hilbert space size, this allows us to scale up the system size L.
We then build the real-valued SSH Hamiltonian in real space by passing the static and dynamic lists, as well as the basis and the data type. After that, we and diagonalise it, storing all eigenenergies and eigenstates.
1 # build real-space Hamiltonian 2 H=hamiltonian(static,dynamic,basis=basis,dtype=np.float64) 3 # diagonalise real-space Hamiltonian 4 E,V=H.eigh() For translation invariant single-particle models, however, the user might prefer to use momentum space, where the Hamiltonian becomes block diagonal, as we showed in the theory section above. This can be achieved using QuSpin's block_tools. The idea behind this tool is simple: the main purpose is to create the full Hamiltonian in block-diagonal form, where the blocks correspond to pre-defined quantum numbers. In our case, we would like to use momentum or kblock's. Note that the unit cell in the SSH model contains two sites, which we encode using the variable a=2. Thus, we can create a list of dictionaries blocks, each element of which defines a single symmetry block. If we combine all blocks, we exhaust the full Hilbert space. All other basis arguments, such as the system size, we store in the variable basis_args. We invite the interested user to check the package documentation for additional optional arguments and functionality of block_tools, cf. App. C. We mention in passing that this procedure is independent of the symmetry, and can be done using all symmetries supported by QuSpin, not only translation. To create the block-diagonal Hamiltonian, we invoke the block_diag_hamiltonian method. It takes both required and optional arguments, and returns the transformation which blockdiagonalises the Hamiltonian (in our case the Fourier transform) and the block-diagonal Hamiltonian object. Required arguments, in order of appearance, are the blocks, the static and dynamic lists, the basis constructor, basis_args, and the data type. Since we anticipate the matrix elements of the momenum-space Hamiltonian to contain the Fourier factors exp(−ik), we know to choose a complex data type. block_diag_hamiltonian also accepts some optional arguments, such as the flags for disabling the automatic built-in symmetry checks. More about this function can be found in the documentation, cf. App. C.
1 # construct block-diagonal Hamiltonian 2 FT,Hblock = block_diag_hamiltonian(blocks,static,dynamic,spinless_fermion_basis_1d, 3 basis_args,np.complex128,get_proj_kwargs=dict(pcon=True)) We can then use all functions and methods of the hamiltonian class to manipulate the blockdiagonal Hamiltonian, for instance the diagonalisation routines: Last, we proceed to calculate the correlation function from Eq. (6). To this end, we shall split the correlator according to the RHS of Eq. (6). Thus, the strategy is to evolve both the Fermi sea |FS(t) and the auxiliary state |nFS(t) in time, and compute the expectation value of the time-independent operator n i (0) in between the two states as a function of time. Keep in mind that we do not need to construct the Fermi sea as a many-body state explicitly, so we rather work with single-particle states.
The first step is to collect all momentum eigenstates into the columns of the array psi. We then build the operators n j=0 and n j=L/2 in real space.
1 ##### prepare the density observables and initial states ##### 2 # grab single-particle states and treat them as initial states 3 psi 0=Vblock 4 # construct operator n_1 = $n_{j= 0}$ Next, we transform these two operators to momentum space using the method rotate_by(). Setting the flag generator=False treats the Fourier transform FT as a unitary transformation, rather than a generator of a unitary transformation.
1 # transform n_j operators to momentum space 2 n_1=n_1.rotate_by(FT,generator=False) 3 n_2=n_2.rotate_by(FT,generator=False) Let us proceed with the time-evolution. We first define the time vector t and the state n_psi 0. We can perform the time evolution in two ways: (i) we calculate the time-evolution operator U using the exp_op class, and apply it to the momentum states psi 0 and n_psi 0. The exp_op class calculates the matrix exponential exp(aH) of an operator H multiplied by a complex number a. One can also conveniently compute a series of matrix exponentials exp(aHt) for every time t by specifying the stating point start, endpoint stop and the number of elements num which define the time array t via t=np.linspace(start,stop,num). Last, by parsing the flag iterate=True we create a python generator -a pre-defined object evaluated only at the time it is called, i.e. not pre-computed, which can save both time and memory. 1 # construct time-evolution operator using exp_op class (sometimes faster) 2 U = exp_op(Hblock,a=-1j,start=t.min(),stop=t.max(),num=len(t),iterate=True) Another way of doing the time evolution, (ii), is to use the evolve() method of the hamiltonian class. The idea here is that every Hamiltonian defines a generator of time translations. This method solves the Schrödinger equation using SciPy's ODE integration routines, see App. C for more details. The required arguments, in order of appearance, are the initial state, the initial time, and the time vector. The evolve() method also supports the option to create the output as a generator using the flag iterate=True. Both ways (i) and (ii) time-evolve all momentum states psi at once, i.e. simultaneously.
To evaluate the correlator, we first preallocate memory by defining the empty array correlators, which will be filled with the correlator values in every single-particle momentum mode |k . Using generators allows us to loop only once over time to obtain the time-evolved states psi_t and n_psi_t. In doing so, we evaluate the expectation value FS(t)|n i (0)|nFS(t) using the matrix_ele() method of the hamiltonian class. The flag diagonal=True makes sure only the diagonal matrix elements are calculated 1 and returned as a one-dimensional array. Finally, we weigh all singe-state correlators by the Fermi-Dirac distribution to obtain the finite-temperature non-equal time correlation function C 0,L/2 (t, β).
The complete code including the lines that produce Fig. 2

Many-Body Localization in the Fermi-Hubbard Model
This example shows how to: • construct Hamiltonians for spinful fermions using the spinful_fermion_basis_1d class, • use quantum_operator class to construct Hamiltonians with varying parameters, • use new basis.index() function to construct Fock states from strings, • use obs_vs_time() functionality to measure observables as a function of time.
A class of exciting new problems in the field of non-equalibrium physics is that of manybody localised (MBL) models. The MBL transition is a dynamical phase transition in the eigenstates of a many-body Hamiltonian. Driven primarily by quenched disorder, the transition is distinguished by the system having ergodic eigenstates in the weak disorder limit and non-ergodic eigenstates in the strong disorder limit. The MBL phase is reminiscent of integrable systems as one can construct quasi-local integrals of motion, but these integrals of motion are much more robust in the sense that they are not sensitive to small perturbations, as is the case in many classes of integrable systems [24,25,39,40,41].
Motivated by recent experiments [42,43,44,45,46,47,48,49,50,51] we explore MBL in the context of fermions using QuSpin. The model we will consider is the Fermi-Hubbard model with quenched random disorder which has the following Hamiltonian: where c iσ and c † iσ are the fermionic creation and annihilation operators on site i for spin σ, respectively. We will work in the sector of 1/4 filling for both up and down spins. The observable of interest is the sublattice imbalance [43,45,49]: where A and B refer to the different sublattices of the chain and N is the particle number. We prepare an initial configuration of fermions of alternating spin on every other site. Evolving this initial state under the Hamiltonian (9), we calculate the time dependence of the imbalance I(t) which provides useful information about the ergodicity of the Hamiltonian H (or the lack thereof). If the Hamiltonian is ergodic, the imbalance should decay to zero in the limit t → ∞, as one would expect due to thermalising dynamics. On the other hand, if the Hamiltonian is MBL then some memory of the initial state will be retained at all times, and therefore the imbalance I(t) will remain non-zero even at infinite times.
Because the Hilbert space dimension grows so quickly with the lattice size for the Fermi-Hubbard Hamiltonian, we only consider the dynamics after a finite amount of time, instead of looking at infinite-time expectation values which require knowledge of the entire basis of the Hamiltonian, which requires more computational resources. Code Analysis-We start out by loading a set of libraries which we need to proceed with the simulation of MBL fermions. 3 from quspin.tools.measurements import obs_vs_time # calculating dynamics 4 import numpy as np # general math functions 5 from numpy.random import uniform,choice # tools for doing random sampling 6 from time import time # tool for calculating computation time While we already encountered most of these libraries and functions, in this example we introduce the quantum_operator class which defines an operator, parametrized by multiple parameters, as opposed to the Hamiltonian which is only parametrized by time. Also, since this example requires us to do many different disorder realisations, we use NumPy's random number library to do random sampling. We load uniform to generate the uniformly distributed random potential as well as choice which we use to estimate the uncertainties of the disorder averages using a bootstrap re-sampling procedure which we explain below. In order to time how long each realization takes, we use the time function from python's time library. After importing all the required libraries and functions we set up the parameters for the simulation including the number of realizations n_real, and the physical couplings J, U, the number of up and down fermions, etc. Next we set up the basis, introducing here the spinful_fermion_basis_1d constructor class.
It works the usual way, except the particle number sector N_f is now a tuple containing the number of up and down fermions. 26 ###### create the basis 27 # build spinful fermions basis 28 basis = spinful_fermion_basis_1d(L,Nf=(N_up,N_down)) The next step in the procedure is to set up the site-coupling and operator lists: Notice here that the "|" character is used to separate the operators which belong to the up (left side of tensor product) and down (right side of tensor product) Hilbert spaces in basis. If no string in present the operator is assumed to be the identity 'I'. The site-coupling list, on the other hand, does not require separating the two sides of the tensor product, as it is assumed that the operator string lines up with the correct site index when the '|' character is removed.
In the last couple of lines defining model (see below), we create a python dictionary object in which we add static operator lists as values, indexed by a particular string known as the corresponding key. This dictionary, which we refer to as operator_dict, is then passed into the quantum_operator class which, for each key in operator_dict, constructs the operator listed in the site-coupling list for that key. When one wants to evaluate this operator for a particular set of parameters, one uses a second dictionary (see params_dict below in lines 74 − 76) with the same keys as operator_dict: the value corresponding to each key in the params_dict multiples the operator stored at that same key in operator_dict. This allows one to parametrize many-body operators in more complicated and general ways. In the present example, we define a key for the Fermi-Hubbard Hamiltonian, and then, as we need to change the disorder strength in between realisations, we also define keys for the local density operator for both up and down spins on each site. By doing this, we can then construct any disordered Hamiltonian by specify the disorder at each site, and changing its value from one realisation to another. 46 # create operator dictionary for quantum_operator class 47 # add key for Hubbard hamiltonian 48 operator_dict=dict(H 0=operator_list_ 0) 49 # add keys for local potential in each site 50 for i in range(L): 51 # add to dictioanry keys h 0,h1,h2,...,hL with local potential operator The quantum_operator class constructs operators in almost an identical manner as a hamiltonian() class with the exception that there is no dynamic list. Next, we construct our initial state with the fermions dispersed over the lattice on every other site. To get the index of the basis state which this initial state corresponds to, one can use the index function of the tensor_basis class. This function takes a string or integer representing the product state for each of the Hilbert spaces and then searches to find the full product state, returning the corresponding index. We can then create an empty array psi_ 0 of dimension the total Hilbert space size, and insert unity at the index corresponding to the initial state. This allows us to easily define the many-body product state used in the rest of the simulation. Now that the operators are all set up, we can proceed with the simulation of the dynamics. First, we define a function which, given a disorder realization of the local potential, calculates the time evolution of I(t). We shall guide the reader through this function step by step. The syntax for this begins as follows: using the tohamiltonian() method of H_dict class, which accepts as argument the dictionary params_dict, which shares the same keys as operator_dict, but whose values are determined by the disorder list which changes from one realisation to another.
Once the Hamiltonian has been constructed we want to time-evolve the initial state with it. To this end, we use the fact that, for time-independent Hamiltonians, the time-evolution operator coincides with the matrix exponential exp(−itH). In QuSpin, a convenient way to define matrix exponentials is offered by the exp_op class. Given an operator A , it calculates exp(atA) for any complex-valued number a. The time grid for t is specified using the optional arguments start, stop and num. If the latter are omitted, default is t = 1. The exp_op objects are either a list of arrays, the elements of which correspond to the operator exp(atA) at the times defined or, if iterate=True a generator for this list is returned. Here, we use exp_op to create a generator list containing the time-evolved states as follows To calculate the expectation value of the imbalance operator in time, we use the obs_vs_time function. Since the usage of this function was extensively discussed in Example 2 and Example 3 of Ref. [14], here we only mention that it accepts a (collection of) state(s) [or a generator] psi_t, a time vector t, and a dictionary dict(I=I), whose values are the observables to calculate the expectation value of. The function returns a dictionary, the keys of which correspond to the keys every observable was defined under. Now we are all set to run the disorder realizations for the different disorder strengths which in principle can be spilt up over multiple simulations, e.g. using joblib [c.f. Example 1 from Ref. [14]] but for completeness we do all of the calculations in one script. Last, we calculate the average and its error bars using bootstrap re-sampling. The idea is that we have a given set of n_real samples from which we can select randomly a set of n_real samples with replacements. This means that sometimes we will get the same individual sample represented multiple times in one of these sets. Then by averaging over this set one can get an estimate of the mean value for the original sample set. This mean value is what is called a bootstrap sample. In principle there are a very large number of possible bootstrap samples so in practice one only takes a small fraction of the total set of bootstrap samples from which you can estimate things like the variance of the mean which is the error of the original sample set. The complete code including the lines that produce Fig. 3

The Bose-Hubbard Model on Translationally Invariant Ladder
This example shows how to: • construct interacting Hamiltonians for bosonic systems, • construct a Hamiltonian on a ladder geometry, • use the block_ops class to evolve a state over several symmetry sectors at once, • measure the entanglement entropy of a state on the ladder.
Physics Setup-In this example we use QuSpin to simulate the dynamics of the Bose-Hubbard model (BHM) on a ladder geometry. The BHM is a minimal model for interacting lattice bosons [36] which is most often experimentally realizable in cold atom experiments [52]. The Hamiltonian is given by where b i and b † i are bosonic creation and annihilation operators on site i, respectively, and ij denotes nearest neighbors on the ladder. We consider a half-filled ladder of length L with N = 2L sites and cylindrical boundary condition, i.e. a periodic boundary condition along the ladder-leg direction. We are interested in the strongly-interacting limit U/J 1, where the mean-field Gross-Pitaevskii theory fails. Hence, we restrict the local Hilbert space to allow at most two particles per site, effectively using a total of three states per site: empty, singly and doubly occupied.
The system is initialised in a random Fock state, and then evolved with the Hamiltonian (11). Since the BHM is non-integrable, we expect that the system eventually thermalizes, so that the long-time occupation becomes roughly uniform over the entire system. Besides the local density, we also measure the entanglement entropy between the two legs of the ladder.
As we consider a translational invariant ladder, the Hilbert space factorizes into subspaces corresponding to the different many-body momentum blocks. This is similar to what was discussed in Sec. 2.2 for the SSH model, but slightly different as we consider translations of the many-body Fock states as opposed to the single particle states [53]. In certain cases, it happens that, even though the Hamiltonian features symmetries, the initial state does not obey them, as is the case in the present example, where the initial state is a random Fock state. Thus, in this section we project the wavefunction to the different symmetry sectors and evolve each of the projections separately under the Hamiltonian for that symmetry sector. After the evolution, each of these symmetry-block wavefunctions is transformed back to the local Fock space basis, and summed up to recover the properly evolved initial state. We can then measure quantities such as the on-site density and the entanglement entropy. Figure 4  show the entanglement density between the two legs as a function of time, while Figure 4(b) shows a heat map of the local density of the bosonic gas. Both quantities show that after a short period of time the gas is completely thermalized. Code Analysis-Let us now show how one can simulate this system using QuSpin. Following the code structure of previous examples, we begin by loading the modules needed for the computation: 1 from __future__ import print_function, division #import python 3 functions 2 from quspin.operators import hamiltonian # Hamiltonians and operators 3 from quspin.basis import boson_basis_1d # bosonic Hilbert space 4 from quspin.tools.block_tools import block_ops # dynamics in symmetry blocks 5 import numpy as np # general math functions 6 import matplotlib.pyplot as plt # plotting library 7 import matplotlib.animation as animation # animating movie of dynamics First, we set up the model parameters defining the length of the ladder L, the total number of sites in the chain N=2*L, as well as the filling factor for the bosons nb, and the maximum number of states per site (i.e. the local on-site Hilbert space dimension) sps. The hopping matrix elements J ⊥ , J ,1 , and J ,2 (see Fig. 1 correspond to the python script variables J_perp, J_par_1, and J_par_2, respectively. The on-site interaction is denoted by U. 9 ##### define model parameters 10 # initial seed for random number generator 11 np.random.seed( 0) # seed is 0 to produce plots from QuSpin2 paper 12  We also define the hopping site-coupling list. In general, QuSpin can set up the Hamiltonian on any graph (thus, including higher-dimensions). By labelling the lattice sites conveniently, we can define the ladder geometry, as shown in Fig. 1 where we used the list_1.extend(list_2) method to concatenate two lists together 2 .
Next, we define the static and dynamic lists, which are needed to construct the Hamiltonian.  Instead of creating a hamiltonian class object in real space, we use the block_ops class to set up the Hamiltonian in a block-diagonal form in momentum space, similar to the SSH model, c.f. Sec. 2.2, but now genuinely many-body. In order to reduce the computational cost, the state is evolved in momentum space and projected to the Fock basis after the calculation. The purpose of block_ops is to provide a simple interface for solving the Schrödinger equation when an initial state does not obey the symmetries of the Hamiltonian it is evolved under. We have seen an example of this in Sec. 2.2 when trying to measure non-equal space-time correlation functions of local operators in a translational invariant system, while in this section we explicitly start out with a state which does not obey translation invariance. To construct the block_ops object we use the follow set of lines, explained below: 39 # create block_ops object 40 blocks=[dict(kblock=kblock) for kblock in range(L)] # blocks to project on to 41 baisis_args = (N,) # boson_basis_1d manditory arguments 42 basis_kwargs = dict(nb=nb,sps=sps,a=2) # boson_basis_1d optional args 43 get_proj_kwargs = dict(pcon=True) # set projection to full particle basis 44 H_block = block_ops(blocks,static,dynamic,boson_basis_1d,baisis_args,np.complex128, 45

basis_kwargs=basis_kwargs,get_proj_kwargs=get_proj_kwargs)
First, we create a list of dictionaries blocks which defines the different symmetry sectors to project the initial state to, before doing the time evolution 3 . The optional arguments basis_args and basis_kwargs apply to every symmetry sector. Last, get_proj_kwargs contains the optional arguments to construct the projectors 4 . This class automatically projects the initial state onto the different symmetry sectors and will evolve the sectors individually in serial or with an additional option in parallel over multiple cpu cores. After the evolution the class reconstructs the state and returns it back to the user. The Hamiltonian and the evolution itself will only be calculated for the sectors which have non-vanishing overlap with the initial state. Furthermore, the class will, by default, construct things on the fly as it needs them to save the user memory. For more information about this class we refer the user to the Documentation, c.f. App. C.

SciPost Physics Submission
Finally, we define the local density operators on the full particle-conserving Hilbert space using the boson_basis_1d class. Having set up the Hamiltonian, we now proceed to the time-evolution part of the problem. We begin by defining the initial random state in the Fock basis. Next we define the times which we would like to solve the Schrödinger equation for. Since the Hamiltonian is time-independent, we use the exp_op class to compute the time-evolution operator as the exponential of the Hamiltonian, see Sec. 2.3. Therefore, we consider linearly spaced time points defined by the variables start,stop, and num. To calculate the states as a function of time, we use the expm function of the block_ops class to first construct the unitary evolution operator as the matrix exponential of the timeindependent Hamiltonian 5 . We define the expm function to have almost identical arguments as that of the exp_op class, but with some major exceptions. For one, because the Hamiltonian factorizes in a block-diagonal form, the evolution over each block can be done separately (e.g. just trivially loop through the blocks sequentially). In some cases, however, e.g. for single particle Hamiltonians [see Sec. 2.2], there are a lot of small blocks, and it actually makes sense to calculate the matrix exponential in block diagonal form which is achieved by setting the optional argument block_diag=True. Another optional argument, n_jobs=int, allows the user to spawn multiple python processes which do the calculations for the different blocks simultaneously for n_jobs>1. On most systems these processes will be distributed over multiple CPUs which can speed up the calculations if there are resources for this available. This also works in conjunction with the block_diag flag where each process creates its own block diagonal matrix for the calculation. Once all the calculations for each block are completed, the results are combined and conveniently projected back to the original local Fock basis automatically. The complete code including the lines that produce Fig. 4

The Gross-Pitaevskii Equation and Nonlinear Time Evolution
This example shows how to: • simulate user-defined time-dependent nonlinear equations of motion using the evolution.evolve() routine, • use imaginary time dynamics to find the lowest energy state.
Physics Setup-The Gross-Pitaevskii wave equation (GPE) has been shown to govern the physics of weakly-interacting bosonic systems. It constitutes the starting point for studying Bose-Einstein condensates, but can also appear in non-linear optics, and represents the natural description of Hamiltonian mechanics in the wave picture. One of its interesting features is that it can exhibit chaotic classical dynamics, a physical manifestation of the presence of a cubic non-linear term.
Here, we study the time-dependent GPE on a one-dimensional lattice: where J is the hopping matrix Whenever U = 0, the system is non-interacting and the GPE reduces to the Heisenberg EOM for the bosonic field operatorψ j (t). Thus, for the purposes of using QuSpin to simulate the GPE, it is instructive to cast Eq. (12) in the following generic form where [ ψ(t)] j = ψ j (t), and • represents the element-wise multiplication The time-dependent single-particle Hamiltonian in real space is represented as an L×L matrix, H sp (t), which comprises the hopping term, and the harmonic trap. We want to initiate the time-evolution of the system at t = 0 in its lowest energy state. To this end, we can define a 'ground state' for the GPE equation, in terms of the configuration which minimises the energy of the system: One way to find the configuration ψ GS , is to solve the GPE in imaginary time (it → τ ), which induces exponential decay in all modes of the system, and singles out the lowest-energy state in the longer run. In doing so, we keep the norm of the solution fixed: Snapshots of the imaginary-time evolution can be seen in Fig. (5), top row. Once we have the initial state ψ GS , we evolve it according to the time-dependent GPE, Eq. (12), and track down the time evolution of the condensate density ρ j (t) = |ψ j (t)| 2 . Fig. 5 (bottom row) shows snapshots of the state as it evolves. Code Analysis-In the following, we demonstrate how one can code the above physics using QuSpin. As usual, we begin by loading the necessary packages: 1 from __future__ import print_function, division 2 from quspin.operators import hamiltonian # Hamiltonians and operators 3 from quspin.basis import boson_basis_1d # Hilbert space boson basis 4 from quspin.tools.evolution import evolve # nonlinear evolution 5 import numpy as np # generic math functions 6 import matplotlib.pyplot as plt # plotting library 7 from six import iteritems # loop over elements of dictionary Next, we define the model parameters. We distinguish between static parameters and dynamic parameters -those involved in the time-dependent trap widening. 9 ##### define model parameters ##### 10 L=3 0 0 # system size 11 # calculate centre of chain 12 if L%2== 0: In order to do time evolution, we code up the trap widening protocol from Eq. (12) in the function ramp. Since we want to make use of QuSpin's time-dependent operator features, the first argument must be the time t, followed by all protocol parameters. These same parameters are then explicitly stored in the variable ramp_args. With this, we are ready to construct the single-particle Hamiltonian H sp (t). The first step is to define the site-coupling lists, and the static and dynamic lists. Note that the dynamic list, which defines the harmonic potential in the single-particle Hamiltonian, contains four elements: apart from the operator string and the corresponding site-coupling list, the third and fourth elements are the time-dependent function ramp and its argument list ramp_args. We emphasize that this order of appearance is crucial. 30  To create the single-particle Hamiltonian, we choose to use the bosonic basis constructor boson_basis_1d specifying the number sector to Nb=1 boson for the entire lattice, and a local Hilbert space of sps=2 states per site (empty and filled). Then we call the hamiltonian constructor to build the single-particle matrix. We can obtain the single-particle ground state without fully diagonalising this matrix, by using the sparse diagonalisation method Hsp.eigsh(). The eigsh() routine accepts the optional flags k=1 and 'which'='SA' which restrict the underlying Lanczos routine to find the first eigenstate starting from the bottom of the spectrum, i.e. the ground state. Any initial value problem requires us to pick an initial state. In the case of imaginary evolution, this state can often be arbitrary, but needs to be in the same symmetry-sector as the true GPE ground state. Here, we choose the ground state of the single-particle Hamiltonian for an initial state, and normalise it to one particle per site. We also define the imaginary time vector tau. This array has to contain sufficiently long times so that we make sure we end up in the long imaginary time limit τ → ∞, as required by Eq. (15). Since imaginary time evolution to report a bug pls visit https://github.com/weinbe58/QuSpin/issues

SciPost Physics Submission
is not unitary, QuSpin will be normalising the vector every τ -step. Thus, one also needs to make sure these steps are small enough to avoid convergence problems with the ODE solver. 52 # define initial state to flow to GS from 53 phi 0=V[:, 0]*np.sqrt(L) # initial state normalised to 1 particle per site 54 # define imaginary time vector 55 tau=np.linspace( 0. 0,35. 0, 71) Performing imaginary time evolution is done using the evolve() method of the evolution submodule. This function accepts an initial state phi 0, initial time tau[ 0], and a time vector tau and solves any user-defined ODE, here GPE_imag_time. The parameters of the ODE are passed using the keyword argument f_params=GPE_params. To ensure the normalisation of the state at each τ -step we use the flag imag_time=True. Real-valued output can be specified by real=True. Last, we request evolve() to create a generator object using the keyword argument iterate=True. Many of the keyword arguments of evolve() are the same as in the H.evolve() method of the hamiltonian class: for instance, one can choose a specific SciPy solver and pass its arguments, or the solver's absolute and relative tolerance. We refer the interested reader to the documentation, cf. App. C. Last, we use our GPE ground state, to time-evolve it in real time according to the trap widening protocol ramp, hard-coded into the single-particle Hamiltonian. We proceed analogously -first we define the real-time GPE and the time vector. In defining the GPE function, we split the ODE into a time-independent static part and a time-dependent dynamic part. The singleparticle Hamiltonian for the former is accessed using the hamiltonian attribute Hsp.static which returns a SciPy sparse array. We can then manually add the non-linear cubic meanfield interaction term. In order to access the time-dependent part of the Hamiltonian, and evaluate it, we loop over the attribute Hsp.dynamic, which is a dictionary, whose keys fun() This example shows how to: • define symmetries for 2D systems, • use the spin_basis_general class to define custom symmetries, • use the obs_vs_time routine with custom user-defined generator to calculate the expectation value of operators as a function of time, • use the new functionality of the basis class to calculate the entanglement entropy.
Physics Setup-In Sec. 2.1, we introduced the spin-1/2 transverse-field Ising (TFI) model and showed how one can map it to an exactly-solvable quadratic fermion Hamiltonian using the Jordan-Wigner transformation. This transformation allows one to obtain an exact, closed-form analytic expression for the eigenstates and eigenenergies of the Hamiltonian (at least when the system obeys translational invariance). The fact that such a transformation exists, is intrinsically related to the model being integrable. Integrable models feature an extensive amount of local integrals of motion -conserved quantities which impose specific selection rules for the transitions between the many-body states of the system. As a result, such models thermalise in a restricted sense, i.e. in conformity with all these integrals of motion, and the long-time behaviour of the system is captured by a generalised Gibbs ensemble [26]. The system thermalises but thermalisation is constrained. On the other hand, non-integrable systems have a far less restricted phase space, and feature unconstrained thermalisation. Their long-time behaviour is, therefore, captured by the Gibbs ensemble. Since thermalising dynamics appears to be intrinsically related to the lack of integrability, if detected, one can use it as an indicator for the absence of a simple, closed form expression for the eigenstates and the eigenenergies of generic quantum many-body systems. Imagine you are given a model, and you want to determine its long-time thermalisation behaviour, from which in general you can infer information about integrability. One way to proceed is to conceive a numerical experiment as follows: if we subject the system to periodic driving at an intermediate driving frequency, this will break energy conservation. Then, in the non-integrable case, the system will keep absorbing energy from the drive until the state reaches infinite-temperature. On the contrary, in the integrable scenario unlimited energy absorption will be inhibited by the existence of additional conservation laws, and the system is likely to get stuck at some energy density. Hence, the long-time dynamics of out-of-equilibrium many-body systems can be used to extract useful information about the complexity of the underlying Hamiltonian.
Below, we take the spin-1/2 transverse field Ising model with periodic boundary conditions, and promote the lattice to be two-dimensional. There is no known simple mapping to a quadratic Hamiltonian in 2d and, therefore, the model is supposedly no longer integrable. Here, we study thermalisation through energy absorption by periodically driving the two spin systems, and demonstrate the difference in the long-time heating behaviour between the two. 1 from __future__ import print_function, division 2 from quspin.operators import hamiltonian, exp_op # operators 3 from quspin.basis import spin_basis_1d, spin_basis_general # spin basis constructor 4 from quspin.tools.measurements import obs_vs_time # calculating dynamics 5 from quspin.tools.Floquet import Floquet_t_vec # period-spaced time vector 6 import numpy as np # general math functions 7 import matplotlib.pyplot as plt # plotting library Then, we define the model parameters. Note that we choose all physical couplings to be unity, so the only parameters are the drive frequency and the drive amplitude, which is set equal to the driving frequency [the value of which is chosen arbitrarily]. The variable L denotes the linear dimension of the two spin-1/2 systems: the 1d chain has a total of L_1d sites, while the 2d lattice has linear dimensions Lx and Ly, and a total number of sites N_2d. All variables indexed by _1d and _2d in the code below refer to the 1d and 2d system, respectively. Before we move on to define the two bases, let us explain how to use the spin_basis_general class two manually define custom symmetries in 2d. QuSpin can handle any unitary symmetry operation Q of multiplicity m Q (Q m Q = 1). We label its eigenvalues exp(−2πiq/m Q ) by an integer q = {0, 1, . . . , m Q − 1}, which is used to set the symmetry block. For instance, if Q = P is parity (reflection), then q = 0, 1 will correspond to the positive and negative parity blocks. If Q = T is translation, then q will label all momentum blocks. In the present case, the symmetries are translations and parity along the x and y lattice directions, and spin inversion.
To this end, we first label all sites with integers, and store them in the variable s. One can think of the 2d system as a snake-like folded (or square-reshaped) s. With this representation each site s can be mapped to a ix and iy coordinate by the mapping: s=ix+Lx*iy. The action of any local symmetry can be programmed on this set of sites. For instance, this helps us define the action of the x and y translation operators T_x and T_y, by shifting all sites in the corresponding direction. Parity (reflection) operations P_x and P_y are also intuitive to define. Last, the spin inversion symmetry is uniquely defined through the mapping s → -(1+s) for every site s it acts on. In this way, the user can define different lattice geometries, and/or sublattice symmetries and QuSpin will do the hard job to reduce the Hilbert space to the corresponding symmetry sector. 16  Let us define the spin bases. For the 1d case, we make use of the spin_basis_1d constructor. We note in passing that one can also study higher-spin systems by using the optional argument S, which accepts a string (integer or half-integer) to specify the spin vector size, see documentation C. Requesting symmetry blocks works as usual, by using the corresponding optional arguments. For the 2d spin system, we use the spin_basis_general class. This constructor allows the user to build basis objects with many user-defined symmetries, based on their action on the lattice sites. Unlike the 1d-basis constructors, that have symmetry blocks with pre-defined variables which take integer values kblock= 0, the general basis constructors accept user-defined block variable names which take a tuple for every symmetry requested: the first entry is the symmetry transformation itself, and the second one -the integer which labels the required symmetry block, e.g. kxblock=(T_x, 0). 26  To set up the site-coupling lists for the two operators in the Hamiltonian we proceed in the usual manner. Using them, we can call the hamiltonian constructor to define the operators − ij S z i S z j and − j S x j . Here we keep the operators separate, in order to do the periodic step-drive evolution, which is why we do not need to define separate static and dynamic lists. 34  In order to do time evolution, we need to define the initial state of our system. In this case, we start from the ground state of the Hamiltonian − ij S z i S z j . We remind the reader that, since we work in a specific symmetry sector, this state may no longer be a product state. To this end, we employ the eigsh() method for sparse hermitian matrices of the hamiltonian class, where we explicitly specify that we are interested in getting a single (k=1) smallest algebraic (i.e. ground state) eigenenergy, and the corresponding eigenstate (a.k.a. the ground state), c.f. Example 0 in Ref. [14] for more details. 49 ###### calculate initial states ###### 50 # calculating bandwidth for non-driven hamiltonian 51 [E_1d_min],psi_1d = Hzz_1d.eigsh(k=1,which="SA") 52 [E_2d_min],psi_2d = Hzz_2d.eigsh(k=1,which="SA") 53 # setting up initial states 54 psi 0_1d = psi_1d.ravel() 55 psi 0_2d = psi_2d.ravel () We are now set to do study the dynamics following the periodic step-drive. Before we go into the details, we note that QuSpin contains a build-in Floquet class under the tools module which can be useful for studying this and other periodically-driven systems, see Example 2 from Ref. [14]. Here, instead, we focus on manually evolving the state. First, we define the number of periods we would like to stroboscopically evolve our system for. Stroboscopic evolution is one where all quantities are evaluated at integer multiple of the driving period. To set up a time vector, which explicitly hits all those points, we use the Floquet_t_vec class, which accepts as arguments the frequency Omega, the number of periods nT, and the number of points per period len_T. The Floquet_t_vec class creates an object which has many useful attributes, including the stroboscopic times and their indices, the period, the starting point, etc. We invite the interested reader to check out the documentation for more information C. 57 ###### time evolution ###### 58 # stroboscopic time vector 59 nT = 2 0 0 # number of periods to evolve to 60 t=Floquet_t_vec(Omega,nT,len_T=1) # t.vals=t, t.i=initial time, t.T=drive period Since the Hamiltonian is piece-wise constant, we can simulate the time evolution by exponentiating the separate terms. Note that, since we choose the driving phase (Floquet gauge) to yield a time-symmetric Hamiltonian, i.e. H(−t) = H(t), this results in evolving the system with the Hamiltonians H zz + AH x , H zz − AH x , H zz + AH x for the durations T /4, T /2, T /4, respectively (think of the phase of the drive as that of a rectilinear cosine drive sgn cos Ωt). To compute the matrix exponential of a static operator, we make use of the exp_op class, where exp(zB) =exp_op(B,a=z) for some complex number z and some operator B.
unitaries within a period to apply them on the state. The generator character of the function means that it will not execute the loops it contains when called for the first time, but rather store information about them, and return the values one by one when prompted to do so later on. This is useful since otherwise we would have to loop over all times to evolve the state first, and then once again to compute the observables. As we can we below, the generator function allows us to get away with a single loop. Finally, we are ready to compute the time-dependent quantities of interest. In order to calculate the expectation ψ(t)|H zz |ψ(t) QuSpin has a routine called obs_vs_time(). It accepts the time-dependent state psi_12_t (or its generator), the time vector t.vals to evaluate the observable at, and a dictionary, which contains all observables of interest (here Hzz_12). The output of obs_vs_time() is a dictionary which contains the results: every observable is being parsed by a unique key (string) (here "E"), under which its expectation value will appear, evaluated at the requested times. Further, if one specifies the optional argument return_state=True, the time-evolved state is also returned under the key "psi_t". 77 ###### compute expectation values of observables ###### 78 # measure Hzz as a function of time 79 Obs_1d_t = obs_vs_time(psi_1d_t,t.vals,dict(E=Hzz_1d),return_state=True) 80 Obs_2d_t = obs_vs_time(psi_2d_t,t.vals,dict(E=Hzz_2d),return_state=True) In fact, obs_vs_time() can also compute the entanglement entropy at every point of time (see documentation or Sec. 2.7). Instead, we decided to show how one can do this using the new functionality of the basis class. Each basis constructor comes with a function method ent_entropy() which evaluates the entanglement entropy of a given state, and may return the reduced density matrix upon request. To compute the entanglement, the user needs to pass the state (here Obs_12_t["psi_t"]), and a subsystem to define the partition for computing the entanglement. The method ent_entropy() can handle vectorised calculations, and will compute the entanglement of the state for each point of time. The output is stored in a dictionary, and the entanglement entropy can be accessed with the key "Sent_A". Finally, to obtain the entanglement entropy density, we also normalise the results by the size of the subsystem of interest. The bosons are subject to an on-site interaction of strength U bb , while the spin-polarised fermion-fermion interaction U ff is effective on nearest-neighbouring sites. The bosonic and fermionic sectors are coupled through an on-site interspecies density-density interaction U bf . We assume unit filling for the bosons and half-filling for the fermions. The BF mixture is initially prepared in the product state |b |f , which is a Mott state for the bosons, and a domain wall for the fermions. A low-frequency periodic drive of amplitude A and frequency Ω couples to the staggered potential in the fermions sector, and pumps energy into the system. We study the growth of the entanglement S ent (t) between the two species, see Fig. 7.
where tr b (·) and tr f (·) are the traces over the boson and fermion sectors, respectively. Code Analysis-Let us explain how to code up the Hamiltonian for the BHM now. As always, we begin by loading the necessary modules for the simulation. New here is the tensor_basis class with the help of which one can construct the basis for a tensor product Hilbert space: 1 from __future__ import print_function, division 2 from quspin.operators import hamiltonian # Hamiltonians and operators 3 from quspin.basis import tensor_basis,spinless_fermion_basis_1d,boson_basis_1d # bases 4 from quspin.tools.measurements import obs_vs_time # calculating dynamics 5 from quspin.tools.Floquet import Floquet_t_vec # period-spaced time vector 6 import numpy as np # general math functions 7 import matplotlib.pyplot as plt # plotting library First, we define the model parameters, and the drive: Next we set up the basis, introducing the tensor_basis constructor class. In its full-fledged generality, tensor_basis takes n basis objects which it uses to construct the matrix elements in the tensor product space: Here we consider the case of two Hilbert spaces, H = H 1 ⊗ H 2 , one for the bosons and one for fermions 6 . One disadvantage of the tensor basis is that it does not allow for the use of symmetries, beyond particle-conservation (magnetisation in the case of spin systems). This is because a tensor basis need not obey the symmetries of the individual bases. The new part comes in specifying the static and dynamic lists. This is where we tell QuSpin that we are dealing with two species tensored together. Notice here that the "|" character is used to separate the operators which belong to the boson (left side of tensor product) and fermion (right side of tensor product) Hilbert spaces in basis. If no operator string in present, the operator is assumed to be the identity 'I'. The site-coupling lists, on the other hand, do not require separating the two sides of the tensor product, as it is assumed that the operator string lines up with the correct site index when the '|' character is removed. For instance, the only term which couples the bosons and the fermions is a density-density interaction. It is In Sec. 2.6, we showed how one can use the method basis.ent_entropy() to compute the time evolution of the entanglement entropy, if we have the time-evolved state. Here, we show a different way of for doing the same. We use the measurements function obs_vs_time(), The output of obs_vs_time() is a dictionary. The results pertaining to the calculation of entanglement are stored under the key "Sent_time". The corresponding value itself is another dictionary, with the output of basis.ent_entropy(), in which the entropy is stored under the key "Sent_A". That second dictionary can also contain other objects, such as the reduced DM evaluated at all time steps, if these are specified in the variable Srdm_args. For more information on that, we refer the user to the documentation. The complete code including the lines that produce Fig. 7

New Horizons for QuSpin
As we demonstrated using various examples, QuSpin is currently capable of simulating a huge class of dynamical quantum systems. Nonetheless, we do not consider it a complete project, as we imagine adding further different functionalities, motivated by the needs of the users. Interested users can put the project on their Github watch list, which will notify them of any new releases.
Being an open project, we invite the interested users to fork the Github QuSpin repository, and actively contribute to the development of the project. We would be more than happy to consider well-documented functions and classes, which build on QuSpin, and include them in further official releases so they can be used by the wider audience. We will also consider patches which allow to combine QuSpin with other open-source packages for studying quantum physics. Just email us, or contact us on Github.
We would also much appreciate it, if users can properly report every single instance of a bug or malfunction in the issues section of the repository, as this is an inevitable part of maturing the code.
Some other new features not mentioned in this article include (See the Documentation C for more details): • OpenMP implimentation of SciPy's expm_multiply function useful for simulating very large systems on supercomputing clusters which have large shared memory nodes (in tools.evolution submodule).
• An efficient OpenMP parallel implementation of csr_matrix-vector multiplication again useful the dynamics of large systems (in tools.misc submodule).
Funding information MB acknowledges support from the Emergent Phenomena in Quantum Systems initiative of the Gordon and Betty Moore Foundation and ERC synergy grant UQUAM.

A Installation Guide in a Few Steps
Detailed installation instructions can be found under:

C Package Documentation
In QuSpin quantum many-body operators are represented as matrices. The computation of these matrices are done through custom code written in Cython. Cython is an optimizing static compiler which takes code written in a syntax similar to Python, and compiles it into a highly efficient C/C++ shared library. These libraries are then easily interfaced with Python, but can run orders of magnitude faster than pure Python code [61]. The matrices are stored in a sparse matrix format using the sparse matrix library of SciPy [62]. This allows QuSpin to easily interface with mature Python packages, such as NumPy, SciPy, any many others. These packages provide reliable state-of-the-art tools for scientific computation as well as support from the Python community to regularly improve and update them [63,64,65,62]. Moreover, we have included specific functionality in QuSpin which uses NumPy and SciPy to do many desired calculations common to ED studies, while making sure the user only has to call a few NumPy or SciPy functions directly.