Electromagnetic trap for polar particles

A novel design for an electromagnetic trap is proposed for confinement of neutral particles having a permanent electric dipole moment. The device uses a combination of a sextupole electric and quadrupole magnetic fields superimposed with a strong constant electric field perfectly aligned along the z-axis. We use the extended dipole model to study the dynamics of particle in this position dependent electromagnetic field. The motion of the centre of mass of the dipole is nonlinearly coupled with its rotation. We find three families of particular solutions for which the orbits are periodic or quasiperiodic. They correspond to the particular cases when the whole energy is accumulated in either translational or rotational degrees of freedom. Stability analysis of chosen particular solutions provides trapping conditions. We present several numerical simulations which illustrate trapping and confinement of an electric dipole in the proposed trap. These simulations were performed for three kinds of polar particles and various fields selections. We hope that the proposed model allows experimental physicists to apply a wide variety of non-destructive methods for manipulation, detection and analysis of trapped particles.


Introduction
Methods of storage and trapping of charged and neutral particles are extremely important for modern ultracold physics and its applications, including high-resolution spectroscopy. For charged particles the strong Lorentz force is used for trapping in harmonic electrostatic and magnetostatic fields [1][2][3][4]. The first electrostatic trap for polar molecules was presented in reference [5]. Stark-decelerated ND 3 molecules were loaded in the quadruple electrostatic field. The electrodes create potential well of depth 0.35 K which is much lower than in hyperbolic Penning traps. Manipulations of polar molecules by electric fields (focussing of molecular beams, deceleration of molecules, and their trapping) are described in [6] in many details.
In reference [7] the trap is designed which exploits the combination of quadruple magnetic field and quadruple electric field for confinement of neutral, ground state OH molecules. The authors use Stark decelerator to slow the beam to 20 m s −1 before it gets into processing chamber. They apply the magnetic field 67 T m −1 and voltage 9 kV. A longitudinal trap depth 400 mK is achieved.
In the recent paper [8], we study the dynamics of an electric dipole in a uniform and static electromagnetic field. We model the dipole as a fixed-length massless rod with two massive point charges at its two edges. We show that the motion of this rigid rotor as a whole heavily depends on the evolution of its electric dipole moment 3-vector. The orbits are unbounded. Thus, to trap the electric dipole effectively, the position dependent fields are necessary.
In the present paper we develop the concept of apparatus for confinement of neutral particles having a permanent electric dipole moment. By confinement we mean a stable configuration in which a particle is localised in a certain limited volume for a long time [2, ch 1.3]. The confinement is achieved by means of a combination of a quadrupole magnetic, sextupole electric and a strong homogeneous electric fields.
The sextupole potential of the shape (z 3 − 3 2 z(x 2 + y 2 )) produces a restoring force along the z-axis, and a repulsive one in the radial direction. Quadrupole magnetic field prevents the radial escaping.
An additional static homogeneous electric field along the z-axis pushes the electric dipole moment vector in the axial direction. If a single confined particle is localised in the plane z = 0 and has no initial axial velocity, and if the dipole moment 3-vector is perfectly aligned along the axial direction, the particle is forced to move on an epitrochoid [2, figure 5.2]. The trajectory of this specifically oriented dipole is just the same as the orbit of a charged particle in a hyperbolic Penning trap, see [1, figure 7], [2, figure 5.5] and [3, figures 3.1 and 3.2]. For a given voltage and magnetic field strength the critical mass, above which the dipole motion is unstable, is proportional to the magnitude of the electric dipole moment. Therefore, the most promising objects are polyatomic molecules and nanoparticles having large permanent dipole moments.
The quadrupole magnetic field can be realised by two orthogonal pairs of coils in anti-Helmholtz configuration. In [9] the superconducting magnet technology is described which produces an ideal quadrupole field of high gradient 250 T m −1 . The authors state that if iron yoke surrounding coils is used, the higher value 300 T m −1 is achieved.
The polar molecule or nanoparticle should be cooled to be trapped. For cooling of polyatomic polar molecules the optoelectrical Sisyphus method is elaborated in reference [10]. Monte Carlo simulations [11] for NH molecules predict that '60 percent of trapped molecules can be cooled from 100 mK to sub-mK temperature in a few seconds.' As technologies are in constant progression (both magnet and cooling), this stimulates the searching for new methods of trapping and confinement low-energy objects with permanent electric dipole moments.
The plan of the paper is the following. In section 2 we briefly review the recent developments in the field of polyatomic molecules and nanoparticles having large permanent dipole moments. In section 3 we derive the Lagrangian governing the dynamics of the electric dipole under the influence of an electromagnetic field. We start with the Lagrangian of two electric charges and derive the function governing the evolution of a linear rigid rotor (a rigid body with one vanishing principle moment of inertia) with zero net charge. Further we design the electromagnetic field that confines a compact dipole in a small volume for long periods of time. In section 4 we derive the Hamiltonian that governs the dynamics of the dipole.
In section 5 we study the non-linear rotation of an electric dipole whose centre of mass is placed at the coordinate origin. In this specific case the dipole rotates as a spherical pendulum. We investigate conical periodic solutions which describe uniform precession of the dipole around the z-axis. They are also called regular precessions. The stability analysis of these solutions was done analytically. It provides the trapping conditions which are expressed in terms of the precession angle, the particle individual characteristics and field strengths.
In section 6 we present the results of numerical simulations. The trajectories of the centre of mass of the dipole and the orbits of arrowhead of dipole moment 3-vector are presented. Section 7 contains the summary and remarks concerning quantum mechanical aspects of the proposed trap. In appendix A we design the shape of electrodes which produce the sextupole potential.

Polar polyatomic molecules and nanoparticles
The basic difference between dynamics of a point charged particles and extended particles is the presence of additional degrees of freedom connected with rotations and vibrations. Even if we limit ourselves to dipole-like particles with only additional rotational degrees of freedom, the dynamics becomes complicated and the aim of this paper is to give its description. Necessity of such analysis was announced by experimental physicists. Recently in [12] the extra fast rotation of optically trapped [13,14] silica nanoparticle has been observed. The rotation is due to transferring spin angular momentum of light to the particle's mechanical angular momentum. For a single 100 nm particle the frequency 1 GHz is achieved. The authors emphasise that the problem of interaction of the centre of mass and the rotational degrees of freedom is of the utmost importance. Obviously, the coupling can be studied in details using a dipole nanoparticle moving in an electromagnetic field. In our opinion, one of the most promising object is the polar DAST nanocrystal [15,16]. The experiments [17] with polar DAST nanocrystals acted upon external electric field demonstrate a large electric dipole moment. In [17, p 390] it was estimated that for 100 nm × 100 nm × 50 nm crystal the moment is 2.8 × 10 7 D and it is directed along the symmetry axis.
In [18] the permanent dipole moment of a cellulose nanocrystal is measured by using the method of transient electric birefringence. Its estimated value is 4440 D. A cellulose nanocrystal (CNC) is a rod of diameter range 3 nm-20 nm. The length is in the range from 100 nm to several microns. The authors believe that the large anomalous dipole moment originates from 'the parallel arrangement of cellulose chains in a non-centrosymmetric crystallographic lattice within each CNC together with the dipole moment borne by each glucosyl monomer'. Shanbhag and Kotov [19] calculated the size-dependent permanent dipole moments 50 D to 100 D of CdS nanocrystals with shape asymmetry by using the PM3 method. The core Cd 84 S 123 nanocrystal of tetrahedral symmetry edge range 3.3 nm-4.8 nm, see [19, figure 2], has no a dipole moment. To obtain nanocrystal with dipole moment, the core cluster is modified chemically by S-H terminal groups, its corners are truncated (a Cd-S-H group is eliminated from apex), and its surface is hydrated. Due to these manipulations the nanocrystal obtains the dipole moment of order 50 D to 110 D, see tables 1-3 in [19]). Authors claim that one of the effect of hydration is the formation of assembles of nanoparticles (nanorods, nanowires) due to dipole-dipole attraction. Their total dipole moments increase linearly. In chapter 20 of the monograph [17] the synthesis of Cd-S-nanoparticle dimers by linking of two nanocrystals with a methylene group is described.
In [20] the large permanent electric dipoles of CdSe nanorods are observed. Their magnitudes range in 50 D-210 D, depending on the length and diameter of a nanoparticle [20, table 1]. The nanorods are coated with organic ligands [17, ch 15].
In the review papers [21,22] the design and synthesis as well as the properties and application of Janus micro-and nanoparticles have been considered. These type of particles have two or more parts with distinct physical properties. In particular, Takei and Shimizu [23] design particles such that positive charges are separated on one hemisphere, and negative charges on the other one. The authors observe the rotation of such spheres in an electric field due to their dipole moments.

Classical theory of the electromagnetic trap for electric dipole
We consider a polar molecule or nanoparticle with zero net charge and a permanent dipole moment. We approximate it as a dipolar rigid rotor consisting of two equal and opposite point charges, q 1 = +q and q 2 = −q, with masses m 1 and m 2 , respectively. They are linked by a massless rigid rod of a fixed length l, and are located in a stationary electromagnetic field. The Lagrangian of the system is where x 1 and x 2 are radii vectors of charges. Φ(x) and A(x) are the scalar and vector potentials of an external electromagnetic field at point x. The charges form a physical dipole which is characterised by its dipole moment d pointing from the negative charge to the positive charge. Its magnitude is d = ql, and d = de, where e is the unit vector. We did not include terms describing mutual interactions between charges q 1 and q 2 . A configuration of the dipole is fully specified by the position of its centre of mass x = (x, y, z), and by the direction of its dipole moment e = (e 1 , e 2 , e 3 ). Thus, the configuration space is R 3 × S 2 , where S 2 denotes the unite sphere.
As mx = m 1 x 1 + m 2 x 2 , and where m = m 1 + m 2 , thus we have and we can write Lagrangian L in terms of (x, e). The length of dipole l is very small, so it is justifiable to expand the potential into power series with respect to l up to a certain order. In further considerations it is necessary to keep terms up to the second order. The truncated Lagrangian has the form where E(x) = −∇Φ(x) and B(x) = ∇ × A(x) are the electric and magnetic field strengths, respectively. By E (x) and B (x) we denote the derivatives of these vector fields. They are 3 × 3 matrices: . We introduce also two parameters: the reduced mass of the dipole μ = m 1 m 2 /m, and parameter which measures the asymmetry of mass distribution.
In the final form of the Lagrangian function (4) we omit terms which are the total time derivative of de · A(x) and the total time derivative of 1 2 dδe · (A (x)e) because they do not change the equations of motion. In the first order approximation (ideal dipole) our expression coincides with that presented in reference [27, equation (6)].
Considering real molecules and particles we have to take into account asymmetry of their mass and charge distribution, which we parametrize by δ. Examples of such asymmetric particles are the water molecule H 2 O and nanoparticle H 54 Cd 80 S 112 which has a shape of deformed tetrahedron [19]. The asymmetry is inherent for Janus particles [21, figure 2], e.g. the hexafluorocyclohexane molecule [28].

Electromagnetic field configuration
Guided by the Penning trap design we assume that the electromagnetic field is axially symmetric and the axis of symmetry coincides with the third axis of the laboratory inertial frame which has the origin at the geometrical centre of the trap. The field has to confine an electric dipole which is modelled as two separated charges of opposite signs. The net force acting on the centre of mass of the dipole caused by a constant electric field vanishes. It is known that an electric dipole can be trapped by a constant magnetic field, however it can occur for very specific initial conditions, see [29][30][31]. A classical electric dipole moves in the presence of an external magnetic field, perpendicular to the dipole's plane of motion, and the integrals of motion should satisfy additional condition. We do not restrict ourselves to this specific configuration. Thus we have to look for coordinate-depended fields.
To bound the dipole in a small region near the centre of the trap both the electric and magnetic fields are necessary. The role of the electric field is to confine the motion of the dipole along the symmetry axis. Unfortunately, it is repulsive in the plane perpendicular to this axis. The magnetic field bends the trajectory of the dipole and enforces it to move around the symmetry axis. It is worth noting, that the described trapping mechanism is similar to that in the Penning trap although the forms of appropriate harmonic fields differ from that in Penning trap. However, to confine the motion of a dipole we have to take into account fact that the polar particle rotates and its rotation influences the motion of its centre of mass. For example, if the dipole rotation axis is parallel to the symmetry axis, then additional force acts along the symmetry axis. It causes a drift of the polar particle out the trap. To prevent this effect we add an appropriate homogeneous electric field which fixes the dipole moment orientation.
For a better understanding of the trapping mechanism we write down the Newton equations of motions for the centre of mass The electromagnetic field is a superposition of harmonic electric and magnetic fields which are considered in Penning traps as imperfections, see [2, ch 6] and [32]. We take the ideal quadrupolar magnetic field [9] B with high gradient B 1 = 250 T m −1 . The magnetic field lines distribution is shown in figure 1. Let us consider ideal nominal situation when the electric dipole moment is parallel to the symmetry axis, i.e., e = (0, 0, 1), and its centre of mass moves in radial plane with velocityẋ = α(x, y, 0), where α is a real parameter. Then, the Lorentz force due to magnetic field alone is The force bends the trajectory because it is perpendicular to the particle's velocity. More importantly, this force is perpendicular to the symmetry axis, so it will not cause the motion along this axis. If e is perpendicular to the axis of symmetry and rotates with angular velocity ω which is approximately constant, i.e., e = (cos(ωt), sin(ωt), 0), then So, in this case the Lorentz force will push the dipole out of the z = 0 plane. According to 6, if the electrostatic field E(x) = −∇Φ(x) acts alone, then the force acting on the dipole is F e = dE (x) · e. To obtain the restoring force along the z-axis we take the sextupole potential  3 where the electrostatic potential Φ 3 is given by equation (10) and E z is the magnitude of the constant homogeneous electric field.
where U 0 is trapping voltage and D is assumed to be the characteristic trap size defined as In this formula a and b are minimal distances from the centre of the trap to endcap electrodes, see appendix A and in particular figures 16 and 17. The corresponding electric field is given in figure 2(a). With this potential the force due to electrostatic field alone acting on the dipole is For the nominal orientation of the dipole e = (0, 0, 1) it becomes So, it is repulsive in the (x, y) plane and it attracts to the origin along z as desired, see figure 3. However, the field given by potential Φ 3 (x) does not satisfy all of our requirements. Namely, it does not guarantee that the axis of the dipole will be close to its nominal position when e is almost parallel to the z axis. To ensure this we add a strong constant electrostatic field along the z axis. Thus, the final form of the potential is The electric field resulting from the potential is displayed in figure 2(b). The sextupole potential Φ 3 (x) can be produced by ideal electrodes which have the shape of the equipotential surfaces, see figure 4(a). The electrodes are shown in figure 4(b): positive are coloured in blue and negative in orange. More details about the shape of electrodes generating this potential are given in appendix A.
Summarising the configuration of our trap is following. The strong constant electrostatic field E z e z is generated by flat electrodes. Between them are placed the coils in anti Helmholtz configuration producing the quadrupole magnetic field B(x) and the hexagonal electrodes generating the sextupole potential Φ 3 (x). The region between the electrodes in which act the electric and magnetic fields is called the processing chamber. Its size is described by the characteristic trap size D.
In our further considerations we normalise the dipole radius vector by the factor D, i.e. x → x/D. Thus, we parameterize the confinement space by dimensionless variables. We introduce also the dimensionless evolution parameter τ = ω c t, where cyclotron frequency depends on the magnetic field constant B 1 , the dipole mass m and magnitude of its moment d: In the rescaled variables the Newton equation (6) can be rewritten in the following form where and In the subsection below we describe the trapping mechanism of a polar particle with perfectly oriented electric dipole moment.

Certain particular solutions
We already mentioned that the desired orientation of the polar particle is such that its dipole moment is parallel to the symmetry axis of the trap. The question arises whether it is possible to keep the orientation of the dipole vector e. More precisely, is it possible that the centre of mass of the dipole moves in coordinate dependent electric field and e = (0, 0, 1)? Clearly, the necessary condition for this kind of motion is that the total torque acting on the dipole vanishes. This torque depends on the electric field E(x) as well as on the magnetic field B(x). Taking into account the axial symmetry of these fields we claim that the total torque acting on the dipole vanishes in two cases (a) the centre of mass moves along the z-axis and the dipole vector is parallel to it; (b) the dipole is symmetric i.e. δ = 0, the centre of mass moves in the (x, y)-plane and the dipole moment vector is parallel to the z-axis.
In the first case we assume that e(t) = (0, 0, 1) and x(t) = (0, 0, z(t)). The function z(t) satisfies the Newton equation (15) which takes the formz Hence, the centre of mass oscillates harmonically with axial frequency ω z = κ/2 about the origin along the z-axis, as the described above mechanism predicts.
Equations (19) are equivalent to the canonical equations They are generated by the following Hamilton function The stability of linear Hamiltonian equation (20) is determined by eigenvalues of the matrix of coefficients.
The characteristic polynomial of this matrix is The necessary condition of stability is that all eigenvalues of the matrix of coefficients are purely imaginary numbers, see e.g. [33, p 31]. If the matrix is diagonalizable, this condition is also sufficient. Substitution (22) gives Requirement that two roots w 1 , w 2 are real and distinct implies that the discriminant of the quadratic equation 1 − κ > 0, that is κ < 1. Under this assumption both roots are positive because w 1 + w 2 = 1 2 (2 − κ) > 0 and w 1 w 2 = κ 2 16 > 0. As result we got the stability condition exactly the same as the trapping condition for the Penning trap, see e.g. [34, equation (3.4)]. Thus κ defined by (17) can be called the trapping parameter. For a given magnetic field strength B 1 the maximum voltage U max follows from the upper stability limit given by condition κ = 1: Both the cases describe idealized situation when the rotation of dipole moment vector is suppressed completely. It is impossible to achieve this in real traps. Exemplary projection of a quasi-periodic trajectory of the centre of mass obtained from numerical simulations for DAST nanocrystal is shown in figure 5. In fact this is the projection of the trajectory in figure 11 (a) onto the (x, y)-plane. For details of these simulations, see section 6. In the next section we describe the real dynamics taking into account all motions of the dipole: its translations and rotations.

Lagrangian and Hamiltonian formulations
In this section we present Lagrange and Hamiltonian functions of the rigid rotor model of the electric dipole in the considered electromagnetic trap. Two coordinate systems for the location of the centre of mass of the dipole are used: the Cartesian and the cylindrical ones. The corresponding Lagrange and Hamilton equations will be studied analytically and numerically in the next two sections. To simplify further considerations we will use the normalized coordinates and time. Now, the parameter (5) transforms into its dimensionless counterpart: We parameterize the components of unit vector e directed along d by spherical coordinates: The truncated Lagrangian (4) takes the form In the formula for L cm symbol R(ϕ) denotes the block diagonal matrix which defines the rotation of the centre of mass variables and their derivatives y T = x, y, z,ẋ,ẏ,ż by the dipole azimuth angle ϕ The 6 × 6 matrix L 6 is where c := cos θ and s := sin θ. The rotational part L d of the total Lagrangian consists of the kinetic energy and the term describing the influence of constant electric field: where The standard Legendre transformation yields the Hamilton function where z T = x, y, z, p x , p y , p z and As both the sextupole electrostatic potential and the quadrupole magnetostatic field are axially symmetric the cylindrical coordinates (ρ, φ, z) relative to geometrical centre between electrodes are a good choice. In this parametrization the dimensionless coordinates of the centre of mass position x are as follows: In these new variables the Lagrangian will be denoted by L. It can be written in the following form with where χ = φ − ϕ, and L d = L d . As Lagrangian depends on the difference χ = φ − ϕ, we can take (ρ, z, χ, ϕ, θ) as the generalised variables. Then angle ϕ is a cyclic variable and its conjugated momentum is a first integral. Using it we reduce the number of degrees of freedom by one. Corresponding Hamilton function takes the form where H d = H d and sin(2θ)(p χ cos χ + ρp ρ sin χ) Remark 1. The only drawback of the introduced coordinates is that they are singular for ρ = 0, or for θ = 0, that is when the centre of mass moves along the z-axis, or when electric dipole moment is aligned along it. For study these cases we can take variables (x, y, z) and modify polar coordinates e 1 = sin ϕ sin θ, e 2 = cos θ, e 3 = cos ϕ sin θ.

Dipole regular precessions and trapping conditions
In this section we distinguish a certain family of exact solutions which describe evolution of the electric dipole vector when its centre of mass is at rest exactly in the centre of the trap. If the dipole moves near the centre of the trap, then its motion can be described by linear variational equations obtained by small perturbations of this exact solution. Stability analysis of specific periodic solutions of this type allows us to derive a trapping condition.

Spherical pendulum solutions
Let us assume that the dipole is symmetric, i.e. δ = 0. Then equations generated by Hamiltonian (33) admit a family of solutions for which x(t) = 0 and p(t) = 0, where p = (p x , p y , p z ). For these solutions the centre of mass of the dipole stays permanently at the origin and its rotations are described by the following equationsθ They are the canonical equations generated by Hamiltonian H d written in (36). These equations coincide with those of spherical pendulum [35, section 16.10]. It is well-known integrable dynamical system and its mathematical properties have been examined in very details. The expressions for exact solutions had been written in the work by Tissot [36], and in [35, section 16.10].
It is worth noting that p ϕ is a first integral of the restricted system (45) but it is not a first integral of Hamilton's equations generated by Hamiltonian (33).
As p ϕ is constant for system (45), we fix its value p ϕ = p and further consider p as a parameter. Let us rewrite the Hamiltonian H d in the form where V 0 (θ) = 1 2η The shape of this potential is shown in figure 6. Its global minimum is the root of the trigonometric equation (48) Figure 6. Graph of the potential V 0 (θ).
Setting x = −cos θ we obtain the quartic equation It is easy to see that it has only one negative root for an arbitrary positive λ. We can express it in terms of elementary functions using the Ferrari method. It yields the only suitable root where the auxiliary function σ(λ) is The minimal energy is E min = V 0 (θ min (λ)) = γ λ sin 2 θ min (λ) − cos θ min (λ) .
If the energy is greater than E min , the polar angle θ oscillates in the infinite potential walls located at points θ = 0 and at θ = π, see figure 6. A non-zero electric field E z > 0 shifts the critical point of the potential from median angle π/2 to the left wall. Moreover, the larger is E z , the deeper is the minimum.
The qualitative properties of the motion depend on the value of energy. For future convenience we rescale the dimensionless time variable τ = τ 2γ/η. We introduce new momenta π θ = p θ / √ 2ηγ, π ϕ = p ϕ / √ 2ηγ, and new Hamiltonian H d = 1 2 H d /γ: Let us fix the value of energy H d = /2. Now π 2 ϕ = λ and π θ =θ. Setting x = −cos θ in the energy integral, after some algebra we obtain the equatioṅ The form of solution of this equation depends on roots of polynomial P( A real motion is possible if all three roots x k are real. Thus, the Cardano discriminant of this polynomial must be non-negative. It can be written in the form where If Δ C > 0, the cubic has three distinct real roots. If Δ C = 0, the roots are real and multiple. The region in ( , λ) plane where Δ C 0 is shown in figure 7 and it is defined by As a stability investigation of first two types of solutions is rather complicated and can be performed only numerically, it will be presented in another paper. For further analysis we take the third type of solutions which have the following form where ω 2 = γ η cos θ min (λ) .
We call it the conic rotation or regular precession.

Linearized system and stability
Let us order the canonical variables in the following way and let X 0 (τ ) denotes a particular solution of Hamilton's equationṡ where H is given by (33) and H (X) denotes the gradient of H(X) written as one column matrix. In the above formula the matrix J is a block diagonal matrix where J n denotes n × n symplectic unit matrix. In order to obtain variational equations we set X = X 0 (τ ) + Y and expand the right hand sides of obtained equations into the power series with respect to Y. Neglecting terms of order higher than one we obtain the following linearized systeṁ where H (X) denotes the Hessian matrix of H(X). For further calculations as X 0 (τ ) we take solution given by (57). For this choice matrix H (X 0 (τ )) has the form H (X 0 (τ )) = S(ωτ ) T HS(ωτ ), where S(ωτ ) is a block diagonal matrix of the form Here I 4 is 4 × 4 identity matrix, and with H 6 is given by (35) and In order to remove the time dependence of A(τ ) we change variational variables Y → Z = S(ωτ )Y. It is equivalent to the passage to the rotating coordinate frame. The linearized system takes the forṁ where we used the identity S(ωτ )JS T (ωτ ) = J. Resulting time-independent matrix B has block diagonal structure B = diag(B 6 , H 4 ), where B 6 = Ω 6 + H 6 with Matrices H 6 and H 4 are defined in (35) and in (66), respectively. Finally, linear system (67) splits into two independent subsystemsŻ for perturbations of angular variables and corresponding momenta Z 1 and for perturbations of the centre of mass and corresponding momenta variables Z 2 , Z = (Z 1 , Z 2 ). The problem of linear stability of the dipole dynamics splits into stability analysis of the rotational dynamics of dipole in linear approximation described by (68) and translational dynamic of the centre of mass determined by linear system (69). The equilibrium Z = 0 of the Hamiltonian system of linear differential equations (67) is stable if and only if all eigenvalues of matrices J 4 H 4 and J 6 B 6 are pure imaginary. In order to study linear stability we analyse characteristic polynomials of these matrices, see e.g. [33, p 31]. Characteristic polynomial of J 4 H 4 has the form We assume that parameters γ, η, ω are positive and 0 θ < π/2 that implies c > 0. The non-zero eigenvalues are given by a pair of complex conjugate purely imaginary values of λ for arbitrary positive values of γ, η, ω and c. But double zero eigenvalue corresponds to two-dimensional Jordan block. Thus the uniform precession of the electric dipole moment is not linearly stable in the Lyapunov sense. However, the linear stability of the motion of the centre of mass does not depend on the stability of rotation. In order to study nonlinear stability we need to take into account mutual interaction between rotations and the motion of the centre of mass of the dipole. We plan to study this problem in our forthcoming paper. The characteristic polynomial of product of matrices J 6 B 6 has the form with coefficients The coefficients of polynomial p 2 (λ) are functions of trapping parameter κ > 0, cosine of polar angle θ denoted by c, and the frequency of dipole rotations ω. The stability conditions of the motion of the centre of mass can be expressed in terms of roots of polynomial p 2 (λ). Substituting λ = iΩ, where Ω ∈ R into p 2 (λ) and introducing the variable where polynomial coefficients are given by (71). To fulfil the stability condition, the roots w 1 , w 2 , and w 3 of the spectral polynomial (72) should be distinct, real, and positive. According to the Vieta formulae, the sum of roots is w 1 + w 2 + w 3 = a 4 , where a 4 > 0 for an arbitrary values of ω and c. The product of roots w 1 w 2 w 3 = a 0 should also be positive. But in fact a 0 > 0 because both factors of a 0 are positive, that is visible if we recall that 0 c 1. The expression w 1 w 2 + w 1 w 3 + w 2 w 3 = a 2 produces condition: a 2 > 0 that gives non-trivial restriction on our parameters. The next condition is due to the assumption that the cubic polynomial p 2 (λ) possesses three distinct real roots. This holds if and only if the Cardano discriminant is positive. We do not write it explicitly because of its very complicated form. Thus the region of stability is defined by two conditions Three dimensional plots of this region are presented in figure 8.

Numerical simulations
In this section we present results of numerical integration of Hamiltonian equations generated by Hamiltonian (42) for DAST and cellulose nanocrystals and for the protein lac repressor. Data for these particles are given in table 1I. In the case of the DAST nanocrystal it is assumed that it is perfectly shaped with δ = 0. Results of these simulations are trajectories of the centre of mass x = (x, y, z) and time evolution of the dipole vector e = (e 1 , e 2 , e 3 ). Calculations are performed for three different values of the initial azimuth angle:     In physical units τ = 1 corresponds to 27.7393 s. The polar particle for τ = 0 is localised at the distance r = 5 × 10 −6 m near the centre of the trap.  If p ρ (0) = 0 and/or p z (0) = 0, then the evolution of the centre of mass changes. It moves in a much larger volume. In the case when only p z (0) = 0, the range of motion is extended along z-axis. If p ρ (0)p z (0) = 0, then extension is in all three dimensions. But even in these cases the orbits look like as three-dimensional Lissajous curves with various frequencies and relative phases depending on initial conditions. In order to get better idea about the depth of the trap it is useful to define the translational temperature T trans as 1 2 mv 2 = kT trans , where v is the velocity of the centre of mass and k is the Boltzmann constant. It changes with time so for a given orbit we give its maximal value. In table 2 we specify also, the potential depth T = mω 2 c D 2 /2k. In figure 11 we show result using DAST I data but with changed initial conditions of p ρ (0) and p z (0). We report here only few results of our numerical experiments. However, all of them demonstrate good stability properties of the proposed trap. Particles in trap can be confined for hours or even days.
For numerical integration we use procedures from GNU Scientific Library [37]. We used a procedure which implements a variable-coefficient linear multistep Adams method in Nordsieck form. The relative precision of integrations was set in the range 10 −10 -10 −14 .

Conclusions and remarks
In this paper we present a new design of electromagnetic trap for a polar particle. We propose a special configuration of electromagnetic field. A neutral particle with a non-zero dipole momentum is trapped in the combination of quadrupole magnetic field and electric field derived from sextupole electrostatic potential. The coordinate dependent fields are superimposed by the static homogeneous electric field, perfectly aligned along the z-axis. They are generated by electrodes and coils located outside the processing chamber where a particle is confined. The benefits of homogeneous electrostatic field are twofold: it suppresses rotation of a dipole moment vector and induces a dipole moment in trapped atoms, molecules, and nanoparticles.
To verify confining efficiency of the trap we analyse dynamics of a dipole in the proposed electromagnetic field using analytical as well as numerical methods. We demonstrate that it is possible to trap and confine it for considerably large intervals of time of order hours of even days.
The dynamics of the dipole is described by the canonical Hamilton's equations. The motion of centre of mass of the dipole is coupled with its rotational motion. The energy is continuously exchanged between the translational and rotational degrees of freedom. This is why the centre of mass orbit is modified by the dipole moment rotations and vice versa.
We found certain particular solutions for which the whole energy is accumulated in dipole moment rotations. For these solutions equations of motion are exactly the same as for the spherical pendulum.
We analyse stability of a periodic solution corresponding to a uniform precession of the dipole around the symmetry axis of the trap. We study the linear perturbation of this motion. Passing to the appropriate rotating coordinate frame we split the equations of motion into two separate subsystems describing perturbation of the motion of the centre of mass and the perturbation of rotational motion of the dipole. Both subsystems are linear differential equations with constant coefficients. In the linear approximation the region of stability is determined in the space of the trapping parameter κ, precession angle θ and angular velocity of precession ω.
Analytical analysis is complemented by results of numerical simulations of Hamiltonian equations of motion for three polar particles with dipole momentum: the cellulose nanorod, the DAST nanocrystal and the protein lac repressor. We show that even for 'large' values of initial translational momenta the mass centre orbits are bounded to a restricted volume. We also propose other physical polar molecules and nanoparticles that can be trapped by means of the proposed electromagnetic fields, and its ruling parameters for trapping are collected in table 2.
Finally, we remark that stability of two families of periodic solutions found in section 3.2 can be used for finding trapping conditions for other regimes of operation of the proposed trap. However, for these solutions the stability can be investigated only numerically. The same remark concerns an arbitrary quasi-periodic spherical pendulum solution described in section 5.1.
When considering classical trajectory, a sextupole electrostatic field confines a neutral dipole along the z-axis (an axis of symmetry of the proposed trap), but repeals it in any direction perpendicular to this axis. A strong quadrupole magnetic field curves the trajectory of the particle in orthogonal plane and ensures confinement. While for polar molecules with ruling parameters collected in table 2 the conditions for trapping are out of reach for current technologies, progress (especially in the magnitude of magnetic gradients applied, and/or size of the dipole moment) enables to envision the trapping of nanoparticles with huge permanent dipole moment, and maybe (with orders of magnitude gain) of massive complex molecules. This is because the voltage has upper limit U max defined by the magnitude of gradient of the quadruple magnetic field. The limits for various samples are given in the second last column of table 2; the formula is written in the last line of table header. The magnitudes vary from μV to mV. Values for DAST II, lac repressor, and AChE enzyme are achievable for current technologies.
To make the trapping more reliable, we propose the smaller voltage presented in the second column of table 2. For the samples mentioned above the magnitudes are near the hundreds of μV.
Unfortunately, for the other samples the values of voltage are out of reach of current technologies. We perform numerical simulations for them in hope of future progress.
The ultimate aim of trapping of polar molecules is to reach the quantum regime where polar molecules are effectively confined in inhomogeneous electrostatic field [38, § 6.4]. A magnetic field is not necessary. In reference [38, § 6.2] the series of experiments is described where the electric fields are used for cooling of beam of polar diatomic molecules CO, YbF, OH, and polyatomic molecules of benzonitrile (C 7 H 5 N).
In reference [7] the neutral ground state OH molecule is confined in the combination of quadruple magnetic field and quadruple electric field. Their magnitudes are not related to each other. In the present paper we design the trap with the combination of non-homogeneous anisotropic electromagnetic fields. We will look further at the quantum dynamics of diatomic polar molecules in our trap in our forthcoming paper. We will calculate the Stark and Zeeman splitting of the energy caused by the interaction between molecular electronic spin and the trap electric and magnetic fields. Since both the non-homogeneous fields vanish at the geometrical centre of the processing chamber, the idealised situation describes a diatomic molecule in a homogeneous electric field. The problem of a quantum rigid rotor is well elaborated [38, § 5.1.2]. In a real trap the rotor is placed in the superposition of infinitesimally small coordinate dependent electric and magnetic fields. We will apply standard perturbation theory to calculate corrections to energy spectrum. This allows us to evaluate the system stability. A detailed description of future calculations for a specific diatomic molecule goes beyond the scope of the present paper. The relevant parameters are the rotational constant, the spin-rotation interaction constant, and the spin-spin coupling [38].
So far as the quantum description of polar polyatomic molecules and nanoparticles, they will be modelled as rigid bodies. Classical Hamiltonian is defined by the tensor of inertia and tensor of charge distribution. The configuration space of a non-degenerate rigid body R 3 × SO(3, R). The quantum mechanics of rigid body develops since 1931 [39,40].
where positively defined parameters a and b are associated with trap's design (see figure 16). On the positive electrode (A2a) the potential takes the value U a = C 1 a 3 + C 2 , while on the negative one (A2b) it is equal to U b = −C 1 b 3 + C 2 . Only difference of the potentials U 0 = C 1 (a 3 + b 3 ) is physically meaningful. It immediately gives the constant C 1 = U 0 /(a 3 + b 3 ). Putting C 2 = 0 we see that, when a voltage U 0 = U a − U b is applied between the positive and negative electrodes, the arrangement creates an electromagnetic potential in the form (10). Let us invert the formulae representing the form of electrode surfaces given by cubic polynomial (A2a). If R = 0, then Z = a where a is the minimal distance to the geometric centre of trap. If 0 < R 2 1/6 a the Cardano's discriminant Δ = 27 −R 6 /2 + a 6 of the depressed cubic is non-negative. There is only one real root which can be represented using hyperbolic functions If R 2 1/6 a, the Cardano's discriminant Δ 0. The cubic polynomial (A2a) has three real roots. We express them in the following form where φ(R) = 1 3 arccos Zeroth root traces the curve which is continuation of segment (a, A) given by hyperbolic functions (A3) and (A4). Rotating this curve around the z-axis we obtain the surface of positive endcup electrode painted blue in figure 17.
The roots Z 1 and Z 2 define the portions of curve which touch each other at point B with coordinates [2 1/6 a, −a/2 1/3 ] (see figure 16). Rotating this curve around the z-axis we obtain the surface of positive ring electrode sketched in black. Resuming, black contours in figure 17 sketch positively charged electrodes (A2a), endcap and ring, which contain the points A and B, respectively. Their surfaces are generated by rotating the curves (A3) and (A5) around the z-axis and painted blue in figure 17.
Substituting −b 3 for a 3 in equations (A4) and (A6) we parameterize the curves in grey in figure 16. Rotation of these curves around the z-axis gives the surfaces of negative electrodes (A2b), endcup and ring painted orange in figure 17.
Thin lines z = ± 3/2ρ in figure 16 represent the cones that are the asymptotes of the endcap electrodes. Ring electrodes have three asymptotes: the cones and radial plane z = 0. Axially symmetric electrodes are used to produce the sextupole potential (10) that gives the electric field E = −∇Φ 3 which is pictured in figure 2(a).