Quantum Control of Exciton Motion in Electric Field

1. Introduction

An exciton is a bound pair of an electron and a hole held together by Coulomb interaction [1]. In semiconductors, excitonic states are typically classified into two types. The first is the Frenkel exciton, characterized by a large binding energy and a short distance between the electron and hole, often localized on a single molecule. The second category is the Wannier-Mott exciton, which exhibits a smaller binding energy and a separation between the electron and hole significantly larger than the lattice constant.

In two-dimensional (2D) materials, the bound states of electrons and holes form 2D excitons. In modern single-layer materials such as transition metal dichalcogenides, due to their ultrathin nature that results in reduced dielectric screening, excitons can exhibit relatively higher binding energies [2] and remain bound even at room temperature. Moreover, 2D excitons are highly tunable under external electric and magnetic fields [3], making them a critical area of research in low-dimensional materials science.

In semiconductors such as GaAs, excitons behave analogously to the proton-electron pairs in a 2D hydrogen atom [4]. The effective Bohr radius of an exciton is given by $a_{\mathrm{ex}}=ϵm{a}_0/m^{*}=11.86$ nm, where $a_0$ is the Bohr radius, $ϵ$ is the dielectric constant of GaAs, m denotes the electron mass, and $m^{*}$ is the reduced mass of the electron-hole pair defined by $1/m^{*}=1/m_e^{*}+1/m_h^{*}.$ Here $m_e^{*}$ and $m_h^{*}$ are the effective masses of the electron and hole, respectively. The energy eigenvalues of exciton are expressed as

$$E_{\mathrm{ex}}=-{\displaystyle \frac{1}{{(n-1/2)}^2}}{\displaystyle \frac{m^{*}}{m}}{\displaystyle \frac{1}{{ϵ}^2}}{\displaystyle \frac{e^2}{2a_0}},$$

(1)

where e is the electron charge and n is the principal quantum number, an integer.

Excitons in semiconductors can be classified as direct or indirect based on whether the electron and hole reside in the same or different quantum well(s) (QW) or material layer(s) [5]. Indirect excitons, characterized by their extended lifetimes ranging from microseconds to milliseconds, facilitate more effective cooling and enable the formation of long-range coherent condensates under applied electric fields at low temperatures [6]. The electric field generated by split-gate electrodes plays a crucial role in controlling the transport of these indirect excitons [7] and enables the switching functionality of optoelectronic transistors [8]. By tuning the voltage applied to the electrodes, the resulting electric field distribution can be precisely modified, allowing for the manipulation of exciton energy and the creation of excitonic potential barriers essential for exciton collection and cooling [9]. Additionally, the electric field of the gate facilitates the formation of interlayer excitons and regulates their transport within excitonic transistors, thereby significantly influencing the performance and efficiency of optoelectronic devices [10]. Excitons trapped in a 2D lattice potential can be manipulated by varying the field direction or strength, enabling control over their transport and phase transitions [5].

Furthermore, excitons offer innovative pathways for developing efficient and low-cost microscale and nanoscale optoelectronic devices. Excitonic solar cells are expected to become a key direction for high-efficiency solar cells in the future [11]. Excitonic micro-LED fabricated on GaN nanowires have achieved record-high efficiencies for submicron scale LED [12], demonstrating the significant potential of excitons in advancing optoelectronic technology. These and others applications of excitons require a control of their position and velocities.

In this paper, we concentrate on direct excitons, where the electron and hole are located in the same quantum well. To achieve control over exciton dynamics, we employ a simple charged stripe-like gate. By modulating charge density on the gate, we produce an electric field with spatially varying time-dependent distribution. This modulation of the electric field enables us to manipulate the quantum degree of freedom [13], specifically the exciton dipole moment [14], thereby controlling its classical motion in this field.

In more detail, when exciton is located in electric field, a displacement occurs between the electron and hole, inducing a dipole moment dependent on the strength and direction of the applied field. As a result, in a nonuniform field a net force acting on the exciton appears and causes its motion. The control of classical motion can be achieved by modifying the gate charge altering the electric field and influencing the dynamics of the exciton via the time-dependent dipole moment. This control should satisfy the desired control tasks, e.g., in terms of the final position and velocity within a specified evolution time and given initial conditions. To achieve this control, we employ an optimal search method for the design of the gate charge producing the time-dependent quantum dipole moment of the exciton resulting in the realization of the control task for its classical motion.

The rest of the paper is organized as follows. In Sec. Section 2, we present the models of the exciton quantum dynamics, its classical motion and of the electric field of the gate. With this knowledge, we introduce in Sec. Section 3 the search problem, develop an optimized search method, and present the search results for various dynamical regimes. Finally, we discuss possible relations to the experiment and formulate the conclusions of this paper in Sec. Section 4.

2. Model and equations of motion

2.1. Dipole Moment

When an exciton, initially positioned at $\overline{x}\left(0\right)$, moves in an external electric field $\mathcal{E}\left(\overline{x}\left(t\right)\right)$ oriented parallel to the x-axis, the system is described by the following Hamiltonian

$$H=H_0+H^{\prime }=H_0-d·[\mathcal{E}\left(\overline{x}\left(t\right)\right)-\mathcal{E}\left(\overline{x}\left(0\right)\right)],$$

(2)

here, $H_0$ denotes the Hamiltonian of the free two-level system $\left\{\right|g_0\rangle ,\left|e_0\rangle \right\}$, with energies $E_g$ and $E_e$, respectively. It can be expressed as $H_0=E_g|g_0\rangle \langle g_0|+E_e|e_0\rangle \langle e_0|=\hslash {\omega }_0{\sigma }_z/2$, where ${\sigma }_z$ is the Pauli matrix, and the atomic transition frequency ${\omega }_0$ is defined by $\hslash {\omega }_0=E_e-E_g$. The dipole moment is given by $d=-e(\widehat{x}-\overline{x}\left(t\right))$ and is oriented parallel to the x-axis. Under the dipole approximation, the interaction Hamiltonian is $H^{\prime }=e(\widehat{x}-\overline{x}\left(0\right))·[\mathcal{E}\left(\overline{x}\left(t\right)\right)-\mathcal{E}\left(\overline{x}\left(0\right)\right)]$.

We expand the dipole moment operator in terms of the eigenstates of $H_0$ as

$$d=d_{ge}|g_0\rangle \langle e_0|+d_{eg}|e_0\rangle \langle g_0|,$$

(3)

where $d_{ge}$ is the matrix element of the dipole moment, and we assume $d_{ge}=d_{eg}=\mu$. At this point, the interaction Hamiltonian in Eq. (2) is given by

$$H^{\prime }=-d·\left[\mathcal{E}\left(\overline{x}\left(t\right)\right)-\mathcal{E}\left(\overline{x}\left(0\right)\right)\right]=-\mu {\sigma }_x·\left[\mathcal{E}\left(\overline{x}\left(t\right)\right)-\mathcal{E}\left(\overline{x}\left(0\right)\right)\right],$$

(4)

where ${\sigma }_x$ is the corresponding Pauli matrix. For this interacting system, we transform the Hamiltonian to the interaction representation, yielding $H_I^{\prime }\left(t\right)=U_0^{†}\left(t\right)H^{\prime }{U}_0\left(t\right)$, where $U_0\left(t\right)=exp\left(-i{H}_{0}t/\hslash \right)=exp\left(-i{\omega }_0{\sigma }_{z}t/2\right)$. All operators, including the Pauli, Hamiltonian, and dipole moment operators, are transformed into the interaction representation. For example, the dipole moment operator becomes $d\left(t\right)=d_{eg}({\sigma }_x+i{\sigma }_y)e^{i{\omega }_{0}t}/2+\mathrm{H}.\mathrm{c}.$, where H.c. denotes the Hermitian conjugate. Assuming the system’s wave function takes the form $|\mathrm{\Psi }\left(t\right)\rangle =a_e\left(t\right)|e_0\rangle +a_g\left(t\right)|g_0\rangle$,where $a_{e,g}\left(t\right)$ are complex coefficients, substituting this into the Schrödinger equation $i\hslash {\partial }_t|\mathrm{\Psi }\left(t\right)\rangle =H_I|\mathrm{\Psi }\left(t\right)\rangle$ yields a set of coupled partial differential equations

$$\begin{array}{cc}\hfill i\hslash {\dot{a}}_e\left(t\right)& =-\mu a_g\left(t\right)e^{i{\omega }_{0}t}[\mathcal{E}\left(\overline{x}\left(t\right)\right)-\mathcal{E}\left(\overline{x}\left(0\right)\right)],\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill i\hslash {\dot{a}}_g\left(t\right)& =-\mu a_e\left(t\right)e^{-i{\omega }_{0}t}[\mathcal{E}\left(\overline{x}\left(t\right)\right)-\mathcal{E}\left(\overline{x}\left(0\right)\right)].\hfill \end{array}$$

(6)

where ${\dot{a}}_g\left(t\right)$ and ${\dot{a}}_g\left(t\right)$ are derivatives with respect to t.

At $t=0$, the total Hamiltonian in the interaction representation is

$$H=\left(\begin{array}{cc}{E}_g& -\mu \mathcal{E}\left(\overline{x}\left(0\right)\right)\\ -\mu \mathcal{E}\left(\overline{x}\left(0\right)\right)& E_e\end{array}\right)=\left(\begin{array}{cc}{E}_g& V\\ V& E_e\end{array}\right).$$

(7)

The eigenvectors of the Hamiltonian given by Eq. (7) are $|{\psi }_{+}\rangle =cos\left(\theta /2\right)|g_0\rangle +sin\left(\theta /2\right)|e_0\rangle$ and $|{\psi }_{-}\rangle =sin\left(\theta /2\right)|g_0\rangle -cos\left(\theta /2\right)|e_0\rangle$, with the mixing angle $\theta \equiv arctan[V/\mathrm{\Delta }]$ and eigenvalues $E_{\pm }=E\pm \sqrt{{(\hslash {\omega }_0/2)}^2+V^2}$ where $E=(E_g+E_e)/2.$ At this point, the expectation value of the dipole moment for the exciton in the ground state $|{\psi }_{-}\rangle$ is given by

$$\langle d\left(0\right)\rangle =-e\langle {\psi }_{-}|\widehat{x}-\overline{x}\left(t\right)|{\psi }_{-}\rangle =sin\left[arctan{\displaystyle \frac{2\mu \mathcal{E}\left(\overline{x}\left(0\right)\right)}{\hslash {\omega }_0}}\right]={\displaystyle \frac{2\mu \mathcal{E}\left(\overline{x}\left(0\right)\right)}{\mathrm{\Delta }\left(\overline{x}\left(0\right)\right)}},$$

(8)

where the gap $\mathrm{\Delta }\left(\overline{x}\left(0\right)\right)\equiv \sqrt{4{\mu }^2{\mathcal{E}}^2\left(\overline{x}\left(0\right)\right)+{\hslash }^2{\omega }_0^2}.$ Equation (8) indicates that the induced dipole moment depends on the external electric field and does not vary linearly with the field. The expectation value of the time-dependent dipole moment in interaction representation is given by

$$\langle d\left(t\right)\rangle =-2Re\left\{d_{eg}{e}^{i{\omega }_{0}t}{a}_e^{*}{a}_g\right\}.$$

(9)

Compared to $H_0$, the term $H^{\prime }$ is small and can be treated as a perturbation. The disturbance induced by the external electric field can be quantified using perturbation theory through energy corrections. For instance, when $t=0$, an exciton initially in the ground state at position a is subjected to the electric field $\mathcal{E}\left(\overline{x}\left(0\right)\right)$. If the external field is nonzero, it will generate an induced dipole moment. The Hamiltonian of the system in this case is given by $H=H_0+H^{\prime }=H_0+\mathcal{E}\left(\overline{x}\right)·(\widehat{x}-\overline{x}\left(0\right))$.

Using time-independent perturbation theory, the first-order correction eigenstates are $|g\rangle =|g_0\rangle -\mathcal{E}\left(\overline{x}\left(0\right)\right)/\hslash {\omega }_0|e_0\rangle$ and $|e\rangle =\mathcal{E}\left(\overline{x}\left(0\right)\right)/\hslash {\omega }_0|g_0\rangle +|e_0\rangle .$ Thus, the expectation value of the dipole moment of the exciton in the ground state $|g\rangle$ at $t=0$ is $\langle d\left(0\right)\rangle =-e\langle g|\widehat{x}-\overline{x}\left(0\right)|g\rangle =2\mathcal{E}\left(\overline{x}\left(0\right)\right)/\hslash {\omega }_0$.

When the external electric field is nonzero, or the initial velocity of the exciton is nonzero, the exciton will move within the field as time evolves. The Hamiltonian governing this motion is given by Eq. (2). By analogy with standard time-dependent perturbation theory, we express $|\mathrm{\Psi }\left(t\right)\rangle$ as an expansion in terms of the previously obtained corrected eigenstates $|\mathrm{\Psi }\left(t\right)\rangle =c_g\left(t\right)e^{-i{E}_{g}t/\hslash }|g\rangle +c_e\left(t\right)e^{-i{E}_{e}t/\hslash }|e\rangle$. Assuming the exciton initially occupies the lower state, the first-order approximation for the coefficients is given by

$$c_g\left(t\right)=1,c_e\left(t\right)=-{\displaystyle \frac{i}{\hslash }}{\int }_0^t{H}_{eg}^{\prime }\left(t^{\prime }\right)e^{i{\omega }_0{t}^{\prime }}d{t}^{\prime },$$

(10)

where $H_{eg}^{\prime }\left(t\right)\equiv \langle e|H^{\prime }|g\rangle$. The induced time-dependent dipole moment is then given by

$$\begin{array}{ccc}& & \langle d\left(t\right)\rangle =-e\langle \mathrm{\Psi }\left(t\right)\left|(\widehat{x}-\overline{x}\left(t\right))\right|\mathrm{\Psi }\left(t\right)\rangle \hfill \\ & & ={\displaystyle \frac{2\mathcal{E}\left(\overline{x}\left(0\right)\right)}{\hslash {\omega }_0}}cos\left({\omega }_{0}t\right)-{\displaystyle \frac{2}{\hslash }}cos\left({\omega }_{0}t\right){\int }_0^{t}sin\left({\omega }_0{t}^{\prime }\right)\mathcal{E}\left(\overline{x}\left(t^{\prime }\right)\right)d{t}^{\prime }\hfill \\ & & +{\displaystyle \frac{2}{\hslash }}sin\left({\omega }_{0}t\right){\int }_0^{t}cos\left({\omega }_0{t}^{\prime }\right)\mathcal{E}\left(\overline{x}\left(t^{\prime }\right)\right)d{t}^{\prime }.\hfill \end{array}$$

(11)

At $t=0$, the time-dependent dipole moment in Eq. (11) is $\langle d\left(0\right)\rangle =2\mathcal{E}\left(\overline{x}\left(0\right)\right)/\hslash {\omega }_0$, which agrees with the result obtained from the time-independent perturbation theory.

The energy of a dipole in an external field is given by $-d\left(t\right)·\mathcal{E}\left(\overline{x}\left(t\right)\right)$. Spatial variations in the potential energy (i.e., the gradient of the potential) induce the force acting on exciton resulting in its acceleration, expressed as $M\ddot{\overline{x}}\left(t\right)=-\nabla U(\overline{x},t)=(d·\nabla )\mathcal{E}$, where $M=m_e^{*}+m_h^{*}$ is the mass of the exciton. By incorporating the dipole moment expression in Eq. (9), the system’s dynamical behavior can be described by

$$M\ddot{\overline{x}}\left(t\right)={\displaystyle \frac{\partial \mathcal{E}\left(\overline{x}\left(t\right)\right)}{\partial \overline{x}\left(t\right)}}\langle d\left(t\right)\rangle .$$

(12)

Henceforth, unless specified, we omit the t-dependence for brevity, e.g. $a_{g,e}\left(t\right)\equiv a_{g,e}$. By combining Newton’s equation of motion, Eq. (12), with the dipole moment expression, Eq. (9), and coupling them with Eqs. (5) and (), we can numerically solve the system under the initial conditions $a_e\left(0\right)=-cos\left(\theta /2\right)$ and $a_g\left(0\right)=sin\left(\theta /2\right)$. This approach allows us to determine the exciton’s trajectory in the external electric field and the induced time-dependent dipole moment.

Under the adiabatic approximation of slow time evolution of the electric field, the dipole moment follows the electric field $\mathcal{E}\left(\overline{x}\left(t\right)\right)$ as described by Eq. (8). The corresponding classical equation of motion in adiabatic limit is given by

$$M\ddot{\overline{x}}\left(t\right)={\displaystyle \frac{\partial \mathcal{E}\left(\overline{x}\left(t\right)\right)}{\partial \overline{x}\left(t\right)}}{\displaystyle \frac{2\mu \mathcal{E}\left(\overline{x}\left(t\right)\right)}{\sqrt{4{\mu }^2\mathcal{E}{\left(\overline{x}\left(t\right)\right)}^2+{\hslash }^2{\omega }_0^2}}}.$$

(13)

The slow motion of the exciton requires the electric field to vary sufficiently slowly, ensuring that the probability of non-adiabatic transitions between instantaneous eigenstates is very small. We present this condition in the form $d\mathrm{\Delta }\left(\overline{x}\left(t\right)\right)/dt\ll {\mathrm{\Delta }}^2\left(\overline{x}\left(t\right)\right),$ where $\mathrm{\Delta }\left(\overline{x}\left(t\right)\right)$ is the instantaneous gap. Here $d\mathrm{\Delta }\left(\overline{x}\left(t\right)\right)/dt=2{\mathrm{\Delta }}^{-1}\left(\overline{x}\left(t\right)\right){\mu }^{2}d{\mathcal{E}}^2\left(\overline{x}\left(t\right)\right)/dt$ with $d{\mathcal{E}}^2\left(\overline{x}\left(t\right)\right)/dt=\partial {\mathcal{E}}^2\left(\overline{x}\left(t\right)\right)/\partial t+\overline{v}\left(t\right)\partial {\mathcal{E}}^2\left(\overline{x}\left(t\right)\right)/\partial x.$ The latter expression implies that the adiabaticity condition depends both on the change of the electric field with time as determined by the charge of the gate and by exciton velocity.

Now we introduce the units. The electron charge $e=-1$, and the reduced Planck constant $\hslash =1.$ The units for frequency, mass and dipole moment are taken as ${\omega }_0=M=\mu =1$. Consequently, the unit of length is defined as $\langle g_0\left|\widehat{x}\right|e_0\rangle$ for the exciton with $\overline{x}=0.$

2.2. Electric field of the gate

As shown in Figure. Figure 1 (a), we model the gate by a uniformly charged stripe infinitely extended in the y-direction $-\infty <y<\infty$ and having a width w along the x-direction. The resulting electric field in Gaussian units at the point $\overline{x}\left(t\right)$ is

$$\mathcal{E}\left(x\right)={\displaystyle \frac{\rho \left(t\right)}{ϵ}}ln{\displaystyle \frac{{(w/2+x)}^2+z_0^2}{{(w/2-x)}^2+z_0^2}}.$$

(14)

Below we define $\overline{\rho }\left(t\right)\equiv \rho \left(t\right)/ϵ.$

Figure 2 illustrates the electric field produced by the gate when the charge density is constant, $\overline{\rho }\left(t\right)=1$, at various positions $x.$ Near the center of the gate, the electric fields from either side almost cancel each other resulting at $\left|x\right|\ll w$ and $\left|x\right|\ll z_0$ in $\mathcal{E}\left(x\right)\approx 2\overline{\rho }wx/(z_0^2+w^2/4).$ At $\left|x\right|\gg w$ and $\left|x\right|\gg z_0$ one obtains $\mathcal{E}\left(x\right)\approx 2\overline{\rho }w/x.$ Figure 2 represents the $x-$dependent electric field demonstrating that its decrease at a large distance is $z_0-$ independent.

3. Control of classical motion

3.1. Search Problem

The motion of an exciton in an external electric field depends on the initial velocity $\overline{v}\left(0\right)$ and position $\overline{x}\left(0\right)$ and responds dynamically to variations in the electric field.

Since we are interested in the quantum control of classical motion, we formulate the conditions for the final position and velocity without specifying the final quantum state as:

$$x_{1,f}<\overline{x}\left(t_f\right)<x_{2,f};v_{1,f}<\overline{v}\left(t_f\right)<v_{2,f},$$

(15)

with the initial conditions $\overline{x}\left(0\right)$ and $\overline{v}\left(0\right).$ Thus, our objective is to control the exciton’s final position and velocity by modifying the external electric field.

For the sake of generality, we mention the relation of this problem to other aspects of quantum control and quantum mechanics. This problem is considerably different from the set of problems studied in the shortcuts to adiabaticity [15] since we consider the effect of the quantum dynamics on the classical motion, not considered in the shortcuts to adiabaticity approaches. However, this problem can be seen as an inverse classical version of the time-of-arrival problem in the fundamentals of quantum mechanics [16]. In addition, it can be related to the problem of quantum backflow [17,18] where the interaction of the quantum motion of the electron and hole can lead to modification of the classical flux of excitons.

3.2. Search Method

Here, we produce a time-dependent electric field by modulating the charge density $\overline{\rho }\left(t\right)$. To control $\overline{\rho }\left(t\right)$, we adopt the first kind Chebyshev polynomial Ansatz $\overline{\rho }\left(t\right)={\sum }_{j=0}^n{c}_j{T}_j\left(t\right)$, subject to the boundary conditions $\overline{\rho }\left(0\right)=\overline{\rho }\left(t_f\right)=0$. The first kind Chebyshev polynomials defined as $T_j\left(x\right)=cos(jarccosx)$ are closely related to trigonometric functions, making it possible to approximate complex functions with relatively simple polynomials. For instance, $T_0\left(x\right)=1$, $T_1\left(x\right)=x$, $T_2\left(x\right)=2x^2-1$, $T_3\left(x\right)=4x^3-3x$, and $T_4\left(x\right)=8x^4-8x^2+1$, illustrate how higher-degree polynomials can be constructed recursively. The coefficients $c_j$ are determined by applying the Ansatz in conjunction with the specified boundary conditions. For $n=4$, the pulse-shaped charge density distribution can be adjusted using three free parameters $c_0,c_1,$ and $c_2$. Imposing $\overline{\rho }\left(0\right)=0$ and $\overline{\rho }\left(t_f\right)=0$ allows $c_3$ and $c_4$ to be expressed in terms of $c_0$, $c_1$, and $c_2$, yielding

$$\overline{\rho }\left(s\right)=c_0{T}_0\left(s\right)+c_1{T}_1\left(s\right)+c_2{T}_2\left(s\right)-(c_1+2c_2)T_3\left(s\right)-(c_0-c_2)T_4\left(s\right),$$

(16)

with $s=t/t_f.$ By adjusting various parameter combinations, such as $c_0$, $c_1$ and $c_2$ in Eq. (16), we can guide the exciton’s motion through the resulting electric field.

To address the search problem defined above, we focus on identifying combinations of three parameters $c_0,c_1,$ and $c_2$ under a certain initial condition. The combinations we are expected to identify should satisfy two final conditions. Firstly, the exciton’s final position falls within the interval $[x_{1,f},x_{2,f}]$. Secondly, the the final velocity falls within $[v_{1,f},v_{2,f}]$. The variables satisfy these conditions are defined as "valid points" while the "invalid points" are the combinations that failed to satisfy the conditions.

This search problem, involving three parameters, can be considered a three-dimensional ($3D$) problem, thus can be mapped as a cube in three dimensional space, as illustrated in Figure 3. The three edges of the cube correspond to the parameters $c_0$, $c_1$, and $c_2$, each with a length of $10.24$, indicating that the the range of the parameters is $[0,10.24]$. Each point within the cube signifies a unique combination of $c_0$, $c_1$, and $c_2$. By discretizing these parameters using the specified step size ${\delta }_c$, we can systematically evaluate whether each point is valid or not within the cube.

To identify the valid points, the most straightforward approach is to directly calculating all points. Although this full-calculation strategy ensures comprehensive exploration, it becomes computationally inefficient for highly complex systems. The inefficiency turns severe when the degrees of freedom is large, or the parameter discretization is high. The computational time may increase exponentially, and a higher portion of the calculating resources might be wasted on irrelevant or inefficient regions. In addition to the direct calculation method, other methods could increase efficiency while losing some precision of identification.

Therefore, the efficient method we expect to develop has to increase the computational efficiency and keep precision. To achieve this, we outline the following steps.

• Step 1: We begin with selecting a relatively larger initial step size ${\delta }_c$ and employing direct calculation. For each point within the cube, the corresponding charge density varies accordingly. Using the Runge-Kutta method, we numerically solve the coupled Eqs. (5), (), (9), and (12). This process yields the exciton’s final position $\overline{x}\left(t_f\right)$ and final velocity $\overline{v}\left(t_f\right)$ at $t=t_f$. Each point satisfied the conditions in Eq. (15) is labeled as valid, others are labeled as invalid.

Here Step 1 only requires relatively less time and resources to get a rough outline of the distribution of valid points. However, a more granular identification of valid points are essential for further analysis. Therefore, the following steps attempt to narrow the search range by only searching for potential points near the valid points under a smaller ${\delta }_c$.

• Step 2: Each valid point obtained from Step 1 is now treated as an initial point to construct a local neighborhood. Toward an initial point ${\mathcal{P}}_0(c_0,c_1,c_2)$, we define its neighborhood ${\mathcal{N}}_1$ by independently adjusting each coordinate $c_i$$(i=0,1,2)$ with a reduced step size ${\delta }_c/2$:

  $${\mathcal{N}}_1({\mathcal{P}}_0,{\delta }_c/2)=\{c_0+{\mathrm{\Delta }}_0{\delta }_c/2,c_1+{\mathrm{\Delta }}_1{\delta }_c/2,c_2+{\mathrm{\Delta }}_2{\delta }_c/2\},$$

  (17)

  where the multipliers ${\mathrm{\Delta }}_0$, ${\mathrm{\Delta }}_1$, and ${\mathrm{\Delta }}_2$ are chosen from the set $\{-1,0,1\}$, which corresponds to decreasing the coordinate by ${\delta }_c/2$, leaving it unchanged, or increasing it by ${\delta }_c/2$, respectively. This approach generates all possible combinations of such adjustments simultaneously for $c_0$, $c_1$, and $c_2$. This yields $3^D=27$ total combinations ($D=3$), one being ${\mathcal{P}}_0$ itself and the remaining $3^D-1$ points surrounding it;
• Step 3: Applying the same numerical method used in Step 1, we evaluate the points in each neighborhood to determine which point remains valid under the reduced step size ${\delta }_c/2$. By focusing on neighborhoods around previously identified valid points, we concentrate computational effort on the regions with potential to find the valid points;
• Step 4: Valid points from Step 3 serve as updated initial points. The step size is further reduced to ${\delta }_c/4$, and the same process is repeated along all $3^D-1$ directions. By iteratively narrowing the step size and concentrating on neighborhoods around valid points, we achieve a more precise and efficient identification of additional valid parameter combinations.

The described search method significantly reduces the calculation of the points and then improves computational efficiency. Intuitively, we define efficiency improvement as the ratio of the number of computational points reduced by the search method to the number of computational points of direct calculation.

Table 1 presents a comparison of the computational efficiency between the search method and direct calculation, utilizing the final conditions $x_{1,f}=8$, $x_{2,f}=8.5$, $v_{1,f}=1.1$, and $v_{2,f}=1.2$ for excitons in an external electric field. On the first round of searching, as mentioned in Step 3, the search method calculated 131 parameter combinations and identified 105 valid points, achieving $78.95\%$ of the valid points found by direct calculation. These 105 points served as the initial points for a subsequent search round, where 1,591 combinations were evaluated, resulting in 799 valid points. In contrast, direct calculation requires evaluating over two million points to identify all 989 valid points. Thus, the search method successfully identified $80.79\%$ of the total valid points while improving computational efficiency by approximately $98.27\%$.

Although the search method significantly enhances computational efficiency, it may overlook valid points, particularly those located further away from the initial valid points. This limitation arises from the neighborhood-based approach, which focuses on evaluating points within a fixed local region around each initial point. As illustrated in Figure 3 (a), green squares represent valid points identified through direct calculation with ${\delta }_c=0.32$, which served as the initial points for subsequent searches. However, the search method evaluates only 26 directions around each initial point, with a reduced step size of ${\delta }_c=0.16$. In the slice of the $(c_0,c_1)$ plane, the blue dots denote the valid points discovered by the first round of searching, while the red rhombus represents valid points found by direct calculation with ${\delta }_c=0.16$. The missing points are located in the lower right region of the slice, indicating areas the search method fails to reach when initialized with the green squares. Consequently, these omissions occur in subsequent search iterations, limiting the method’s overall effectiveness.

To address this limitation, it becomes necessary to expand the search range. To enhance the coverage, the neighborhood ${\mathcal{N}}_2$ is extended in all directions, allowing the search to encompass ${(D+2)}^3-1$ combinations with a reduced step size of ${\delta }_c/2$, similar to Eq. (17) as:

$${\mathcal{N}}_2({\mathcal{P}}_0,{\delta }_c/2)=\{c_0+{\mathrm{\Delta }}_0{\delta }_c/2,c_1+{\mathrm{\Delta }}_1{\delta }_c/2,c_2+{\mathrm{\Delta }}_2{\delta }_c/2\}.$$

(18)

In this definition, the multipliers ${\mathrm{\Delta }}_0$, ${\mathrm{\Delta }}_1$, and ${\mathrm{\Delta }}_2$ are chosen independently from the set $\{\pm 1,0,\pm 2\}$, indicating that each coordinate $c_i$ is adjusted by integer multiples of ${\delta }_c/2$. This expanded neighborhood increases the discretized search points in all directions within the cube, thereby mitigating the risk of missing valid points located further away from the initial points. Consequently, Steps 3 and 4 of the search method are adjusted to accommodate the expanded scope.

In Step 3, the points within the extended neighborhood ${\mathcal{N}}_2({\mathcal{P}}_0,{\delta }_c/2)$ are evaluated using the same numerical methods outlined in Step 1. This extension allows the identification of valid points located further from the initial points obtained in Step 1, thereby improving the overall coverage of the search space. In Step 4, the newly identified valid points from Step 3 are used as updated initial points. The step size is then reduced to ${\delta }_c/4$, and the search process is repeated within the further refined neighborhood ${\mathcal{N}}_2({\mathcal{P}}_0,{\delta }_c/4)$, ensuring a more thorough exploration of the parameter space.

3.3. Search Results

We maintain the initial and final conditions of excitons in an external electric field and apply the modified search method to identify valid points anew. As shown in Table 2, after the first round of searching (with a step size of ${\delta }_c/2$), the modified search method identified 129 valid points. In the second round of searching (with a step size of ${\delta }_c/4$), a total of 984 valid points were identified. Compared to the original search method, this refined neighborhood approach increased the proportion of valid points found in the two search rounds by $19.52\%$ and $18.70\%$, achieving $98.47\%$ and $99.49\%$ of the total valid points, respectively.

Figure 3 (b) illustrates the $(c_0,c_1)$ plane at $c_2=1.6$ obtained using the modified search method. While the original method, with its limited search range, omits many valid points, especially those further from the initial green squares, the expanded search successfully identifies these previously missed points, including those in the lower right region of the same slice.

This expanded neighborhood approach significantly enhances the effectiveness of the search method by addressing its primary limitation, ensuring a higher proportion of valid points are identified while maintaining computational efficiency. This modified search method effectively solves such optimization problems and can significantly conserve computational resources.

Figure 4 shows how the exciton’s position, velocity, dipole moment, and external electric field evolve over time under a stripe-like gate with modulated charge density. The blue and red solid lines correspond to two different sets of charge density parameters: $c_0=0.32$, $c_1=0.96$, $c_2=1.6$ and $c_0=0.32$, $c_1=1.92$, $c_2=0.48$, respectively. Figure 4 (a) and (b) present the evolution of $\overline{x}\left(t\right)$ and $\overline{v}\left(t\right)$, both displaying smooth increases. Figure 4 (c) plots the dipole moment $\langle d\left(t\right)\rangle$, which grows with the electric field initially and around the expected value of 1 when ${\omega }_0{t}_f\gg 1$ in the adiabatic regime. This value corresponds to the ground state $|{\psi }_{-}\rangle =\left(\right|g_0\rangle -|e_0\rangle )/\sqrt{2}$, where the diagonal energy terms are relatively small. The dipole moment oscillates around this value due to quantum effects. Figure 4 (d) illustrates the external electric field $\mathcal{E}\left(\overline{x}\left(t\right)\right)$, demonstrating the robustness of the optimized electric field in precisely controlling the exciton’s motion.

Next, we consider a relatively short $t_f=3$, where the adiabatic picture is not applicable. When $t_f=3$, for two selected sets of parameters $c_0,c_1,c_2$, the time evolution of the exciton’s position, velocity, dipole moment, and the external electric field is presented in Figure 5. The velocity $v\left(t\right)$ in Figure 5 (b) oscillates more strongly, indicating less smooth variation over shorter times. In Figure 5 (c), the dipole moment $\langle d\left(t\right)\rangle$ initially follows the electric field. However, due to the rapid variation of the field, which violates the adiabatic condition, $\langle d\left(t\right)\rangle$ exhibits significant oscillations and no longer strictly follows the field evolution, resulting in a noticeable deviation from the adiabatic evolution. Consequently, at the final time $t=t_f$, the dipole moment does not decrease in tandem with the diminishing electric field, starkly contrasting with the adiabatic evolution scenario.

Figure 6 compares the results of the original and modified search methods on the same $(c_0,c_2)$ plane at $c_1=6.4$. The original search method, shown in Figure 6 (a), is constrained by a limited search range and consequently identifies fewer valid points, particularly missing those further from the initial points. In contrast, the modified search method, as depicted in Figure 6 (b), significantly enhances coverage by expanding the search range, thereby successfully locating additional valid points, including those initially missed. As summarized in Table 3, the search method identifies 977 valid points in the first search round, corresponding to $99.19\%$ of the total 985 valid points. In the second search round, the method successfully identifies 8,178 valid points, accounting for $99.70\%$ of the total 8,203 valid points. This demonstrates the improved efficiency and accuracy of the modified search method in identifying valid parameter combinations within a reduced timeframe.

4. Conclusions and outlook

We proposed an algorithm of quantum control of the classical motion of excitons in two-dimensional semiconductors. This approach is based on the design of time- and position-dependent electric fields produced by a stripe-like gate. This electric field produces a time-dependent dipole moment resulting in a net force acting on the exciton and producing its classical acceleration. The ways to improve the efficiency of the search algorithm are presented and analyzed.

The proposed approach can be used for the control of exciton energy transfer processes and optics of semiconducting quantum wells by various local static and time-dependent electric fields. In general, it can be applied to the analysis of dynamics of two-level systems such as spins of single electrons and excitons in the presence of spin-orbit coupling and system magnetization [19,20].