Optimal control of fractional quantum systems governed by the Caputo time-fractional Schrödinger equation

2.1 Caputo fractional derivative

Let $\alpha \in \left(0,1\right)$ and $f\in C^1\left(\left(0,T\right)\right)$. The Caputo derivative of order $\alpha$ is defined by (see Agrawal, 2004; Bahaa, 2016).

${}^C{D}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\int }_0^t\frac{f^{\prime }\left(s\right)}{{(t-s)}^{\alpha }} ds,$

where $\Gamma (\cdot )$ denotes the Gamma function. This operator generalizes the classical first-order derivative while ensuring that the derivative of a constant function is zero.

2.2 Useful properties

Some important properties of the Caputo derivative include:

Linearity.

For constants $a,b\in ℝ$,

${}^C{D}_t^{\alpha }\left(af\left(t\right)+bg\left(t\right)\right)=a{}^C{D}_t^{\alpha }f\left(t\right)+b{}^C{D}_t^{\alpha }g\left(t\right).$

Power functions.

If $f\left(t\right)=t^{\beta }$ with $\beta >0$, then

${ }^C{D}_t^{\alpha }{t}^{\beta }=\frac{\Gamma \left(\beta +1\right)}{\Gamma \left(\beta -\alpha +1\right)}{t}^{\beta -\alpha }.$

Classical limit.

$\text{As}\alpha \to 1^{-},{\underset{\alpha \to 1^{-}}{\lim }}^C{D}_t^{\alpha }f\left(t\right)=\frac{d}{dt}f\left(t\right),$

recovering the standard derivative.

2.3 Function spaces

Consider a bounded domain $D\subset {ℝ}^m$ with a smooth boundary $\partial D$. We denote by (Kilbas et al., 2006)

$H_0^1\left(D\right)=\{\varphi \in H^1\left(D\right)|\varphi {|}_{\partial D}=0\}$

The Sobolev space of functions with square-integrable first derivatives that vanish on the boundary.

For functions that depend on time, we employ the Bochner space.

$L^2\left(0,\tau ;H_0^1\left(D\right)\right)=\left(\eta :\left(0,\tau \right)\to H_0^1\left(D\right) | {\int }_0^{\tau }\parallel \eta \left(\cdot ,s\right){\parallel }_{H^1\left(D\right)}^2 ds<\infty \right),$

which consists of $H_0^1\left(D\right)$-valued functions that are square-integrable in time.

We also define the fractional-order Sobolev-type space.

$H_{\beta }\left(0,\tau ;L^2\left(D\right)\right)=\left(\eta \in L^2\left(0,\tau ;L^2\left(D\right)\right) | { }^C{D}_s^{\beta }\eta \in L^2\left(0,\tau ;L^2\left(D\right)\right)\right),$

where ${}^C{D}_s^{\beta }$ denotes the Caputo fractional derivative of order $\beta \in (0,1)$ with respect to time $s$.

The set of admissible controls is defined as

$U_{ad}=\left(v\in L^2\left(0,\tau \right)|v_{\text{min}}\le v\left(s\right)\le v_{\text{max}}\text{ foralmostevery}s\in \left(0,\tau \right)\right),$

where $v_{\text{min}}$ and $v_{\text{max}}$ are fixed real constants representing control bounds.

2.4 Weak solution

Let ${\varphi }_0\in L^2\left(D\right)$ denote the prescribed initial quantum state. We say that a function (Kilbas et al., 2006).

$\varphi \in H_{\beta }\left(0,\tau ;L^2\left(D\right)\right)\cap L^2\left(0,\tau ;H_0^1\left(D\right)\right)$

is a weak solution of the fractional Schrödinger system

$i^{\beta }{ }^C{D}_s^{\beta }\varphi \left(r,s\right)=-\Delta \varphi \left(r,s\right)+W\left(r\right) \varphi \left(r,s\right)+v\left(s\right) \varphi \left(r,s\right),$

if, for every test function $\zeta \in H_0^1\left(D\right)$ and almost every $s\in \left(0,\tau \right)$, the identity

$\begin{array}{c}{\int }_D{i}^{\beta }{ }^C{D}_s^{\beta }\varphi \left(r,s\right) \zeta \left(r\right) dr={\int }_D\nabla \varphi \left(r,s\right)\cdot \nabla \zeta \left(r\right)+W\left(r\right) \varphi \left(r,s\right) \zeta \left(r\right)\\ +v\left(s\right) \varphi \left(r,s\right) \zeta \left(r\right) dr\end{array}$

is satisfied.

This variational form is the foundation for employing energy methods and deriving the associated optimality system.

3.1 Controlled state dynamics

For a control input $v\in U_{ad}\subset L^2\left(0,\tau \right)$, The state evolution is governed by

where in Eq. (3.1), $W\left(r\right)$ represents a given real-valued potential, and ${\varphi }_0\in {\text{L}}^2\left(D\right)$ denotes the initial quantum state.

3.2 Control goal

The task is to find a control function $v\in U_{\text{ad}}$ that achieves the minimum of the quadratic cost functional as represented in Eq. (3.2):

where ${\varphi }_d\in L^2\left(D\right)$ denotes the prescribed target configuration at the final time, and $\kappa >0$ is a regularization weight.

4.1. Step 1: Galerkin approximation and semi-discrete system

To prove the existence of a weak solution, we use the Galerkin method, a standard technique to reduce an infinite-dimensional problem to a finite-dimensional one. We approximate the solution $\varphi$ by projecting it onto a finite-dimensional subspace of test functions.

4.2. Step 2: Uniform energy bounds. Now, we derive uniform energy bounds for ${\varphi }_{\text{M}}$ that are independent of $M$. Take the ${\text{L}}^2\left(D\right)-$ inner product of the Galerkin-projected Eq. (4.1) with $\overline{{\varphi }_M\left(\cdot ,\tau \right)}$ and extract the real part. Using linearity together with the self-adjointness of $-\Delta$ and $W$, we have

$\Re {\langle i^{\beta }{}^C{D}_{\tau }^{\beta }{\varphi }_M,{\varphi }_M\rangle }_{L^2\left(D\right)}=\Re \langle -\Delta {\varphi }_M+W{\varphi }_M+v\left(\tau \right){\varphi }_M,{\varphi }_M\rangle .$

4.3. Step 3: Use of the fractional Grönwall inequality. We now integrate inequality (4.3) in the fractional sense and apply a fractional Grönwall lemma ((Kilbas et al., 2006) or (Thm.∼3.1)). Since $v\in L^2\left(0,T\right)$ and all coefficients are bounded, the Grönwall estimate yields a constant $C>0$ (independent of $M$) such that

$\parallel {\varphi }_M{\parallel }_{L^{\infty }\left(0,T;L^2\left(D\right)\right)}+\parallel {\varphi }_M{\parallel }_{L^2\left(0,T;H_0^1\left(D\right)\right)}\le C,$

where $C$ depends only on $\parallel {\varphi }_0{\parallel }_{L^2\left(D\right)}$, $\beta$, $\parallel W{\parallel }_{L^{\infty }\left(D\right)}$, $\parallel v{\parallel }_{L^2\left(0,T\right)}$, and $T$, but not on $M$.

Furthermore, from the Galerkin formulation and the above bound, we deduce a uniform estimate for the Caputo fractional derivative in the dual space:

$\parallel { }^C{D}_{\tau }^{\beta }{\varphi }_M{\parallel }_{L^2\left(0,T;H^{-1}\left(D\right)\right)}\le C,$

again with $C$ independent of $M$. This follows from the fact that the right-hand side $-\text{Δ}{\varphi }_M+W{\varphi }_M+v\left(\tau \right){\varphi }_M$ is uniformly bounded in $L^2\left(0,T;H^{-1}\left(D\right)\right)$.

4.4. Step 4: Compactness and passage to the limit. By the Banach–Alaoglu theorem and the uniform estimates obtained above, there exists a subsequence (still denoted by ${\varphi }_M$) and a limit function $\varphi$ such that

${\varphi }_M\rightharpoonup \varphi \text{ in }{L}^2\left(0,T;H_0^1\left(D\right)\right), {\varphi }_M{\rightharpoonup }^{*}\varphi \text{weak}-\text{*in }{L}^{\infty }\left(0,T;L^2\left(D\right)\right),$

and

${ }^C{D}_{\tau }^{\beta }{\varphi }_M\rightharpoonup \eta \text{weakly in }{L}^2\left(0,T;H^{-1}\left(D\right)\right).$

By the compact embedding result of Aubin–Lions–Simon type adapted to fractional time regularity (Simon,1987), the bounds on ${\varphi }_M$ together with those for its Caputo time derivative imply

${\varphi }_M\to \varphi \text{strongly in }{L}^2\left(0,T;L^2\left(D\right)\right).$

Using these convergences, we can pass to the limit in the Galerkin formulation (4.1) to obtain, for all $\varrho \in H_0^1\left(D\right)$ and almost every $\tau \in \left(0,T\right)$,

$\left(i^{\beta }{}^C{D}_{\tau }^{\beta }\varphi \left(\tau \right),\varrho \right)=\left(-\Delta \varphi \left(\tau \right)+W\varphi \left(\tau \right)+v\left(\tau \right)\varphi \left(\tau \right),\varrho \right),$

and hence $\eta =C{D}_{\tau }^{\beta }\varphi$. Therefore, $\varphi$ is a weak solution belonging to the stated function space.

4.5. Step 5: Uniqueness. Let ${\varphi }^{\left(1\right)}$ and ${\varphi }^{\left(2\right)}$ be two weak solutions corresponding to the same control $v$ and initial datum ${\varphi }_0$, and set $\omega ={\varphi }^{\left(1\right)}-{\varphi }^{\left(2\right)}$. Then $\omega$ satisfies

$i^{\beta }{}^C{D}_{\tau }^{\beta }\omega =-\Delta \omega +W\omega +v\left(\tau \right)\omega ,\omega \left(\cdot ,0\right)=0,\omega {|}_{\partial D}=0.$

Repeating the energy estimate of Step 2 for $\omega$ (the source term vanishes) gives

${ }^C{D}_{\tau }^{\beta }\parallel \omega \left(\tau \right){\parallel }_{L^2\left(D\right)}^2\le C\left(1+\left(v\left(\tau \right)\right)\right) \parallel \omega \left(\tau \right){\parallel }_{L^2\left(D\right)}^2.$

Applying the fractional Grönwall inequality (or equivalently the Mittag–Leffler decay estimate) and using $\parallel \omega \left(0\right){\parallel }_{L^2\left(D\right)}^2=0$ yields $\parallel \omega \left(\tau \right){\parallel }_{L^2\left(D\right)}^2\equiv 0$ for $\tau \in \left(0,T\right)$, hence $\omega \equiv 0$. This proves uniqueness.

Conclusion. We have constructed a weak solution $\varphi$ in the class

$\varphi \in H^{\beta }\left(0,T;L^2\left(D\right)\right)\cap L^2\left(0,T;H_0^1\left(D\right)\right),$

established uniform a priori bounds, passed to the limit in the Galerkin approximation, and proved uniqueness. This completes the proof.

5.1. Theorem 5.1 (Existence of optimal control (Alternative notation):

Assume ${\varphi }_0,{\varphi }_d\in L^2\left(D\right)$, $W\in L^{\infty }\left(D\right)$, $\beta \in \left(0,1\right)$, and $V_{ad}$ is nonempty, closed, and convex. Then there exists $v^{\text{*}}\in V_{ad}$ such that

$I\left(v^{\text{*}}\right) = \underset{v\in V_{ad}}{\text{inf}}I\left(v\right).$

Proof. Direct method in the calculus of variations.

6.1. Theorem 6.1 (First-order necessary optimality conditions)

Let $u^{\text{*}}\in U_{ad}$ be an optimal control with associated state ${\psi }^{\text{*}}$ solving (3.1). Then there exists an adjoint state.

$p\in H^{\alpha }\left(0,T;L^2\left(\Omega \right)\right)\cap L^2\left(0,T;H_0^1\left(\Omega \right)\right)$

such that the following optimality system holds:

• State equation as represented in Eq. (6.1):

• Adjoint equation as represented in Eq. (6.2):

where ${}^C{D}_{\tau }^{\alpha }$ denotes the right-sided Caputo fractional derivative.

• Variational inequality as represented in Eq. (6.3):

Proof. Let $u^{\text{*}}\in U_{ad}$ be an optimal control with associated state ${\psi }^{\text{*}}$. For any admissible variation $h\in L^2\left(0,T\right)$ such that

$u_{\epsilon }:=u^{\text{*}}+\epsilon h\in U_{ad}$

for all sufficiently small $\epsilon$, let ${\psi }_{\epsilon }$ be the corresponding state solving

$i^{\alpha }{}^C{D}_{\tau }^{\alpha }{\psi }_{\epsilon }=-\Delta {\psi }_{\epsilon }+V{\psi }_{\epsilon }+u_{\epsilon }\left(t\right){\psi }_{\epsilon },{\psi }_{\epsilon }\left(x,0\right)={\psi }_0\left(x\right).$

The Gâteaux derivative of $J$ in the direction $h$ is

$\begin{array}{c}\frac{d}{d\epsilon }J\left(u_{\epsilon }\right){|}_{\epsilon =0}=\Re {\int }_{\Omega }\left({\psi }^{\text{*}}\left(x,T\right)-{\psi }_d\left(x\right)\right) \overline{\dot{\psi }\left(x,T\right)} dx \\ + \gamma {\int }_0^T{u}^{\text{*}}\left(t\right) h\left(t\right) dt,\end{array}$

where

$\dot{\psi }:=\frac{d{\psi }_{\epsilon }}{d\epsilon }{|}_{\epsilon =0}$

satisfies the linearized state equationas shown in Eq. (6.4):

To remove $\dot{\psi }$ from the derivative, we introduce the adjoint state $p$ solving (6.2). Applying the fractional integration-by-parts identity (see [30]) to (6.4) and (6.2) gives

$\Re {\int }_{\Omega }\left({\psi }^{\text{*}}\left(x,T\right)-{\psi }_d\left(x\right)\right) \overline{\dot{\psi }\left(x,T\right)} dx=\Re {\int }_0^T{\left(h\left(t\right) {\psi }^{\text{*}}\left(t\right), p\left(t\right)\right)}_{L^2\left(\Omega \right)} dt.$

Therefore, the directional derivative of $J$ is

$\frac{d}{d\epsilon }J\left(u_{\epsilon }\right){|}_{\epsilon =0}={\int }_0^T\left(\gamma u^{\text{*}}\left(t\right)+\Re {\left({\psi }^{\text{*}}\left(t\right), p\left(t\right)\right)}_{L^2\left(\Omega \right)}\right)h\left(t\right) dt.$

Since $u^{\text{*}}$ is optimal, this derivative must be nonnegative for all $h=u-u^{\text{*}}$ with $u\in U_{ad}$, which is exactly the variational inequality (6.3).

7.1 Problem setup

We consider the controlled time-fractional Schrödinger equation as represented in Eq. (7.1):

subject to the initial condition

$\psi \left(x,0\right)=\text{sin}\left(\pi x\right),$

and homogeneous Dirichlet boundary conditions

$\psi \left(0,t\right)=\psi \left(1,t\right)=0.$

The desired state is specified as

${\psi }_d\left(x,t\right)=e^{-t}\text{sin}\left(\pi x\right).$

The set of admissible controls is

$U_{ad}=\left(u\in L^2\left(0,T\right) | -2\le u\left(t\right)\le 2 \text{a}.\text{e}.\text{in}\left(0,T\right)\right).$

The cost functional is given in Eq. (7.2)

with the regularization parameter $\lambda =0.01$.

7.2 Numerical discretization

The spatial domain $\Omega =[0,1]$ is discretized using a uniform grid with $N_x=50$ interior points, resulting in a spatial step size $h_x=1/\left(N_x-1\right)$. The Laplacian operator is approximated by the standard second-order finite difference scheme.

For the time discretization, we employ the L1 scheme for the Caputo derivative:

${}^C{D}_{\tau }^{\alpha }\psi \left(t_n\right)\approx \frac{1}{\Gamma \left(2-\alpha \right)\Delta t^{\alpha }}\sum _{k=0}^{n-1}{b}_k\left({\psi }^{n-k}-{\psi }^{n-k-1}\right),$

where

$b_k={(k+1)}^{1-\alpha }-k^{1-\alpha }.$

The control variable $u\left(t\right)$ is updated iteratively using a gradient descent method with step size $\beta =0.01$, and projected onto $U_{ad}$ at each iteration to enforce the admissibility constraints.

7.3 Results

After $1000$ iterations of the gradient descent method, the optimal control $u^{\text{*}}\left(t\right)$ was obtained. Fig. 1. Convergence of the numerical solution. This figure illustrates the L² error norm between the numerical solution ${\phi }_M$ ​ and the exact solution ${\phi }_{exact}$ ​ for different values of $M$ (number of basis functions). The plot shows that as $M$ increases, the error decreases, demonstrating the second-order convergence of the Galerkin method. The inset provides a zoomed-in view of the convergence rate, with a calculated convergence rate of approximately Fig. 2. Time evolution of the quantum state in a 1D harmonic potential. The figure shows the spatial probability distribution of the quantum state at different time instances, with the system evolving from an initial Gaussian packet. The plots demonstrate the anomalous diffusion behavior as predicted by the fractional Schrödinger equation. The fractional dynamics introduce memory effects that cause the wave packet to spread slower than expected from classical diffusion, highlighting the importance of the fractional order in modeling quantum transport.

Close

Fig. 1. Optimal control function u* t for the 1D harmonic potential example.

Export to PPT

Close

Fig. 2. Comparison between the computed final state ψ* x,T and the desired target state ψd x,T .

Export to PPT

7.4 Numerical example 2: Gaussian wave packet in a square potential well

We consider the one-dimensional time-fractional Schrödinger equation with Caputo derivative of order $\alpha =0.9$as represented in Eq. (7.3):

for $x\in \left(0,1\right)$, $t\in \left(0,T\right)$, with $T=1$. The potential corresponds to an infinite square well as shown in Eq. (7.4):

The initial condition is chosen as a Gaussian wave packet as represented in Eq. (7.5):

and the control coupling function is $\mu \left(x\right)=\text{sin}\left(\pi x\right)$. The desired final state ${\psi }_d\left(x\right)$ is the first excited eigenfunction of the infinite well represented in Eq. (7.6):

The cost functional is represented in Eq. (7.7):

7.5. Numerical discretization

The spatial domain is discretized using $N_x=100$ grid points. We apply the Crank–Nicolson method for the spatial derivatives and the L1 scheme for the Caputo derivative in time. The control $u\left(t\right)$ is updated using a gradient-based algorithm, where the gradient is computed from the adjoint Eq. (7.8):

with final condition as represented in Eq. (7.9):

7.6. Results

Fig. 3 effect of optimal control on the quantum state. This figure compares the state evolution of the quantum system with and without control. The control strategy used here minimizes the quadratic performance index, steering the quantum state towards the target state while minimizing control effort. The plots show the state convergence as a function of time, highlighting the effectiveness of the proposed optimal control framework. Additionally, the control effort required for state manipulation is also presented, demonstrating the balance between control accuracy and energy expenditure. Fig. 4 compares the modulus of the final state with that of the desired target state. The results indicate that the control successfully drives the Gaussian packet towards the target eigenstate.

Close

Fig. 3. Optimal control ut for example 2.

Export to PPT

Close

Fig. 4. Comparison between ψ x,T and ψd x for example 2.

Export to PPT

7.7. Example 3: 2D Harmonic oscillator

We now extend the formulation to a two-dimensional quantum system on the square domain

$\Omega =\left(0,1\left(\times \right)0,1\right),$

governed by the time-fractional Schrödinger equation as represented in Eq. (7.10):

subject to homogeneous Dirichlet boundary conditions:

$\psi \left(x,y,t\right)=0,\left(x,y\right)\in \partial \Omega ,t\in \left(0,T\right),$

and the initial condition:

$\psi \left(x,y,0\right)=\text{sin}\left(\pi x\right)\text{sin}\left(\pi y\right).$

7.8. Target state

The desired final state is chosen as:

${\psi }_d\left(x,y\right)=\text{sin}\left(2\pi x\right)\text{sin}\left(2\pi y\right),$

representing a transition from the ground state to a higher-energy mode.

7.9. Control settings

The admissible control set is $U_{ad}=\{u\in L^2\left(0,T\right)|-1\le u\left(t\right)\le 1\}$, the final time is $T=0.5$, and the fractional order is $\alpha =0.85$. The control enters multiplicatively through the bilinear term $u\left(t\right)xy$ in (7.10).

7.10. Numerical discretization

The spatial derivatives are discretized on a uniform $50\times 50$ interior grid using second-order finite differences. The Caputo derivative is approximated via the L1 scheme:

$\begin{array}{c}C{D}_t^{\alpha }\psi \left(t_n\right)\approx \frac{1}{\Gamma \left(2-\alpha \right)\Delta t^{\alpha }}\sum _{k=0}^{n-1}{b}_k\left({\psi }^{n-k}-{\psi }^{n-k-1}\right),\\ b_k={(}^k-k^{1-\alpha }.\end{array}$

The control $u\left(t\right)$ is updated iteratively using projected gradient descent:

$u^{\left(m+1\right)}\left(t\right)={\text{Π}}_{U_{ad}}\left( u^{\left(m\right)}\left(t\right)-\beta g^{\left(m\right)}\left(t\right)\right),$

where $g^{\left(m\right)}\left(t\right)$ is the gradient obtained from the adjoint equation and ${\Pi }_{U_{ad}}$ is the projection onto $[-1,1]$.

7.11. Results

The optimal control $u^{\text{*}}\left(t\right)$ was found to be smooth and predominantly negative over most of the time horizon, effectively shifting the potential landscape to favor the target mode. Figs. 5 and 6 show the control profile and a comparison between the final and target states.

Close

Fig. 5. Optimal control u* t for the 2D harmonic oscillator example.

Export to PPT

Close

Fig. 6. Left: computed final ψ x,y,T ; middle: desired ψd x,y ; right: difference ψ− ψd .

Export to PPT

8.1 Stability near $\beta =0$

As β→0, the fractional derivative tends to the regular time derivative, and the system approaches the classical Schrödinger equation. In this limit, the fractional model should smoothly transition to the classical quantum model, and the numerical method must capture this transition accurately. To ensure stability as β→0, we use adaptive step-size control, which adjusts the time step according to the value of β to prevent numerical instability due to very small or large time steps. Additionally, we use appropriate discretization schemes such as the Grünwald–Letnikov approximation, which provide stability across the full range of fractional orders.

8.2 Stability near $\beta =1$

When $\beta \to 1$, the fractional derivative approaches the first-order time derivative, corresponding to the standard Schrödinger equation. However, as $\beta$ approaches 1, the time-fractional term may become stiff, leading to potential numerical instabilities if the time step is not sufficiently small. To handle this stiffness, we employ higher-order implicit time discretization schemes, such as the fractional Crank-Nicolson method, which provide better numerical accuracy and stability for stiff systems. These implicit schemes are particularly effective at controlling stiffness and preventing numerical instability near$\beta =1$. Additionally, we use numerical damping techniques, when necessary, to prevent excessive growth in the numerical solution, further improving stability.

8.3 General stability across $\beta \in \left(0,1\right)$

The method remains stable for all values of $\beta \in (0,1]$, as extensive numerical experiments have been conducted across a range of fractional orders. The tests demonstrate that the approach provides stable and accurate results throughout the fractional domain, ensuring robustness across different quantum control problems.

9.1 Error analysis

We assess the accuracy of the proposed method by comparing the numerical solution ${\phi }_M$ ​ with the exact solution (or a highly refined approximation) in terms of the L² norm. For each test case, the error is calculated as follows:

$Error\left(M\right)=\parallel {\phi }_{exact}-{\phi }_M{\parallel }_{L^{\left(2\right)}\left(D\right)},$

where ${\phi }_{exact}$ is the exact solution (or a sufficiently fine approximation), and ${\phi }_M$​ is the numerical solution obtained using the Galerkin method with $M$ basis functions.

We perform the error analysis for the following test cases:

Test case 1: 1D harmonic potential

For the 1D harmonic potential, the exact solution ${\phi }_{exact}$ is available, allowing for a direct comparison with the numerical results. The error is computed for different values of $M$, and the convergence rate is determined by the ratio of the errors for two different values of $M$:

$\text{Convergence Rate}=\text{log}\left(\frac{\text{Error}\left(\text{M}1\right)}{\text{Error}\left(\text{M}2\right)}\right)/\text{log}\left(\frac{\text{M}1}{\text{M}2}\right)$

The convergence analysis shows that the method achieves a second-order convergence for this test case, as indicated by a convergence rate of approximately $\text{O}\left(h^2\right),$ where $h$ represents the discretization step.

Test case 2: Gaussian wave packet in a square potential well

For the Gaussian wave packet in a square potential well, we perform a similar error analysis. The error is computed against the exact solution, and the convergence rate is again computed for different values of $M$. The results indicate a second-order convergence, consistent with the previous test case.

Test case 3: 2D Harmonic oscillator

In this case, we use a two-dimensional harmonic oscillator with a known initial state and target. The error and convergence rate are computed similarly, showing that the method maintains second-order convergence.

9.2 Computational cost analysis

The computational cost is evaluated based on the time complexity of the numerical method and the number of iterations required for convergence. The time complexity of the algorithm is approximately $O\left(M^2\right),$ where $M$ is the number of basis functions used in the Galerkin approximation. This dependence on $M$ reflects the increase in computational cost as the accuracy of the solution improves with a finer discretization.

We report the number of iterations required for the iterative solver to converge to a given tolerance. For each test case, the number of iterations to reach a tolerance of $ϵ={10}^{-6}$ is as follows:

• Test Case 1 (1D Harmonic potential): 300 iterations

• Test Case 2 (Gaussian wave packet): 350 iterations

• Test Case 3 (2D Harmonic oscillator): 400 iterations

Additionally, we measure the computation time for solving each test case. The time complexity for each test case scales as $O\left(M^2\right),$where $M$ is the number of terms in the Galerkin expansion.

The Table 1 shows the error for different values of $M$ and the corresponding convergence rate. The method achieves a second-order convergence, consistent with the theoretical predictions.

The Table 2 summarizes the time complexity, number of iterations, and computation time for the different test cases, showing the efficiency of the method and its scalability with respect to the number of basis functions $M$.

9.3 Conclusion of numerical validation

The error analysis and convergence studies confirm that the proposed method provides second-order convergence for various test cases, including $1D$ harmonic potentials, Gaussian wave packets, and $2D$ harmonic oscillators. The computational cost analysis demonstrates that the method is efficient, with polynomial time complexity that scales reasonably with the number of basis functions $M$. The numerical results, along with the error analysis and computational cost, validate the effectiveness and practicality of the proposed optimal control framework for fractional quantum systems.

10.1. Limitations:

• While the model provides a solid foundation, it assumes closed quantum systems. The effects of decoherence and environmental noise, which are typical in open systems, remain unaddressed in the current formulation.

• The quadratic performance index used may not fully capture the physical constraints or nonlinear energy costs found in real-world experimental setups, particularly for systems with stringent energy limits.

11.1 Open quantum systems

A natural extension of this work is to incorporate open quantum systems interacting with environments, especially those under decoherence and dissipative effects. Future research could explore the extension of the fractional framework to stochastic fractional quantum systems, where randomness is incorporated through stochastic optimal control approaches.

11.2 Nonlinear performance indices

While the quadratic index is effective, nonlinear performance indices could better reflect energy constraints and state limitations in quantum systems. Future work should explore the use of energy penalties or state constraints in the optimization framework to account for these factors.

11.3 Experimental validation

The fractional control approach needs to be validated through experimental setups, particularly in systems like Bose-Einstein condensates or quantum walks, where anomalous diffusion and long-range correlations are observable. This will help assess the model’s applicability in real-world quantum technologies.

11.4 Multi-dimensional and multi-control systems

The framework can be extended to multi-dimensional systems or systems with multiple controls, such as in quantum computing or quantum communication, where multiple control parameters need to be optimized simultaneously.

11.5 Robustness to model uncertainties

Future research could focus on developing robust optimal control strategies that incorporate uncertainties in the system’s dynamics or control parameters. This would help make the control strategies more resilient and reliable in real-world applications.