# Creating Hamiltonian Simulation Operator in Q#

I am trying to create a unitary operator $$U = \sum^{T - 1}_{k=0}$$ $$|k\rangle$$ $$\langle k |$$ $$\otimes$$ $$e^{i A k}$$ in Q#, where A is a Hermitian matrix. For the beginning, I just want A to be a combination of 2 Pauli matrices, say $$A = X + 2Z$$. This is the code that I have, where input is of type Qubit[3] and register of Qubit[2]:

   ...
let SIZE_OF_MATRIX = 2;
let unitaryGenerator = (2 ^ SIZE_OF_MATRIX, ConstructU);
let registerLE = LittleEndian(register);

MultiplexOperationsFromGenerator(unitaryGenerator, registerLE, input);
...

function ConstructU (j : Int) : (Qubit[] => Unit is Adj + Ctl) {
let generatorSystem = GeneratorSystem(2, MapToGeneratorIndex);

let evolutionGenerator = EvolutionGenerator(PauliEvolutionSet(), generatorSystem);

let unitaryOperator = TrotterStep(evolutionGenerator, 1, - IntAsDouble(j));

return unitaryOperator;
}

// The purpose of this function is to map each part of the Hamiltonian generator
// to a generator index.
// Initially, we want to test the matrix A = X + 2Z.
function MapToGeneratorIndex (index : Int) : GeneratorIndex {
// We only have 2 terms, hence index can only be 0 or 1

if (index == 0) {
// Here we just want X
return GeneratorIndex(([1], [1.0]), [0]);
}

elif (index == 1) {
// Here we want 2Z
return GeneratorIndex(([3], [2.0]), [0]);
}

// TODO: throw an error
return GeneratorIndex( ([1000], [1000.0]), [0]);
}


Does anyone know what am I doing incorrectly? I am not getting the result that I am expecting. I know the code is messy but I am just trying to make it work for a basic 2x2 matrix first.

Thanks for the help!

• From an initial read, my thought is that the problem is that you've told MultiplexOperationsFromGenerator to expect four unitary operations. Fst(unitaryGenerator) should specify the number of operations returned by Snd(unitaryGenerator). In this case, if I understand your definition of MapToGeneratorIndex, that means Fst(unitaryGenerator) should be 2 instead of 2 ^ SIZE_OF_MATRIX. – Chris Granade May 18 at 17:13
• In the problem I am trying to solve, T (the upper limit of the sum), is defined as 2 ^ SIZE_OF_MATRIX, so I am indeed expecting 4 unitary operations. How many matrices are used in the Hamiltonian should be specified in the GeneratorSystem, which I did in the first line of the ConstructU function, and they are 2 indeed. – Martin May 18 at 17:39
• Ah, thanks for clarifying. Could you provide more detail on what result you're getting, then, to help diagnose? Thanks! – Chris Granade May 18 at 17:41
• No problem! When I print out the input register using the DumpRegister function right after I apply U, the resulting state is 0 both in the |0> and |1> components for some reason. – Martin May 19 at 13:48

Check your arguments to TrotterStep(). You are passing -IntAsDouble(j) as the trotterStepSize argument, but I believe this should just be a constant, since it indicates the simulation time for each Trotter step.
let unitaryOperator = TrotterStep(evolutionGenerator, 1, 0.1);

where I arbitrarily chose a value of 0.1 for trotterStepSize, I get some behavior that I believe is more like what you are expecting.
• But then where would you use the k from $U = \sum^{T - 1}_{k=0}$ $|k\rangle$ $\langle k |$ $\otimes$ $e^{i A k}$ ? (in the code that j which you mentioned was actually this k, sorry for the confusion). I thought the simulation time depends on that variable, since you will basically have 4 unitaries of the form e^(iAk) and k changes (when you use MultiplexOperationsFromGenerator). – Martin May 27 at 10:38
• Ok, thanks, now I understand a little better what you are trying to do. So I also see your original code doing something non-trivial if I just replace IntAsDouble(j) with something like IntAsDouble(j+1) or even IntAsDouble(j-1). It isn't immediately obvious to me why this is happening, but my guess is that there is some subtle bug in your implementation where the unitary is only actually getting applied for the j==0 case. Do you need to ensure that both registers are in superposition before applying the unitary? e.g.: ApplyToEach(H, input); ApplyToEach(H, register); – Ryan Shaffer May 27 at 23:12