A system composed of a harmonic oscillator coupled to a two-level atom is one of the quantum systems, which can be completely solved. Although this system is simple, it is never a easy work for the quantum calculations, especially when the system consists of many such simple constituent parts. In this paper, we present a programming method, by which the calculation tasks for the matrix representation of the Hamiltonian of system can be automatically fulfilled. Coupled-cavity array systems are used to demonstrate our programming method. Some quantum properties of these systems are also discussed.

## I. INTRODUCTION

Many experiments showing quantum phenomena have been conducted in different ways. Considerable effort for realizing photonic devices has been made. Various methods to study this subject have been proposed, such as single-photon switching,^{1–11} single-photon transistors^{12–15} and so on. When working with these photonic devices, we need to deal with the interaction between photons and atoms or other particles.^{16–18} For the photons to be under good conditions, i.e., a great chance for the light involved in the physical processes to interact with the atoms or other particles, researchers usually turn to optical cavities^{19,20} for help. For the sake of simplicity, when we analyze such systems in quantum mechanics, they can be abstracted as a quantum oscillator coupling to a two-level atom. This model is simple but very important and has been studied by many researchers.^{6,7,21}

Although the abstracted model is quite simple, no analytical solution has been reached. The most common way in dealing with the problem of such a model is to apply the rotating-wave approximation (RWA) method,^{22,23} which was first introduced by Jaynes and Cummings in 1963.^{24} The RWA method has been used in a large amount of works. It helps to solve such problems, yet the calculation tasks are still very heavy. In view of this, we come up with a method using the computer, together with the RWA method, to efficiently solve the quantum computing problems. In the process of dealing with the problems concerning finite dimensional quantum systems, which can be spanned by a set of discrete bases, we managed to calculate the matrix representation of Hamiltonian using a programming method (see supplementary material). In this paper, we present our methods and provide some examples.

## II. METHOD

### A. One-dimensional coupled-cavity array systems

First, we consider a coupled-cavity array system consisting of a finite number of cavities. For simplicity, each cavity is embedded with a two-level system, and only one single field is considered. For convenience, we first study a system containing only two cavities, just as illustrated in Ref. 21. For such a system, its Hamiltonian can be written as follows ($\u210f=1$):

where $\omega c$ and $\omega a$ represent the cavity and atom frequencies, respectively; *A* is the inter-cavity coupling constant and *g* is the coupling strength between atom and cavity, both being in the same cavity. $\xe2i\u2020$ ($\xe2i$) is the creation (annihilation) operator for the field in the *i*th cavity. In such a two-level system, $|gi$ and $|ei$ represent the ground and excited states, respectively, where the suffix *i* indicates which cavity they are in. Thus, the atomic lowering and raising operators in the *i*th cavity can be represented as $|gi\u2009ei|$ and $|ei\u2009gi|$, respectively.

Although the Hamiltonian of this form is clear and concise, if no state vector is given, the Hamiltonian would have only some literal meaning. In practical calculations, dealing with operators acting on vectors in some certain vector spaces is what we expect. In this case, if we restrict the total excitation number to 1, the basis states may be expressed as follows:

Thus the matrix representation of the Hamiltonian Eq. (1) in bases $|\varphi 1\u2009,|\varphi 2\u2009,|\varphi 3\u2009,|\varphi 4$ can be expressed as

The process for performing these calculations is very easy. However, if the number of cavities is increased to even greater, such as 10 cavities, and the total excitations are increased to 2, then the number of basis states would be 200. The work for writing out all the basis states along with calculating the matrix representation would be too heavy for manual calculation. For this reason, we introduce our method for the calculation work using computers, in the hope of improving work efficiency. Systems that can be expressed within finite dimensional Hilbert spaces are considered below. Two steps are involved in our method: (i) calculating the basis states, for example, $|\varphi 1\u2009,|\varphi 2\u2009,\u2026,|\varphi n$, which compose the Hilbert space; (ii) calculating the matrix representation of the Hamiltonian according to the basis states.

To construct a set of complete orthogonal basis states, a straightforward method is to initially enumerate all the possible states and then select the ones that satisfy certain conditions. For example, if we restrict the number of total excitations in Eq. (1) to 2, then the possible values for the number of photons in each cavity are 0, 1, and 2, whereas the possible states of atom are 0 and 1, where we have denoted states $|g$ and $|e$ as $|0$ and $|1$ (i.e. $|g\u21d2|0$ and $|e\u21d2|1$), respectively, for the sake of programming. By now it should be clear that numbers in different positions have different meanings. For example, if a number *n* appears in the position for photonic state, then *n* photons exist; similarly, if a number 0 (1) appears in the position for atomic state, then it means the atom is in ground (excited) state. In this way, a state such as $|\psi =|2111\u2297|2212$ is one of our candidates until the restriction that the number of total excitations be two is applied to it. As a result, this state $|\psi $ should of course be eliminated. A total of $3\xd72\xd73\xd72=36$ possible candidates exist; we can easily list them out. Furthermore, for some software, dealing with the list of our candidates as a whole is recommended to be more efficient. However, when processing the high-dimensional Hilbert space, storing all the possible candidates will consume quite a large amount of resources and would then cause the computer to crash or to behave abnormally.

To enumerate and check all the possible candidate states in a convenient way, we introduce the mixed radix number systems in our method. Following the above example, three possible values exist for the number of photons and two possible values for the energy levels of atom in one cavity. Thus, the mixed radix can be denoted as {3, 2, 3, 2} and all the possible candidate states in this base format should run from 0_{3}0_{2}0_{3}0_{2} = 0_{10} to 2_{3}1_{2}2_{3}1_{2} = 35_{10}, where the subscripts indicate the bases for the digits. For convenience, their relationship can be expressed in a general way as

where *n* is a decimal number, the sequence {*b*1, *b*2, *b*3, *b*4} is the mixed radix, and the list of numbers {*d*1, *d*2, *d*3, *d*4} is the representation of *n* in this mixed radix number system. The length of the list being 4 is determined by the number of particles involved in our system, which consists of two cavities, each containing photons and an atom. In programming, we check the list {*d*1, *d*2, *d*3, *d*4}, in accordance with the restriction that the total excitations are two and decide if this list should be used to construct our basis states for the Hilbert space. For example, as 3_{10} = 0_{3}0_{2}1_{3}1_{2} and 0 + 0 + 1 + 1 = 2, we thus use this sequence of numbers, which denotes a state $|01g1\u2297|12e2$ indicating that there is one photon, and the atom is in the excited state in the second cavity. When the number increases to 4_{10} = 0_{3}0_{2}2_{3}0_{2}, which indicates two photons in the second cavity, we also use this sequence {0, 0, 2, 0}. By contrast, the number 0_{3}0_{2}2_{3}1_{2}, which is 5_{10}, indicates two photons and one excited atom exist in the second cavity. It fails to satisfy our restriction. We should just increase the value of *n* and do the check in accordance with our restriction.

For different systems, the mixed radix bases and the restrictions used to select the basis states may be different. We should have no problem extending our method to more complicated systems, for example, systems with more cavities, more atoms, more total excitations, and more complicated energy levels for the atoms.

Using the method illustrated above, we would obtain a set of basis states making up the Hilbert space for our quantum system. Given the basis states, we can generate the matrix representation of the Hamiltonian *H*, whose element looks like

*C*_{k} reflects the operations that the Hamiltonian operator needs to perform when converting a vector $|\varphi j$ into another vector $|\varphi k$ or just the same vector $|\varphi j$. One operation of the Hamiltonian, acting on a basis state, can do the following things:

Creating one photon in one cavity, while annihilating one photon in an adjacent cavity and vice versa;

Raising the atom level from the ground state $|g$ to the excited state $|e$, while annihilating a photon within the same cavity and vice versa;

Keeping the photons and atoms in their eigenmodes unchanged.

Each of the nontrivial elements, say *H*_{i,j}, represents one (just only one) of the above three operations.

To describe our method in a general and convenient way, we consider a system containing three cavities (each contains a two-level atom) and restrict the number of excitations to be two. In our programming method, to obtain *C*_{k} as expressed in Eq. (5), we first compare the differences, regarding the photons and atom in the same cavity as a whole, between the two vectors $|\varphi i$ and $|\varphi j$. Four different situations exist.

First, more than two differences exist. For example,

where the subscripts indicate the sequence numbers of cavities. In this situation, just performing one operation of the Hamiltonian cannot change a vector, $|\varphi i$ ($|\varphi j$), into another, $|\varphi j$ ($|\varphi i$). Thus, we have

Second, 0 difference exists, implying that the two vectors $|\varphi i$ and $|\varphi j$ are the same. To keep the state unchanged, the operation of the Hamiltonian keeps the photons and atoms in their eigenmodes. For example, for the following states,

the corresponding elements for the matrix representation are

Third, one difference exists. In this case, the Hamiltonian changes both the photon and atom states in one of the cavities, while leaving the states for all the other cavities unchanged. Performing the second operation only once can fulfill this task. For example,

These two states differ from each other by the states in the third cavity. In order for the Hamiltonian to change the state $|\varphi i$ into state $|\varphi j$, it should annihilate one photon and raise the ground state of atom up to the excited state. For the change from $|\varphi j$ to $|\varphi i$, it should lower the excited state to the ground state and create a photon. We have

For states

the matrix elements are

Finally, two differences exist, i.e., for two-state vectors, after performing one operation of the Hamiltonian, two cavities differ from the previous ones. Performing the first operation once is deemed appropriate for this task, which results in a nontrivial matrix element. However, things are slightly complicated here, and we need to make further study. The different situations are presented as follows:

Both the first and second operations are involved. Although only one of the first operation is involved, however, the interaction between atom and photons may still change the states in each cavity. For example,

the transform from $|\varphi i$ to $|\varphi j$ can be processed in the following way: first hopping a photon from cavity 3 to cavity 2, and then annihilating a photon while raising the ground state of the atom to the excited state in the second cavity. Thus, as more than one operations of the Hamiltonian is processed, we will have

for these situations. Cases such as the following also exist:

with the total excitations of cavity 2 in state $|\varphi i$ being 0, whereas the total excitations of cavity 2 in state $|\varphi j$ being 2, indicating that at least two of the first operations are involved. Besides the operations of communicating photons between cavities 2 and 3, the second operation is also involved. Now that more than one operation are involved, we have

More than one of the first operation are involved. For example,

In this situation, at least two photon hopping processes are needed. Thus, we have

In cases such as

we also have

After eliminating the above nonsignificant bare state pairs, some selected ones that can transform from/to one another by performing the first operation only once remain. For example,

$|\varphi i$ can be transformed into $|\varphi j$ by hopping a photon from cavity 3 to cavity 2, whereas $|\varphi j$ can be transformed into $|\varphi i$ by hopping a photon from cavity 2 to cavity 3. Thus, for this case, we obtain

For situations such as

we have

The above is our programming method. Now for our three-cavity system, following the processes illustrated above, running from 0_{10} = 0_{3}0_{2}0_{3}0_{2}0_{3}0_{2} to 215_{10} = 2_{3}1_{2}2_{3}1_{2}2_{3}1_{2} and selecting the appropriate list of numbers to make up a complete set of basis states, for example, ${|\varphi 1\u2009,|\varphi 2,\u2026,|\varphi n}$, by which the Hilbert space can be composed, we would get the matrix representation of the Hamiltonian as follows:

In this expression, we introduce the energy detuning, which is expressed as $\Delta =\omega a\u2212\omega c$, to make it fit the width of layout without changing the eigenvectors and the order of the eigenvalues. After getting the matrix representation of the Hamiltonian, we can calculate the eigensystem of the quantum system, which is not a problem nowadays because our computers (even the plain personal laptops) are quite good at processing things such as this, and get most of the information about the system.

### B. High-dimensional coupled-cavity systems

In the above description about our method, we take the one-dimensional coupled cavities as examples. However, extending our method to two- or three-dimensional systems is easy. Things to be aware of are the adjacent conditions. We take a two-dimensional system for the detailed introduction. Figure 1 shows nine cavities: five nonempty cavities each containing a two-level atom, and four empty cavities. If any two cavities are indexed as $m,n$ and $p,q$, then the conditions of adjacent cavities can be expressed as

We consider the subspace with total excitations being 2. The mixed radix of the system, counting the cavities from left to right and from top to bottom, can be expressed as {3, 2, **3**, **1**, 3, 2, **3**, **1**, 3, 2, **3**, **1**, 3, 2, **3**, **1**, 3, 2}, where we used the normal digits to describe the nonempty cavities, whereas the bold digits are used to describe the empty cavities. Following the processes illustrated above, we can get a set of basis states composing this Hilbert space. When calculating the matrix element, these new adjacent conditions should be taken into consideration. Obtaining the matrix representation of this Hamiltonian is as easy as in the last example.

## III. TYPICAL EXAMPLES

Next, we use our method to solve some more “practical” problems. Let us consider a one-dimensional coupled resonator waveguide, as depicted in Refs. 6,7. We can see that the total excitation is 1, and some two-level atoms are embedded in the middle cavities of the cavity array, whereas the left and right end up being empty cavities. In this paper, we first study the time evolution of a certain state, say $|\psi (0)$. Once we get the system state $|\psi (t)$ at time *t*, we can give the probability distribution of the excitation as a function of time; that is to say, we can determine at any time whether the excitation exists as a photon in the *i*th cavity or as an excited atom in the *j*th cavity. However, in cases where many cavities exist, these distributions are slightly complex for us to read. To derive results which are concise and of great significance, we divide the cavity array into three parts: the left part consisting of the left *n*_{L} empty cavities; the middle parts, the middle *n*_{M} cavities (each of which is embedded with a two-level atom); and the right part, *n*_{R} empty cavities. If we can turn off the interaction between the photon and atom, that cavity would be empty for photon. By contrast, if we turn the interaction on, that cavity would be nonempty for photon. Thus, we can adjust the numbers *n*_{L}, *n*_{M}, and *n*_{R} with this method. In the following, we mainly discuss the properties of the system based on the probability distributions for finding the excitation to be in these three parts.

Let us denote the bases of the system as $|\varphi i$, $i=1,\u2026,m$. Using our method discussed above, we can easily obtain the eigensystems of the quantum systems when the numbers of empty and nonempty cavities are different. Let us denote the eigenstates of the Hamiltonian as $|Ep$, thus the completeness condition satisfies

If the state at the beginning is illustrated as

then the state at time *t* can be expressed as

where *E*_{p} is the eigenvalue for state $|Ep$. If the initial state is one of the eigenstates of the system, the state will not change over time, except for a global phase factor. Thus the probability distributions of the excitation in the three parts would be constants. Let us start with a simple system, where *n*_{L} = *n*_{R} = 1 and *n*_{M} = 2. We choose the eigen states as the initial states. The probability distributions, say *p*_{L}, *p*_{M} and *p*_{R}, as functions of time are shown in Fig. 2. Clearly, they are all constants over time.

We are not interested in the eigenstate cases. In this paper, we assume that, at the very beginning, no photon exists in any of the cavities, and the atoms in this coupled-cavity array are all in their ground states. Then, suddenly, a photon is injected into the leftmost cavity. If we number the cavities from left to right as $1,2,\u2026,n$, such an initial state can be expressed as $|1,g1\u2297|0,g2\u2297\u2026\u2297|0,gn$, where *n* = *n*_{L} + *n*_{M} + *n*_{R}, or just $|1\u2297|0\u2297\u2026\u2297|0\ufe38nL\u22121\u2297|0,g\u2297\u2026\u2297|0,g\ufe38nM\u2297|0\u2297\u2026\u2297|0\ufe38nR$. In most of the cases, such a state would not be one of the eigenstates of the system but rather be a linear combination of all the eigenstates, such as Eq. (29), with coefficients being $\u27e8Ep|\psi (0)\u27e9$. As time goes on, the coefficients $eiEpt\u27e8Ep|\psi (0)\u27e9$ vary for different eigenvalues; that is to say, the relative phases of the evolutionary state vary with time. The changes in the relative phases, which is resulted from different eigenvalues of different eigenstates, lead to the photon hop processes between the adjacent cavities. The length of the cavity array, the number of the atoms, the detuning between atom energy level and photon frequencies, and the hopping strength contribute a lot to the state evolution.

We first consider changing the number of empty cavities, and keeping other factors unchanged. The parameters are *n*_{M} = 1, $\Delta =0$, and *A* = 0.01*g*, and we assign 1, 2, 3, 4 to *n*_{L} and *n*_{R}, where *n*_{L} = *n*_{R}. The results are shown in Fig. 3. The probability distribution curves exhibit the sinusoidal/cosine oscillations over time. For the longer-length cavity array systems, it takes more time for the excitation to be transmitted to the other end of the cavity array. Thus the oscillation periods become longer as we increase the length, by adding more empty cavities, of the cavity array. Throughout the process, the excitation as a photon state exists in the left or the right empty cavities. The nonempty cavity seems inactive in the evolution process.

Next, we increase the number of atoms while keeping all the other parameters unchanged to conduct further study. We obtain the following results as shown in Figs. 4 and 5. By increasing the number of nonempty cavities, we see that when we increase the number of nonempty cavities rather than the empty cavities, the oscillation periods increase more significantly. Generally, a greater number of cavities involved in the system means more time required for the oscillation to go through a period.

In the above discussions, the excitation always exists as a photon state in the empty cavities. The length of the cavity array seems to be a trivial factor for finding the excitation to exist as an excited atom state. Thus, we turn to the detuning, that is, $\Delta $, for answers. The results are shown in Fig. 6. As the detuning increases, the probability for us to find the excitation to be in the nonempty cavities increases obviously. When $\Delta =0$, the interaction between the photon and atom is very intense. As a result, the photon is blocked from being transmitted to the other end. If we increase the detuning, the interactions within the system seem to be less intense, which leads to the transmission of the photon throughout the system. As expected, the oscillating period gets shorter as we keep on increasing the detuning. The length of the cavity array still plays an important role in the evolution processes. The situation is just as we illustrated above. If we increase the length of the cavity array, the oscillating period will get longer. This result is shown in Fig. 7.

Finally, we consider the hopping strength. The results are shown in Fig. 8. It makes the photon hopping processes between adjacent cavities more intense as we increase the hopping strength. As a result, it speeds up the time evolution process of the system.

In the above examples, our programs completed the computational tasks quite well. Now we are still working to improve them to be more useful. There have been many programming projects for solving the quantum computing problems. The most famous one among these projects must be QuTiP, the Quantum Toolbox in Python,^{25,26} which is very powerful. In dealing with quantum computing problems, the exponentially increasing dimensionality of the Hilbert space often severely limits the size of system that can be analytically calculated. In this regard, our programming method just consider some certain subpart of the whole Hilbert space. A similar approach also exists in QuTiP. In order to make the problem simple and without losing the generality, just as illustrated in the above, a subspace in each physical system containing a total of one excitation is considered. Regarding the quantum dynamics of physical systems in these smaller spaces, we obtain the analytical result for each curve. However, it should be clear that as the quantum system grows too large, where even the subspace is too big for us to obtain the analytical solution, we should turn to the numerical method, just as illustrated in Refs. 25,26, for help.

## IV. CONCLUSION

In this paper, the quantum systems consisting of harmonic oscillators and two-level atoms are considered. We first introduced our method for calculating the matrix representation of Hamiltonian in finite dimensional quantum systems using a programming method. To implement this method, two main procedures are described in detail.

For the first process, we introduced two ways to build a Hilbert space. The first way is to enumerate all possible candidates first and then select the ones that meet the requirements. It consumes a lot of memory, which is a disaster for physical systems with large dimensions. Thus, we thought of another way, that is, checking the candidates one by one and deciding whether they meet the requirements. In this way, all the chosen ones would make up a set of bases for the Hilbert space. The mixed radix number systems are introduced for efficiency.

For the second process, we introduced three important operations, describing how a Hamiltonian would change a state vector into another vector, and described in detail the construction of a matrix representation of a Hamiltonian. In order for the off-diagonal elements of the matrix to be nontrivial, the most important thing that we should emphasize is that the excitation exchanging requirement must be satisfied.

Finally, to test our method, we gave some examples about cavities and discussed some of their properties. Nonetheless, this method is not limited to this model only. This method can be applied to more complicated models.

## SUPPLEMENTARY MATERIAL

See supplementary material for more details about the partial source code of program for calculating the excitation number restricted state space, not reported in the main text. In this program (programmed by Mathematica software), the state space for a system containing five cavities, each consists of one two-level atom and some photons, are calculated.

## ACKNOWLEDGMENTS

We acknowledge the support by the National key research and development program of China (2017YFA0303700) and the National Natural Science Foundation of China (11534006).