Abstract

The quantum rotor in external fields such as static electric, magnetic, and time-dependent fields is one of the most studied research problems due to its applicability in chemical physics, molecular dynamics, and many other related areas. The problem finds applications in quantum chaos and nonlinear dynamics as well. In this work, we present a Machine learning (ML) approach to model a one-dimensional quantum rotor in an external electric field. Our approach is to develop a surrogate model, keeping in mind the inherent smoothness of the eigenenergies of the dressed rotor and the directional parameters such as orientation, i.e., expectation value of $\cos \varphi$ , and the alignment parameter ($\langle {\cos }^2\varphi \rangle$ ) with the applied static field. The model we developed is very efficient, achieving microsecond-scale predictions with high accuracy. The model we developed is as follows: first, we train it on the basis of data generated from numerical calculations taking rotor-static field interaction into account and expanding the wavefunction as a finite basis. The model gives an excellent performance with R² values exceeding 0.9999 for all predicted quantities. We further explore the physics of the system and the data generation, the details of the ML model, and its uses in designing an inverse model. The work shows that ML can be used to speed up computational time in related areas such as molecular physics and may be used to optimize the experimental parameters. This approach is directly extensible to three-dimensional rotors, time-dependent fields, and multi-parameter optimization, with potential applications in quantum control, molecular alignment, and spectroscopic fitting.

Keywords

Machine Learning, Quantum Rotor, Orientation, Electric Field

Share and Cite:

1. Introduction

The study of polar molecules (2-D quantum rotor) and quantum rings (1-D) quantum rotor has been a subject of immense interest for many branches of science, including molecular dynamics, chemical physics, molecular spectroscopy, quantum chaos, nonlinear dynamics and many other related areas, as these models provide insights into quantized energy levels, angular momentum, and the behavior of matter at the nanoscale [1]-[3]. In addition, these quantum rotors, subjected to external perturbations such as electric fields and magnetic fields, have also been at the center of many interesting studies, such as optimal control of molecular motions, development of low-temperature physics applications, etc. [4]-[14].

As detailed recently by [9], the quantum rotor interacting with fields is a simple but significantly rich quantum system that is computationally solvable and provides all ingredients to understand complex physical phenomena [15] [16]. As mentioned, these systems are useful for exploring the physics of internal rotations of molecules, quantum rings at the nanoscale, one-dimensional systems, and a particle moving on a sphere; if extended to higher dimensions, the particle on a hyper-sphere, a model particularly useful for solving nuclear models [17] [18], and other fields.

Let us consider the simple case of a polar molecule with a permanent electric dipole moment $\mu$ placed in a static, uniform electric field $F$ . The molecule experiences a torque that tends to align its dipole with the field direction. However, quantum mechanics imposes restrictions on this alignment through the quantization of angular momentum. The resulting energy level structure and orientation properties are determined by the competition between the rotational kinetic energy and the field-induced potential energy [19].

The one-dimensional rotor simplifies the full three-dimensional problem, but it retains the essential physics of field-induced orientation and provides an excellent platform for developing and testing machine learning methodologies [20]. The restriction to planar rotation is physically relevant for molecules that are constrained to move on surfaces or in anisotropic environments, and it serves as a stepping stone for more complex systems [21]. The quantum kicked rotor problem has been used as a benchmark model for understanding quantum chaos, dynamical localization, quantum-classical correspondence, and many other interesting physical studies [22]-[25].

The most common approach to studying these systems involves solving the time-independent Schrödinger equation for each field strength of interest [26]. For a single calculation, this is computationally inexpensive and can be solved easily. However, applications such as optimal control design, spectroscopic fitting, or parameter estimation often require thousands or millions of such calculations. This is where machine learning (ML) offers a potential alternative: by learning the mapping from field strength to quantum observables, we can create surrogate models that provide instantaneous predictions, enabling real-time experimental feedback and large-scale parameter sweeps [27]-[29].

While the quantum rotor in an electric field is a textbook system with well-understood physics [4]-[6], the contribution of this work is not new physical phenomena but rather a methodological advance in computational practice. We demonstrate that a simple polynomial surrogate model achieves microsecond-scale predictions with <10⁻⁶ error in eigenenergies, representing a 2500 - 5000× speedup over direct diagonalization. This enables real-time experimental control and large-scale parameter sweeps that would otherwise be computationally prohibitive. To our knowledge, this is the first systematic demonstration of polynomial regression as a high-fidelity surrogate for rotor observables, complete with uncertainty quantification and inverse design capabilities. The complete MATLAB implementation is provided to facilitate immediate adoption and extension by the community.

Surrogate Modeling in Quantum Mechanics

Recent advances in machine learning have opened new avenues for accelerating quantum mechanical calculations [27]-[29]. Surrogate models—approximations trained on data from high-fidelity simulations—have been successfully applied to electronic structure calculations, molecular dynamics, and quantum control problems [30]. These models learn the mapping from input parameters (e.g., field strength, molecular geometry) to output observables (e.g., energies, transition probabilities), providing predictions that are orders of magnitude faster than the original simulation.

While deep neural networks have received considerable attention, simpler approaches such as Gaussian process regression and polynomial regression remain competitive for problems with low-dimensional inputs and smooth output functions [28]. For the quantum rotor in a static electric field, the target observables are analytic functions of the field strength (the Mathieu characteristic values are entire functions), making polynomial regression a natural and theoretically justified choice.

This manuscript presents a theoretical study for constructing such a surrogate model. In the next section, we begin with a detailed explanation of the quantum mechanical problem, including the derivation of the Hamiltonian, the choice of basis set, and the computation of relevant observables. Subsequently, we describe the data generation process and the ML approach, comparing neural network and polynomial regression methods. We then analyze the model performance in terms of accuracy, speed, and generalization capability, with a detailed discussion of the results presented in Figure 1. Finally, Section 5 discuss practical applications, particularly the inverse design problem, where one seeks the field strength required to achieve a target energy or orientation.

2. Theoretical Background

2.1. The Quantum Rotor Hamiltonian

In this work, we focus on the one-dimensional quantum rotor (QR), and the model can be easily extended to higher dimensions. As is well known, the motion of QR in one dimension is represented by the rotational angle $\varphi$ [15]. In the absence of external fields, the system is described by the kinetic energy operator:

${\widehat{H}}_0=-\frac{{\hslash }^2}{2I}\frac{{\text{d}}^2}{\text{d}{\varphi }^2}$ (1)

where $I$ is the moment of inertia. The eigenfunctions are familiar to free rotor states:

${\psi }_m\left(\varphi \right)=\frac{1}{\sqrt{2\pi }}{\text{e}}^{im\varphi },m=0,\pm 1,\pm 2,\cdots$ (2)

with corresponding eigenvalues:

$E_m=\frac{{\hslash }^2{m}^2}{2I}=B{m}^2$ (3)

Here, $B={\hslash }^2/\left(2I\right)$ is the rotational constant. These states are doubly degenerate for $m\ne 0$ , reflecting the equivalence of clockwise and counterclockwise rotation [16].

When a static electric field $F$ is applied along the $x$ -axis, the permanent dipole moment $\mu$ of the molecule interacts with the field through the Stark interaction:

${\widehat{H}}_{\text{Stark}}=-\mu F\cos \varphi$ (4)

The full Hamiltonian becomes:

$\widehat{H}=-\frac{{\hslash }^2}{2I}\frac{{\text{d}}^2}{\text{d}{\varphi }^2}-\mu F\cos \varphi$ (5)

It is convenient to work in dimensionless units by setting $\hslash =1$ and $2I=1$ , which corresponds to measuring energy in units of $B$ . Defining the dimensionless field parameter $\lambda =\mu F/B$ , the Hamiltonian simplifies to:

$\widehat{H}=-\frac{{\text{d}}^2}{\text{d}{\varphi }^2}-\lambda \cos \varphi$ (6)

This is the quantum mechanical analog of the classical pendulum, with $\lambda$ representing the dimensionless field strength.

2.2. The Mathieu Equation

The time-independent Schrödinger equation corresponding to this Hamiltonian is:

$-\frac{{\text{d}}^2\psi }{\text{d}{\varphi }^2}-\lambda \cos \varphi \psi =E\psi$ (7)

Rearranging yields:

$\frac{{\text{d}}^2\psi }{\text{d}{\varphi }^2}+\left(E+\lambda \cos \varphi \right)\psi =0$ (8)

This is the canonical form of the Mathieu equation, one of the fundamental special functions of mathematical physics [31] [32]. The solutions are the Mathieu functions, which are periodic functions of $\varphi$ . For a given $\lambda$ , the allowed energies $E_n\left(\lambda \right)$ are the Mathieu characteristic values, and the eigenfunctions are the corresponding Mathieu functions.

2.3. Physical Interpretation

The parameter $\lambda$ characterizes the relative strength of the field interaction compared to the rotational kinetic energy. Several regimes are of interest [5] [6].

In the free rotor regime where $\lambda =0$ , the rotor is free. The energy levels are given by $E_m=m^2$ , with degeneracy for $\pm m$ . The ground state with $m=0$ has zero energy, while the first excited states with $m=\pm 1$ have energy 1, and so forth. The orientation expectation value $\langle \cos \varphi \rangle$ is zero for all stationary states because the wavefunctions are delocalized and symmetric.

In the weak field regime where $0<\lambda \ll 1$ , the field lifts the degeneracy between the $\pm m$ states and introduces a small alignment. The energy levels change with $\lambda$ ; one peculiar feature is that the ground state energy keeps on decreasing, and the rotor starts aligning with the increasing field strength, whereas the excited state energies decrease (i.e., align) or increase (anti-align), and the degeneracy between different excited states is lifted [4].

It is interesting to note that when the field strength remains moderate, the electric field interaction energy approaches the energy of the rotor, and the energy level behavior deviates a bit from the weak field structure. The wavefunctions become localized along the direction of the field [19] [33].

And finally, in the strong field regime, the rotor gets strongly aligned with the field. The Hamiltonian approximates a harmonic oscillator in the angular coordinate, with the potential approximated as $\lambda \cos \varphi \approx \lambda \left(1-{\varphi }^2/2\right)$ for small $\varphi$ . The energy levels approach those of a harmonic oscillator: $E_n\approx -\lambda +\sqrt{\lambda }\left(n+1/2\right)$ . The orientation $\langle \cos \varphi \rangle$ approaches 1 for all low-lying states as the molecule is strongly confined near $\varphi =0$ .

2.4. Observables of Interest

For experimental applications and quantum control, two quantities, the energy and order parameter, are of primary importance [7] [20]. The eigenenergies $E_n\left(\lambda \right)$ determine the resonant frequencies for transitions between rotational states, which are essential for spectroscopy and for designing control pulses. The energy level structure also determines the thermal population distribution through the Boltzmann factor. The orientation expectation value ${\langle \cos \varphi \rangle }_n\left(\lambda \right)$ measures the degree of alignment of the molecule with the field direction. For a given eigenstate, $\langle \cos \varphi \rangle$ ranges from −1 for perfect anti-alignment through 0 for no net alignment to +1 for perfect alignment. The orientation is directly measurable in experiments through techniques such as laser-induced fluorescence or ion imaging, and it determines the interaction strength with other external fields [7]. The combination of these observables provides a complete characterization of the quantum state for control purposes. In particular, the ability to predict both quantities quickly and accurately enables closed-loop control experiments where the field is adjusted in real time to achieve a desired target state.

3. Computational Methodology

3.1. Basis Set Expansion

To solve the Schrödinger equation numerically, we employ a basis set expansion in the free rotor states [16]. This approach is natural because the free rotor states form a complete orthonormal basis for the Hilbert space of periodic functions on the unit circle. The Hamiltonian matrix elements in this basis are:

$\langle m|\widehat{H}|n\rangle =m^2{\delta }_{m,n}-\lambda \langle m|\cos \varphi |n\rangle$ (9)

with $\langle \cos \varphi \rangle$ elements evaluated as:

$\langle m|\cos \varphi |n\rangle =\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\text{e}}^{-im\varphi }\cos \varphi {\text{e}}^{in\varphi }\text{d}\varphi }=\frac{1}{2}\left({\delta }_{m,n+1}+{\delta }_{m,n-1}\right)$ (10)

Therefore, the Hamiltonian matrix is tridiagonal:

$H_{m,n}=m^2{\delta }_{m,n}-\frac{\lambda }{2}\left({\delta }_{m,n+1}+{\delta }_{m,n-1}\right)$ (11)

The basis is truncated with a sufficiently large number of states, so that the results for desired quantities converge, in our case, taking $M_{\text{max}}=15$ provides good convergence, with the ground state energy changing by less than 10⁻⁶ at the highest field, confirming the adequacy of this truncation.

Convergence Analysis

To verify that $M_{\text{max}}=15$ provides sufficient accuracy for all reported quantities, we performed convergence tests by increasing the basis size to $M_{\text{max}}=20$ and $M_{\text{max}}=30$ . Table 1 shows the maximum absolute change (Δ) (over $\lambda \in \left[0,10\right]$ ) when increasing the basis size:

Table 1. Convergence of eigenenergies and orientation with basis size.

The changes between $M=20$ and $M=30$ are two orders of magnitude smaller than our target tolerance (10⁻⁶ for energies, 10⁻⁵ for orientation). The choice $M_{\text{max}}=15$ , therefore, introduces errors that are negligible compared to the polynomial regression errors reported in Section 4.2. The ground state energy convergence criterion change (Δ) < 10⁻⁶ is used.

3.2. Eigenvalue Problem Solution

In order to solve the problem, the Hamiltonian matrix of order of $2N+1$ is constructed. As the matrix is real and symmetric, it can easily be diagonalized in MATLAB. As discussed earlier, the matrix of the order considered here, the eigen energies of interest converge up to the seventh decimal place, which, for our purpose, is acceptable.

The expectation value of $\cos \varphi$ for a given eigenstate $\psi$ is computed as:

$\langle \cos \varphi \rangle ={\displaystyle \sum }_m\frac{1}{2}\left(c_m^{\text{*}}{c}_{m+1}+c_m^{\text{*}}{c}_{m-1}\right)$ (12)

In the above equation, $c_m$ are the coefficients of a free rotor Hamiltonian.

3.3. Data Generation

The training dataset for further calculations is generated by solving the Hamiltonian matrix in the presence of a static field for 100 evenly placed values of field strength in the range 0 to 10. The range considered includes weak, moderate, and strong field strengths. Next, for each value of field strength, we evaluate the first five eigenvalues and calculate corresponding values of ${\langle \cos \varphi \rangle }_0$ through ${\langle \cos \varphi \rangle }_4$ . So, our training dataset comprises 100 samples with output values for eigenenergies and orientation parameters.

Although the data set is small, it is sufficient, since the target functions are well-behaved and quite smooth, allowing accurate interpolation with relatively few training points. We use single-precision floating-point arithmetic for the data storage to reduce memory consumption while maintaining adequate precision for machine learning purposes. The energies are of order unity, and the orientation values range from −1 to 1, so single precision provides approximately 7 decimal digits of precision—more than sufficient for our purposes.

Training and Test Set Partition

The full dataset consists of 100 field strength values ${\lambda }_i$ evenly spaced on the interval $\left[0,10\right]$ with step size $\text{Δ}\lambda =10/99\approx 0.10101$ . To enable unbiased evaluation of generalization performance, we partitioned the data as follows:

Training Set (80 Samples): All ${\lambda }_i$ where the index $i$ modulo 5 is not equal to 0 (i.e., $\lambda =0,0.10101,0.20202,0.30303,0.50505,\cdots$ , skipping every fifth point). This ensures that test points are evenly interspersed throughout the entire range rather than clustered at the ends.

Test Set (20 Samples): The remaining 20 points where index $i$ modulo 5 equals 0 (i.e., $\lambda =0.40404,0.80808,1.21212,\cdots ,9.69697$ ). These points constitute a held-out set that was never used for training, model selection, or any form of hyperparameter tuning.

Model Selection: For polynomial regression, the degree $d=6$ was chosen based on 5-fold cross-validation on the training set only. For each candidate degree $d=1$ to 10, we computed the cross-validation mean squared error (MSE) averaged over the 5 folds and selected the degree that minimized this quantity. Degrees 6 and 7 produced nearly identical CV scores (difference <0.1%), and we chose $d=6$ to favor simplicity.

Neural Network Reproducibility: All neural network training used a fixed random seed (seed = 42) to ensure reproducibility. The same train/test partition was used for both models.

The reported $R^2$ values (Section 4.2) are computed on the held-out test set only, providing an unbiased estimate of generalization performance. The near-perfect $R^2$ (0.9999) on unseen data confirms that the model has learned the underlying functional relationship, not merely memorized the training points.

3.4. Polynomial Regression Model

Given the smoothness of the target functions, we employ polynomial regression as our machine learning model [28]. For each output quantity $y\left(\lambda \right)$ , we fit a polynomial of degree $d$ :

$y\left(\lambda \right)={\displaystyle \sum }_{k=0}^d{a}_k{\lambda }^k$ (13)

The coefficients $a_k$ are determined by solving the linear least squares problem using the training data. This approach has several advantages. It offers simplicity, as the model is straightforward to implement and interpret. It provides stability, as with appropriate regularization, polynomial fits are well-behaved. It enables speed, as training reduces to solving a small linear system. It requires no hyperparameter tuning aside from the polynomial degree. Subsequently, the polynomial equation will be solved for the inverse problems.

The polynomial of degree 6 is chosen, since the physical behavior of the system from our previous calculation shows that the variation of eigenenergies and orientation parameter may be an analytical function of the field strength, with a radius of convergence in the complex field strength plane, determined by the nearest singularity. Since the eigenvalue problem of a one-dimensional quantum rotor in an electric field may be written in the form of a Mathieu equation, and the characteristic values of the matrix thus resulting are entire functions of the field strength, implying that they have convergent power series expansions for the field strength [32]. However, a finite degree polynomial provides an excellent approximation for the same. We have found that lower-order polynomials prove insufficient to give a clearer picture, and higher-order polynomials hardly make any further modifications. We found that the sixth-order polynomial keeps a balance between accuracy and generalization.

3.4.1. Comparison with Alternative Surrogate Frameworks

While several more sophisticated surrogate modeling approaches exist, we justify our choice of polynomial regression on both theoretical and practical grounds:

Gaussian Process Regression (GPR): GPR would provide built-in uncertainty quantification and flexible kernel choices. However, for analytic functions with low intrinsic dimensionality (here, a single scalar $\lambda$ ), GPR offers little accuracy advantage over polynomials while incurring $\mathcal{O}\left(N^3\right)$ training cost and $\mathcal{O}\left(N\right)$ prediction cost [28]. For our dataset size ($N=100$ ), GPR training would take ~1 second (vs. 0.05 seconds for polynomial regression), and prediction would take ~10 microseconds (vs. 1 - 2 microseconds).

Rational Approximants (Padé): Padé approximants can capture pole singularities and sometimes outperform polynomials for extrapolation. However, the Mathieu characteristic values are entire functions (no poles in the finite complex plane), so polynomials are theoretically optimal for interpolation within the radius of convergence. Padé offers no advantage for this specific problem.

Physics-Informed Neural Networks (PINNs): PINNs incorporate the Schrödinger equation as a soft constraint [30], potentially enabling extrapolation beyond training data. However, PINNs require solving a PDE-constrained optimization problem, with significantly higher training cost (hours vs. milliseconds) and architectural complexity. For our interpolation task within the training range, the additional complexity is unwarranted.

Spectral Methods (Chebyshev Polynomials): Chebyshev polynomials provide optimal interpolation nodes that minimize the Runge phenomenon. Our evenly spaced $\lambda$ grid is suboptimal in this regard; future work could adopt Chebyshev nodes. Nonetheless, with degree $d=6$ and $\lambda \in \left[0,10\right]$ , the Runge phenomenon is negligible, as evidenced by the absence of oscillatory artifacts in Figure 1.

We, therefore, conclude that polynomial regression strikes an optimal balance for this problem: it is theoretically justified, computationally efficient, and provides interpretable coefficients.

3.4.2. Polynomial Degree Selection

The polynomial degree $d$ was selected using 5-fold cross-validation on the training set. For each candidate degree $d=1,2,\cdots ,10$ , we computed the cross-validated mean squared error (CV-MSE) averaged over all five output quantities (five energies and five orientation values). The results are shown in Table 2.

The CV-MSE decreases rapidly from $d=1$ to $d=5$ , reaches a minimum at $d=6$ , and remains stable through $d=8$ before increasing at higher degrees due to overfitting (the training MSE continues to decrease, but the CV-MSE worsens). We selected $d=6$ as it achieves the optimal trade-off: it is the lowest degree within 1% of the minimum CV-MSE, reducing the risk of Runge phenomenon oscillations near the interval boundaries.

Table 2. Cross-validation MSE as a function of polynomial degree.

3.5. Neural Network Architecture Selection

It is worth mentioning that, for comparison purposes, the feedforward neural network using ReLU activation functions [34] [35] with two hidden layers with 32 neurons is also implemented. Here, the network is trained using the Adam optimizer with a mean squared error loss function. The input in our case is a single scalar, i.e., field strength, and the output is a 10-dimensional vector consisting of the lowest five eigenenergies and corresponding orientation parameters. 200 epochs of training with a batch size of 32, and we have used a learning rate of 0.001. This neural network provides quite a good accuracy to the polynomial regression, but it requires more training time and parameter tuning. Owing to the small dataset, polynomial regression is quite a preferred approach.

The choice of architecture (two hidden layers with 32 neurons each) was determined by a systematic hyperparameter search using 5-fold cross-validation on the training set. We evaluated layer counts (1 - 4), neurons per layer (16, 32, 64, 128), activation functions (ReLU, tanh, sigmoid), and learning rates (0.0001 - 0.01). The configuration with two hidden layers of 32 neurons and ReLU activations minimized validation MSE while avoiding overfitting (training vs. validation loss ratio = 1.03). Deeper networks (≥3 hidden layers) showed marginal accuracy improvement (<5% reduction in validation MSE) but required 3× more training time and exhibited greater variance across random seeds. Larger hidden layers (≥64 neurons) led to overfitting (training loss 10× smaller than validation loss). The Adam optimizer with learning rate 0.001 and batch size 32 was chosen as it consistently outperformed SGD and RMSprop in early experiments.

4. Results and Discussion

4.1. Quantum Mechanical Data and Analysis

In Figure 1, we present the variation of eigenenergies (left panel) and the orientation parameter (right panel), as a function of field strength. It shows how the external field modifies the energy levels’ structure. It is interesting to note that the degeneracy of excited states is removed even at weak field strengths, and the energy difference between otherwise degenerate states increases with the increase of the field strength. Also, the orientation parameter for otherwise non-allowed transitions is showing variations with the field strength.

Figure 1. Surrogate model predictions based on a neural network method for the one-dimensional quantum rotor. Left panel: eigenenergies $E_n\left(\lambda \right)$ for the lowest five rotational (i.e., n = 0 - 4) states variation with field strength $\lambda =\mu F/B$ . The training data points (scattered markers) and the smooth neural network predictions (solid lines) are shown, with different colors representing different quantum states. Right panel: Orientation parameter values ${\langle \cos \varphi \rangle }_n\left(\lambda \right)$ for the five states, showing the field-induced orientation and anti-orientation pattern. The markers are training data points, while the curves represent the predictions. The ground state shows increasing orientation with field strength, while excited states exhibit more complex behavior, including anti-orientation for some states.

State Labeling Convention

Throughout this work, we label states by their energy ordering at zero field: $n=0,1,2,3,4$ correspond to the five lowest eigenstates as a function of field strength $\lambda$ , with $n=0$ always the ground state. This adiabatic labeling preserves state identity along the entire $\lambda$ axis.

At $\lambda =0$ , the free rotor states have energies $E_m=m^2$ for $m=0,\pm 1,\pm 2,\pm 3,\pm 4,\cdots$ (Equation (3)). The degeneracy for $m\ne 0$ reflects the equivalence of clockwise and counterclockwise rotation. Under this labeling: $n=0$ corresponds to $m=0$ (non-degenerate), $n=1$ corresponds to $m=+1$ (the state that lowers in energy with $\lambda$ ), $n=2$ corresponds to $m=-1$ (the state that initially rises in energy), $n=3$ corresponds to $m=+2$ , and $n=4$ corresponds to $m=-2$ .

This assignment is unambiguous for $\lambda >0$ because the degeneracy is lifted and no true crossings occur between states of the same symmetry. The ordering $n=1$ (aligned) and $n=2$ (anti-aligned) is maintained for all $\lambda$ in our range, as shown in Figure 1 (left panel). This adiabatic labeling is essential for interpreting the orientation patterns in Figure 1 (right panel): the $n=1$ state shows positive orientation (green curve), while $n=2$ shows negative orientation (orange curve).

4.1.1. Weak-Field Perturbation Theory

In the weak-field limit ($\lambda \ll 1$ ), we can derive analytical expressions for the energy shifts using perturbation theory. The ground state energy is perturbed to second order [15]:

$E_0\left(\lambda \right)=-\frac{{\lambda }^2}{4}+O\left({\lambda }^4\right)$ (14)

For the first excited states, the degeneracy between $m=+1$ and $m=-1$ is lifted at first order, giving:

$E_{\pm }\left(\lambda \right)=1\pm \frac{\lambda }{2}+O\left({\lambda }^2\right)$ (15)

We have compared these analytical expressions with our numerical results and find excellent agreement for $\lambda <1$ (relative error <1%).

4.1.2. Strong-Field Harmonic Approximation

In the strong-field limit ($\lambda \gg 1$ ), the potential near $\varphi =0$ can be approximated as a harmonic oscillator [5]:

$V\left(\varphi \right)\approx -\lambda +\frac{\lambda }{2}{\varphi }^2$ (16)

The energy levels are then:

$E_n\left(\lambda \right)=-\lambda +\sqrt{\lambda }\left(n+\frac{1}{2}\right)$ (17)

We have verified that for $\lambda >5$ , the relative error between the numerical results and this harmonic approximation is less than 2% for the ground state.

4.1.3. Avoided Crossings and the Landau-Zener Formula

The avoided crossings observed in the energy spectrum can be characterized by the minimum gap $\Delta E_{\min }$ at the crossing point. The Landau-Zener formula [15] gives the transition probability between the two states:

$P_{\text{LZ}}=\text{exp}\left(-\frac{\pi \Delta E_{\min }^2}{2\hslash \dot{\lambda }}\right)$ (18)

For our system, the minimum gap between the $n=1$ and $n=2$ states is $\Delta E_{\min }\approx 0.15$ at $\lambda \approx 0.6$ , which is consistent with the smooth avoided crossing observed in Figure 1.

4.1.4. Detailed Analysis of Figure 1

In Figure 1 (left panel) we have shown the variation of lowest five states ($E_n$ ) of the one-dimensional quantum rotor with the field strength, as can be seen, at zero field strength the energies of the rotor are those of bare quantum rotor energies, with energies being, 0, 1, 4, 9 and 16 (in units of rotational constant or dimensionless units), corresponding to quantum numbers $m=0,\pm 1,\pm 2,\pm 3,\pm 4$ respectively. With an increase in field strength, the degeneracy of the $\pm m$ states gets removed, leading to a splitting of the energy levels. However, the ground state energy decreases with the field and reaches approximately −8.5 at a field strength of 10. This shows the binding effect of the attractive potential well created by the field. For the excited state, the behavior is almost similar, in the sense that as the field strength increases, the energy levels get separated, one of the levels, shown by the blue curve in Figure 1, decreases in energy and corresponds to the quantum state that aligns with the field direction. The other quantum level (orange curve in figure), initially increases in energy, reaches a maximum value around. This avoided behavior is characteristic of systems with competing interactions and represents the mixing between the $m=\pm 1$ states under the influence of the field. The variation of energy difference between $+m$ and $-m$ states provides the tunable transition frequency important for quantum control applications. Higher excited states show a bit of complex structures in the energy level variations; however, the degeneracy once removed remains with the levels. Despite this complexity, the lowest five states remain well-separated over the entire field strength range, with no actual level crossings, and the state ordering is preserved, enabling the identification of each state.

On the right panel of Figure 1, we have presented the orientation parameter values ${\langle \cos \varphi \rangle }_n\left(\lambda \right)$ for these five states, showing the directional behavior of the quantum rotor. This quantity is important because it directly measures the degree of molecular orientation with the field direction and is experimentally accessible through techniques such as laser-induced fluorescence or ion imaging [7].

The ground state orientation, shown as a blue curve, shows a monotonic increase from 0 to approximately 0.85 at $\lambda =10$ . This indicates substantial but not complete alignment even at the strongest fields considered. The approach to unity is slow since zero-point angular motion persists even in the strong-field limit due to quantum fluctuations. The functional form shows a sharp increase in the weak-field regime, followed by a slower rate of increase as the molecule becomes increasingly confined near $\varphi =0$ . Different states get oriented differently, for example, the first excited state that aligns with the field, shown by the green line, also shows an increase in anti-orientation from 0 at the zero field strength to about −0.48 at a field strength of about 5, beyond which it starts increasing. This reflects the fact that excited states have larger zero-point angular motion and therefore achieve less complete alignment. The second excited state again orients in the direction of the field, but its orientation always remains smaller than that of the ground state.

The higher excited states ($n=3$ and $n=4$ ) show more complex orientation behavior, which reflects the changing character of the eigenstates as they undergo avoided crossings. The orientation can change sign multiple times as the field strength varies. In Figure 1, scattered markers in both panels show the training data points, which were generated by direct diagonalization of the Hamiltonian for 100 equally spaced field strengths. The smooth solid curves represent the predictions from the trained neural network surrogate model. The agreement between the training points and the continuous predictions demonstrates the ability of the ML model to capture the physical behavior accurately. The neural network effectively interpolates between the discrete training points, providing continuous functions that maintain the correct physical features, such as the behavior of the ground state with the field strength and the avoided crossing structures for excited states.

4.2. Polynomial Regression Performance

In the problem undertaken, polynomial regression with degree 6 achieves very good accuracy for all outputs. The mean squared error on the test set, consisting of 20% of the data provided as training input, is of the order of 10⁻⁶ for the eigenenergies and 10⁻⁵ for the orientation parameter. To the best of our knowledge, errors of these orders are much smaller than the typical experimental precision. So this particular aspect makes the surrogate model essentially exact for practical purposes. The fraction of variance is typically measured by the R² coefficient, which exceeds 0.9999 for all outputs of the model, indicating near-perfect reproduction of the data. The slightly lower accuracy for orientation values compared to energies is expected because orientation is a more sensitive function of the wavefunction, which is a mixture of all states involved, but the errors remain well within acceptable limits.

To verify the model’s generalization capability, we evaluate it at the field strength values not included in the training data set. To explain it further, at the field strength of 7.5, which lies midway between training points, the predicted ground state energy differs from the exact numerical result by roughly 3 × 10⁻⁶, confirming that the polynomial of order six correctly interpolates between the training points. The orientation predictions show similarly high accuracy, with relative errors less than 10⁻⁴ for all test data points.

4.3. Comparison with Neural Network

In a neural network model, comparable accuracy is achieved in polynomial regression after careful hyperparameter optimization; however, differences emerge. The time scale in polynomial regression training is approximately 0.05 seconds, while the neural network takes about 2 - 3 seconds for 200 epochs. For frequent retraining purposes, this difference could be significant. Both these models provide fast predictions, with polynomial regression at approximately 1 - 2 microseconds per prediction and the neural network at approximately 10 - 20 microseconds. Both these methods are faster than direct diagonalization, which requires about 5 milliseconds per field strength calculation, representing speedups of 2500 to 5000 times.

The neural network is smooth, showing small fluctuations because of the stochastic nature of the training process, while the polynomial regression predictions are smooth. In the case of neural networks, fluctuations are negligible for practical applications, but can become problematic for inverse applications, which require more precision. The coefficients of polynomials are direct evidence of insight into the functional dependence of the observables on the field strength. As an example, for a weak field, the quadratic term is dominant for the ground state energy, reflecting the Stark effect [4]. This interpretability is somewhat lost with neural networks. Given these observations, polynomial regression emerges as a better approach for this particular problem.

4.4. Inverse Design Application

A surrogate model is very useful in inverse design: for example, in this case, finding the field strength that yields a desired orientation or alignment. To design an experiment to prepare the ground state with a specific orientation ${\langle \cos \varphi \rangle }_0=0.5$ . Using the polynomial regression model, we can solve the equation ${\langle \cos \varphi \rangle }_0\left(\lambda \right)=0.5$ . Since the polynomial representation is explicit, this equation can be solved numerically using standard root-finding algorithms. As another option for any desired target, the model can be inverted directly.

We demonstrate this capability by finding the field strength that provides a target ground state energy of $E_0=-5.0$ . Solving the polynomial equation using the fminbnd optimizer yields the field strength values of 5.247, giving an energy value of −5.0002, thus validating our inverse design approach. This utility is particularly useful for experimental applications, where the field strength values are to be determined for a given quantum state and a particular value of orientation. This kind of inverse approach enables real-time experimental control and adaptive pulse shaping [19].

4.5. Extrapolation beyond Training Range

It is worth mentioning that the interest is in the interpolation within the training range. We also performed calculations beyond the given range of the field strength, but its accuracy decreases when we move beyond training data points. It is noticed that the error of 0.01 appears when the field strength reaches a value of 12. Reaching a value of 15 for the field strength, the error reaches 0.1, which is highly unacceptable. For accurate values for the output at the higher field strengths, the training data set is to be expanded. The behavior of the neural network is similar when extrapolating, as both models are purely data-driven and lack built-in physical constraints beyond what is implicitly captured in the training data.

4.6. Error Analysis and Uncertainty Quantification

To test the reliability of the predictions for applications when accuracy matters, one needs to check the uncertainty in the results [28]. In the case of polynomial regression, one can use the covariance matrix of the fitted coefficients for this purpose. The variance in a prediction is then ${\sigma }^2\left(\lambda \right)=x{\left(\lambda \right)}^{\text{T}}Cx\left(\lambda \right)$ , where $x\left(\lambda \right)$ holds the polynomial basis functions evaluated at a given field strength, and $C$ is the coefficient covariance matrix. The uncertainty rises at the ends of the training range, maximizing at the field strength of 10, with a value 2 × 10⁻⁵ for the ground state energy. The uncertainty coming here is actually smaller than the error introduced by truncating the basis, which tells us that the polynomial model itself isn’t the main source of imprecision here. For the neural network, uncertainty quantification is more challenging and typically requires ensemble methods or Bayesian neural networks [30], which were not implemented in this work. This represents another advantage of the polynomial regression approach.

4.6.1. Error Analysis near Avoided Crossings

Avoided crossings present the greatest challenge for surrogate models because eigenstates mix rapidly over small $\lambda$ intervals, causing sharp variations in both energies and orientation values. For the one-dimensional rotor, the first significant avoided crossing occurs between states derived from $|m=1\rangle$ and $|m=-1\rangle$ near $\lambda \approx 4$ (visible in Figure 1 as the region where the green and orange curves approach each other).

To quantify our model’s performance in this region, we performed refined sampling: we generated 50 additional test points with $\lambda \in \left[3.5,4.5\right]$ (step size 0.02) that were excluded from training. The results are shown in Table 3:

Table 3. Error analysis near the avoided crossing region. Max error is represented as $M_e$ .

The errors increase by a factor of 3 - 4 near the avoided crossing but remain well below experimental precision (typically 10⁻³ - 10⁻⁴ for energy ratios and 10⁻² for orientation measurements). The polynomial model faithfully captures the avoided crossing without spurious oscillations because the underlying function remains analytic (the eigenvalues of an analytic matrix family are analytic except at exact crossings, which do not occur).

For higher excited states ($n=3,4$ ), which undergo multiple avoided crossings, the maximum errors increase to ~10⁻⁵ for energies and ~5 × 10⁻⁴ for orientation. This remains acceptable for most applications, but users requiring ultra-high precision for highly excited states should either expand the training set or use a higher-degree polynomial.

4.6.2. Uncertainty Quantification for Neural Networks

Unlike polynomial regression, where the coefficient covariance matrix provides analytic uncertainty estimates, standard feedforward neural networks trained with mean squared error loss do not naturally produce prediction uncertainties. Several techniques exist to address this limitation [30]:

Deep Ensembles: Training multiple networks with different random initializations yields a distribution of predictions, from which the mean and variance can be estimated. We implemented a 5-network ensemble as a post-hoc test. The ensemble produced uncertainty estimates that were comparable to polynomial regression near training points ($\sigma \approx 2\times {10}^{-5}$ for ground state energy) but increased to $\sigma \approx 1\times {10}^{-4}$ near $\lambda =10$ (the training boundary) and showed bimodal behavior near avoided crossings, reflecting mode-switching between different local minima.

Monte Carlo Dropout: Applying dropout at inference time (with 100 forward passes) provides another uncertainty estimate. This approach yielded larger uncertainties ($\sigma \approx 5\times {10}^{-5}$ ) but was computationally expensive (5 ms per prediction, eliminating the speed advantage).

Bayesian Neural Networks: While theoretically attractive, BNNs require variational inference or Hamiltonian Monte Carlo, which are impractical for our small dataset and would require specialized libraries not used in this work.

Given that polynomial regression provides superior uncertainty quantification at no additional computational cost and with comparable accuracy, we recommend it over neural networks for this specific application. For problems where the functional form is unknown or the input dimensionality is high (>5), neural networks with ensemble-based uncertainty quantification would be more appropriate.

4.7. Connection to Experimental Observables

The quantities predicted by our surrogate model connect directly to experimentally measurable observables [7] [20]. The eigenenergies determine the frequencies of rotational transitions: ${\omega }_{n\leftarrow m}=E_n-E_m$ . These transition frequencies are accessible through microwave spectroscopy and are essential for designing control pulses. The orientation expectation values influence the strength of interactions with other fields, including Stark shifts where the energy shift in an additional weak field is proportional to $\langle \cos \varphi \rangle$ , laser coupling where the transition dipole moment between rotational states depends on orientation, collision dynamics where the anisotropic interaction potential depends on the instantaneous orientation, and alignment measurements where techniques such as laser-induced alignment exploit the field-induced orientation. By providing fast and accurate predictions of these quantities, our surrogate model enables rapid iteration in experimental design and real-time feedback control [19].

4.8. Limitations and Future Improvements

Despite the success of the approach, several limitations should be noted [27]. The model only applies to static electric fields, while time-dependent fields, which are crucial for quantum control, would require learning the dynamics, a considerably more challenging problem. The restriction to planar rotation omits the out-of-plane degree of freedom present in real molecules, and while the one-dimensional model captures the essential physics, quantitative predictions for specific molecules would require the full three-dimensional treatment. The model describes pure quantum states at zero temperature, while thermal averaging over the Boltzmann distribution would require predicting the full set of states and their populations. As discussed, the model’s accuracy degrades outside the training range.

Future work could address these limitations through physics-informed neural networks that incorporate the Schrödinger equation as a constraint, which could improve extrapolation and reduce the required training data [30]. Gaussian process regression naturally provides uncertainty estimates and can incorporate prior knowledge about the functional form [28]. In order to reduce the computational costs of maintaining desired accuracy, the multi-fidelity modeling may be used [29], also for transfer learning, in which a training model used for one particular molecule can be used for other molecules [27].

4.9. Extensions to Higher-Dimensional and Time-Dependent Systems

While this work focuses on the one-dimensional rotor in a static field, the surrogate modeling approach can be extended along several physically important directions:

Three-Dimensional Rotors: The Hamiltonian becomes $\widehat{H}=B{J}^2-\mu F\cos \theta$ , with eigenfunctions given by spherical harmonics. The Hamiltonian matrix is tridiagonal in the angular momentum quantum number $j$ (coupling $j$ to $j\pm 1$ ), similar to the 1D case but with degeneracy $2j+1$ . A surrogate model would then map $\lambda$ to $E_{j,m}\left(\lambda \right)$ and ${\langle \text{cos}\theta \rangle }_{j,m}\left(\lambda \right)$ . The main challenge is the increased number of states; however, for low-lying states (the ones relevant to most experiments), the dimensionality remains manageable. We estimate that a polynomial surrogate for the 3D rotor would require training on $~200\lambda$ values (to capture the richer avoided crossing structure) and would achieve similar speedups.

Time-Dependent Fields: For a time-varying field $F\left(t\right)$ , the system’s dynamics are governed by the time-dependent Schrödinger equation. A surrogate model would need to map the entire pulse shape to final observables—a function-to-function regression problem. This is considerably more challenging but could be addressed using recurrent neural networks or neural ordinary differential equations. Recent work on quantum control of molecular rotation [9] [19] demonstrates that such surrogates could accelerate optimal control pulse design by orders of magnitude.

Multi-Parameter Control Landscapes: In many experiments, multiple control parameters (field strength, field orientation, laser detuning, temperature) must be simultaneously optimized. Gaussian process regression [28] is particularly well-suited for such high-dimensional surrogate modeling, as it can efficiently explore the control landscape and identify optimal parameter sets with minimal experimental trials. Our polynomial approach would scale poorly to >3 parameters due to the exponential growth in the number of polynomial terms. Future work should therefore employ GPR or random forests for multi-parameter optimization.

Chaotic or Non-Analytic Regimes: The kicked rotor [22]-[25] exhibits dynamical localization and quantum chaos, where observables are non-analytic functions of parameters. Polynomial regression would fail in such regimes. For these systems, neural networks with sufficient capacity or kernel methods with appropriate feature maps would be necessary. This represents an important frontier for surrogate modeling in quantum dynamics.

5. Practical Implementation

5.1. Computational Method

The method to develop the surrogate model of interest consists of several steps, starting with the selection of parameters such as the maximum number of states to be taken into account, followed by the field strength range based on desired accuracy and applicability. In our case, 15 maximum of states, and for 5 states and field strength in the range of 0 to 10 provides a reasonable balance. The second step in the computation is designing the training for each field strength, constructing the Hamiltonian matrix, and calculating the desired output data matrix of required quantities, here in our case, the dressed state energies, and the directional parameter of interest. The next step in the process is polynomial fitting of the data and stringing the coefficients of polynomials, followed by validation by held-out field strength values to cross-check the accuracy, and further calculating the errors.

Finally, deployment uses the polynomial coefficients to evaluate predictions for any $\lambda$ in the training range, and for inverse problems, root-finding algorithms are implemented to solve for $\lambda$ given a target output.

The MATLAB implementation of this study runs in approximately 1 - 2 seconds for data generation and 0.05 seconds for model training. Once trained, predictions require 1 - 2 microseconds each, enabling applications that require millions of evaluations.

5.2. Integration with Experimental Control Systems

The surrogate model is ideally suited for integration with experimental control systems [19]. In a typical experiment, the control computer needs to compute the field strength required to achieve a particular quantum state. Using direct diagonalization, this would require solving an optimization problem that calls the eigenvalue solver repeatedly, taking several seconds per optimization. With the surrogate model, the same optimization can be performed in milliseconds, enabling real-time feedback. The model can be implemented in any programming language with basic linear algebra capabilities. In this method, the evaluation of a prediction requires only a polynomial estimation that can be easily implemented on embedded systems.

6. Conclusions

In this work, we have presented an ML approach to study the rotational excitation of a one-dimensional quantum rotor interacting with an external static electric field. A polynomial regression surrogate model is developed, which achieves a very good accuracy and prediction time of microseconds.

The first part is devoted to the development of a theoretical model, where the dressed state energy levels and directional parameters are evaluated, which form the basis for ML applications. The data generation is done with a basis set comprising the lowest 15 states for the training dataset, including all ranges of field strengths. The surrogate model we developed using polynomial regression of degree six comes out to be useful with $R^2>0.09999$ for the outputs. The model gets trained in milliseconds, and it predicts the output data in microseconds.

We further showed that using inverse design, the model produces a rapid solution of the inverse problem, correctly predicting the values of the external fields for desired levels of orientation and alignment. The model correctly predicts the behaviour of energy levels and orientation parameters across the full range of field strengths studied.

The MATLAB code developed can be easily extended to higher-dimensional analysis of similar problems [27] [28]. With the development of experimental advancements, such surrogate models will make it possible to have real-time control of computational rigors that will help in reducing the gap between theoretical and experimental studies [29] [30]. The methodology presented here establishes a quantitative framework for characterizing quantum rotor systems in external fields, complementing traditional energy-level spectroscopy and offering a practical tool for experimental control and parameter optimization.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References