# Example: Dynamics in the Bose-Hubbard Model¶

We want to simulate the Bose-Hubbard model, and look at some basic dynamics. So first we write our Hamiltonian that defines this model.

(1)

Equation (1) is the general expression for the Hamiltonian of the Bose-Hubbard model. For this exercise we work in the canonical ensemble (fixed particle number), and we ignore the chemical potential and write the Bose-Hubbard model as Eq. (2).

(2)

Here we have our creation and annihilation operators, , and our number operator, , as usual. Additionally, our second term is scaled by the parameter that describes site interactions, .

First we must include the necessary python libraries. The two that are needed are the MPSPyLib, which includes how python processes the MPS matrix inputs and format the user's Hamiltonian into the proper format for the Fortran compiler to run. The post-process flag is an argument in the main function and can be passed as command line argument as --PostProcess=F or --PostProcess=F. We have imported two different libraries, and use the handles mps and np to refer to functions within each library. Thus any time a function is preceded with one of these handles, (i.e. Ops.BuildBoseOperators())

 1 2 3 4 5 import MPSPyLib as mps import numpy as np from sys import version_info import matplotlib.pyplot as plt import sys 

Next we need to create which operators we want to observe in our Hamiltonian. So for the Bose-Hubbard model, we want bosonic operators, so we call Ops.BuildBoseOperators(), with a maximum of 5 particles allowed on a given lattice site. This is the term inside the parentheses. We then want build a custom operator: we are choosing to look at our interaction term, which will be defined as . The following lines correspond to creating this operator, and setting up our system.

 20 21 22 23 24  # Build operators Operators = mps.BuildBoseOperators(5) Operators['interaction'] = 0.5 * (np.dot(Operators['nbtotal'], Operators['nbtotal']) - Operators['nbtotal']) 

Then we build our Hamiltonian. We need to add two terms (or two sums) to our Hamiltonian. So we first define a variable that will act as our Hamiltonian, which in this case will be H.

 26 27 28 29  # Define Hamiltonian MPO H = mps.MPO(Operators) H.AddMPOTerm('bond', ['bdagger', 'b'], hparam='t', weight=-1.0) H.AddMPOTerm('site', 'interaction', hparam='U', weight=1.0) 

Note that we have two different types of MPO terms present in our Hamiltonian (here we have 'bond' and 'site' terms). This is because the summations for each term in the Hamiltonian are different. The bond summation considers nearest-neighbor interactions, and has different site indices, , for each operator. MPS takes each operator fed to the function (in this case the 'bdagger', 'b' expression) and assigns indices in order of the list. The other type of summation is simply over all of the sites. This type of MPO term expects only one operator (or a single expression of operators).

Next we define what observables we want to measure in our ground state. We can look at the number operator at each site

 31 32 33 34  # Ground state observables myObservables = mps.Observables(Operators) # Site terms myObservables.AddObservable('site', 'nbtotal', 'n') 

Then we want to define our correlation operators, when looks at neighboring sites. We define two observables, and . This is done by the following commands:

 35 36 37  # Correlation functions myObservables.AddObservable('corr', ['nbtotal','nbtotal'], 'nn') myObservables.AddObservable('corr', ['bdagger','b'], 'spdm') 

Next we set the convergence parameters for our simulation. This is important to set the convergence for finding our ground state to a reasonable tolerance to get an accurate ground state, but should be optimized along with simulation run times.

 39 40 41 42 43 44 45  # Convergence parameters myConv = mps.MPSConvParam(max_bond_dimension=20, max_num_sweeps=6) myConv.AddModifiedConvergenceParameters(0, ['max_bond_dimension', 'local_tol'], [50, 1E-14]) myKrylovConv = mps.KrylovConvParam(max_num_lanczos_iter=20, lanczos_tol=1E-6) 

We have modified our convergence parameters to check that our max bond dimension in our system is no larger than 50 and our local tolerance is . We set our Krylov-Convergence parameters to do a maximum of 20 Lanczos Iterations with a tolerance of for each calculation. It should be noted that computation time is largely impacted by these convergence parameters. If tolerances are set too small for a given system, computation time can be greatly increased.

In this simulation we also define a second set of observables to show that we can look at different observables for different portions of our simulation. Here we want to look at the number operator for each site as our system is evolving. So we define a new observable dictionary:

 47 48 49 50  # Dynamics time evolution observables dynObservables = mps.Observables(Operators) # Site terms dynObservables.AddObservable('site', 'nbtotal', 'n') 

The next step is to define values for our Hamiltonian we created, assign what location and name our simulation output will be run with, and define what type of simulation we want to run. We first define our parameters and in our Hamiltonian. We also define the number of lattice sites.

 52 53 54 55 56 57  # Specify constants and parameter lists U = 10.0 t = 1.0 staticsParameters = [] L = 6 tlist = [5.0, 20.0] 

Then we want to run a simulation over two different variations, where the value of U is changed at different rates. This will define our quenching function which will allow us to ramp our function value . This is done through:

 60 61 62 63 64 65  # Quench function ramping down def Ufuncdown(t, tau=tau): return 10.0 + 2.0 * (1.0 - 10.0) * t / tau # Quench function ramping back up def Ufuncup(t, tau=tau): return 1.0 + 2.0 * (10.0 - 1.0) * (t - 0.5 * tau) / tau 

We then define our quenches, which will change the value of one of our parameters, in this case .

 71 72 73 74 75  Quenches = mps.QuenchList(H) Quenches.AddQuench(['U'], 0.5 * tau, min(0.5 * tau / 100.0, 0.1), [Ufuncdown], ConvergenceParameters=myKrylovConv) Quenches.AddQuench(['U'], 0.5 * tau, min(0.5 * tau / 100.0, 0.1), [Ufuncup], ConvergenceParameters=myKrylovConv) 

We note that Dynamics.AddQuench() takes the Hamiltonian, the variable in our Hamiltonian which we want to ramp, the range of values for which we will apply this function, assigns a function which calculates the new value, and the set of convergence parameters that we want to abide by.

Then we create the dictionary that will be used for our simulation. This is done through the following lines of code:

 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93  staticsParameters.append({ 'simtype' : 'Finite', # Directories 'job_ID' : 'Bose_Hubbard_mod', 'unique_ID' : 'tau_'+str(tau), 'Write_Directory' : 'TMP/', 'Output_Directory' : 'OUTPUTS/', # System size and Hamiltonian parameters 'L' : L, 't' : t, 'U' : U, # Specification of symmetries and good quantum numbers 'Abelian_generators' : ['nbtotal'], 'Abelian_quantum_numbers' : [L], # Convergence parameters 'verbose' : 1, 'logfile' : True, 'MPSObservables' : myObservables, 'MPSConvergenceParameters' : myConv, 'Quenches' : Quenches, 'DynamicsObservables' : dynObservables 

Here we have three different sections of assignments. We have our simulation references, 'job\_ID' which is what each simulation file will be named, with a unique ending given by the value assigned to the key 'unique\_ID'. Then we assign where these files should be saved. We want our compilation files to be saved in the local folder TMP and all of our outputs to the local folder OUTPUTS.

Next we define our system. We need to call every unassigned variable in our Hamiltonian, and also specify our system size. We do this through the key assignments for 'L', our system size, 't' and 'U' the variables discussed in our Hamiltonian in Eq. (2).

Next we have the section which we can define any symmetries or quantum numbers for our system. These are set by the keys 'Abelian_generators' and 'Abelian_quantum_numbers'.

And finally we describe which convergence parameters we want to use for each section of our analysis, and which observables we want to measure for each portion of our analysis. This is done through the lines 89-92.

Next we have the primary function calls that convert this python dictionary to the fortran files that will be used to actually run the MPS algorithms. This is done by calling the dictionary we just defined, and passing the operators, the Hamiltonian, and the option to post process data for our simulation. Then we call the function tools.runMPS() to run the simulation.

  95 96 97 98 99 100 101 102  # Write Fortran-readable main files MainFiles = mps.WriteFiles(staticsParameters, Operators, H, PostProcess=PostProcess) # Run the simulations and quit if not just post processing if(not PostProcess): mps.runMPS(MainFiles, RunDir='./') 

Then to check that this simulation works correctly, we can use the following code write our time data for version of our simulation. Here we have run two simulations with and . This can be done through the following commands:

 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120  # Postprocessing and plotting # --------------------------- Outputs=mps.ReadStaticObservables(staticsParameters) DynOutputs=mps.ReadDynamicObservables(staticsParameters) ii = 0 le_list = [[], []] for t in tlist: myfile='outmod'+str(t)+'.dat' hfile=open(myfile,'w') for p in DynOutputs[ii]: hfile.write('%30.15E'%(p['time']) + '%30.15E'%(p['U']) + '%30.15E'%(p['Loschmidt_Echo'].real) + '%30.15E'%(p['Loschmidt_Echo'].imag) + '%30.15E'%(p['bond_dimension']) + '\n') le_list[ii].append(abs(p['Loschmidt_Echo'])**2) hfile.close() 

In the last part we collected already the information of the Losschmidt echo over time. The remaining lines plot the overlap with the initial state over time for the two different quench times:

 126 127 128 129 130 131 132 133 134 135 136 137 138 139  plt.rc('text', usetex=True) plt.rc('font', family='serif') plt.plot(np.linspace(0, 1, len(le_list[0])), le_list[0], 'r-', label="tq=5") plt.plot(np.linspace(0, 1, len(le_list[1])), le_list[1], 'g-', label="tq=20") plt.xlabel(r"\textbf{Time in percent} " r"$t / t_q$", fontsize=16) plt.ylabel(r"\textbf{Loschmidt Echo} " r"$| \langle \psi(t) | \psi(0) \rangle |^2$", fontsize=16) plt.legend(loc="lower left") plt.show()