Nonlinear dynamics and energy transfer for two rotating dipoles in an external field: A three-dimensional analysis

TL;DR

This paper analyzes the complex nonlinear dynamics and energy transfer mechanisms of two coupled rotating dipoles in an external electric field, revealing transitions between energy regimes and chaotic behavior.

Contribution

It provides a three-dimensional classical analysis of dipole interactions under external fields, highlighting energy transfer and chaos in a novel setting.

Findings

Abrupt transition between energy equipartition and non-equipartition regimes

Chaotic behavior dominates in strong electric fields

Energy transfer mechanisms are independent of initial rotation directions

Abstract

We investigate the structure and the nonlinear dynamics of two rigid polar rotors coupled through the dipole-dipole interaction in an external homogeneous electric field. In the field-free stable head-tail configuration, an excess energy is provided to one of the dipoles, and we explore the resulting three-dimensional classical dynamics. This dynamics is characterized in terms of the kinetic energy transfer between the dipoles, their orientation along the electric field, as well as their chaotic behavior. The field-free energy transfer mechanism shows an abrupt transition between equipartition and non-equipartition regimes, which is independent of the initial direction of rotation due to the existence of an infinite set of equivalent manifolds. The field-dressed dynamics is highly complex and strongly depends on the electric field strength and on the initial conditions. In the strong field regime, the energy equipartition and chaotic behavior dominate the dynamics.

Equations

$$V_d=\frac{1}{4\pi {ϵ}_0}\frac{({\mathbf{d}}_1\cdot {\mathbf{d}}_2)r^2-3({\mathbf{d}}_1\cdot \mathbf{r})({\mathbf{d}}_2\cdot \mathbf{r})}{r^5},$$

$$V({\mathbf{l}}_1,{\mathbf{l}}_2,t)=-q\mathcal{E}(t)(z_1+z_2)+\frac{q^2}{4\pi {ϵ}_0{a}_l^3}[{\mathbf{l}}_1\cdot {\mathbf{l}}_2-3x_1{x}_2],$$

$$f(t)=\{\begin{array}{ccc}\frac{t}{t_1}& \text{mbox}if& 0\le t<t_1\\ 1& \text{mbox}if& t\ge t_1\end{array}$$

$$E\equiv H=\sum_{i=1}^{2}\frac{1}{2\mu}\bigg{[}P_{x_{i}}^{2}+P_{y_{i}}^{2}+P_{z_{i}}^{2}\bigg{]}+V(\mathbf{l}_{1},\mathbf{l}_{2},t),$$

$${\mathbf{l}}_i^{\prime }=(x_i^{\prime },y_i^{\prime },z_i^{\prime })=(x_i/l,y_i/l,z_i/l),i=1,2,$$

$$E^{\prime}\equiv H^{\prime}=\frac{H}{\chi}=\sum_{i=1}^{2}\frac{1}{2}\bigg{[}P_{x_{i}}^{\prime 2}+P_{y_{i}}^{\prime 2}+P_{z_{i}}^{\prime 2}\bigg{]}+{\cal V}_{1}(\mathbf{l}^{\prime}_{1},\mathbf{l}^{\prime}_{2},t^{\prime}).$$

$${\mathcal{V}}_1({\mathbf{l}}_1^{\prime },{\mathbf{l}}_2^{\prime },t^{\prime })=-\beta f(t^{\prime })(z_1^{\prime }+z_2^{\prime })+({\mathbf{l}}_1^{\prime }\cdot {\mathbf{l}}_2^{\prime }-3x_1^{\prime }{x}_2^{\prime }).$$

$$x_i=\cos {\theta }_i,y_i=\sin {\theta }_i\cos {\varphi }_i,z_i=\sin {\theta }_i\sin {\varphi }_i,i=1,2$$

$$H=\sum_{i=1}^{2}\frac{1}{2}\bigg{[}P_{\theta_{i}}^{2}+\frac{P_{\phi_{i}}^{2}}{\sin^{2}\theta_{i}}\bigg{]}+{\cal V}_{2}(\theta_{1},\phi_{1},\theta_{2},\phi_{2};t),$$

$${\mathcal{V}}_2({\theta }_1,{\varphi }_1,{\theta }_2,{\varphi }_2;t)$$

$${\dot{\theta }}_1$$

$${\dot{P}}_{{\theta }_1}$$

$${\dot{P}}_{{\theta }_2}$$

$${\dot{P}}_{{\varphi }_1}$$

$${\dot{P}}_{{\varphi }_2}$$

$$\mathcal{M}=\{({\theta }_1,P_{{\theta }_1},{\theta }_2,P_{{\theta }_2})|{\varphi }_1-{\varphi }_2=0,$$

$$\text{mbox}and{P}_{{\varphi }_1}=P_{{\varphi }_2}=0;\beta =0\}.$$

$${\mathcal{H}}_{\mathcal{M}}=\frac{P_{{\theta }_1}^2+P_{{\theta }_2}^2}{2}+\sin {\theta }_1\sin {\theta }_2-2\cos {\theta }_1\cos {\theta }_2,$$

$${\mathcal{M}}_{XZ}=\{({\theta }_1,P_{{\theta }_1},{\theta }_2,P_{{\theta }_2})|{\varphi }_1={\varphi }_2=0,$$

$$\text{mbox}and{P}_{{\varphi }_1}=P_{{\varphi }_2}=0\}.$$

$$\mathcal{N}=\{({\varphi }_1,P_{{\varphi }_1},{\varphi }_2,P_{{\varphi }_2})|{\theta }_1={\theta }_2=\pi /2,P_{{\theta }_1}=P_{{\theta }_2}=0\}.$$

$${\mathcal{H}}_{\mathcal{N}}=\frac{P_{{\varphi }_1}^2}{2}+\frac{P_{{\varphi }_2}^2}{2}-\beta f(t)(\sin {\varphi }_1+\sin {\varphi }_2)+\cos {\varphi }_{12}.$$

$${\mathcal{V}}_1({\mathbf{l}}_1,{\mathbf{l}}_2)=-\beta (z_1+z_2)+[{\mathbf{l}}_1\cdot {\mathbf{l}}_2-3x_1{x}_2]$$

$${\mathcal{V}}_L({\mathbf{l}}_1,{\mathbf{l}}_2,{\lambda }_1,{\lambda }_2)$$

$$x_2+{\lambda }_1{x}_1=0,x_1+{\lambda }_2{x}_2=0,$$

$$y_2-2{\lambda }_1{y}_1=0,y_1-2{\lambda }_2{y}_2=0,$$

$$z_2-2{\lambda }_1{z}_1-\beta =0,z_1-2{\lambda }_2{z}_2-\beta =0,$$

$$1-|{\mathbf{l}}_1{|}^2=0,1-|{\mathbf{l}}_2{|}^2=0.$$

$$x_1(0)$$

$$P_{y_1}(0)$$

$$P_{x_1}(0)$$

$${\hat{K}}_i=\frac{⟨K_i⟩}{⟨K_1⟩+⟨K_2⟩},$$

$$⟨K_i⟩=\frac{1}{2t_f}{\int }_0^{t_f}[P_{xi}^2(t)+P_{yi}^2(t)+P_{zi}^2(t)]dt,$$

$$⟨z_i⟩=\frac{1}{t_f}{\int }_0^{t_f}{z}_i(t)dt,i=1,2.$$

$$\frac{d\mathbf{r}}{dt}=\mathbf{f}(\mathbf{r},t),$$

$$\frac{d\delta \mathbf{r}}{dt}=\frac{\partial \mathbf{f}(\mathbf{r},t)}{\partial \mathbf{r}}\delta \mathbf{r}.$$

Full text

I Introduction

The experimental availability of ultracold dipolar gases represents a strong motivation for the investigation of the physical phenomena related to the long-range and anisotropic dipole-dipole interaction pfau12 . By tuning the dipole-dipole interaction, the dipolar gas properties can be significantly modified giving rise to a rich variety of novel applications, such as, the creation of novel quantum liquids pfau16 ; ferlaino16 ; pfau16_2 , exotic spin dynamics naylor16 ; rey19 , and the control of ultracold chemical reactions ye10 ; ye11 ; ye13 , thermalization lev18 or energy exchange bakker99 ; tolbert00 ; noel15 ; wenger19 .

The mechanism of energy transport mediated by the dipole-dipole interaction is closely related to the nonlinear behavior of many-body systems and has attracted special attention in classical dynamics  michl ; sim ; jonge ; rotors2D ; pre2018 ; estevez18 . Most of the nonlinear dynamical studies on coupled rotating dipoles restrict their motion to planar rotations. As a consequence, each dipole is described by one angle and the number of degrees of freedom of the system is equal to the number of dipoles. A natural extension of these studies is to allow the dipoles to perform three-dimensional (3D) overall rotations, and, to the best of our knowledge, the corresponding literature is very scarce. The 3D classical dynamics has been explored for isolated rotors, such as diatomic or symmetric top molecules, exposed to combinations of external electric fields A516 ; A286

In our previous work rotors2D , we studied the classical dynamics of two polar rotors, coupled by the dipole-dipole interaction, in the presence of an external electric field in a planar invariant manifold. By restricting the motion of the dipoles to the invariant manifold, the system has two degrees of freedom. The energy transfer mechanism between the dipoles for varying field strength has been analyzed in terms of the phase space structure of the system. Here, we extend this previous study by exploring the complete 3D classical dynamics allowing the dipoles to rotate in any direction in space. Since the rotors motion is not restricted to a planar manifold as in Ref. rotors2D , each dipole is described by two angles, and we encounter a Hamiltonian system with four degrees of freedom. We assume that the dipoles are initially in the stable head-tail configuration with fixed spatial positions. This stable configuration is perturbed by adding a certain excess energy to one of the dipoles, which starts rotating from the head-tail configuration axis, while at the same time the electric field is turned on. The follow up dynamics is investigated in terms of the energy exchange mechanism and the orientations of the rotors induced by the electric field. We also explore the chaoticity of the system by using a fast chaos indicator.

Due to the existence of the infinite set of equivalent manifolds, our field-free results show that the energy transfer mechanism is independent of the direction for which the dipole starts to rotate, which is in complete agreement with our previous results rotors2D . The field-dressed dynamics is highly complex and strongly depends on the strengths of the dipole and the electric field interactions as well as on the initial conditions. In the very weak field regime, when the dipole interaction is dominant, the energy transfer mechanism resemble the field-free case, although the border between the non-equipartition and equipartition regimes becomes more diffuse. As the electric field increases, the energy equipartition regime dominates the dynamics, although we still encounter regions of non-equipartition energy for certain initial conditions. For strong fields, the two dipoles are significantly oriented along the electric field axis. When the electric field is turned on, the system becomes non-integrable so that there appears chaotic motion. For low and intermediate values of the field, the degree of chaoticity increases as the interaction with the electric field becomes more dominant. Surprisingly, for strong electric fields the chaoticity of the system remains very pronounced which is quite unexpected because the gradual increase of the field would eventually lead the system to its integrable limit. Indeed, we verified that, for very large values of the electric field, a slow crossover to an integrable phase space takes place.

This work is organized as follows. In Sec. II we discuss the classical Hamiltonian of the system, its symmetries and invariant manifolds. The underlying equilibrium points are also presented. In Sec. III we explore the time evolution of the energy transport between the dipoles and their orientations with varying initial conditions and for several the electric field strengths. The chaoticity of this systems is discussed in Sec. IV. Our conclusions are provided in Sec. V. In the Appendix, an exhaustive analysis of the existence, stability and bifurcations of the equilibria is given.

II Classical Hamiltonian, symmetries and invariant manifolds

II.1 Classical Hamiltonian

The potential energy $V_{d}$ between two dipoles with dipole moments ${\bf d}_{1}$ and ${\bf d}_{2}$ due to the mutual dipole-dipole interaction (DDI) is given by A753

[TABLE]

with ${\bf r}$ being their relative position. In our case, we consider two identical rigid rotors having electric dipole moments ${\bf d}_{1}=q\ {\bf l_{1}}$ and ${\bf d}_{2}=q\ {\bf l_{2}}$, with $d=|{\bf d}_{1}|=|{\bf d}_{2}|=q\ l$, being $q$ and $l$ the charge and length of the dipoles.

Besides the DDI, the dipoles are in the presence of an external homogeneous time-dependent electric field ${\cal E}(t)$ parallel to the Laboratory Fixed Frame (LFF $XYZ$) $Z$-axis. Assuming that the positions of the rotors are fixed along the LFF $X$-axis and separated by a constant distance $a_{l}$ (see Fig.1), the total interaction potential reads as follows

[TABLE]

where the vectors $\{\mathbf{l}_{i}=(x_{i},y_{i},z_{i})||\ ,|\mathbf{l}_{i}|=l,\ i=1,2\}$ determine the orientation of the dipoles. The first term in (2) stands for the interaction with the external field and the second one for the DDI. Here, we assume that the homogeneous electric field has a turn-on modeled by a linear ramp $f(t)$, so that ${\cal E}(t)=E_{s}f(t)$ with $f(t)$

[TABLE]

This linear ramp was used in Ref. chandre and it also mimics the turn-on of laser pulses  A182 ; A747 .

If $\mu$ is the reduced mass of the dipoles, the dynamics of the system is governed by the classical Hamiltonian (the energy $E$)

[TABLE]

subject to the holonomic constraints $\{\mathbf{l}_{i}=(x_{i},y_{i},z_{i})||\ ,|\mathbf{l}_{i}|=l,\ i=1,2\}$. At this point, and without loss of generality, we use a dimensionless version of the Hamiltonian (4) by expressing the energy in units of the parameter $\chi=d^{2}/4\pi\epsilon_{0}a_{l}^{3}$ that controls the DDI. To do that, we introduce the dimensionless coordinates

[TABLE]

and the dimensionless time $t^{\prime}=t/t_{d}$, where $t_{d}=\sqrt{I/\chi}$ is the new unit of time being $I=\mu\ l^{2}$ the moment of inertia of the dipoles. After applying these transformations to (4), we arrive at the following dimensionless Hamiltonian

[TABLE]

where the potential ${\cal V}_{1}(\mathbf{l}^{\prime}_{1},\mathbf{l}^{\prime}_{2},t^{\prime})$ reads as

[TABLE]

In this way, the dynamics depends on the energy $E^{\prime}\equiv H^{\prime}=H/\chi$ and on the new (dimensionless) electric field parameter $\beta=dE_{s}/\chi$, which is the ratio between the strengths of the electric and the dipole-dipole interactions. For the sake of simplicity, we omit in the following the primes in Hamiltonian (5)

The aforementioned holonomic constraints between the Cartesian coordinates of the dipoles reduce the number of degrees of freedom from 6$+1/2$ to 4$+1/2$. By taking the LFF $X$-axis as the polar axis, the transformation between the Cartesian and the Euler angles $(\theta_{1},\phi_{1},\theta_{2},\phi_{2})$ of each rotor reads (see Fig.1)

[TABLE]

and the Hamiltonian (5) converts to

[TABLE]

where the interaction potential (2) is

[TABLE]

In spherical coordinates, the Hamiltonian (7) defines a $(4+1/2)$-degree-of-freedom dynamical system with coordinates $(\theta_{1},\phi_{1},\theta_{2},\phi_{2})$ and the corresponding momenta $(P_{\theta_{1}},P_{\phi_{1}},P_{\theta_{2}},P_{\phi_{2}})$.

II.2 Symmetries and invariant manifolds

Since the rotors are identical, Hamiltonians (5) and (7) possess a exchange symmetry. Besides this discrete symmetry, the field-free system, i. e., considering only the dipole interaction and $\beta=0$, is also invariant under rotations around the common LFF $X$-axis. Besides the energy, this continuous symmetry implies that the $X$-component of the total angular momentum $L_{X}=P_{\phi_{1}}+P_{\phi_{2}}$ is conserved.

The Hamiltonian equations of motion arising from (7) read as follows:

[TABLE]

For $\beta=0$, the infinite set of manifolds ${\cal M}$ of codimension four given by

[TABLE]

are invariant under the dynamics. On each of these manifolds, the Hamiltonian (7) for $\beta=0$ reduces to the two degrees of freedom Hamiltonian

[TABLE]

and the rotational motion of the dipoles is restricted to a given common polar plane of constant azimuthal inclination $\phi_{1}-\phi_{2}=0$ where the polar angles $(\theta_{1},\theta_{2})$ vary in the interval $[-\pi,\pi)$. The existence of the manifolds ${\cal M}$ is associated to the aforementioned rotational invariance of the Hamiltonian (5) around the LFF $X$-axis. For $\beta\neq 0$, the electric field breaks this rotational symmetry, and $L_{X}=P_{\phi_{1}}+P_{\phi_{2}}$ is no longer an integral of the motion. It is worth noticing that the presence of the electric field reduces the (infinite) invariant manifolds ${\cal M}$ to a single one along the direction $\phi_{1}=\phi_{2}=0$, namely

[TABLE]

The planar dynamics and the energy transfer on the manifold ${\cal M}_{XZ}$ has already been studied in Ref. rotors2D .

The equations of motion II.2 provide an additional invariant manifold ${\cal N}$ of codimension four,

[TABLE]

When the system is moving on this manifold ${\cal N}$, the dynamics is governed by the $(2+1/2)$-degree-of-freedom Hamiltonian

[TABLE]

On the manifold ${\cal N}$, the rotational motion of the dipoles is restricted to the parallel $y_{1}-z_{1}$ and $y_{2}-z_{2}$ planes, respectively. This configuration was already considered in Ref. jonge for the general case of a chain of $N\geq 2$ dipoles and zero electric field.

II.3 Equilibrium configurations

The equilibrium points of a dynamical system provide useful information about its behavior. For $t\geq t_{1}$, the electric field parameter reaches its maximal value $\beta$ with $f(t)=1$, and, using the Cartesian Hamiltonian (5), the equilibrium points are the critical points of the potential

[TABLE]

which is ${\cal V}_{1}(\mathbf{l}_{1},\mathbf{l}_{2},t)$ ((6)) with $t\geq t_{1}$, under the constraints $\{\mathbf{l}_{i}=(x_{i},y_{i},z_{i})||\,|\mathbf{l}_{i}|^{2}=1,\ i=1,2\}$, together with the conditions $P_{x_{i}}=P_{y_{i}}=P_{z_{i}}=0$. Thence, by introducing the Lagrange multipliers $\lambda_{1}$ and $\lambda_{2}$, the critical points of (13) are the extrema of the Lagrange function ${\cal V}_{L}$,

[TABLE]

Thus, the critical points are roots of the system of equations $\vec{\nabla}_{x_{i},y_{i},z_{i},\lambda_{1},\lambda_{2}}{\cal V}_{L}=0$, given by

[TABLE]

In the presence of an electric field ($\beta\neq 0$), the roots of (15) provide the following six sets of (isolated) critical points:

• i)

Two critical points $P_{1}$: $z_{1}=z_{2}=\beta/3$, $y_{1}=y_{2}=0$, $x_{1}=x_{2}=\pm\sqrt{1-\beta^{2}/9}$. These points exist when $\beta<3$ and their energy is $E_{1}=-(6+\beta^{2})/3$. In the field-free case, i. e., for $\beta=0$, they correspond to the well-known stable head-tail configurations of the two dipoles along the common LFF $X$-axis.

• ii)

Two critical points $P_{2}$: $z_{1}=-z_{2}=\pm 1$ and $x_{1}=x_{2}=y_{1}=y_{2}=0$. The energy of these points is $E_{2}=-1$ and they exist for values $\beta>0$.

• iii)

The critical point $P_{3}$: $z_{1}=z_{2}=1$ and $x_{1}=x_{2}=y_{1}=y_{2}=0$. Its energy is $E_{3}=1-2\beta$, and it exists for $\beta>0$.

• iv)

The critical point $P_{4}$: $z_{1}=z_{2}=-1$ and $x_{1}=x_{2}=y_{1}=y_{2}=0$. The energy of this point is $E_{4}=1+2\beta$ and it exists for $\beta>0$.

• v)

Two critical points $P_{5}$: $z_{1}=z_{2}=-\beta$, $y_{1}=y_{2}=0$, $x_{1}=-x_{2}=\pm\sqrt{1-\beta^{2}}$. These points exist when $0\leq\beta<1$ and their energy is $E_{5}=(2+\beta^{2})$. In the absence of an electric field $\beta=0$, they correspond to the well-known unstable tail-tail or head-head configurations of the two dipoles along the common LFF $X$-axis.

• vi)

Two critical points $P_{6}$: $z_{1}=z_{2}=\beta/2$, $y_{1}=-y_{2}=\pm\sqrt{1-\beta^{2}/4}$, $x_{1}=x_{2}=0$. These points exist when $0\leq\beta<2$ and their energy is $E_{6}=-(2+\beta^{2})/2$.

In contrast, for the field-free case $\beta=0$, the number of critical points reduces to $P_{1}$ and $P_{5}$, and to two degenerate circles of equilibria $D_{1}$ and $D_{2}$:

• •

The set $D_{1}$ is given by $z_{1}=z_{2}=\cos\alpha$, $y_{1}=y_{2}=\sin\alpha$, $x_{1}=x_{2}=0$, with $\alpha=[0,2\pi)$. The energy of this circle of stationary points is $E_{3}=1$, and, when $\beta=0$, the former isolated critical points $P_{3}$ and $P_{4}$ for $\beta\neq 0$ are included in $D_{1}$.

• •

The set $D_{2}$ is given by $z_{1}=-z_{2}=\cos\alpha$, $y_{1}=-y_{2}=\sin\alpha$, $x_{1}=x_{2}=0$, with $\alpha=[0,2\pi)$. The energy of these equilibria is $E_{4}=-1$, and, when $\beta=0$, the former isolated equilibria $P_{2}$ and $P_{6}$ for $\beta\neq 0$ are included in $D_{2}$.

A detailed study of the stability and existence of the critical points as well as their bifurcations as $\beta$ is varied is provided in the Appendix.

Let us explain the differences between the equilibria of the full 3D system and when the dynamics is restricted to the planar manifold ${\cal M}_{XZ}$ rotors2D . For $\beta=0$, there are two isolated equilibria $P_{1,5}$ and the degenerate circles of equilibria $D_{1,2}$; whereas on the manifold ${\cal M}_{XZ}$, there are five isolated equilibria $P_{1,2,3,4,5}$. For $\beta\neq 0$, all the equilibria are isolated, such that in the 3D system there appears an additional equilibrium $P_{6}$ located off the manifold ${\cal M}_{XZ}$.

III Energy transfer processes

In this section, we explore the classical dynamics of the two dipoles in an external electric field. We assume that initially at $t=0$, the two dipoles are at rest, with zero kinetic energy, in the stable head-tail configuration along the LFF $X$-axis given by the equilibrium points $P_{1}$. In this state, the system has the minimal energy $E_{1}=-2$. From this initial configuration, the field is turned on at $t=0$ by the ramp-up function (3), and a certain excess energy $\delta K$ is given to dipole one. Therefore, taking into account the holonomic constraints $\{\mathbf{l}_{i}=(x_{i},y_{i},z_{i})||\ |\mathbf{l}_{i}|^{2}=1,i=1,2\}$, the system leaves the stable equilibrium configuration in such a way that, at $t=0$, the initial conditions of the dipoles are

[TABLE]

The angle $\alpha\in[-\pi,\pi)$ allows us to consider all possible directions of rotations perpendicular to the head-tail axis, and all initial conditions of its momentum. Using the initial conditions III, the equations of motion given by the Hamiltonian (5) are numerically integrated up to a final time $t_{f}$ by means of the so-called Störmer-Verlet algorithm Andersen . This numerical algorithm is a symplectic integrator that preserves the holonomic constraints of the system. The final integration time is fixed to $t_{f}=400$. Our numerical tests have shown that this stopping or final time is appropriate for a proper characterization of the outcomes. In particular, for a LiCs molecule in its ground state trapped in an optical lattice with $a_{l}=450$ nm, the time unit is $t_{d}=\sqrt{I/\chi}\approx 6.5$ ns, and the final time $t_{f}=400\approx 2500$ ns. The ramp-up time is fixed to $t_{1}=2$, that roughly corresponds to $12$ ns and that can be achieved in current experiments with realistic field strengths.

For these initial conditions, we explore the dynamics of this system for different ratios of the electric field interaction and the dipole-dipole one, specifically $\beta=0,0.1,1,10,100,1000$. To do so, we compute the normalized time-average of the kinetic energy of each dipole, $\widehat{K}_{i}$, given by

[TABLE]

with

[TABLE]

where $t_{f}$ is the final time. In the presence of the electric field, the axial symmetry no longer exists and the dipoles tend to orient along the electric field direction. We characterize their orientations by computing the time-average of the Cartesian coordinate $z_{i}$

[TABLE]

Computationally, these time-averaged quantities $\widehat{K}_{i}$ and $\langle z_{i}\rangle$ are easily calculated as a function of the kinetic energy $\delta K$ given to dipole one and of the angle $\alpha$ with which this dipole begins to rotate. However, from a realistic perspective, this angle is difficult to be externally controlled. In this way, we also compute the averages of $\widehat{K}_{i}$ and $\langle z_{i}\rangle$ over all possible realizations of this angle $\alpha$. These $\alpha$-averaged quantities $\langle\widehat{K}_{i}\rangle_{\alpha}$ and $\langle z_{i}\rangle_{\alpha}$ provide insight, as a function of the excess energy $\delta K$, into the global extent of the energy equipartition and the orientation of the dipoles, respectively.

Let us start analyzing the field-free system, i. e., $\beta=0$. The dynamics depends on $\delta K$, which is the excess kinetic energy given to the system, and on the angle $\alpha$. The normalized time-averaged kinetic energy of the dipole one is presented in Fig. 2 for $0\leq\delta K\leq 15$ and $\alpha\in[-\pi,\pi)$. According to this color map, the kinetic energy transfer between the dipoles does not depend on $\alpha$, that is, on the direction in which the first dipole begins to rotate. In the initial state both dipoles are aligned along the symmetry axis of the field-free system, the LFF $X$-axis, which is the generatrix of the invariant manifolds ${\cal M}$. Since the excess kinetic energy is added to one of the dipoles, the system will be always moving on one of these invariant manifolds ${\cal M}$. As the dynamics on these manifolds ${\cal M}$ is equivalent, due to the rotational symmetry around the LFF $X$-axis, the direction in which dipole one starts to rotate has no impact on the global dynamics. Thus, the normalized time-averaged kinetic energy is independent of $\alpha$, as it is shown in Fig. 2.

For an excess energy $\delta K$ smaller than the critical value $\delta K_{c}\approx 6$, Fig. 2 shows that the system always reaches the equipartition energy regime, $\langle K_{1}\rangle$ is very close to $\langle K_{2}\rangle$, and a continuous energy flow between the rotors occurs. For $\delta K\approx 6$, this equipartition regime abruptly breaks, so that most of the kinetic energy remains always in dipole one for $\delta K\gtrsim 6$. As a consequence, the equipartition energy regime inside the invariant manifolds ${\cal M}$ applies only for low values of the excess energy $\delta K<\delta K_{c}$. This feature was detected in Ref. jonge . The authors provided a dynamical explanation of this phenomenon in rotors2D when they studied the energy transfer only in the invariant manifold ${\cal M}_{XZ}$ given by II.2. Let us emphasize that taking into account all the invariant manifolds ${\cal M}$, Fig. 2 proves that the system exhibits the same feature in all of them, as they are dynamically equivalent due to the field-free axial symmetry.

For a weak electric field, the electric field interaction can be considered as a perturbation to the dipole-dipole one. Indeed, for $\beta=0.1$, the normalized time-averaged kinetic energy of dipole one $\widehat{K}_{1}$ in Fig. 3(a) shows qualitatively a similar behaviour as the field-free results of Fig. 2, except for a region around the critical value $\delta K_{c}=6$. Thence, for $\delta K\lesssim 6$, the system reaches the equipartition regime, while for $\delta K\gtrsim 6$ most of the kinetic energy remains stored in dipole one. It is worth noticing that, in the neighborhood of $\delta K_{c}\approx 6$, the time-averaged $\widehat{K}_{1}$ shows sudden (irregular) variations, which indicate that in this region the system is very sensitive to the initial conditions, i. e., to the values of $\delta K$ and $\alpha$. The evolution of $\langle z_{1}\rangle$ depicted in Fig. 3(b) indicates that this weak electric field $\beta=0.1$ does not cause any significant orientation in dipole one. In contrast, for $\delta K\gtrsim 6$, the color map of Fig. 3(c) shows regions of initial conditions leading to a moderate orientation of the dipole two along the electric field with $\langle z_{2}\rangle\approx 0.5$, see the light brown regions of Fig. 3(c). A correlation between the behaviors of $\widehat{K}_{i}$ and $\langle z_{i}\rangle$ arises when the $\alpha$-averaged $\langle\widehat{K}_{i}\rangle_{\alpha}$ and $\langle z_{i}\rangle_{\alpha}$ are analyzed in Fig. 3(d). While the system is in the equipartition region, i. e., $\delta K\lesssim 6$, the orientation of both dipoles is negligible $\langle z_{i}\rangle_{\alpha}\approx 0$. However, in the region $\delta K\gtrsim 6$ where most of the kinetic energy is in dipole one, the $\alpha$-averaged orientation $\langle z_{2}\rangle_{\alpha}$ of dipole two begins to increase monotonically with $\delta K$, whereas the dipole one remains non-oriented $\langle z_{1}\rangle_{\alpha}\approx 0$. This behavior is expected because in the non-equipartition regime most of the kinetic energy is stored in dipole one, which swings very fast compared to dipole two. Thence, dipole two is more likely to be oriented along the field than dipole one.

Now, we analyze the dynamics when the dipole-dipole and the electric field interactions are of the same order of magnitude, specifically we take $\beta=1$. The corresponding results are shown in Fig. 4. Roughly speaking, for $\delta K\lesssim 8$ the system eventually relaxes to equipartition, see Fig. 4(a). Whereas, for $\delta K\gtrsim 8$ the system progressively moves away from equipartition, and as $\delta K$ increases, most of the kinetic energy remains in the initially excited dipole one. For $\delta K\gtrsim 8$, it is worth noticing the remarkable complex behavior observed in the evolution of $\widehat{K}_{1}$ in Fig. 4(a), i. e., the dynamics shows a high sensitivity to the initial conditions $\delta K$ and $\alpha$. This feature will be addressed in the next section. For $\delta K\lesssim 8$, there is no remarkable difference in the orientations of the dipoles in Fig. 4(b)-(c). We only highlight the slight orientation of both dipoles for small values of $\delta K$, see the light brown regions in Figs. 4(b)-(c) for $\delta K\lesssim 2$. By further increasing the excess energy, $\delta K\gtrsim 8$, the orientation of each dipole evolves differently. While the orientation of dipole one slightly decreases for increasing $\delta K$, see Fig. 4(b), there are regions in Fig. 4(c) where the orientation of the dipole two presents a significant increase. Furthermore, Figs. 4(a) and (c) show similar patterns, which confirm a clear correlation between the non-equipartition regime with most of the kinetic energy in dipole one, and dipole two having a larger orientation.

The evolutions of the $\alpha$-averaged $\langle\widehat{K}_{1}\rangle_{\alpha}$ and $\langle z_{i}\rangle_{\alpha}$ depicted in Fig. 4(d) show that, indeed, there is equipartition up to $\delta K\lesssim 8$. For this situation, $\langle z_{1}\rangle_{\alpha}$ and $\langle z_{2}\rangle_{\alpha}$ are equal and take the maximal orientation for small excess energy values, and they have very similar values and decrease monotonically for increasing $\delta K$. When the equipartition regime is lost for $\delta K\gtrsim 8$, the evolutions of the $\alpha$-averaged $\langle\widehat{K}_{1}\rangle_{\alpha}$ and $\langle\widehat{K}_{2}\rangle_{\alpha}$, red and blue solid lines in Fig. 4(d), respectively, indicate a gradual but non-monotonic growth of the kinetic energy of dipole one at expense of the one of dipole two. At the same time, $\langle z_{1}\rangle_{\alpha}$ continues decreasing for increasing $\delta K$, whereas, analogously to $\langle\widehat{K}_{1}\rangle_{\alpha}$, $\langle z_{2}\rangle_{\alpha}$ non-monotonically increases.

Now, we increase the electric field strength up to $\beta=10$, i. e., the electric field interaction is one order of magnitude larger than the dipole-dipole one, the corresponding results are presented in Fig. 5. In the color map of Fig. 5(a), we encounter that most of the initial conditions lead the systems to equipartition for $\delta K\lesssim 7$. Analogously to the $\beta=1$ dynamics, a more complex behavior appears for $\delta K\gtrsim 7$. Except for a region around $\alpha=\pi/2$, the energy transfer mechanism is dominated by a non-equipartition regime, and most of the initial excess energy $\delta K$ remains in dipole one. Although in this region around $\alpha=\pi/2$ the equipartition regime is dominant, there exist regions of initial conditions that lead to non-equipartitioning.

The time-averaged inclinations $\langle z_{1}\rangle$ and $\langle z_{2}\rangle$ shown in Figs. 5(b)-(c) present similar patterns as those observed for $\widehat{K}_{1}$ in Fig. 5(a). In the region dominated by energy equipartition $\delta K\lesssim 7$, the maps of Figs. 5(b)-(c) indicate that the dipoles are significantly oriented along the electric field axis. For $\delta K\gtrsim 7$, the orientation of dipole one in Fig. 5(b) decreases as $\delta K$ increases, and around $\alpha=-\pi/2$, there are initial conditions leading even to antiorientation, which correspond to this dipole having most of the energy in Fig. 5 (a). In contrast, the orientation of dipole two is globally enhanced in the non-equipartition region $\delta K\gtrsim 7$, see Fig. 5(c). Only for initial conditions in the region around $\alpha=\pi/2$ both dipoles have similar orientations slightly directed towards the electric field. In other words, when the dipole one begins to rotate in the direction of the electric field and equipartition is dominant, the dipoles end up having enough energy to prevent their orientation. That is, for initial conditions around $\alpha=\pi/2$, the non-equipartition regime existing in the absence of the electric field $\beta=0$ for $\delta K\gtrsim 6$, see Fig. 2, is mostly broken when the field is switched on. In contrast, for initial conditions out of that region, the system tends to remain in the same non-equipartitioning regime exhibited in the absence of the electric field, with most of the energy in dipole one. Thus, the electric field interaction is still not able to relax the system out of the non-equipartition regime, and, therefore, only the second dipole is significantly oriented by the electric field.

Again, the $\delta K$-evolution of the $\alpha$-averaged quantities in Fig. 5(d) confirms the behavior of the time-averaged ones presented in Figs. 5(a)-(c). During the equipartition regime for $\delta K\lesssim 7$ the dipoles have similar orientation. This orientation is maximal for small values of $\delta K$ and decreases as $\delta K$ increases. For $\delta K\gtrsim 7$, $\langle z_{1}\rangle_{\alpha}$ continues decreasing because most of the energy is stored in dipole one, whereas $\langle z_{2}\rangle_{\alpha}$ grows slightly for increasing $\delta K$.

By further increasing the electric field to $\beta=100$, the dynamics is dominated by the interaction with the electric field. In Fig. 6(a), the normalized time-averaged kinetic energy $\widehat{K}_{1}$ shows a complex behavior strongly depending on the initial excess energy $\delta K$ and on the angle $\alpha$. The map of Fig. 6(a) shows wide regions of energy equipartition. Interestingly, there is a blue-colored (brown-colored) area centered around $\alpha=\pi/2$ ($\alpha=-\pi/2$) where most of the kinetic energy is stored in dipole two (one). This means that when dipole one starts to rotate on an axis perpendicular to its $x_{1}z_{1}$ plane, the equipartition regime is not reached. For this strong electric field, a substantial orientation of the dipoles is expected, which is confirmed by the large values of $\langle z_{1}\rangle$ and $\langle z_{2}\rangle$ in Fig. 6(b)-(c). Regardless of this orientation, Figs. 6(b) and (c) present similar structures as Fig. 6(a), as was previously observed for $\beta=1$ and $10$. Thus, the non-equipartition regions in Fig. 6(a) with most of the kinetic energy located in one of the dipoles correspond to the lighter and darker red colored regions of Figs. 6(b)-(c), where the orientation is minimal and maximal, respectively. The $\alpha$-averaged $\langle\widehat{K}_{1}\rangle_{\alpha}$ and $\langle\widehat{K}_{2}\rangle_{\alpha}$ of Fig. 6(d) confirm that the system is always near equipartition, with a small positive energy balance for dipole one. Due to the strong electric field $\beta=100$, $\langle\widehat{K}_{1}\rangle_{\alpha}$ and $\langle\widehat{K}_{1}\rangle_{\alpha}$ show just a slight variation as $\delta K$ increases. A similar situation is found for the $\alpha$-averaged $\langle z_{i}\rangle_{\alpha}$, where the large orientation induced by the field is only slightly counteracted for increasing initial energy $\delta K$. Again, the less energetic dipole two is more oriented than the more energetic dipole one.

When the electric field interaction is three order of magnitude larger than the dipole-dipole one i. e., $\beta=1000$, most of the initial conditions leave the system close to the equipartition regime see Fig. 7(a). Moreover, the normalized time-averaged kinetic energy $\widehat{K}_{1}$ and the orientations show similar behaviours as those for the $\beta=100$ case, compare Figs. 6(a)-(b)-(c) and Figs. 7(a)-(b)-(c). We again encounter the non-equipartition regions around $\alpha=\pi/2$ and $-\pi/2$, where most of the kinetic energy is stored in dipole two and one, respectively. Thus, the energy equipartition regime is closely related to a smaller orientation of dipole two and viceversa. However, due to the strong electric field the two dipoles are significantly oriented along the electric field direction, see Figs. 7(b)-(c). The $\alpha$-averaged $\langle\widehat{K}_{i}\rangle_{\alpha}$ and $\langle z_{i}\rangle_{\alpha}$ of Fig. 7 assert that the system is always close to the equipartition regime, whereas the strong field orients the dipoles with $\langle z_{i}\rangle_{\alpha}>0.8$ with $i=1,2$.

IV Regular and chaotic dynamics

In this section, we study the chaoticity of the dynamics and its possible relation to the time-averaged energy transfer. The chaotic character of an orbit in a dynamical system is related to its sensitivity to the corresponding initial conditions. As a numerical measure of this sensitivity, we use the so-called, Fast Lyapunov Indicator (FLI) Froeschle1 ; Froeschle2 ; Fouchard to determine the degree of chaoticity of an orbit. Given a $n-$dimensional flow of a dynamical system

[TABLE]

the time evolution of the variational vector $\delta{\bf r}(t)$ is provided by the variational equation

[TABLE]

For a given orbit with initial conditions ${\bf r(0)}$ and $\delta{\bf r(0)}$, the numerical integration of equations (20) and (21) up to a final time $t_{f}$ yields the value of the FLI of that orbit defined by

[TABLE]

The variational vector $\delta{\bf r}$ increases linearly with time for regular periodic and quasiperiodic orbits, whereas it increases exponentially for chaotic orbits.

We have calculated FLI for varying initial conditions $\delta K$ and $\alpha$. To do so, we have integrated the equations of motion and their first variational equations in Cartesian coordinates by means of the Störmer-Verlet algorithm Andersen using the same initial conditions as in the previous section. Now, we stop the integration if the FLI reaches the saturation value nine, that characterizes an orbit as chaotic, or if the integration time reaches the stopping value $t_{f}=3000$. Our numerical tests have shown that this stopping time and saturation limit are adequate to correctly characterize any orbit.

The FLI color maps are presented in Fig. 8 for the electric field strengths $\beta=0,0.1,1,10,100$ and $1000$, the initial kinetic energy excess $\delta{\cal K}\in[0,15]$ and the angle $\alpha\in[-\pi,\pi)$. It is worth noting that, along the horizontal lines $\alpha=0$ and $\alpha=\pm\pi$ in Fig. 8, the system evolves on the invariant manifold ${\cal M}$ given by the planes $y_{1}=y_{2}=0$, i. e., both dipoles are restricted to rotate in the common $x_{1}z_{1}\equiv x_{2}z_{2}$ plane.

In the absence of the electric field, $\beta=0$, we observe in Fig. 8(a) that the FLI values do not depend on the angle $\alpha$, that is, on the direction in which the first dipole begins to rotate. As explained above, this is because regardless of the value of $\alpha$, the system will always move in one of the equivalent invariant manifolds ${\cal M}$. Moreover, for zero electric field, the FLI value increases with the kinetic energy excess $\delta{\cal K}$, although the dynamics is regular for all initial conditions. By switching on the electric field, most of the dynamics becomes chaotic even for weak electric fields such as $\beta=0.1$, see Fig. 8(b). For electric field values up to $\beta=10$, the regular behavior only persists for low values of the kinetic energy excess $\delta{\cal K}<1$ and particularly around $\alpha=0$ and $\alpha=\pm\pi$, see Figs. 8(b-d). For stronger electric fields, the chaotic behavior of the system still persists in wide regions of the FLI maps Fig. 8(e-f) for $\beta=$100 and 1000 and for the range of $\delta K$ values shown. This is somehow an unexpected result because for increasing field strength the dynamics is gradually dominated by the interaction of the dipoles with the field. This would lead the system to a (quasi) integrable regime, where most of the dynamics would be regular. However, the dipolar interaction, although being small compared to the electric field interaction, it is still able to cause a substantial volume of chaotic motion. In this sense, we have calculated the FLI maps in the range $0<\delta K\leq 15$ for very large values of the field parameter, namely for $\beta=10^{5}$ and $10^{6}$ (see Fig.9). Indeed, we observe in the color maps of Fig.9. that, for increasing values of $\beta$, the regions of initial conditions leading to regular motion grow in size while the regions of initial conditions leading to chaotic orbits shrink. However, this global tendency to regular motion is very slow, which confirms the relevant role that the DDI is playing in the dynamics even for very large electric field values.

A comparison of the results obtained for the time-averaged energy transfer, see Figs. 3-7, with those for the degree of chaoticity presented in Fig. 8 does not show any significant correlation. This may be due to the fact that regular and chaotic orbits can achieve the same time-averaged kinetic energy. We have additionally computed another chaos indicator, the so called Smaller Alignment Index (SALI) Skokos1 ; Skokos2 , and the corresponding results are in perfect agreement with the FLI ones presented here.

V Conclusions

We have theoretically investigated two interacting classical dipoles fixed in space, which are described as rigid rotors, in the presence of a homogenous electric field. The two dipoles are initially at rest in the stable head-tail configuration. Initially, one of the dipoles is pushed out of this stable configuration by giving it a specific velocity in a direction perpendicular to the head-tail axis, and at the same time the electric field is switched on. The following classical dynamics is explored in terms of the energy transfer mechanisms between the two dipoles and their orientations along the electric field axis.

The field-free dynamics of this system was previously investigated in the invariant manifold ${\cal M}_{XZ}$ II.2 in Ref. rotors2D . Here, we have shown that the dynamics is independent of the initial rotation angle of the dipole because it is restricted to one of the invariant manifolds of the infinite set ${\cal M}$. Depending on the initial excess energy given to one of the dipoles, the system falls to either an energy equipartition regime or a non-equipartition one, and the dynamics is regular.

The classical field-dressed dynamics strongly depends on the electric field strength and on the initial conditions. For weak external fields, the dynamics is still dominated by the dipole-dipole interaction and the energy transfer dynamics resembles the field-free dynamics with the two dipoles having a small orientation. By increasing the electric field, the interaction with this field dominates the classical dynamics. The size of the non-equipartition region appearing for large values of the initial excess energy decreases as the electric field increases, and the system tends to the energy equipartition regime. In addition, the orientation of the two dipoles along the electric field direction increases. We note that when the initial momentum is almost parallel of antiparallel to the electric field direction, we observe that even for very strong electric fields, the non-equipartition behavior dominates the dynamics. Finally, we find that, even for large electric field strengths and for our considered excess energies, the system shows a highly chaotic behavior. This is a remarkable feature because for increasing electric field strength it would be expected that the system gradually tends to a (quasi) integrable state where the dynamics is fully dominated by the interaction of the dipoles with the field and where the mutual dipole interaction would be considered as a small perturbation. However, as we have observed for very large values of the electric field, the DDI is still able to generate significant regions of chaotic motion.

A natural continuation of this work would be the investigation of the energy transfer and the possible collective phenomena of systems with many 3D-dipoles including in particular, the study of a linear chain of 3D-dipoles.

Acknowledgements.

M.I. and J.P.S. acknowledge financial support by the Spanish Project No. MTM 2017-88137-C2-2-P (MINECO). R.G.F. gratefully acknowledges financial support by the Spanish Project No. FIS2017-89349-P (MINECO), and by the Andalusian research group FQM-207. This study has been partially financed by the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), Ref. SOMM17/6105/UGR.

Appendix A Stability, existence and bifurcations of the equilibrium points

Here, we discuss the existence, stability and energy of the equilibria of the systems formed by two dipoles in an external electric field. These equilibria are summarized in Table 1. We also analyze the bifurcations appearing between these equilibria as the electric field strength is increased.

A.1 Stability and existence of the equilibria for $\beta=0$

For zero electric field, $\beta=0$, the only isolated equilibria are $P_{1}$ and $P_{5}$. Their nature is deduced by applying the method of Lagrange multipliers to the potential energy surface ${\cal V}_{1}(\mathbf{l}_{1},\mathbf{l}_{2},t=0)$  Eq. (6) with $\beta=0$, imposing the constraints $\{\mathbf{l}_{i}=(x_{i},y_{i},z_{i})||\,|\mathbf{l}_{i}|^{2}=1\}_{i=1}^{2}$. There exist the two stable equilibria $P_{1}$ with minimal energy $E_{1}=-2$, which are minima, and the two unstable equilibria $P_{5}$ with maximal energy $E_{5}=2$, which are maxima.

A.2 Stability and existence of the equilibria for $\beta\neq 0$ and $t\geq t_{1}$

For the general case $\beta\neq 0$, instead of using Cartesian variables, the analysis of the stability of the equilibria is carried out in spherical variables $\{(\theta_{i},P_{\theta_{i}}),(\phi_{i},P_{\phi_{i}})\}_{i=1}^{2}$. The equilibrium points of the Hamiltonian flux II.2 for $t\geq t_{1}$ are the critical points of the potential ${\cal V}_{2}$ [see Eq. (8)],

[TABLE]

together with the conditions $P_{\theta_{i}}=P_{\phi_{i}}=0$. Obviously, the potential ${\cal V}_{2}(\theta_{1},\phi_{1},\theta_{2},\phi_{2},t\geq t_{1})$ presents the same critical points as ${\cal V}_{1}(\mathbf{l}_{1},\mathbf{l}_{2})$ in Eq. (13). Thence, after expressing in spherical coordinates the (isolated) critical points $\{P_{i}\}^{6}_{i=1}$ detailed in subsection II.3, their stability can be inferred from the nature of the corresponding eigenvalues of the Hessian matrix associated to ${\cal V}_{2}(\theta_{1},\phi_{1},\theta_{2},\phi_{2},t\geq t_{1})$.

A.2.1 The equilibria $P_{1}$

For the two equilibria $P_{1}$, the eigenvalues of the Hessian matrix are:

[TABLE]

Because the equilibria $P_{1}$ only exist when $0<\beta<3$, its four eigenvalues are positive, which indicate that $P_{1}$ are minima. In other words, when they exist, equilibria $P_{1}$ are stable.

A.2.2 The equilibria $P_{2}$

The eigenvalues of the Hessian matrix of the two equilibria $P_{2}$ are:

[TABLE]

These eigenvalues indicate that $P_{2}$ are always rank-two saddle points because they have two positive and two negative eigenvalues. Then, equilibria $P_{2}$ are always unstable.

A.2.3 The equilibrium $P_{3}$

The eigenvalues of the Hessian matrix of the equilibrium $P_{3}$ are:

[TABLE]

The first and second eigenvalues are always positive. For $0<\beta<2$, the equilibrium $P_{3}$ is a rank-two saddle point because it has two positive and two negative eigenvalues. In the interval $2<\beta<3$, $P_{3}$ is a unstable rank-one saddle point (the third eigenvalue is positive). For $\beta>3$, all the eigenvalues are positive, e.g., the equilibrium $P_{3}$ is a (stable) minimum.

A.2.4 The equilibrium $P_{4}$

The eigenvalues of the Hessian matrix of the equilibrium $P_{4}$ are:

[TABLE]

The first, second and third eigenvalues are always negative, while the fourth one changes from positive to negative for $\beta>1$. Then, $P_{4}$ is an unstable rank-three saddle point when $0<\beta<1$, and a (unstable) maximum for $\beta>1$.

A.2.5 The equilibria $P_{5}$

The eigenvalues of the Hessian matrix of the two equilibria $P_{5}$ are:

[TABLE]

Because equilibria $P_{5}$ only exist for $\beta<1$, all the eigenvalues are negative, such that $P_{5}$ are always (unstable) maxima.

A.2.6 The equilibria $P_{6}$

The eigenvalues of the Hessian matrix of the two equilibria $P_{6}$ are:

[TABLE]

The first and second eigenvalues are positive, the third one is negative and the fourth one is positive for $0<\beta<2$. Because equilibria $P_{6}$ only exist for $0<\beta<2$, they are unstable rank-one saddle points.

A.3 Bifurcations

The bifurcations between the critical points are presented in Fig. 10 by the evolution of the equilibria energies as a function of the electric field parameter $\beta$. For the field-free case $\beta=0$, besides the (isolated) equilibria $P_{1}$ and $P_{5}$, which are respectively the head-tail ground state (minimum) configurations and the head-head and tail-tail maxima, the system presents the aforementioned two circles $D_{1,2}$ of degenerate equilibria, see subsection II.3. For $\beta>0$, the circles of equilibria blow-up, and from $D_{1}$ and $D_{2}$ emerge the two isolated equilibria $P_{3,4}$ and $P_{2,6}$, respectively. In the interval $0<\beta<1$, the study of the stability of the equilibria shows that $P_{1}$ and $P_{5}$ are the absolute minima and maxima of the system, see Fig. 10, while $P_{2,3,4,6}$ are saddle points of rank-one, rank-two, rank-three and rank-one, respectively. As $\beta$ increases towards $\beta=1$, the equilibria $P_{4}$ and $P_{5}$ approach each other see Fig. 10, such that, at $\beta=1$, a first Pitchfork bifurcation takes place: the two equilibria collide and only the equilibrium $P_{4}$ survives afterwards, becoming the equilibrium of maximal energy. This is the expected configuration of maximal energy where both dipoles are located along the negative $z_{1}$ and $z_{2}$ axes, i. e., oriented antiparallel to electric field direction.

In the interval $0<\beta<3$, the equilibria $P_{1,3,6}$ approach each other (see Fig. 10), and at $\beta=2$, equilibria $P_{3}$ and $P_{6}$ coincide. From this second pitchfork bifurcation, only equilibrium $P_{3}$ survives, becoming a saddle point of rank-one. Finally, at $\beta=3$, a third Pitchfork bifurcation between the equilibria $P_{1}$ and $P_{3}$ occurs such that, for $\beta>3$, only $P_{3}$ survives becoming the equilibrium of minimal energy. In this configuration of minimal energy, both dipoles are oriented along the field, i. e., along the positive $z_{1,2}$ axes. By further increasing $\beta$, $\beta>3$, the dynamics is gradually dominated by the interaction with the electric field, and the landscape of the potential energy surface ${\cal V}_{1}(\mathbf{l}_{1},\mathbf{l}_{2})$ resembles the one obtained by neglecting the dipole-dipole interaction.