Driven nonlinear dynamics of two coupled exchange-only qubits

Inspired by creation of a fast exchange-only qubit (Medford et al., Phys. Rev. Lett., 111, 050501 (2013)), we develop a theory describing the nonlinear dynamics of two such qubits that are capacitively coupled, when one of them is driven resonantly at a frequency equal to its level splitting. We include conditions of strong driving, where the Rabi frequency is a significant fraction of the level splitting, and we consider situations where the splitting for the second qubit may be the same or different than the first. We demonstrate that coupling between qubits can be detected by reading the response of the second qubit, even when the coupling between them is only of about $1\%$ of their level splittings, and calculate entanglement between qubits. Patterns of nonlinear dynamics of coupled qubits and their entanglement are strongly dependent on the geometry of the system, and the specific mechanism of inter-qubit coupling deeply influences dynamics of both qubits. In particular, we describe the development of irregular dynamics in a two-qubit system, explore approaches for inhibiting it, and demonstrate existence of an optimal range of coupling strength maintaining stability during the operational time.


I. INTRODUCTION
Electrically gated quantum dots provide a promising platform for realizing qubits serving as building blocks of a quantum computer [1][2][3][4][5]. Two-level systems created from one or a few electrons in quantum dots have been experimentally isolated, initialized, and coherently manipulated. Such qubits were fabricated from single [6][7][8][9][10], double [11][12][13][14], and triple quantum dot [15,16] structures. Reduction of decoherence is achieved by applying dynamical decoupling techniques [17][18][19]. All techniques based on a single quantum dot or double quantum dots require, for performing two-axis rotations of the electron spin (or a pseudospin) on the Bloch sphere, either high-frequency magnetic fields or magnetic-field gradients (from micromagnets or dynamical nuclear polarization), or spin-orbit coupling. Employing triple-dot qubits allows one to perform two-axis rotations by using the Heisenberg exchange only and completely by electrical means, which allows one to achieve fast performance by applying voltages to the gates.
Such an approach based on coded qubits with the total electron spin S ¼ 1=2 and its projection S z ¼ 1=2 was proposed by DiVincenzo et al. [20] and realized experimentally by Laird et al. [15]. More recently, two-axis rotations in such exchange-only qubits and their tomographic description were achieved [21,22]. A specific property of these qubits, requiring nonstandard techniques for their operation, stems from the fact that their natural rotation axes are not mutually perpendicular but intersect at an angle 2π=3 [see Fig. 1(a)]. This is also the case for the three-electron double-dot qubits of Ref. [13]. For single-qubit operations, this problem was resolved in Ref. [24] by tuning the qubit into the resonant exchange regime, which is defined in Sec. II below. Operations performed on such a resonant exchange (RX) qubit are very fast, at a scale of a few nanoseconds, with the Rabi nutation frequency comparable to the qubit-level splitting. With the relaxation times T 1 ∼ 40 μs and T 2 ∼ 1 μs (even without applying CPMG pulse sequences) [24], this qubit is a highly promising candidate for developing few-qubit systems. We take advantage of this separation of time scales between the operational and relaxation times and concentrate on the coupled dynamics of two qubits. For a reasonable magnitude of interqubit coupling, the duration of two-qubit operations is within the window of coherent operation.
Similarly to single-and double-dot qubits [9,12,25], the next step is achieving two-qubit entanglement, and the first proposals for capacitively [26] and exchange [27] coupled double RX qubits have already been made. Establishing entanglement between two qubits is a demanding task, and in this paper, we propose and pursue in depth a protocol based on resonant driving of one of the qubits (qubit A) at its level splitting and measuring the response of the second qubit (qubit B) to the electrical signal produced by qubit A. Therefore, qubit B serves as a detector of coupling between the two qubits. In the discussion below, we examine several models of capacitive coupling between qubits. This coupling, in the absence of external mechanisms, keeps both qubits inside their logical spaces.
First, our analysis shows that not only the magnitude of the signal induced in qubit B but also a gross pattern of its Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published articles title, journal citation, and DOI. dynamics depend critically on the specific geometry of interqubit coupling and its strength. There is an optimal range of the coupling strength because at stronger coupling (albeit rather modest) dynamics of the two-qubit system might become irregular. Second, the backaction of qubit B onto qubit A, an essential part of entanglement, profoundly influences its own dynamics. Third, the crossed-axis geometry of both qubits, while not critical for the single-qubit operation, significantly affects the coupling between qubits and two-qubit dynamics. In particular, a specific choice of the mechanism of interqubit coupling, Figs. 1(b)-1(d), is critical for the spectrum of frequencies dominating the dynamics of coupled qubits and avoiding early switching to the irregular regime. While our analysis was performed for capacitive coupling, it has implications also for exchange coupling between qubits.
In what follows, we apply our protocol to three different schemes of capacitive coupling between two exchangeonly qubits and demonstrate spectacular differences in their dynamics.

II. SINGLE EXCHANGE-ONLY QUBIT
An exchange-only qubit is a three-electron triple quantum dot with Zeeman-split-energy levels operated in the regime where the total electronic spin S ¼ 1=2 and its projection onto the magnetic field is S z ¼ 1=2 [15,20]. Its resonant exchange modification is operated in the parameter range where the two-dimensional qubit space is well separated from the other states of the system, including all S ¼ 3=2 states and S ¼ 1=2, S z ¼ −1=2 states [24].
In what follows, the three dots of each of the qubits are designated as L, M, and R, for the left, middle, and right dots. As outlined in more detail in Appendix A, basic properties of the qubit can be described in the framework of a Hubbard model whose basis includes two states, (201) and (102), with inhomogeneous charge distributions, and two "neutral" (111) states; the numbers show populations of ðL; M; RÞ dots. An electron pair in a doubly occupied dot is in a singlet state since triplet states are essentially higher in energy. The effective Schroedinger equation for this Hubbard model is for a four-spinor ðv L ; v 0 ; v 1 ; v R Þ T defined by Eq. (A5) of Appendix A; here, T stands for transpose.
Exchange-only qubits operate predominantly in the (111) region, with (201) and (102) serving mostly for initialization and projective measurements. However, an admixture of these states is critical for operation of the qubit and the interqubit coupling. In the center of the (111) region, the ground state of the qubit is j0i ≈ 1 ffiffi where ðσ z ; σ x Þ are Pauli matrices of the qubit pseudospin (for brevity, we use the term spin in what follows) acting in the qubit space. The exchange integrals are where the tunneling matrix elements t L and t R connect the (111) sector to (201) and (102) states, respectively, U > 0 is a single-site repulsion energy, and V and −V are potentials applied to the left and right dots. The V terms in the denominators reflect the effect of the j201i and j102i components of the spinor. Equations (1) and (2) are applicable when t L =U, t R =U ≪ 1; see Appendix A for details. In our discussion, we set ℏ ¼ 1, and we refer to J z as the qubit-level splitting. The resonant-exchange qubit operates in the regime t L ¼ t R near the point V ¼ 0, where J x ¼ 0 and J z is stationary, dJ z =dV ¼ 0. For simplicity, we disregard the dependence of t L and t R on V. As seen from Eq. (1), matrices ðσ x ; σ z Þ perform single-qubit rotations about the axes orthogonal in the qubit space.

III. INTERQUBIT COUPLING
Assuming the qubits are operated in a regime when tunneling between qubits A and B is negligible, the qubit-qubit interaction is capacitive through the electrostatic charges. The dots in qubits A and B are arranged in such a way that exchange coupling inside each of the qubits exists only between adjacent dots, but we consider three geometries, as shown in Figs. 1(a)-1(c). They differ in the capacitive coupling between qubits, and in all cases, we consider only the capacitive coupling between two adjacent dots belonging to different qubits. This provides us with three models that show rather different two-qubit dynamics, as shown in Sec. IV below. In each of the geometries, we use symbols L, M, and R to numerate dots inside qubits.
We begin with a linear array of Fig. 1(b), where the rightmost dot of A is proximate to the leftmost dot of B. The excess electron densities δn α at the edges of these qubits are equal to jv A R j 2 and jv B L j 2 and can be conveniently expressed in terms of the amplitudes (v 0 , v 1 ) of the corresponding qubits using the first and fourth rows of Eq. (A6). In terms of the unit vectorsê L andê R of Fig. 1 1=2Þ by their Cartesian coordinates in the xz plane, the excess electron densities on the edge dots are Hence, charges on the outermost dots are proportional to squares of projections of spin wave functions onto the tilted axes. Even though the single-qubit gate operations are immune to the nonorthogonality of ðê L ;ê R Þ axes, the electric charges, and therefore the interaction between qubits, are affected by it. The qubit-qubit interaction expressed in terms of the wave-function amplitudes of qubits A and B when it is dominated by excess charges on the right dot of A and the left dot of B [see Fig. 1 Here, κ is the dielectric constant, and R is the distance between neighboring end dots of qubits A and B. Having in mind a complicated geometry of the system and patterns of screening, R can also be considered as an effective coupling constant that should be measured experimentally.
Equation (4) includes both renormalization corrections to single-qubit Hamiltonians and interqubit correlations. Omitting the former contributions as described in Appendix B, we arrive at the interaction Hamiltonian in terms of the qubit-spin vector matrices We note thatĤ int includes σ x and σ z matrices of both qubits.
In the geometry of Fig. 1 There is one more geometry [see Fig. 1 where two qubits are coupled through their M dots. Its effective Hamiltonian isĤ The implementation of the CPhase gate using the Hamiltonian in Eq. (7) can be achieved by following standard protocols [27]. Three models of Eqs. (5-7) are based on the assumption of capacitive coupling only between the two closest dots of A and B qubits. The experimental fact that sensor dots respond predominantly to charges on the qubit dots closest to them can serve as a justification for using these models. From the theoretical perspective, this regularity can be ascribed to screening of the long-range Coulomb interaction by metallic gates. We expect that the drastic difference in spin dynamics of these models unveiled in Sec. IV below will serve to choose optimal geometries of two-qubit systems.

IV. DYNAMICS AND ENTANGLEMENT
We consider a protocol when both qubits are tuned to with an additional oscillating term J A x =2 ¼ ϵ cosðωt þ φÞ produced by an ac voltage applied between the outermost dots of qubit A beginning at time t ¼ 0. Qubit B responds only to a capacitive signal produced by qubit A, and therefore, dynamics of qubit B is a direct indication of the coupling between A and B.
To get an outlook of the signal produced by A, we choose ω ¼ J A z . If the ac perturbation is small, we can solve the problem analytically in the rotating wave approxi- . Subsequently, we transform back to the laboratory frame, and all the following equations are presented in this frame.
In the zeroth order in J AB , the solution for qubit A starting at the north pole of the Bloch sphere at the instant : The components of the spin S A of qubit A on its Bloch sphere are equal to the mean values of the components of the σ A matrix, S A ¼ hσ A i. Projections of the spin along the ðê L ;ê R Þ axes give the electron densities entering the interaction Hamiltonian according to Eqs. (5-7), Calculation based on Eq. (8) results in S A y ðtÞ ¼ sin ϵt cos ðωt þ φÞ; Now, we evaluate the response of qubit B due to the coupling J AB at first order. Comparison of Eq. (12) with Eqs. (5-7) demonstrates a critical effect of the interqubit coupling mechanism onto the electric signal driving qubit B. In the geometry of Fig. 1(b), the spectrum of the signal includes three frequencies: two Raman frequencies ω AE ϵ coming from Eqs. (10,11) and the Rabi frequency ϵ coming from Eq. (12). In the geometries of Figs. 1(c) and 1(d) only the Rabi frequency is present. The dynamics of qubit A at frequencies that differ from the driving frequency ω stems from the fact that qubits are anharmonic systems, and it was already observed in double-dot qubits [29,30].
It is a special property of the geometry of Fig. 1(d) that qubits are coupled through their M dots, and therefore, the ac perturbation experienced by qubit B is symmetric rather than antisymmetric in its L and R dots. When qubit B is tuned to the point J B x ¼ 0, its effective Hamiltonian iŝ where the time-dependent term is a perturbation exerted by qubit A in the RWA (with ϵ B ∼ J AB ), and for simplicity, we have chosen the phase φ B ¼ 0. BecauseĤ B q commutes with σ B z , the z projection of the qubit spin S B is conserved, S B z ¼ const. Solving the Schroedinger equation with the initial condition S B x ð0Þ ¼ 1, S B y ¼ 0, and applying Eq. (9), we arrive at Therefore, S B precesses about the z axis with a phase modulated by the interqubit coupling J AB .
In the experiments of Ref. [24], the Rabi frequency ϵ was large and comparable to the qubit-level splitting J A z . Under such conditions, the accuracy of RWA is reduced, and Eqs. (10)(11)(12) provide only a qualitative outlook onto the dynamics of qubit A. Therefore, in our simulations presented below, we solve for the dynamics of qubit A exactly. These simulations are also necessary to account for the backaction of qubit B onto qubit A, which provides the entanglement mechanism.
The magnitude of J AB is highly sensitive to the ratio t=U; we assume t L ¼ t R ¼ t. If to estimate t=U from J A z ≈ 2t 2 =U with t ≈ 17 μev [24], then t=U ≈ 4 × 10 −2 , with the qubit-level splitting ≈ 1.4 μev. J AB estimated with κ ≈ 10 and R ≈ 200 nm is about J AB ≈ 10 −3 J A z . However, different estimates result in larger values of t=U [31], and our simulations were performed for J AB =J A z of about 10 −2 , which seems realistic. In practice, however, the coupling J AB should be obtained from experiment. An important purpose of the protocols proposed in this paper is to enable a measurement of this coupling.
The RX qubit Hamiltonian for quantum dots in GaAs includes, besides the dominant Heisenberg exchange, two additional contributions from spin-orbit and hyperfine interactions [26,[32][33][34] as in the experiments of Ref. [24], and a similar inequality for J AB , we disregard dephasing and spin relaxation produced by these perturbations, which is a subject of future work. Therefore, our results describe coherent dynamics of two entangled qubits. Under these conditions, a two-qubit system that starts at t ¼ 0 in a product state evolves as a pure two-qubit state jψ AB i. Entanglement between qubits manifests itself through the Schmidt decomposition of jψ AB i in terms of single-qubit states [3], with c i ≥ 0, P i c 2 i ¼ 1. Entanglement can be quantified in terms of reduced density matrices ρ A ¼ Tr B ρ AB and ρ B ¼ Tr A ρ AB , where ρ AB ¼ jψ AB ihψ AB j, and partial traces are taken over the states of qubits B and A, respectively. The von Neumann entropy associated with the reduced density matrices is a measure of entanglement in the system [35,36], The entanglement E takes its maximum value E ¼ 1 Below, we present simulation results for the geometries of Figs. 1(b) -1(d). We perform our simulation in the qubit space where the interaction Hamitonian is given in terms of Pauli matrices. The total Hamiltonian governing the qubit evolution isĤ q þĤ int . The ode15s function on MATLAB was used to evaluate the results. It is a variable-order solver that employs the Gear's method for solving differential equations and is suited for solving stiff problems.
Dynamics of spin S B reflects the strength of the coupling between qubits and intricacies of their joint dynamics, while E reflects the degree of quantum entanglement.
A. Geometry of Fig. 1

(b): R-L coupling
In this section, we present the results of simulation for two qubits of Fig. 1(b) with the capacitively coupled R dot of qubit A and the L dot of qubit B. Exact dynamics of both qubits and the effect of interqubit coupling are consistently taken into account. Qubit A is driven by a resonant voltage ϵ cos ωt with the frequency ω ¼ J A z applied between its L and R dots at t ¼ 0. At t ¼ 0, both dots are at the north poles of their Bloch spheres. Dimensionless time τ ¼ J A z t=2π is measured in periods of the free precession of qubit A. In Fig. 2(a), color reflects Rabi nutation of qubit B driven by its coupling to qubit A with J AB ¼ 0.02J A z . The ratio J B z =J A z is plotted along the vertical axis, and four resonances are seen. Three of these occur when J B z is equal to J A z and J A z AE ϵ, while the fourth one occurs when J B z equals the Rabi frequency ϵ. While only J A z AE ϵ bands of the upper triplet are predicted by Eqs. (12,11), the strongest resonance corresponds to J B z ¼ J A z ; it corresponds to the resonant transfer of excitation between A and B [37]. The modulation of S B z may be understood as a two-step process, where, first, the oscillatory voltage at frequency ω mixes the states S A z ¼ 1 and S A z ¼ −1 of qubit A, and then the term J AB σ A x σ B x leads to oscillations in S B z . Figure 2(b) gives a similar plot but for the Rabi nutation of qubit A. It shows pronounced anomalies near all four frequencies discussed above. They provide direct demonstration of coupled dynamics of the qubits.
Insets to Fig. 2(a) demonstrate a π=2 nutation of qubit B from the north pole to the equatorial plane and back to the north pole, and its fast precession about the z axis. Observation of the nutation, by techniques of Ref. [24], should prove the coupled dynamics of qubits and allow one to evaluate J AB . Figure 2(c) plotted for J AB ¼ 0.04J A z shows, in addition to four strong resonances, a number of weaker features which are due to anharmonicity of both qubits and are not seen for J AB ¼ 0.02J A z . They increase quickly with J AB , and for J AB ¼ 0.2J A z [ Fig. 2(d)], patterns are highly irregular. The period in τ varies significantly with changes in J B z =J A z . All panels of Fig. 2 were plotted under the condition of exact resonance, ω ¼ J A z . To figure out whether this condition played an essential role in developing an irregular pattern of Fig. 2, we recalculated the figure for ω ¼ 1.1J A z and found no drastic changes to it. While details changed, the basic pattern did not. This stability may be attributed to the width of the resonance of about ϵ, which was rather large, equal to ϵ ¼ 0.3J A z in our simulations. It is worth mentioning that a large Rabi frequency ϵ ¼ 0.3J A z , while maintaining fast operation of qubit A, contributes to the multifrequency irregular dynamics. Fig. 2(a) into a single feature [see Fig. 3(a)]. It is accompanied by a weaker feature at J B z ≈ 0.04J A z , which we attribute to slow nutation of qubit B controlled by its J AB coupling to qubit A. Despite of weak driving, entanglement reaches the value E ≈ 0.95 at τ ≈ 50 (for Fig. 3(b), the J B z ≈ J A z resonance disappears while the low J B z feature grows stronger, propagates to the region of higher J B z ∼ 0.2J A z , and experiences pronounced oscillations of the Rabi nutation of B with a period Δτ ≈ 5. It may also be noted that the time scale of entanglement generation is unaffected by the decrease in Rabi frequency. Therefore, coupled nonlinear dynamics of two qubits shows high multiformity dependent on the specific choice of parameter values.
While irregular behavior is not surprising given the interplay of a number of frequencies and nonlinearity of the system, the fact that such a behavior sets in [ Fig. 2(d)] at the values of the coupling constant as small as J AB ∼ 0.04J A z and at the time scale of a single nutation of qubit B has important experimental implications. It suggests that establishing a well-controlled entanglement of two qubits imposes restrictions on the coupling strength not only from below but also from above. With increasing J AB , development of volatile behavior forestalls the increase in the operation speed. We note that early set-in of volatile behavior is related to the presence of three frequencies in the responses δn α of the L and R dots, and the geometries of Secs. IVB and IVC show lesser volatility.
While the above discussion focuses on coupled dynamics of spins S A and S B , insets to Figs. 2(b)-2(d) display a purely quantum quantity, entanglement E between qubits. Remarkably, while nutation of B in Fig. 2(c), for ϵ ¼ 0.3J A z , reaches its maximum at τ M ≈ 14, the entanglement E is represented by a two-humped curve with a maxima on both sides of τ M and a deep minimum near τ M . The same is true for insets to Figs. 2(a) and 2(b) with τ M ≈ 32. Such a behavior can be understood at a qualitative level because the nutation of B is entirely due to its coupling to A, while the nutation of A is primarily due to the driving. Hence, it is nutation of B that reflects the entanglement. Near the first maximum of irregular nutation, the dynamics of S B z is slow, and this indicates that S A and S B are nearly decoupled; hence, entanglement E is small. On the contrary, a fast change of S B z indicates strong correlations between S A and S B ; hence, E is large. The inset to Fig. 2(d) demonstrates strong, while irregular, entanglement with a pronounced peak E ≈ 1 at τ ≈ 1.5, near the first maximum of the nutation of S B .
We note that despite quite a moderate coupling strength, J AB ∼ 0.04J A z , strong nutation of B and entanglement E develop at a scale of time τ ∼ 10, still before the effects of dephasing and decoherence are expected to manifest themselves. Therefore, we envision the existence of a considerable time scale where nonlinear dynamics dominates over dephasing. Fig. 1

(c): M-L coupling
Coupled dynamics of qubits A and B in the geometry of Fig. 1(c) was calculated similarly to Sec. IVA but with a Hamiltonian of Eq. (6) rather than of Eq. (5). Because coupling of qubit B to qubit A is represented by a single matrix σ A z , the input from A is dominated by Rabi frequency ϵ, as follows from Eq. (12). Correspondingly, the response of qubit B at J AB ¼ 0.01J A z is dominated by a single frequency ϵ ¼ 0.3J A z , as seen from Fig. 4(a); satellites originating from nonlinearities are visible but very weak. Figure 4(c) demonstrates a π=2 nutation of qubit B and its precession. Figure 4(b) shows that patterns of the nutation of qubit B remain basically regular even for J AB ¼ 0.1J A z , with four oscillations occurring in the illustrated time period; however, detailed patterns show more features, as seen from Fig. 4(d). This demonstrates that purification of the signal coming from qubit A preserves stability of entangled two-qubit dynamics and suppresses transition into the irregular regime. Purification can be achieved by a proper choice of geometry or by using cavities similar to Josephson qubits [39].  3 (color online). Rabi nutation of two coupled qubits in the geometry of Fig. 1(b), where both qubits are initialized in their j0i states. Qubit A is driven at its level splitting J A z with a Rabi frequency ϵ ¼ 0.03J A z , and qubit B is coupled to A with a coupling constant J AB . Time τ is in periods of the free precession of A. The color scale is the value of S B z between −1 and 1.  Fig. 1(d). External resonant driving is applied to qubit A, and interqubit coupling is of σ A z σ B z type and described by the Hamiltonian of Eq. (7). This system possesses one integral of motion, σ B z , and therefore, the only dynamics of S B is precession about the z axis. Figure 5 displays precession of qubit B with the initial conditions S B x ðt ¼ 0Þ ¼ 1 and S A z ðt ¼ 0Þ ¼ 1. It is seen that with the exclusion of small values of x is very regular. Gross features of the column structure of the plot are described by the phase modulation of Eq. (13), while a gradual change of the color of the columns should be attributed to backaction. Fast precession of S B x is an intrinsic property of qubit B controlled by J B z , but the columnar structure of Fig. 5 [distinctly reflected in the S B x ðτÞ plot of the left inset] originates from the interqubit coupling and is controlled by J AB . This observation deserves more detailed comments. As has been mentioned below Eq. (15), E ¼ 1 is only possible for c 1 ¼ c 2 ¼ 1= ffiffi ffi 2 p ; hence, the reduced density matrices ρ A ¼ ρ B ¼ σ 0 =2, where σ 0 is a unit matrix. Therefore, S B ¼ Trðρ B σ B Þ ¼ 0 because traces of Pauli matrices vanish, and similarly, S A ¼ 0. For qubit B, S z B is conserved and equals identically to 0 due to initial conditions, S B x ≈ 0 at the maximum of E, as seen in the left inset to Fig. 5, and it has been checked that S B y ≈ 0 there (not shown). We checked that the same is true for S A .
The results of this section, especially the data of Fig. 4(d), prove the efficiency of the proposed protocol for producing entanglement of two qubits and evaluating coupling J AB between them. They also provide a warning against irregular dynamics in a two-qubit system, Fig. 2(d). It is typical of (i) the regime of fast operation, with comparable values of Rabi frequency and qubit-level splitting, and (ii) an interaction Hamiltonian with multiple noncommuting terms. In such a regime, efficient control of a twoqubit system bounds the magnitude of J AB not only from below but also from above.

V. DEPENDENCE ON INITIAL CONDITIONS
In Sec. IV, we investigated the effect of geometry on coupled two-qubit dynamics. The next problem is sensitivity of the dynamics to initial conditions. It is important both from the practical and from the conceptual point of view. First, it is difficult to control with high accuracy the initial state of a two-qubit system. Second, there is an open question of whether dynamics of such a system is stable or chaotic at a large time scale. Indeed, dynamics of a two-qubit system is described by four nonlinear differential equations for two polar angles, θ A and θ B , and two azimuths, ϕ A and ϕ B , on their Bloch spheres. Because, generically, solutions of systems of nonlinear differential equations, beginning from three equations, are prone to deterministic chaos, such a behavior is also expected for two coupled qubits. To be specific, we have in mind a deterministic classical Lorenzian chaos [40] rather than quantum chaos typical of mesoscopic systems [41,42] or the quantum chaotic border on the number of qubits [43]. Deterministic chaos manifests itself as a high sensitivity of solutions to small changes in the initial conditions.
To this end, we solved dynamic equations for a set of initial conditions as applied to the parameters of Fig. 2(d), where the dynamics was most irregular. The results, for ϵ ¼ 0.3J A z , are presented in Fig. 6 for five values of θ 0 B ¼ θ B ðt ¼ 0Þ, the polar angle of qubit B at t ¼ 0. In all cases, ϕ B ðt ¼ 0Þ ¼ 0 and θ A ðt ¼ 0Þ ¼ 0. It is seen that while each of the curves shows a rich pattern, with a number of minima and maxima, these pattern are rather similar for all five curves plotted for values of θ 0 B ranging from θ 0 B ¼ 0 to θ 0 B ¼ 48 0 . In particular, positions of all major features and even their magnitudes (in units of their initial values) remain very close. Therefore, we conclude that, at least at this time scale (which is important for experiment), the complexity of dynamics is dominated by nonlinearities of both qubits and their crosstalk and shows no signatures of chaotic behavior.

VI. CONCLUSIONS
We proposed a protocol for producing entanglement between two qubits and studied it analytically and numerically as applied to three geometries of capacitively coupled exchange-only qubits. The protocol is based on driving qubit A at its level splitting J A z and reading the response of qubit B. We found that the protocol generates entanglement when the coupling J AB between qubits is only about 1% of J A z . We have also found that the patterns of the entanglement are highly sensitive to the double-qubit geometry and the mechanism of interqubit coupling. For fast qubits with Rabi frequency ϵ comparable to J A z , dynamics shows early switching to irregular behavior if qubit B experiences a multifrequency signal from the harmonically driven qubit A. Because anharmonicity is inherent in two-level systems, we compare various coupling schemes and select geometries allowing the inhibition of early transition to irregular dynamics. To be specific, the geometry of Fig. 1(c) is more promising than the geometry of Fig. 1(b). We note that the requirement of stable entanglement of two qubits (i) may impose not only a lower but also an upper bound on the interqubit coupling constant J AB and (ii) is more easily satisfied when the interqubit coupling is dominated by a single frequency. Therefore, using cavities can facilitate entanglement not only through (i) enhancing the interqubit coupling and (ii) establishing a connection between spatially separated qubits, but also through (iii) purification of the frequency spectrum of interqubit coupling. While our calculations were performed for capacitively coupled qubits, basic conclusions are also applicable to exchange coupled qubits because both charge and exchange densities are bilinear in the wave-function amplitudes.
j111i ¼ jψ rl 111 i: (A4) We choose qubit states as j0i ≈ j111i and j1i ≈ j111i. In the basis ðj201i; j111i; j111i; j102iÞ, the total wave function is In the above basis, the eigenvalue equation takes the form where components of the vector v ¼ ðv L ; v 0 ; v 1 ; v R Þ T are amplitudes of the states j201i, j111i, j111i, and j102i, respectively. This Hamiltonian can be projected into the 2 × 2 qubit space by eliminating the ðv L ; v R Þ components The qubit is operated under the conditions t L , t R ≪ U, and therefore, E dependence of J L ðEÞ and J R ðEÞ is important mostly during the preparation and measurement cycles when jVj ∼ U [15,21]. In the central region of jVj ≪ U, the energy is small, E ∼ t 2 =U ≪ U, and can be omitted in the denominators of J L ðEÞ and J R ðEÞ. Then, after shifting the origin of the energy in Eq. (A7)by J L þ J R , we arrive at Eqs. (1) and (2).

APPENDIX B: DERIVATION OF THE INTERACTION HAMILTONIAN
For capacitively coupled qubits, the interaction Hamiltonian between excess charges on nearest-neighbor dots of the two qubits originate from the Coulomb interaction given by Eq. (4). The excess charge on each dot can be expanded in terms of the wave-function amplitudes. For a symmetric system (t L ¼ t R ¼ t), the charges on the dots