QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems part I: spin chains

We present a new open-source Python package for exact diagonalization and quantum dynamics of spin(-photon) chains, called QuSpin, supporting the use of various symmetries in 1-dimension and (imaginary) time evolution for chains up to 32 sites in length. The package is well-suited to study, among others, quantum quenches at finite and infinite times, the Eigenstate Thermalisation hypothesis, many-body localisation and other dynamical phase transitions, periodically-driven (Floquet) systems, adiabatic and counter-diabatic ramps, and spin-photon interactions. Moreover, QuSpin's user-friendly interface can easily be used in combination with other Python packages which makes it amenable to a high-level customisation. We explain how to use QuSpin using four detailed examples: (i) Standard exact diagonalisation of XXZ chain (ii) adiabatic ramping of parameters in the many-body localised XXZ model, (iii) heating in the periodically-driven transverse-field Ising model in a parallel field, and (iv) quantised light-atom interactions: recovering the periodically-driven atom in the semi-classical limit of a static Hamiltonian.

• A major representative feature of QuSpin is the construction of spin Hamiltonians containing arbitrary (possibly non-local in space) many-body operators. One example is the four-spin operator O = j σ z j σ + j+1 σ − j+2 σ z j+3 + h.c.. Such multi-spin operators are often times generated by the nested commutators typically appearing in higher-order terms of perturbative expansions, such as the Schrieffer-Wolff transformation [19,20,21] and the inverse-frequency expansion [22,23]. Sometimes they appear in the study of exactly solvable reverse-engineered models.
• Another important feature is the availability to use symmetries which, if present in a given model, give rise to conservation laws leading to selection rules between the manybody states. As a result, the Hilbert space reduces to a tensor product of the Hilbert spaces corresponding to the underlying symmetry blocks. Consequently, the presence of symmetries effectively diminishes the relevant Hilbert space dimension which, in turn, allows one to study larger systems. Currently, QuSpin supports the following spin chain symmetries: -total magnetisation (particle number in the case of hard-core bosons) -parity (i.e. reflection w.r.t. the middle of the chain) -spin inversion (on the entire chain but also individually for sublattices A and B) -the joint application of parity and spin inversion (present e.g. when studying staggered or linear external potentials) -translation symmetry -all physically meaningful combinations of the above We shall see in Sec. 2, constructing Hamiltonians with given symmetries is done by specifying the desired combination of symmetry blocks.
• As of present date, ED methods represent one of the most reliable ways to safely study the dynamics of a generic quantum many-body system. In this respect, it is important to emphasise that with QuSpin the user can build arbitrary time-dependent Hamiltonians. The package contains built-in routines to calculate the real (and imaginary) time evolution of any quantum state under a user-defined time-dependent Hamiltonian based on SciPy's integration tool for ordinary differential equations [24].
• Besides spin chains, QuSpin also allows the user to couple an arbitrary interacting spin chain to a single photon mode (i.e. quantum harmonic oscillator). In this case, the total magnetisation symmetry is replaced by the combined total photon and spin number conservation. Such an example is discussed in Sec. 2.4.
To set up any Hamiltonian, we need to calculate the basis of the Hilbert space it is defined on, see line 13 below. Note that, since we work with spin operators here, it is required to pass the flag pauli=False; failure to do so will result in a Hamiltonian defined in terms of the Pauli spin matrices σ µ i = 2S µ i . One can also display the basis using the command print basis. 11 ##### set up Heisenberg Hamiltonian in an enternal z-field ##### 12 # compute spin-1/2 basis 13 basis = spin_basis_1d(L,pauli=False) The XXZ Hamiltonian H from Eq. (1) obeys certain symmetries. In particular, one can specify a magnetisation sector (a.k.a. filling) using the basis optional argument N_up=int, where int∈ [ 0, L] is any integer to specify the number of up-spins, see line 14. However, magnetisation is not the only integral of motion -the model also conserves parity, i.e. reflection w.r.t. the middle of the chain. The parity operator has eigenvalues ±1 and thus further divides the Hilbert space into two new subspaces. To restrict the Hamiltonian to either one of them, we use the basis optional argument pblock=±1. Since parity and magnetisation commute, it is also possible to request them simultaneously, see line 15. We stress that each one of the lines 13-15 is sufficient to build the basis on its own and we only show them all here for clarity.
In QuSpin, many-body operators JS µ 1 i 1 . . . S µn in are defined by a string of letters µ 1 , . . . µ n , representing the operator types, µ i ∈["+","-","x","y","z"], together with a site-coupling list [J, i 1 , . . . , i n ] which holds the coupling and the indices for the sites i that each spin operator acts at on the lattice. Setting up the spin-spin operators in the XXZ model goes as follows. First, we need to define the site-coupling lists J_zz, J_xy and h_z. To uniquely specify a twospin interaction, we need (i) the coupling, and (ii) -the labels of the sites the two operators act on. QuSpin uses Python's indexing convention meaning that the first lattice site is always i = 0, and the last one: i = L − 1. For example, for the zz-interaction, the coupling is Jzz, while the two sites are the nearest neighbours i,i+1. Hence, the list [Jzz,i,i+1] defines the bond operator J zz (0)S µ i S ν i+1 (we specify µ and ν in the next step). To define the total interaction energy J zz (0) L−2 i=0 S µ i S ν i+1 , all we need is to loop over the L − 2 bonds of the open chain 12 . In the same spirit one can define boundary or single-site operators, such as h_z. It is also possible to set up multi-spin operators, as we show in Sec. 2.3. 16 # define operators with OBC using site-coupling lists 17  The above lines of code specify the coupling but not yet which spin operators are being coupled. i.e. we have not yet fixed µ and ν. To do this, we need to create a static and/or dynamic operator list. As the name suggests, static lists define time-independent operators. Given the site-coupling list J_xy from above, it is easy to define the operator J xy /2 L−2 i=0 S + i S − i+1 by specifying the spin operator type in the same order as the site indices appear in the corresponding site-coupling list: [["+-",J_xy]]. In other words, the order "+-" corresponds directly to the site-index order "i,i+1". Similarly, one should set up the hermitian conjugate In the end, one can concatenate these operator lists to produce the static part of the Hamiltonian. If the Hamiltonian has time dependence, it is defined using dynamic lists. Since we are dealing with a static problem in this section, we set the dynamic to an empty list. In the following three sections, we show how to set up non-trivial time-dependent Hamiltonians.

dynamic=[]
Once the static and dynamic lists are set up, building up the corresponding Hamiltonian is a one-liner. In QuSpin, this is done using the hamiltonian constructor class, see line 24 below. The first required argument is the static list, while the second one -the dynamic list. These two arguments must appear in this order. Another argument is the basis, which carries the necessary information about symmetries. Yet whether a given Hamiltonian has these symmetries, depends on the operators defined in the static and dynamic lists. The hamiltonian class performs an automatic check on the Hamiltonian for hermiticity, the presence of magnetisation conservation and other requested symmetries, and raises an error if these checks fail. 23  Often times, one does not need to fully diagonalise H, but a part of the spectrum suffices. For instance, if one is interested in the many-body bandwidth of the model, it can be computed from the smallest and largest eigenvalues. This can be done efficiently using the eigsh attribute (eigenvalues of a sparse hermitian matrix), see line 32 below. The optional argument k=2 ensures that only two eigenstates are calculated. To determine which ones, the The entire code is available in Example Code 1.

Adiabatic Control of Parameters in Many-Body Localised Phases
Physics Setup-Strongly disordered many-body models have recently enjoyed a lot of attention in the theoretical condensed matter community. It has been shown that, beyond a critical disorder strength, these models undergo a dynamical phase transition from an delocalised ergodic (thermalising) phase to a many-body localised (MBL), i.e. non-conducting, non-thermalising phase, in which the system violates the Eigenstate Thermalisation hypothesis [25,26,27,28,29,30,31].
In our first QuSpin example, we show how one can study the adiabatic control of model parameters in many-body localised phases. It was recently argued that the adiabatic theorem does not apply to disordered systems [32]. On the other hand, controlling the system parameters in MBL phases is of crucial experimental [33,34,35,36,37] significance. Thus, the question as to whether there exists an adiabatic window for some, possibly intermediate, ramp speeds (as is the case for periodically-driven systems [38]), is of particular and increasing importance.
Let us consider the XXZ open chain in a disordered z-field with the time-dependent Hamiltonian where J xy is the spin-spin interaction in the xy-plane, disorder is modelled by a uniformly distributed random field h j ∈ [−h 0 , h 0 ] of strength h 0 along the z-direction, and the spin-spin interaction along the z-direction -J zz (t) -is the adiabatically-modulated (ramp) parameter.
In the following, we set J zz (0) = 1 as the energy unit. It has been demonstrated that this model exhibits a transition to an MBL phase [39].  system is in a many-body localised phase, while for h 0 = h ETH = 0.1 the system is in the ergodic (ETH) delocalised phase. We now choose the ramp protocol J zz (t) = (1/2 + vt)J zz (0) with the ramp speed v, and evolve the system with the Hamiltonian H(t) from t i = 0 to 5 t f = (2v) −1 . We choose the initial state |ψ i = |ψ(t i ) from the middle of the spectrum of H(t i ) to ensure typicality; more specifically we choose |ψ i to be that eigenstate of H(t i ) whose energy is closest to the middle of the spectrum of H(t i ), where the density of states, and thus the thermodynamic entropy, is largest.
To determine whether the system can adiabatically follow the ramp, we use two different indicators: (i) we evolve the state up to time t f and project it onto the eigenstates of H(t f ). The corresponding diagonal entropy density: in the basis {|n f } of H(t f ) at small enough ramp speeds v, is a measure of the delocalisation of the time-evolved state ψ(t f ) among the energy eigenstates of H(t f ). If, for instance, after the ramp the system still occupies a single eigenstate |ñ f , then s d = 0.
The second measure of adiabaticity we use is (ii) the entanglement entropy density of subsystem A, defined to contain the left half of the chain and |A| = L/2. We denoted the reduced density matrix of subsystem A by ρ A , and A c is the complement of A. Figure 1 shows the entropies vs. ramp speed in the MBL and ETH phases. The interesting underlying physics is, however, beyond the purpose of this paper.
Since we want to produce many realisations of the data and average over disorder, we specify the simulations parameters: n_real is the number of disorder realisations, while n_jobs is the joblib parallelisation parameter which determines how many Python processes to run simultaneously 6 . 9 ##### define simulation parameters ##### 10 n_real=2 0 # number of disorder realisations 11 n_jobs=2 # number of spawned processes used for parallelisation Next, we define the physical model parameters. The time-dependent disordered Hamiltonian consists of two parts: the time-dependent XXZ model which is disorder-free, and the disorder field whose values differ from one realisation to another. We focus on the XXZ part first. Let us code up the driving protocol J zz (t) = (1/2+vt)J zz (0). As already explained, our goal is to obtain the disorder-averaged entropies as a function of the ramp speed v. Hence, for each disorder realisation, we need to evolve the initial state many times, each corresponding to a different ramp speed. However, calculating the Hamiltonian every single time is not particularly efficient from the point of view of simulation runtime. We thus want to set up a family of Hamiltonians {v : H(t; v)} at once, and we shall employ Python's features to do so. This will require that the drive speed v is not a parameter of the function ramp, see line 29, but is declared beforehand as a global variable. Once, ramp has been defined, reassigning v dynamically induces the corresponding change in ramp without any need to re-define ramp itself. We shall comment on how this works later on in the code. 26 ##### set up Heisenberg Hamiltonian with linearly varying zz-interaction ##### 27 # define linear ramp function 28 v = 1. 0 # declare ramp speed variable 29 def ramp(t): To set up the static part of the Hamiltonian, we follow the same steps as in Sec. 2.1. Since the Hamiltonian (2) conserves the total magnetisation, the overlap betweens states of different magnetisation sectors vanishes trivially, and we can reach larger system sizes by working in a fixed magnetisation sector. A natural choice is the zero-magnetisation sector which contains the ground state. Parity is broken by the disorder field, so we leave it out.
The time-dependent part of the Hamiltonian is defined using dynamic lists. Similar to their static counterparts, one needs to define an operator string, say "zz" to declare the specific operator from a site-coupling list. Apart from the site-coupling list J_zz, however, a dynamic list also requires a time-dependent function and its arguments, see line 34 below 7 . In the linearly driven XXZ-Hamiltonian we are setting up here, the function arguments ramp_args is an empty list, see line 31 above. The careful reader might have noticed that there is a certain freedom in coding the coupling of the time-dependent term, J zz (t) = (1/2 + vt)J zz (0): here we choose to include the constant Jzz_ 0 in the zz site-coupling list and hence this factor is absent in the definition of the ramp function. Building the Hamiltonian is straightforward, as we explained in Sec. 2.1. To produce the entropies vs. ramp speed data over many disorder realisations, we define the function realization which returns a two NumPy arrays for the MBL and ETH phases respectively. Each array contains values of the values of the diagonal entropy density s d , and the values of the entanglement entropy density s ent for each velocity, as columns of the array. We now walk the reader step by step through the definition of realization. The first argument is the vector of ramp speeds, vs, required for the dynamics. The second argument is the time-dependent XXZ Hamiltonian H_XXZ to which we shall be adding a disordered zfield for each disorder realisation. The third argument is the spin basis which is required to calculate s ent . The fourth (last) argument is the realisation number, which is only used to print a message about the duration of the single realisation run.  In order to properly be able to use H XXZ (t; v) as a family of Hamiltonians in v (we shall see exactly how this works in a moment), we explicitly declare the variable v global. 50 global v # declare ramp speed v a global variable Since the problem involves disorder, we have to generate multiple disorder realisations. In this case, it is recommended to reset the pseudo-random number generator before any random numbers have been drawn. This is because the code spawns multiple python processes to do the disorder realization in parallel. Therefore, if the pseudo-random number generator is seeded before the new processes are spawned, all the parallel jobs will produce the same disorder realizations. 52 seed() # the random number needs to be seeded for each parallel process Next, we set up the full disordered time-dependent Hamiltonian of the problem The random field h j differs from one realisation to another and is, therefore, defined inside the realization function. Recall that we want to compare the localised with the delocalised regimes, corresponding to the disorder strengths h MBL and h ETH , respectively. For each lattice site i we draw a random number unscaled_fields[i] uniformly in the interval [−1, 1], and store it in the vector unscaled_fields, see line 55 below. Building the external z-field proceeds in exactly the same way as before: (i) we calculate the site-coupling list, line 57, (ii) we designate that the operator is along the z-axis by defining a static operator list, line 59, and (iii) we use the already computed spin basis to construct the operator matrix with the hamiltonian class, lines 61-62. QuSpin has the option to disable the default checks on hermiticity, magnetisation (particle number) conservation, and symmetries using the auxiliary dictionary no_checks passed straight to hamiltonian as keyword arguments. This can allow the user to define non-hermitian operators. Last, in lines 64-65, we define the MBL and ETH time-dependent Hamiltonians, corresponding to the two disorder strengths h ETH and h M BL . Let us first focus on the MBL phase. We want the initial state to be as close as possible to an infinite-temperature state within the given symmetry sector. To this end, we can first calculate the minimum and maximum energy, Emin and Emax of the spectrum of H MBL (t = 0). Then, by taking the 'centre-of-mass' we obtain a number, E_inf_temp, which represents the infinite-temperature energy up to finite-size effects, line 72. Note that the **eigsh_args is a standard Python way of reading off the arguments by name from a dictionary. 67 ### ramp in MBL phase ### 68 v=1. 0 # reset ramp speed to unity 69 # calculate the energy at infinite temperature for initial MBL Hamiltonian 70 eigsh_args={"k":2,"which":"BE","maxiter":1E4,"return_eigenvectors":False} 71 Emin,Emax=H_MBL.eigsh(time= 0. 0,**eigsh_args) 72

E_inf_temp=(Emax+Emin)/2. 0
The initial state psi_ 0 is then that eigenstate of H MBL (t = 0), whose energy is closest to E_inf_temp, using optional argument sigma=E_inf_temp. It remains to discuss the helper function _do_ramp. Its job is to evolve the initial state psi_ 0 with the hamiltonian object H and to calculate the entropies at the end of the ramp.  Once we have the state at the end of the ramp, we can obtain the entropies as follows. Calculating s ent is done using the measurements function ent_entropy which we imported in line 3. It requires the quantum state (here the pure state psi), and the basis the state is stored in 8 . Optionally, one can specify the site indices which define the subsystem retained after the partial trace using the argument chain_subsys. Note that ent_entropy returns a dictionary, in which the value of the entanglement entropy is stored under the key "Sent".
The function ent_entropy has a many further features, described in the documentation, see App. D. This concludes the definition of _do_ramp.
Back to the function realization, we have already seen how to obtain the entropies in the MBL phase. We now do the same thing in the delocalised ETH phase. The code is the same as the MBL one: We can now display how long the single iteration took 8 The basis is required since the subsystem may not share the same symmetries as the entire chain. The complete code including the lines that produce Fig. 1 is available in Example Code 2.

Heating in Periodically Driven Spin Chains
Physics Setup-As a second example, we now show how one can easily study heating in the periodically-driven transverse-field Ising model with a parallel field [40,41,42]. This model is non-integrable even without the time-dependent driving protocol. The time-periodic Hamiltonian is defined as a two-step protocol as follows: Unlike the previous example, here we consider a closed spin chain with a periodic boundary (i.e. a ring). The spin-spin interaction strength is denoted by J, the transverse field -by g, and the parallel field -by h. The period of the drive is T and, although the periodic step protocol contains infinitely many Fourier harmonics, we shall refer to Ω = 2π/T as the frequency of the drive. Since the Hamiltonian is periodic, H(t + T ) = H(t), Floquet's theorem applies and postulates that the dynamics of the system at times lT , integer multiple of the driving period (a.k.a. stroboscopic times), is governed by the time-independent Floquet Hamiltonian H F . In other words, the evolution operator is stroboscopically given by While the Floquet Hamiltonian for this system cannot be calculated analytically, a suitable approximation can be found at high drive frequencies by means of the van Vleck inversefrequency expansion [22,43]. However, this expansion is known to calculate the effective Floquet Hamiltonian H eff in a different basis than the original stroboscopic one: which requires the additional calculation of the so-called Kick operator K eff (0) to 'rotate' to the original basis.
In the inverse-frequency expansion, we expand both H eff and K eff (0) in powers of the inverse frequency. Let us label these approximate objects by the superscript (n) , suggesting that the corresponding operators are of order O(Ω −n ): Using the short-hand notation one can show that, for this problem, all odd-order terms in the van Vleck expansion vanish [see App. G of Ref. [44]] while the first few even-order ones are given by It was recently argued based on the aforementioned Floquet theorem that, in a closed periodically driven system, stroboscopic dynamics is sufficient to completely quantify heating [45], to report a bug pls visit https://github.com/weinbe58/QuSpin/issues Regimes of slow and fast heating can then be easily detected by looking at the energy density E absorbed by the system from the drive and the entanglement entropy of a subsystem. We call this subsystem A and define it to contain L/2 consecutive chain sites 11 : where the partial trace in the definition of the reduced density matrix (DM) ρ A is over the complement of A, denoted A c , and L A = L/2 denotes the length of subsystem A. Since heating can be exponentially slow at high frequencies [46,47,48,49], one might be interested in calculating also the infinite-time quantities where ρ A is the infinite-time reduced DM of subsystem A, and ρ F d is the DM of the Diagonal ensemble [50] in the exact Floquet basis {|n F : H F |n F = E F |n F }: We note in passing that in general s rdm = lim N →∞ N −1 N l=0 s ent (lT ) due to interference terms, although the two may happen to be close.
In Fig. 2 we show the time evolution of E(lT ) and s ent (lT ) as a function of the number of driving cycles l for a given drive frequency, together with their infinite-time values.
Code Analysis-Let us now discuss the QuSpin code for this problem in detail. First we load the required classes, methods and functions required for the computation: 1 from quspin.operators import hamiltonian # Hamiltonians and operators 2 from quspin.basis import spin_basis_1d # Hilbert space spin basis 3 from quspin.tools.measurements import obs_vs_time, diag_ensemble # t_dep measurements 4 from quspin.tools.Floquet import Floquet, Floquet_t_vec # Floquet Hamiltonian 5 import numpy as np # generic math functions After that, we define the model parameters:  Next, we define the basis, similar to the example in Sec. 2.2. One can convince oneself that the Hamiltonian in Eq. (5) possesses two symmetries at all times t which are, therefore, also inherited by the Floquet Hamiltonian. These are translation invariance and parity (i.e. reflection w.r.t. the centre of the chain). To incorporate them, one needs to specify the desired block for each symmetry: kblock=int selects the many-body states of total momentum 2π/L*int, while pblock=±1 sets the parity sector. For all total momenta different from 0 and π, the translation operator does not commute with parity, in which case semi-momentum states producing a real Hamiltonian are the natural choice [51]. The optional argument a=1 specifies the number of sites per unit cell 12 .
19 # compute basis in the 0-total momentum and +1-parity sector 20 basis=spin_basis_1d(L=L,a=1,kblock= 0,pblock=1) The definition of the site-coupling lists proceeds similarly to the MBL example above. It is interesting to note how the periodic boundary condition is encoded in line 25 using the to report a bug pls visit https://github.com/weinbe58/QuSpin/issues To program the full Hamiltonian H(t), we use the second line of Eq. (5). The time-independent part is defined using the static operator list. For the time-dependent part, we need to pass the function drive and its arguments drive_args, defined in lines 15-18, to all operators the drive couples to. In fact, QuSpin is smart enough to automatically sum up all operators multiplied by the same time-dependent function in any dynamic list created. Note that since we are dealing with a Hamiltonian defined by Pauli matrices and not the spin-1/2 operators, we drop the optional argument pauli for the hamiltonian class, since by default it is set to pauli=True. 26  The following lines define the approximate van Vleck Floquet Hamiltonian, cf. Eq. (8). Of particular interest is line 37 where we define the site-coupling list for the three-spin operator "zxz". Apart from the coupling J**2*g, we now need to specify the three site indices i,(i+1)%L,(i+2)%L for each of the operators "zxz", respectively. In a similar fashion, one can define any multi-spin operator. In order to rotate the state from the van Vleck to the stroboscopic (Floquet-Magnus) picture, we also have to calculate the kick operator at time t = 0. While the procedure is the same as above, note that K eff (0) has imaginary matrix elements, whence the variable dtype=np.complex128 is used (in fact this is the default dtype optional argument that the hamiltonian class assumes if one does not pass this argument explicitly). If the user tries to force define a real-valued Hamiltonian which, however, has complex matrix elements, QuSpin will raise an error. 54    Now that we have concluded the initialisation of the approximate Floquet Hamiltonian, it is time to discuss how to study the dynamics of the system. We start by defining a vector of times t, particularly suitable for the study of periodically driven systems. We initialise this time vector as an object of the Floquet_t_vec class. The arguments we need are the drive frequency Omega, the number of periods (here 1 0 0), and the number of time points per period len_T (here set to 1). Once initialised, t has many useful attributes, such as the time values t.vals, the drive period t.T, the stroboscopic times t.strobo.vals, or their indices t.strobo.inds. The Floquet_t_vec class has further useful properties, described in the documentation in App. D. To calculate the exact stroboscopic Floquet Hamiltonian H F , one can conveniently make use of the Floquet class. Currently, it supports three different ways of obtaining the Floquet Hamiltonian: (i) passing an arbitrary time-periodic hamiltonian object it will evolve each Fock state for one period to obtain the evolution operator U (T ). This calculation can be parallelised using the Python module joblib, activated by setting the optional argument n_jobs. (ii) one can pass a list of static hamiltonian objects, accompanied by a list of time steps to apply each of these Hamiltonians at. In this case, the Floquet class will make use of the matrix exponential to find U (T ). Instead, here we choose, (iii), to use a single dynamic hamiltonian object H(t), accompanied by a list of times {t i } to evaluate it at, and a list of time steps {δt i } to compute the time-ordered matrix exponential as i exp(−iH(t i )δt i ). The Finally, we can calculate the time-dependence of the energy density E(t) and the entanglement entropy density s ent (t). This is done using the measurements function obs_vs_time. If one evolves with a constant Hamiltonian (which is effectively the case for stroboscopic time evolution), QuSpin offers two different but equivalent options, that we now discuss. (i) As a first required argument of obs_vs_time one passes a tuple (psi_i,E,V) with the initial state, the spectrum, and the eigenbasis of the Hamiltonian to do the evolution with. The second argument is the time vector (here t.vals), and the third one -a dictionary with the operator one would like to measure (here the approximate energy density HF_ 02/L, see line 83 below. If the observable is time-dependent, obs_vs_time will evaluate it at the appropriate times: ψ(t)|O(t)|ψ(t) . To obtain the entanglement entropy, obs_vs_time calls the measurements function ent_entropy, whose arguments are passed using the variable Sent_args. ent_entropy requires the basis, and optionally -the subsystem chain_subsys which would otherwise be set to the first L/2 sites of the chain. To learn more about how to obtain the reduced density matrix or other features of ent_entropy, consult the documentation, App. D. The other way to calculate a time-dependent observable (ii) is more generic and works for arbitrary time-dependent Hamiltonians. It makes use of Schrödinger evolution to find the time-dependent state using the evolve method of the hamiltonian class. While we introduced evolve in Sec. 2.2, here we explain an important feature: if the optional argument iterate=True is passed, then QuSpin will not do the calculation of the state immediately; instead -it will create a generator object. This generator object will calculate the time dependent state one by one upon request. By doing so one can avoid the causal loop over the times t.vals to first find the state, and then looping once more over time to evaluate observables. The evolve method typically works for larger system sizes, than the ones that allow full ED. One can then simply pass the generator psi_t into obs_vs_time instead of the initial tuple. Last, we compute the infinite-time value of the energy density, the entropy of the infinite-time reduced density matrix, as well as the Floquet diagonal entropy. They are, in fact, closely related to the expectation values of the Diagonal ensemble of the initial state in the Floquet basis [45]. The measurements tool contains the function diag_ensemble specifically designed for this purpose. The required arguments are the system size L, the initial state psi_i, as well as the Floquet spectrum EF and states VF. The optional arguments are packed in the auxiliary dictionary DE_args, and contain the observable Obs, the diagonal entropy Sd_Renyi, and the entanglement entropy of the reduced DM Srdm_Renyi with its arguments Srdm_args. The additional label _Renyi is used since in general one can also compute the Renyi entropy with parameter α, if desired. The function diag_ensemble will automatically return the densities of the requested quantities, unless the flag densities=False is specified. It has more features which allow one to calculate the temporal and quantum fluctuations of an observable at infinite times (i.e. in the Diagonal ensemble), and return the diagonal density matrix. Moreover, it can do additional averages of all diagonal ensemble quantities over a user-specified energy distribution, which may prove useful in calculating thermal expectations at infinite times, cf. App. D. 93 ##### calculate diagonal ensemble measurements 94 DE_args = {"Obs":HF_ 02,"Sd_Renyi":True,"Srdm_Renyi":True,"Srdm_args":Sent_args} 95 DE = diag_ensemble(L,psi_i,EF,VF,**DE_args) 96 The complete code including the lines that produce Fig. 2 is available in Example Code 3.

Quantised Light-Atom Interactions in the Semi-classical Limit: Recovering the Periodically Driven Atom
Physics Setup-The last example we show deals with the quantisation of the (monochromatic) electromagnetic (EM) field. For the purpose of our study, we take a two-level atom (i.e. a singlesite spin chain) and couple it to a single photon mode (i.e. a quantum harmonic oscillator). The Hamiltonian reads where the operator a † creates a photon in the mode, and the atom is modelled by a twolevel system described by the Pauli spin operators σ x,y,z . The photon frequency is Ω, N ph is the average number of photons in the mode, A -the coupling between the EM field E = N −1 ph a † + a and the dipole operator σ x , and ∆ measures the energy difference between the two atomic states.
An interesting question to ask is under what conditions the atom can be described 16 by the time-periodic semi-classical Hamiltonian: Curiously, despite its simple form, one cannot solve in a closed form for the dynamics generated by the semi-classical Hamiltonian H sc (t).
To address the above question, we prepare the system such that the atom is in its ground state, while we put the photon mode in a coherent state with mean number of photons N ph , as required to by the semi-classical regime [52]: We then calculate the exact dynamics generated by the spin-photon Hamiltonian H, measure the Pauli spin matrix σ z which represents the energy of the atom, σ x -the 'dipole' operator, and the photon number n = a † a: and compare these to the semi-classical expectation values Figure 3 a shows a comparison between the quantum and the semi-classical time evolution of all observables O as defined above. As expected, we find a reasonable agreement with the deviation at longer times increasing, depending on the number of photons used, and the drive frequency.
Code Analysis-We used the following compact QuSpin code to produce these results. First we load the required classes, methods and functions to do the calculation: 1 from quspin.basis import spin_basis_1d,photon_basis # Hilbert space bases 2 from quspin.operators import hamiltonian # Hamiltonian and observables 3 from quspin.tools.measurements import obs_vs_time # t_dep measurements 4 from quspin.tools.Floquet import Floquet,Floquet_t_vec # Floquet Hamiltonian 5 from quspin.basis.photon import coherent_state # HO coherent state 6 import numpy as np # generic math functions Next, we define the model parameters as follows: 8 ##### define model parameters ##### 9 Nph_tot=6 0 # maximum photon occupation 10 Nph=Nph_tot/2 # mean number of photons in initial coherent state 11 Omega=3.5 # drive frequency 12 A= 0.8 # spin-photon coupling strength (drive amplitude) 16 Strictly speaking the Hamiltonian Hsc(t) describes the spin dynamics in the rotating frame of the photon, defined by a → ae −iΩt ; however, all three observables of interest: a † a, and σ y,z are invariant under this transformation. To set up the spin-photon Hamiltonian, we first build the site-coupling lists. The ph_energy list does not require the specification of a lattice site index, since the latter is not defined for the photon sector. The at_energy list, on the other hand, requires the input of the lattice site for the σ z operator: since we consider a single two-level system or, equivalently -a singlesite chain, this index is 0. The spin-photon coupling lists absorb and emit also require the site index which refers to the corresponding Pauli matrices: in this model -0 again due to dimensional constraints. 16  Next, we define the initial state as a product state, see Eq. (15). Notice that in the QuSpin spin_basis_1d basis convention the state | ↓ = (1, 0) t . This is because the spin basis states are coded using their bit representations and the state of all spins pointing down is assigned the integer 0. To define the oscillator (a.k.a. photon) coherent state with mean photon number N ph , we use the function coherent_state: its first argument is the eigenvalue of the annihilation operator a, while the second argument is the total number of oscillator states 17 . The next step is to define a vector of stroboscopic times, using the class Floquet_t_vec. Unlike in Sec. 2.3, here we are also interested in the non-stroboscopic times in between the perfect periods lT . Thus, we omit the optional argument len_T making use of the default value set to len_T=1 0 0, meaning that there are now 100 time points within each period. We now time evolve the initial state both in the atom-photon, and the semi-classical cases using the hamiltonian class method evolve, as before. Once again, we define the solution psi_t as a generator expression using the optional argument iterate=True. 53 # evolve atom-photon state with Hamiltonian H 54 psi_t=H.evolve(psi_i,t.i,t.vals,iterate=True,rtol=1E-9,atol=1E-9) Finally, we calculate the time-dependent expectation values using the measurements tool function obs_vs_time. Its arguments are the time-dependent state psi_t, the vector of times t.vals, and a dictionary of all observables of interest, and were discussed in Sec. 2.3. obs_vs_time returns a dictionary with all time-dependent expectations stored under the same keys they were passed. They can be accessed as shown in lines 75 and 78, respectively. 73 # in atom-photon Hilbert space 74 Obs_t = obs_vs_time(psi_t,t.vals,{"n":n,"sz":sz,"sy":sy}) 75 O_n, O_sz, O_sy = Obs_t["n"], Obs_t["sz"], Obs_t["sy"] 76 # in the semi-classical Hilbert space 77 Obs_sc_t = obs_vs_time(psi_sc_t,t.vals,{"sz_sc":sz_sc,"sy_sc":sy_sc}) 78 O_sz_sc, O_sy_sc = Obs_sc_t["sz_sc"], Obs_sc_t ["sy_sc"] The complete code including the lines that produce Fig. 3 is available in Example Code 4.

Future Perspectives for QuSpin
We have demonstrated that the QuSpin functionality allows the user to do many different kinds of ED calculations. In one spatial dimension, one also has the option of using a wide range of available symmetries. The reader might have noticed that, provided the study does not require the use of symmetries (e.g. a fully disordered 2D model), one can specify the site-coupling lists such as to build higher-dimensional Hamiltonians, using the spin_basis_1d class. Setting up higher-dimensional Hamiltonians with symmetries is possible in limited cases, too, when they can be uniquely mapped to one-dimensional systems.
In addition to the features we have discussed in this paper, there are many other functions defined inQuSpin which are all listed in the Documentation (Appendix D). Some of the more interesting ones include the tensor_basis class which constructs a new basis object out of two other basis objects, thus implementing the tensor product. This can be employed, e.g., to study interacting ladders with hard-core bosons. Another class which is useful for state-of-theart calculations is the HamiltonianOperator class. It does the matrix vector product without actually storing the matrix elements which significantly reduces the amount of memory needed to do this operation. This is particularly suited for diagonalising very large spin chains using eigsh, as it only requires on the order of a hundred calls of the matrix vector product to solve for a few eigenvalues and eigenvectors (for a specific example, see the documentation in App. D). A recent addition to the tools module is the block_tools module. This module contains a class which projects an initial state in the full Hilbert space to a set of user provided symmetry sectors and then evolves each block in parallel (possibly over many CPU core's) with a single function call. This is useful in cases where the intial state may not obey the given symmetries of the Hamiltonian used to evolve the state, for example when calculating nonequal time correlation functions. Finally the hamiltonian class is not just limited to matrices generated in our code from the operator strings. In general, this class also takes arbitrary matrices as inputs for both static and dynamic operators; therefore, one can use all of the packages functionality for any user-defined matrix.
We have set up the code to make it easily generalisable to different types of systems. We are currently working towards adding the one-dimensional symmetries for spinless and spinful fermions, but we are hoping to eventually add higher spins and bosons, too. Farther into the future we may implement a number of two dimensional lattices as well as their symmetries. We also welcome anyone who is interested in contributing to this project to reach out to the authors with any questions they may have about the package organisation. All modifications can be proposed through the pull request system on github.com.
We would much appreciate it if the users could report bugs using the issues forum in the QuSpin online repository. use miniconda.

A.1 Mac OS X/Linux
To install Anaconda/miniconda all one has to do is execute the installation script with administrative privilege. To do this, open up the terminal and go to the folder containing the downloaded installation file and execute the following command: $ sudo bash <installation_file> You will be prompted to enter your password. Follow the prompts of the installation. We recommend that you allow the installer to prepend the installation directory to your PATH variable which will make sure this installation of Python will be called when executing a Python script in the terminal. If this is not done then you will have to do this manually in your bash profile file: $ export PATH="path_to/anaconda/bin:$PATH" Installing via Anaconda.-Once you have Anaconda/miniconda installed, all you have to do to install QuSpin is to execute the following command into the terminal: $ conda install -c weinbe58 quspin If asked to install new packages just say 'yes'. To keep the code up-to-date, just run this command regularly.
Installing Manually.-Installing the package manually is not recommended unless the above method failed. Note that you must have the Python packages NumPy, SciPy, and Joblib installed before installing QuSpin. Once all the prerequisite packages are installed, one can download the source code from github and then extract the code to whichever directory one desires. Open the terminal and go to the top level directory of the source code and execute: $ python setup.py install --record install_file.txt This will compile the source code and copy it to the installation directory of Python recording the installation location to install_file.txt. To update the code, you must first completely remove the current version installed and then install the new code. The install_file.txt can be used to remove the package by running: $ cat install_file.txt | xargs rm -rf.

A.2 Windows
To install Anaconda/miniconda on Windows, download the installer and execute it to install the program. Once Anaconda/miniconda is installed open the conda terminal and do one of the following to install the package: Installing via Anaconda.-Once you have Anaconda/miniconda installed all you have to do to install QuSpin is to execute the following command into the terminal: > conda install -c weinbe58 quspin If asked to install new packages just say 'yes'. To update the code just run this command regularly.
Installing Manually.-Installing the package manually is not recommended unless the above method failed. Note that you must have NumPy, SciPy, and Joblib installed before installing QuSpin. Once all the prerequisite packages are installed, one can download the source code from github and then extract the code to whichever directory one desires. Open the terminal and go to the top level directory of the source code and then execute: > python setup.py install --record install_file.txt This will compile the source code and copy it to the installation directory of Python and record the installation location to install_file.txt. To update the code you must first completely remove the current version installed and then install the new code.

B Basic Use of Command Line to Run Python
In this appendix we will review how to use the command line for Windows and OS X/Linux to navigate your computer's folders/directories and run the Python scripts.

B.1 Mac OS X/Linux
Some basic commands: • change directory: $ cd < path_to_directory > • list files in current directory: $ ls list files in another directory: $ ls < path_to_directory > • make new directory: $ mkdir <path>/< directory_name > • copy file: $ cp < path >/< file_name > < new_path >/< new_file_name > • move file or change file name: $ mv < path >/< file_name > < new_path >/< new_file_name > • remove file: $ rm < path_to_file >/< file_name > Unix also has an auto complete feature if one hits the TAB key. It will complete a word or stop when it matches more than one file/folder name. The current directory is denoted by "." and the directory above is "..". Now, to execute a Python script all one has to do is open your terminal and navigate to the directory which contains the python script. To execute the script just use the following command: Windows also has a auto complete feature using the TAB key but instead of stopping when there multiple files/folders with the same name, it will complete it with the first file alphabetically. The current directory is denoted by "." and the directory above is "..".

B.3 Execute Python Script (any operating system)
To execute a Python script all one has to do is open up a terminal and navigate to the directory which contains the Python script. Python can be recognised by the extension .py. To execute the script just use the following command:

D 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 [53]. The matrices are stored in a sparse matrix format using the sparse matrix library of SciPy [24]. 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 [54,55,56,24]. 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. The complete up-to-date documentation for the package is available online under: https://github.com/weinbe58/QuSpin/#quspin