Quantum Dot Phase Transition Simulation with Hybrid Quantum Annealing via Metropolis-Adjusted Stochastic Gradient Langevin Dynamics

We report a hybrid quantum-classical simulation approach for simulating the optical phase transition observed experimentally in the ultrahigh-density type-II InAs quantum dot array. A hybrid simulation scheme, which contains stochastic gradient Langevin dynamics (a well-known Bayesian machine learning algorithm for big data) along with adiabatic quantum annealing, is developed to reproduce the experimentally observed phase transition. By implementing the simulation scheme into a quantum circuit, we successfully veriﬁed the phase transition observed in the experiment. Our work demonstrates for the ﬁrst time the feasibility of hybridizing quantum computation with classical Langevin dynamics for the analysis of carrier dynamics and quantum phase transition of the quantum dot.


Introduction
Semiconductor nanostructures have drawn significant research attention for decades due to their extraordinary lowdimensional physical properties. Among them, quantum dots (QD) are among the most noticed and researched materials and are proposed in various applications ranging from QD laser, QD memory device and QD single-photon source, and QD quantum bits for quantum computation and coherent quantum information processing [1][2][3]. e spatial confinement in QD causes the reduction of the average dielectric constant, resulting in a weak dielectric screening. e diluted screening effect in QD enables strong Coulomb interactions between the electron and hole and thus facilitates the formation of strongly paired electron-hole exciton upon optical excitation [4]. e excitonic property of QD dominates its optical and optoelectronic response and plays a significant role in quantum computation and coherent information processing [5]. Exciton Rabi oscillation in a single QD has been reported and opens excellent opportunities for coherent manipulation of the two-states quantum system [2,3]. However, a successful quantum system requires a great deal of scalability and controllability for massive quantum computation. It is technically challenging for the bottom-up strategy to fabricate a large-scale quantum system by fine-tuning and manipulating a single QD within a coherent time. On the other hand, a quantum system utilizing an ensemble of QD macroatoms (or quantum well) has been proposed in theory, which possesses the advantage to overcome these issues while enabling massive quantum computation [1]. Manipulating the quantum information process in such a large quantum structure requires deep insights into the dynamics regarding the evolution of quantum states. However, the carrier dynamics involved in the ensemble of macroatoms are much harder to clarify by a purely classical approach due to the increased quantum structure size. Recently, it has also been shown that vice-versa of the realization of the quantum computational device will enable more efficient quantum mechanics simulation over the classical computer for calculations such as electronic energy [6], reaction rate [7] as well as phase transition [8]. Meanwhile, due to the incompetence of the current quantum-computing hardware, the so-called hybrid quantum-classical algorithm specially designed for NISQ (Noisy Intermediate-Scale Quantum device) is gaining a lot of attentions. Various famous NISQ type hybrid quantum-classical algorithm such as QAOA (Quantum approximate optimization algorithm) and VQE (Variational Quantum Eigensolver) have been proposed [9].
is work proposes a hybrid quantum-classical simulation approach for simulating the phase transition in type-II InAs/GaAs QD, a typical quantum structure in III-V materials. Carrier dynamics analysis has been a hot topic from both experimental and theoretical points of view for various types of III-V nanostructure such as a nanowire [10][11][12] and quantum well tube [13][14][15]. Carrier dynamics of conventional InAs/GaAs QD have also been extensively studied experimentally and theoretically [16][17][18][19][20]. However, the excitation energy-dependent carrier dynamics and the accompanied phase transition in ultrahigh-density up to ∼10 12 cm − 2 type-II InAs/GaAs QD are poorly understood due to simulation complexity and fabrication difficulties. In this work, first, we show the experimental observation of phase transition between two well-defined excitonic states in ultrahigh-density type-II QD arrays. e temperaturedriven phase transition was further formulated using manybody electron-hole Hamiltonian schematically written as H � (H lp + H cp ) + (H es + H env ). Here, the H lp denotes the local potential within each QD, H cp defines the three Coulomb potentials through the interactions among electron-electron, hole-hole, and electron-hole), H cp defines the transition between the quantum states under external potential excitation, and H env defines the energy dissipation between the quantum system and the environments and is the main cause of quantum decoherence. At last, a hybrid quantum-classical simulation involving a time-dependent Hamiltonian was constructed to perform the adiabatic quantum annealing using a quantum gatebased quantum circuit along with the stochastic gradient Langevin dynamics.

Quantum Dot Growth and Phase Transition
e in-plane ultrahigh-density ((0.5 ∼ 1 × 10 12 cm − 2 ) InAs QDs layer was grown on a GaAsSb buffer layer and an InAsSb wetting layer (WL) on GaAs (001) substrates by molecular beam epitaxy (MBE). Detailed information regarding the experimental conditions and growth procedure can be found in literature [21]. Photoluminescence (PL) properties of this type-II heterostructure of high-density InAs QDs/GaAsSb were measured by using a time-resolved micro-PL system including an InGaAs avalanche photodiode (APD) and an InGaAs diode array. e monochromatic excitation light was selected from a white fiber-laser beam (6 W). e diameter of an excitation laser beam was about 50 μm. e sample temperature was controlled from 15 K to 300 K. Figure 1 shows the sketch of proposed carrier excitation, transportation, and recombination paths (Figure 1(a)) as well as the temperature-dependent PL experimental results of the InAs QD layer on GaAsSb under the excitation energy of 1.44 eV (Figure 1(b)). As experimentally demonstrated in our previous report [10], a critical excitation energy of 1.496 eV exists. Above or below the critical energy results in different types of PL results dominated by carrier recombination from either the ground state (GS) or excitation state (ES) of InAs QD. We have also confirmed that in the case of excitation energy below 1.496 eV, the ES carrier is not due to the state filling effect for the InAs QD under high excitation power density but rather the carrier relaxation from the InAs QD/WL/GaAsSb layer. Figure 1(a) sketches the various paths for the carrier excitation, transportation, and ES-dominating recombination under the excitation energy of 1.44 eV. It should be noted that the recombination of the carrier from ES or GS shows a typical type-II feature where there exists a spatial deviation between the electron in the QD and the hole in the WL/GaAsSb layer. Meanwhile, two ES states, ES1 and ES2, are involved in the carrier recombination and give rise to the corresponding PL features. Figure 1(b) shows the temperature dependence of the PL spectroscopy, where the three PL components: GS, ES1, and ES2 evolve in a dramatically different manner with the increase of the temperature. e GS-featured PL peak energy shifted to the longer wavelength region with the rise of temperature based on the well-known Varshni relation. On the other hand, a prominent feature was observed for ES1 and ES2 peak energy shift with temperature. It can be clearly seen that there exists a phase transition between the ES1 and ES2 PL energy peaks. e ES1 peak emerges at the onset temperature of 15 K while vanishes at around 80 K, followed by the emergence of the ES2 peak. It is worth noting here that the phase transition is not due to the multimodal distribution in QD size. It is well understood that the PL peak of large-sized QD will become dominant with the increase of temperature due to the carrier escape from QDs of smaller size families via tunneling to adjacent QDs. As a result, the PL peak energy shifts to the lower energy region at higher temperatures, opposite to the phase transition between ES1 and ES2 observed here [22]. Qualitatively, these thermal phase transitions were interpreted as the result of phonon perturbation involved relaxation, and the change of the phonon relaxation provides the change of exciton formation.

Theoretical Simulation Schemes
In the following work, we provide a rigorous theoretical simulation of the thermal dynamic behavior of this phase transition using the exciton-based quantum annealing circuit model. e central idea of the proposed model is that we introduce the excitonic distribution operators ES n , where n denotes the nth layer of high-density InAs QDs/GaAsSb ensemble (Figure 2(a)). e two types of PL peaks: ES1 and ES2 energy states with the two eigenvalues ES n � 0 and ES n � 1 correspond to the temperature-dependent excitonic peak energies in the n th layer of QD ensemble, which forms the single qubit basis as |0 and |1. e whole Hilbert space H is then spanned by the basis as |ES � n n |ES n , (ES n � 0, 1). Restricted to this computational space, the full many-body Here, ε n denotes the individual excitonic energy in the n th layer of QD ensemble array and ε nn′ specifies the Coulomb coupling effect between the QD layer of n and n ′ . e effective Hamiltonian expressed in formula (1) shares the same structure as the one proposed in the literature [23] and the one used in most of the NMR-based quantum-computing schemes [24]. In this work, to fit the model (2) into the experimental results shown in Figure 1(b), we focus only on one layer of QD ensemble, i.e., single-qubit system, thus formula (1) is reduced to the following simplest form: H � εES.
As mentioned before, two eigenvalues: ES � 0, 1 of the excitonic distribution operators ES correspond to the two PL peaks: ES1 and ES2. In this work, instead of calculating the eigenvalue directly, the excitonic distribution is represented as the distribution of exciton radius, i.e., the distance between the electron and hole based on the following consideration: e ultrahigh density of QD array studied in this    Advances in Condensed Matter Physics work could be treated as a pseudo-two-dimensional quantum well structure due to the strong coupling between neighboring QDs. Moreover, meanwhile, the exciton binding energy of quantum well can be approximated as Eexciton ≈ − (ε r r eh ) − 1 based on formulas (28), (29), and (31) in reference [25], where ε r is the dielectric constant and r eh is the mean electron-hole distance. Experimentally, the exciton binding energy can be estimated from the PL peak energy using the following relation: E exciton � E G (orE HOMO− LUMO ) − E PL . us, in the condition where the bandgap E G or HOMO and LUMO energy levels stay constant, the PL peak energy directly reflects the exciton binding energy and vice versa. In addition, we treat the distinct evolution of PL peaks of ES1 and ES2 with the change of temperature as a temperature-dependent phase transition problem. e dynamic evolution of electron-hole pair distribution morphology governs the phase separation. In order to map electron-hole pair distribution, we adopted a Ising model-based morphology generation method widely used in the field of small molecule blends in the organic photovoltaic field [26,27].
A hybrid quantum-classical version of the Metropolisadjusted Langevin algorithm was employed in this work to investigate the temperature-dependent dynamics and the phase transition. At first, we will give a general picture of the simulation scheme designed in this work. Figure 2(b) shows the exciton model constructed based on the experimental results. In this model, the electron is localized in each QD while the hole can freely move in the in-plane band of the InAs wetting layer, thus possessing a large Bohr radius. e proposed model is consistent with the experimentally verified carrier dynamics model shown in reference [28]. Due to this irregular arrangement, the exciton in this work can only be treated as a hybrid type of the Frenkel and Wannier exciton [29]. Figure 2(c) shows the coupling distance distribution among electrons and holes. In order to facilitate simulation efficiency and save simulation cost, we define an effective coupling radius R eff and only the exciton within R eff are used to create the Ising morphology.
e Ising model is constructed as follows: Here, q i , q j are either fixed electrons or mobile holes; ε is the dielectric constant; and R ij is the distance between particle q, q j .
Moreover, in order to investigate the temperature-induced quantum phase transition, we built a quantized version of the Ising model by converting the particle q i , q j to σ Z i and σ Z j as well as introducing a transverse operator σ x i and time-varied excitation temperature T(t) (the transverse operator can be viewed as the phonon perturbation mentioned in section 2 to trigger the phase transition between (ES1 and ES2)). e total quantum Ising Hamiltonian is then represented as the following quantum annealing model: In this work, we adopt a more special quantum annealing model: so-called adiabatic quantum annealing model [30], which takes the form as follows: where s is a normalized temperature parameter taking the value within 0∼1, instead of varying the quantum state, the Hamiltonian shape is changed overtime to reach the ground state eigenenergy. e adiabatic quantum annealing model can be solved by a classical computer using Suzuki-Trotter decomposition [20], the so-called quantum Monte Carlo approach, or solved by quantum annealer hardware such as the one developed by D-Wave [31]. In this work, we adopted the quantum circuit-based quantum annealing approach. However, since the quantum annealer hardware is not available for us at this moment, thus an in-house developed quantum annealer simulator is used as an alternative. We also claim that using quantum annealing to find the ground state energy of formula (3) is not the only choice. For small size problems, we could quickly obtain a rigorous and exact solution for formula (3). e aim of employing quantum annealing in our work is purposed on the case where the computation time of the problem under investigation increases exponentially, thus suffering the so-called curse of dimension.
A hybrid simulation scheme (Langevin dynamics + quantum annealer) is developed to fulfill the purpose of the simulation of the dynamic morphology evolution of the excitons. Simply speaking, the kinetic update of the hole positions is conducted by following the stochastic gradient Langevin dynamics [32] over the Ising energy landscape. Stochastic gradient Langevin dynamics is a well-known Bayesian machine learning algorithm for learning big data. Langevin dynamics is usually expressed in the following way by introducing a stochastic process into Newton's second law: Here, m represents mass, a is the acceleration, F represents the force, v is velocity, and ε is stochastic variable with the mean being zero and variance being some timedependent parameter. Since mass is usually small, the left side of formula (4) can be assumed as zero. Under this assumption, along with slight reformulation, formula (4) can be further written as follows: Formula (5) can be further simplified as the follows by assuming the force F is the gradient of potential U: 4 Advances in Condensed Matter Physics e abovementioned is the core contents of the Langevin dynamics. e stochastic gradient Langevin dynamics introduce the so-called minibatch-sized stochastic gradient when calculating ∇U, which significantly reduce the computation time by sampling over a minibatch with sample size |S| being far smaller than the whole sample size N. Recently, to enhance further the converge rate, Metropolis-adjusted Langevin algorithms have been developed [33,34]. In this work, we adopted the method from reference [34] by adding a Metropolis-Hasting step to stochastic gradient Langevin dynamics. us, after the update of the holes, ∅ ij in formula (2) is calculated and serve as the input of formula (10) for the adiabatic quantum annealer. e acceptance/rejection process is finally implemented based on the Metropolis algorithm using the output energy from quantum annealing. After the Metropolis process is converged by monitoring the update of hole position at each specified temperature, the histogram over the sampled distance R ij between the hole-electron pairs is used to identify the quantum state of exciton ES1 and ES2, as shown in Figure 2(e). Since the energy is inverse proportional to the distance R ij , high energy state ES2 is supposed to be located at the right side of low energy state ES1.
Next, we will give a brief discussion about the quantum annealing circuit used to simulate the ground state energy of Hamiltonian in formula (7) at each specified temperature. From the time-dependent Schrödinger equation, the wavefunction evolves based on the following unitary operator: Following the Suzuki-Trotter decomposition, the unitary operator can be approximated as follows: where

Advances in Condensed Matter Physics
In principle, the number of qubits should equal the number of particles under investigation. For instance, if 100 particles (40 electrons ⊕ 60 holes) are used in the simulation, we thus need 100 qubits to prepare the quantum circuit, which is not possible currently due to the hardware and memory constraints for both classical and quantum computers. We, therefore, constructed the following approximation for a simulation to be performed within the reachability of available computation resources. (i) since we adopted a stochastic gradient Langevin dynamics, an update of the hole location is performed sequentially; in other words, we have set the minibatch size of |S| � 1 for performing the stochastic gradient Langevin dynamics simulation. Under this condition, the required qubit will be reduced to one hole plus the coupling electrons. In the 100 particle examples, the qubit will be reduced to 41 (40 electrons ⊕ one hole). (ii) We could further reduce the number of qubits by treating the electrons localized in each quantum dot as one single ensemble electron. e same assumption is also applied to the other holes because the location change of holes is extremely small within each update. Under this assumption, we further treat the one ensemble electron and one ensemble hole as one ensemble exciton. erefore, the final qubit needed to implement the

Algorithm: Phase Transition Simulation by Hybrid Quantum Classical Algorithm
13 end if 14 15 Store the value of the exciton distance r that is within the effective range adiabatic quantum annealing circuit is reduced to two qubits (one ensemble exciton ⊕ one hole). Based on the approximation mentioned above, we introduce the following approximation formula: where i stands for the i th hole; υ stands for the single ensemble exciton; ∅ ij is the Coulomb coupling between each particle; ∅ i is the marginalized Coulomb coupling over the j th particle inside the ensemble exciton. Figure 3 shows the quantum circuit based on the assumption used in formula (10) to implement the unitary operator U zz (s) and U x (s).
e coupling strength ∅ ij can be calculated using formula (2) and is converted to the phase angle θ as the input parameter for unitary operator U 3 (θ) and U 1 (θ) shown in the upper part of Figure 3. Iterating k times of the single unit showing in the lower part of Figure 3 corresponds to the summation of the coupling energy between i th and j th particle in the exponential function of formula (7). In this work, the number of iterating times k are optimized to by considering both the accuracy and computation cost, and k � 300 is used for all the simulations. e pseudocode of the whole simulation algorithm is given in Figure 4. Figure 5 shows the simulation results based on the adiabatic quantum annealing simulator described in Figure 3. Here, we have fixed the total number of electrons and holes as 100 while tuning the ratio between electron and hole. Figure 5 shows the ratio of 4 : 6, i.e., 40 electrons at fixed positions and 60 holes updating the position based on the stochastic gradient Langevin dynamics and Metropolis acceptance probability. For better clarity, the final hole positions at each temperature (left side in Figure 5) and the histogram of electron-hole distance distribution (right side in Figure 4) for four typical temperatures: 15 K, 80 K, 120 K, and 150 K were presented. Since the x-axis is defined as the distance between the electron and hole, the peak located on the left side represents higher exciton binding energy following the equation: E exciton ≈ (− ε r r eh ) − 1 . Meanwhile, in order to extract the ES1 and ES2 states, double peak Gaussian fitting were conducted, shown as the bold solid lines for each temperature. It can be seen that there exists a phase transition at the temperature between 80 K and 120 K, at which the lower energy ES1 state vanishes while the ES2 state emerges. is tendency corresponds well with one of the experimental results shown in Figure 1(b), indicating the validity of our simulation results. Figure 6 shows the convergence behavior of the potential expectation value of the operator U ZZ , which represents the coupling energy between the one ensemble exciton and one hole. Before we explain the simulation results, we first briefly explain how the expectation value is calculated. e following calculation is performed at each temperature to monitor the potential energy convergence.

Results and Discussion
Here, the coefficients c 1 c 2 c 3 c 4 of wavefunction, Ψ is obtained analytically from the quantum annealing simulator. erefore, formula (11) can be further simplified as follows:  Figure 6: Convergence behavior of hybrid of adiabatic quantum annealing simulator with Metropolis-adjusted stochastic gradient Langevin dynamics; the expectation potential value for operator U zz is simulated at four typical temperatures: 15K, 80K, 120K, and 150 K. 8 Advances in Condensed Matter Physics We must point out that in a real quantum computer, the coefficients of Ψ are usually not available and the expectation value of U ZZ can then be evaluated using formula (13) by measuring the probability P of the two qubits (qb 1 , qb 2 ) being either "0" or "1".
e potential energy drops for the four temperatures under investigation as the number of iterations increases. Although most convergence display a quite similar behavior, the one studied under 80 K indicates an insufficient number of iterations. We have tried to increase the iteration steps, but the resultant exciton distribution becomes dramatically different from the one shown in Figure 5. As having been mentioned before, a temperature of 80 K corresponds approximately to the critical temperature for the phase transition thus, we recon it is difficult for the system to reach an equilibrium between the two phases. For the temperature of 15 K and 120 K, the convergence curve shows stair-like reduction, indicating the effectiveness of the hybrid quantumclassical algorithm overcoming the saddle point (or local minimum). For the temperature of 150 K, the convergence curve shows an early drop onset and a monotonic decrease with the increase of iteration due to the additional thermal driving force. We have also conducted simulations by varying the electron/hole ratio parameters. Figure 7 shows the simulation results under three different electron/hole ratios. It can be easily seen here that an optimal electron/hole ratio exists to reproduce the experimental results. Compared to the ratio of 1 : 9 and 7 : 3, the ratio of 4 : 6 generates the closest results to the experimental ones. Although many uncertain factors are  altering the simulation results, a roughly balanced number of electrons and holes support the exciton model proposed in this work. At last, we must point out that the current simulation has a strong dependence on the random number generator seed. Since this is the first trial of hybridizing quantum annealing simulator with the Metropolis-adjusted stochastic gradient Langevin dynamics, there are still many rooms to improve the quality of the simulation results presented in the current form in future. Parameters such as the iteration number for the Suzuki-Trotter were not optimized but simply determined based on our experience due to the constraints of simulation times. Due to a similar reason, the temperature interval is also adjusted based on experience. We have put forward our code on GitHub with an open-access link so the results could be easily verified and reproduced.

Conclusions
In conclusion, we reported a theoretical analysis of the optical phase transition observed experimentally in the ultrahighdensity quantum dot array investigated by the exciton-based quantum annealing circuit. A hybrid simulation scheme, which contains stochastic gradient Langevin dynamics along with adiabatic quantum annealing, is developed to reproduce the experimentally observed phase transition. By implementing the simulation scheme into the quantum circuit, we successfully verified the phase transition observed in the experiment. Our work demonstrates for the first time the feasibility of hybridizing the quantum computation with classical Langevin dynamics for the analysis of carrier dynamics and quantum phase transition of the quantum dot. In the current work, we focused only on the in-plane coupling of single layers of QD array and built a single qubit model, and it will be our next step to construct and investigate the multiqubit system and quantum dynamics by including the vertical coupling effect in the QD array.

Data Availability
All the simulation code for reproducing the simulation results can be accessed from the following link: https:// github.com/KShiba24/QDPL_simulator.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.