We describe the computational model of quantum energy transfer in JCH like chains based on qubit representation of quantum states. The chain consists of optical cavities, each of which contains one two level atom interacting with field inside the cavity. Photons can jump between cavities that creates the effect of conductivity in the chain. Mathematical model of the conductivity is based on the qubit representation of quantum states in which any qubit has the fixed physical sense. This allows to include the dephasing noise, input and output intensity and analyse the contra intuitive quantum effects like dephasing assisted transport (DAT) and quantum bottleneck. The program realization of the model is based on the package Mathematica and utilizes such advantages of object oriented programming as module structure of the program. We represent the results of numerical simulation: the optimal output and input intensities, DAT and bottleneck.
Computer simulation is the important tool for the investigation of open quantum systems
(see [
It is convenient to measure the memory in terms of bits; for the quantum states it will be qubit memory. For the classical simulation with n qubits we must store 2n bits, and this value grows so fast that for the simulation on personal computers we are limited by the value n = 12, for supercomputers the values n ≈ 30 are only accessible now and up to n = 60 in the foreseeable future. The supercomputer simulation is effective if only we can select carefully very limited computational task; to find the statement of such a task we should fulfil preliminary simulation on the more flexible personal computers, for which the preparation and tuning of the program is substantially easier.
In addition, even the usage of a supercomputer for the simulation of quantum dynamics requires the peculiar algorithms substantially optimizing the standard technique of linear algebra that can be found only in the framework of one ten of qubits, where the Hamiltonian matrix is in principle visible.
In this work we propose the qubit model, in which any qubit has the fixed physical sense. For example, if we encode the state of the single two level atom in the optical cavity in three qubits, we can reserve the first qubit for its state: |0i1 — ground, and |1i2 — excited, and the rest two qubits — for the state of the field, such that |00i2,3, |01i2,3, |10i2,3, |11i2,3 encode vacuum state, one, two and three photon states correspondingly. This approach makes possible to measure exactly the memory that we reserve for quantum states and manipulate it as it is convenient for the different forms of the model.
We will describe the qubit form of Jaynes-Cummings-Hubbard model (see [
The cavity stores the photons of the fixed frequency ωc close to the atomic frequency that
creates the condition of valuable coupling between atoms and field. Fock states of the photons
can be also used as the logical states for the quantum computer; this scheme has been proposed
in the works [
We measure the effectiveness of the transport by the conductivity of excitations, that is
the time in which the excitation travels from the beginning of the chain and reaches the end.
We fix this time by the additional qubit called a sink, which interaction with the last cavity
is irreversible. We also consider the noise, which affects conductivity in the form of dephasing
that brings the interesting and practically important effect of increasing conductivity with the
growth of the noise intensity: dephasing assisted transport — DAT (see the works [
The conventional mathematical representation of such open quantum systems is Markov
master equation in the form of Kossakowski-Lindblad (see [
The influence of noise is very important for the conductivity in JCH chains — in artificial
devices as well as in living organisms (see [
The conductivity of a single excitation can be formulated in the form of quantum walks
( [
DAT effect was investigated in details in the works [
Here we use excitation conductivity as the field of application of programming complex based on qubit approach. In the next section we explain the mathematical model, then we describe the computer program, and show the result of simulation made on working stations. In the last section we briefly discuss the results and perspectives.
We give a definition, starting with the most detailed model ”excitations, photons, phonons.” We consider a linear chain consisting of identical optical cavities with two level atoms inside. Each cavity holds the photons with frequency ωc. Fock state with n photons inside i-th cavity we denote by |nifii. Inside each i-th cavity there is one two-level atom, which eigenstates: ground and excited we denote by |0iat i, |1iat i correspondingly. The difference between frequencies of these states (detuning) d = ωc − ωa is small in comparison with each of them: d ≪ ωc that allows to use rotating wave approximation for the atom-photon interaction at the large enough time frame, and Jaynes -Cummings model (JC). In addition, each atom is placed in a bath of thermal phonons that have dephasing effect on its excitation. In this paper all phonons have the same frequency ωp, for which DAT effect is maximal. In fact, this frequency is equal to the difference between eigen frequencies of the Hamiltonian of JCH model (see below) that are responsible for the transfer of excitations (see the work [12]).
The transfer of energy from the cavity to cavity occurs through the flight of photons between 1 the cavities where − h δi,j is the amplitude of the flight from the cavity j to the cavity i per unit of time, so that the amplitudes of the opposing flights are mutually conjugate, that is δi,j = δ¯j,i.
The Hamiltonian of our model is thus obtained by adding to the Jaynes-Cummings Hubbard Hamiltonian HJCH the term of exciton-phonon interaction Hint ep: H =
HJCH + Hint ep, HJCH = hωa P σi+σi + hωc P ai+ai + P(γai+σi + γ¯aiσi+) + Pi6=j (δi,j ai+aj + δ¯j, iaj+ai), i i i Hint ep = g P(bi+ + bi)σi+σi.
i (1) Here ai+, ai denote operators of creation and annihilation of a photon in i-th cavity, bi+, bi+- denote creation and annihilation of a phonon in i-th cavity, σi+, σi- denote creation and annihilation of the excited state of the corresponding atom (excitation), γ is the amplitude of photon emitting by excited atom in the unit of time. We assume that the photon jump is possible between the neighbouring cavities only where it occurs with the same amplitude so that all δi,j = 0 for |i − j| 6= 1, and for |i − j| = 1, δi,j = δ. Constant g is the intensity of exciton-phonon interactions (square root of Huang-Phys factors).
In the simplified model, which we call ”excitation-phonon” photons are ignored, and excitations are transmitted from one atom to the other; its Hamiltonian Hep has the form Hep = He + Hint ep, He = hωa X σi+σi + X(μi,j σi+σj + μ¯j,iσj+σi),
i i6=j The inflow of the energy can be viewed either as a constant occurrence of the excitation on the first node, taking the energy income from the outer bath, or simply set the initial state of the first atom as excited, regardless of photons. In the first case the bath inflows can create excitation irreversibly, or we can simply consider the initial state of a strongly excited field in the first cavity, which interacts with the first atom.
There are two ways for the simulation of dephasing as well. Either we consider it as the interaction with the explicit phonons, or as the irreversible process with only excitations and without explicit phonons.
Combining all of the above methods of studying the conductivity, we get all kinds of particular computing models. In any case, for the irreversible models master equation of KossakowskiLindblad should be applied: (2) (3) Uδt = e − i · δt H ¯h ,
1 ρ˜t+δt = δt X(LiρLi∗ − 2 (Li∗Liρ + ρLi∗Li)),
i
which contains the unitary dynamics Uδt as the particular case. Here Lindblad operators Li describe the irreversible part of the process. For example, the irreversible runoff of excitations to the sink from the last node (end) is expressed by the operator Lsink = gsink|0iend|1isinkh0|sinkh1|end, where the positive coefficient gsink expresses the intensity of runoff, dephasing in i-th node resulted from the interaction with the implicit phonons is reflected by the operator Ldeph, i = σi+σi, photon leak through the walls of the i-th cavity — by Ldet, i = ai, etc. We consider the resulting quantum effects of excitation conductivity. We determine the conductivity by the degree of filling sink in a certain time.
1). Quantum bottleneck. Reduction of conductivity with increasing intensity of the runoff. Graph of conductivity depending on the rate of the runoff is not monotonically increasing as 0 would be in the case of classical conductivity, but has a local maximum at a finite rate gsink. The nature of this counter-intuitive effect is purely quantum.
It can be explained by considering the dynamics of the population of the initial state according to the time when the runoff intensity is large. If there is no runoff, the dynamics will be similar to the graph of the cosine, which has a maximum point at the initial time. The runoff constantly destroys the excitation at the end of the chain that makes the density matrix close to the initial state in which a reduction of the excitonic population at the first atom decreases with time very weak. This process competes with the obvious classic excitonic population decrease with increasing runoff.
If the runoff is strong enough, it constantly holds the density matrix (and hence population of excitations in the first atom) in the initial state in which there is practically no reduction of excitonic population at all that results in ” freezing ” of the population, like in quantum Zeno effect, in which the effect of ” freezing ” comes from frequent measurements of the quantum state. In the case of the classical conductivity instead of cosine fall of excitonic population we would have a linear decrease, and there would be no bottleneck.
2). Dephasing assisted transport — DAT. A plot of the conductivity on the intensity
of the noise (the number of phonons of the resonant frequency, or intensity of the interaction
between excitations and phonons, expressed in coefficients of Lindblad for dephasing) has a local
maximum at a non-zero point. This means that there is some non-zero noise intensity at which
the conductivity is maximal. There are quantum processes, connected with the conductivity,
the effectiveness of which is enhanced with the increase of noise, for example, mixing at random
walks on graphs ( [
DAT as a purely quantum effect. Dephasing is direct suppression of the off-diagonal elements of the density matrix. Quantum dynamics is obtained by adding the amplitude states, the evolution of which leads to the same point of the classical space. Dephasing shifts the phase of the amplitude to a certain area. If this area would be narrower than the natural spread of the phase of states in the absence of dephasing, this influence will make the interference more constructive, thus increasing conductivity. This effect can be illustrated by the example of road traffic. If possibilities of cars are the same and allow rapid acceleration and braking, then synchronous mode of motion (phase conservation) will support high bandwidth of the road. But if the possibilities of vehicles are very different, then to increase bandwidth it is needed to the contrary, smooth mode of braking and acceleration, i.e. dephasing.
tion:
Kossakowski–Lindlbad equation in diagonal form is used for calculating non-unitary evolui¯h ∂ρ ∂t
N2−1 1 = [H, ρ] + i X γi(AiρAi∗ − 2 (Ai∗Aiρ + ρAi∗Ai)) i (4) Unitary evolution is optimized by diagonalizing the Hamiltonian.
In both cases, the evolution operator is determined by a Hamiltonian H with a set of nonnormalized matrices Li = √γiAi.
For general (non-unitary) case, the evolution is computed by the equations (3) directy.
These evolution equations, as well as the measurement implementation, are independent of the used model. Various models could be created by constructing a Hamiltonian H and a set of Lindblad matrices Li according to that model rules and semantics.
In order to allow simple operator-based constructing of Hamiltonians and Lindblad matrices
we build the models by binding one a several qubits (a qubit group) to a specific state parameter,
for example to the excitation of the first atom (see [
A general chain model was constructed. This model implements a chain of cavities, each of which has similar multi-level excitation, photon field, and phonons. It has energy limiting support and supports attaching an input source and an output sink.
• bi is the phonon state for atom i, 1 — the corresponding phonon is present, 0 — the corresponding phonon is absent; • ai is the exciton state for atom i, 1 — excited state, 0 — not excited; • pi is the photon state for atom i, 1 — there is a photon, 0 — no photon. The chain model has the following parameters: • Natoms — the number of atoms, • k — photon tunnelling rate between neighbouring cavities, • μ — photon–atom interaction strength, • g — phonon–atom interaction strength, • ωa — frequency of one excitation (atomic transition frequency), • ωp — frequency of one photon (cavity frequency), • ωg — frequency of one phonon, • in — input rate, replenishment coefficient for the first excitation a1, • out — output rate, sink coefficient for transferring excitation from aNatoms to s. Hamiltonian of the system equals to Jaynes–Cummings–Hubbard Hamiltonian: H = X ωppi+pi− + X ωaai+ai− + X ωbbi+bi− + X(kpi++1pi− + k∗pi+pi−+1) + i i i i X(μpi−ai+ + μ∗pi+ai−) + (g + g∗) X((bi− + bi+)ai+ai−).
i i (5)
The sink s is attached to the last exciton in the chain. Input and output is performed using Lindblad operators: + Lin = in ∗ p1
+ − Lout = out ∗ s pNatoms (6)
In all the follow-up graphs and results, it is taken ωa = ωp = 0.1 and ωg = 0.01.
Two dephasing models were implemented: • Unitary-based: with explicit phonons as shown on 2.
Hamiltonian of the system has the (g + g∗) P((bi− + bi+)ai+ai−) part.
i
• Lindblad-based: there are no explicit photons in the system. Hamiltonian of the system
does not have the corresponding part, Lindblad-like dephasing operators Di = pi+pi− are
used instead [
With the addition of Lindblad-like dephasing operators, evaluation equation 3 takes the following form:
1 1 ρ˜t+δt = δt X(LiρLi∗ − 2 (Li∗Liρ + ρLi∗Li)) + δt X di(Di∗DiρDi∗Di − 2 (Di∗Diρ + ρDi∗Di)), i i ρt+δt = Uδ∗tρ˜t+δtUδt. (7)
Here di are the coefficients of the corresponding dephasing operators. In this model, all of di are equal to g.
The above-described model was implemented using the Wolfram Mathematica computer algebra system, and was split into several modules, and the most important of them are QuantumSystem, QuantumExperiment, QuantumM easures, and the models modules (one module for each model). For this work the corresponding model module was called QuantumM odelJ CHqU ni. Those four modules would be referenced as System, Experiment, M easures, and M odelJ CH below.
Figure 3 shows the relationships between the modules in the system, with Experiment.nb being the main notebook file with all the parameters and output. Dotted lines represent internal dependencies, solid lines represent the data flow.
Some micro-optimization details are discarded to simplify the description.
As shown above, the model could be fully described using a Hamiltonian H and set of Lindblad-like dephasing operators with corresponding coefficients: {(Li, li)} and {(Dj, dj)}. The model is internally represented as a list of three elements: • The description of the model and its parameters. For M odelJ CH, stores the number of possible particles in the system. Not used in computation, could be of any type; • Hamiltonian H of the system — a square two-dimensional list of complex numbers; • A list of Lindblad-based operators, each of which is by itself a list of three elements: – The type of the operator, being either ”Lindblad” or ”Dephasing” — a string; – The coefficient li or di, depending on the type — a complex number. This one is optional, if the list has only two items, the coefficient is taken as being equal to 1; – The matrix Li or Di, depending on the type — a square two-dimensional list of complex numbers with the same size as H.
The M odelJ CH module constructs such an list using the supplied parameters of the model and the helper functions from the System module.
The Experiment module contains everything related to the evolution logic, the most used exported function is QuantumExperiment.
QuantumExperiment method accepts one mandatory argument system — the model representation in the form mentioned in takes the model representation as described in 4.1, and some named optional arguments (OptionsP attern): initialState, evolution, maxT , dT , plotStep, measure, criteria, return.
initialState describes the initial state of the system, it is either a state vector number from the basis (from 1 to N , where N ∗ N is the size of the Hamiltonian), a state vector of length N , or a density matrix (two-dimensional list of size N ∗ N ). There is no future difference between passing in a number of some state vector or the corresponding state vector itself (for example, 2 corresponds to {0, 1, 0, ..., 0}).
evolution is the string that selects the numeric method for evolution, one of ”unitary”, ”eigen”, ”measurement”, and ”lindblad”, or a method that returns an evolution function itself. Those methods will be described below. maxT maxT is the total time of emulation, dT is the step time (δt). In an ordinary case, dT is the maximum number emulation steps. plotStep is the number of emulation steps between the measurements. For an optimized evolution method that implements an accurate analytical solution (at the moment only the ”eigen” method for the unitary evolution), plotStep has the effect of dT multiplier, where the effective dTe is equal to dT ∗ plotStep, and the maximum maxT number of emulation steps is equal to .
dT ∗ plotStep measure is either a string ”diagonal” or ”entropy” or a measurement function, that accepts either one (current density matrix R) or two (current density matrix R and a density matrix Rprev of the previous step) arguments.
criteria accepts a function of one argument — the current density matrix R. If criteria is specified, the evolution stops on the step when the passed function returns T rue.
return is a string that describes what the method should return: ”measures” — the list of measurements results, ”plot” — either two- or three-dimensional plot of the measurements results (depends on whether the measurement returns a scalar or a vector), ”last” — the final density matrix R, ”lastM easure” — the measurement result on the final density matrix, ”time” — the total evolution time (for use with criteria), or ”time, lastM easure” — a list of two elements: the total evolution time and the measurement result on the final matrix.
The evolution, system, dT , and initialState arguments are passed to the internal QuantumEvolution method.
This method is designed to return an iteration function that performs the next step of the evolution and returns a density matrix.
If evolution is not a string, but is a function — the method returns evolution[system, dT, initialState].
For all the built-in function, the output matrix is normalized with M = (M + M ∗)/2 and M = M/T race[M ] to partially eliminate the effect of inaccurate number representation.
evolution = ”unitary” returns the most simple method Rt+δt = U · Rt · U ∗, where U = − i · δt H e ¯h . Lindblad-based part is discarded.
evolution = ”eigen” returns a method similar to ”unitary”, but does not use an iterative algorithm and instead computes everything from the start. If the initialState is not a density matrix, this method could be expressed as:
Rt = Vt · Vt∗
Vt = ER · DiagonalM atrix[EU n] · ES
Where Vt is the current state vector, EU is a vector of EUj = e− i ·¯hδt λj , {λj } is the vector of eigenvalues of H, ES is the set of eigenvectors of H, V0 is the initial state vector, ER = V0·ES−1, n is the number of step from 1.
The evolution of a density matrix has a similar, but has more computational operations. evolution = ”measurement” implements the following equations:
Rat = Rt−δt + Rt−δt ∗ C + C · Rt + X XAj · Rt−δt · XBj
j
Rt = U · Rat · U ∗ Matrices Aj , Bj , and C could be easily computed from the sets {Li, li} and {Di, di}: evolution = ”lindblad” implements the following equations:
Ali = Li
∗ Bli = dT · li · Li Adi = Di∗ · Di Bdi = dT · di · Di∗ · Di
1 C = − 2 dT (X li · Li∗ · Li + X dj · Dj∗ · Dj ) i j
St = St−δt + ANt−δt · (X XAj · Rt−δt · XBj ) · ANt∗−δt
j ANt = AN d · ANt−δt APt = AP d · APt−δt Rt = APt · St · APt∗
AN0 = AP0 = IdentityM atrix[N ], AP d = e AN d = AP d−1, S0 = R0.
−i · δt ¯h
H+C , (11) (12) space of the matrix −i · t H + C each time instead of an iterative algorithm. ¯h
−i · t H+C This results in ANt = e ¯h , APt = ANt−1, St = ANt · Rt · ANt∗.
There is also a ”lindblad−eigen” method that calculates ANt and APt using the eigenvector
This model implements the two pre-defined measures, ”diagonal” and ”entropy”, as well as a set of some helper functions.
The entropy is calculated as P λi · log(λi) where λi is the i-th eigenvalue of the supplied i density matrix ρ.
Some of the helper functions are made to select the diagonal of the partial trace (to observe the behaviour of a subsystem), to calculate the partial entropy of a subsystem or the mutual entropy (for example the mutual entropy between subsystems or the mutual entropy between the photons and the excitons in the system), to measure the distance between two density matrices (to calculate the error rate or to plot a simple graph of how close did we get to the target state over time), to calculate the distance of partial trace matrices (ignoring some qubits).
The System module also contains methods to construct density matrices from the state vector or the basis state number (used by the Experiment module), and some helper functions to ease the creation of Hamiltonians and Lindblad/Dephasing matrices in the model.
For example, QuantumHamiltonian is intended to ease the creation of models that are based on qubit representation: this method accepts the total number of qubits in the system and a list of rules, each one of which consists of: • The positions of the qubits that this rule deals with, • The initial state, • The target state, • The amplitude of this rule.
For example, two rules {{1, 2}, {0, 1}, {1, 0}, k} and {{1, 2}, {1, 0}, {0, 1}, k∗} could be used to represent the interaction k · a2+a1− + k∗ · a1+a2−, and one rule {{5}, {1}, {1}, ω} could be used + − to represent ω · a5 a5 . A similar method QuantumLindblad constructs the Lindblad part of the model — the set of matrices Li.
There are also helper methods that ease the implementation of limiting the total energy of the system (and to pass the energy limit as the parameter to the model), for example the system of 10 excitons with energy limited to 1 has only 10 basis states and Hamiltonian of size 10 ∗ 10, not 1024∗1024. Based on those helper methods, QuantumHamiltonian and QuantumLindblad implement an additional three-argument version, with the third argument being the definition of the subspace to which the model is limited. That argument has to be set in order to directly construct the Hamiltonian and Lindblad matrices for the reduced system — that is considerably faster than constructing the matrices for the whole system and then selecting a subspace.
For all numeric experiments with input enabled the total energy was not explicitly limited, except for the natural limits imposed by the model. For numeric experiments with input disabled, the initial state has one photon and the total energy is limited to 1 (because it could not raise higher than that value).
μ = 0.2, k = 0.8, time = 150 μ = 0.8, k = 0.2, time = 60
For high enough times, the absolute maximum over all possible output rates is reached when g = 0. But non-zero values of g can greatly improve the conductivity in cases when the output rate does not match with the optimal output rate. That is shown on figure 5. It could also be seen that the effect does not depend on g being equal to some exact value in some cases, and there is a broad range of possible values for g that improves the conductivity.
For very low times (or low target sink values, depending on the stop criterion) these graphs are different and could show other short-lived effects, but we are inspecting target sink values close to 1.0.
In the first part of figure 5 (μ = 0.2, k = 0.8) it could be seen that non-zero values for g improve the conductivity for the case out < outOpt. In the second part (μ = 0.8, k = 0.2) it could be seen that non-zero values for g improve the conductivity for the case out > outOpt.
The exact range of possible output rate values for which the conductivity could be improved by adding dephasing depends on the parameters of the chain.
μ = 0.2, k = 0.8, time = 150 μ = 0.8, k = 0.2, time = 60
For the unitary dephasing model, the main observed effects are the same as for the Lindbladbased dephasing model in 5.2 — non-zero values of g can improve the conductivity in cases of inoptimal values of the output rate, though the overall effect of unitary dephasing (with the supplied parameters) is milder.
Figure 6 depicts the same experiment as figure 5, but with the unitary dephasing model.
It is also worth noting that setups that were optimal to reach low target sink values (for example 0.3 or 0.4) are not always optimal to reach target sink values close to 1.
A qubit-based approach to modelling the evolution of a quantum system is described.
A description of a reusable computational model is provided, mentioning some common optimizations while still keeping the model generic — an optimized system in Hilbert space of dimension 10 with an energy limit of 1 and a system in Hilbert space of dimension 210 system with the photon pumping are defined by the exact same code, and the main part of the computational system is not bound to any specific system layout at all.
Quantum bottleneck effect was explained and reproduced over both input and output rates.
Two models of dephasing showed similar results in the long term (when the numeric experiments were run long enough for the sink to reach a value close to 1).
For both dephasing models non-zero dephasing rates did not give a positive effect on the conductivity in the case of optimal output rate, but gave a large positive effect in some cases when the output rate was not optimal and when the conductivity is capped by the quantum bottleneck effect. 12. Plenio M. B., Huelga S. F. Vibration, Quanta and Biology // Contemp. Phys. 2013. Vol. 54,
N. 4. P. 181–207.