Derivation of equations of the model of the dynamic behaviour of the three-dimensional atmospheric cloud of electrically charged ice crystals under the influence of electrostatic forces

The purpose of this research is to create a theoretical model describing the dynamics of a cloud of electrically charged ice crystals in a three-dimensional space. The linear motion of the crystals was declined and considerations were limited to rotation of these crystals under the influence of electrostatic forces. The crystals were modeled as thin circular discs to simulate hexagonal plate ice crystals. This has reduced problem to the problem of two degrees of freedom (that allow for rotation around its center of gravity) per crystal. The basic assumption of the model is the crystals in the cloud are synchronized, which means that the differences in rotation of adjacent crystals are small. With this assumption, it was possible to derive continuous differential equations describing the dynamic behavior of such a medium on a macro scale from equilibrium of moments (from electrostatic forces) acting on the single crystal and discrete equations of motion of the cloud. Finally, the system of differential equations is second order with respect to time and space with trigonometric variable coefficients is derived. Examples of solutions of the system of equations are shown and compared with both solutions of similar theoretical problems available in the literature and optical effects in the clouds described by observers. It can be used in the future to better understand the dynamic optical phenomena occurring in the atmosphere.


Analysis of the problem
The problem of modelling the dynamic behavior of electrically charged ice crystals is relatively new issue. At the beginning of the twenty-first century, Foster and Hallett in a number of works, e.g., [1,2] on the basis of laboratory experience, analyze the distribution of ice crystals in an electric field. This research was inspired by the fact that in many earlier works [3] it was found that clouds containing ice crystals, similarly to clouds made of water droplets, may accumulate electric charges. Later, in [4] a hypothesis explaining an atmospheric optical phenomenon commonly called the miracle of the sun as a result of the passing of sunlight through a synchronously vibrating electrically charged cloud of ice crystals was presented. This hypothesis was then extended to other similar phenomena by other authors [5]. In the following years, the theoretical model of vibrations of synchronous clouds of electrically charged ice crystals was extended. The mechanical model presented in that paper [4] was one-dimensional and had many limitations, including small angles of crystal rotation and point distribution of electrical loads. Then in [6] the authors has presented a non-linear two-dimensional model of such a medium. Although at the same time a lot of 3D models describing the performance of ice crystals in clouds in relation to various physical phenomena were created [7,8], to date no threedimensional theoretical model of synchronous vibrations of an electrically charged cloud of ice crystals has been formed, which would be an extension of the model proposed in the paper [4] and which could be the basis for explaining the observed optical effects in the phenomenon called the ''miracle of the sun''.

Purpose and scope of work
The purpose of this work is to create a threedimensional, continuous mathematical model describing the dynamic behavior of a cloud of electrically charged ice crystals. Based on assumptions, described in detail in the next chapter, the number of degrees of freedom of such medium is reduced to two. The work involved the topic of non-linear vibrations occurring in both small and large angles. The final equations of motion of the medium were obtained using the finite difference method and the expansion of certain functions into the Taylor series. Derived nonlinear system of two differential equations of the second order with respect to time and space have been solved for the simplest case of crystal vibrations in one plane. he obtained results were compared with both theoretical solutions of 2D issues available in the literature and observations of the real optical phenomenon in the atmosphere.

Assumptions of theoretical model
Three-dimensional space be given in which a cloud of electrically charged ice crystals covering the X area is considered. The following assumptions were done: 1. Cloud covers an area X & R 3 , is homogenous and its crystals are identical.
2. The crystals have round discs shape with radius r and height h, where h ( R (Fig. 1). 3. The crystals are evenly distributed in all three directions (Fig. 1). 4. Linear motion of the crystals in all directions is negligible (we assume that their centers of gravity do not move towards each other). 5. The rotation of the crystal around an axis perpendicular to the surface of the crystal is negligible. 6. The electric charge is evenly distributed on the upper and lower surfaces of each crystal. 7. There is synchronization between the crystals, wich means that the differences in rotation angles between neighboring crystals are relatively small.
Note that assumption 2 simplifies the problem to modeling plate crystals [6]. For simplicity by using a circle we approximate a naturally occurring form of a regular hexagon. Assumptions 4 and 5 reduce the spatial problem in which each crystal would have 6 degrees of freedom, to the simplified problem with 2 degrees of freedom. Further considerations of the model will show that it is enough to analyse in detail the influence of the nearest 8 crystals (marked in blue in Fig. 1) and all further crystals, if we group them properly, provide the influence of identical character and value decreasing with the distance. Therefore, our considerations may be developed into the whole R 3 space and infinite number of crystals through the transition to the border with the distance between the central crystal and the group of analysed crystals interacting with it. A detailed analysis of the relationship between the value of moments and the number of crystals taken into account will be carried out in chapter 2.7. The last assumption leads to a continuous and homogenous model of such a cloud. The physical basis of such assumption is clear and caused by the tendency of the system to organising itself by minimizing the sum of the potential energy of the crystals. Similar assumptions are the basis of many physical phenomena in the field of liquid crystals and nanomaterials, where, under the influence of the external electromagnetic field, the orientation of the particles is aligned [9, 10].

Vector analysis of 3D model geometry
It was considered the case of two crystals located in the space R 3 . One of them named as central and denoted by the index ''0'', the second is named as interacting and denoted by the index ''j''. Let's assume the global coordinate system in the center of the central crystal ( Fig. 1). We define local cartesian coordinates O 0 X 0 Y 0 Z 0 and O 00 X 00 Y 00 Z 00 related to the central and interacting crystal, respectively. For simplicity, we require that the axes O 0 X 0 and O 00 X 00 respectively be the lines of the steepest descent on the crystals defined in the global coordinate system (Fig. 1). Such a definition causes local coordinate systems not to be permanently assigned to crystals, but they are variables in time.
Lets introduce two-dimensional polar coordinates related to both crystals O 0 r 0 x 0 and O 00 r j x j . The angles x 0 and x j will be measured clockwise from the axes O 0 X 0 and O 00 X 00 respectively. Let us consider an infinitesimal parts of the two plates with dimensions dr 0 Â dx 0 and dr j Â dx j respectively. Let us denote the position vectors of these infinitesimal parts by W 0 and W j . Thus, in the local systems O 0 X 0 Y 0 Z 0 and O 00 X 00 Y 00 Z 00 related to the corresponding crystals these vectors have the coordinates: Let us denote by R Y ðH 0 Þ and R Y ðH j Þ rotation matrix by an angle H 0 and H j about O 0 Y 0 and O 00 Y 00 axis respectively, and by R Z ðN 0 Þ and R Z ðN j Þ rotation matrix by an angle N 0 and N j about O 0 Z 0 and O 00 Z 00 axis respectively: where i j; 0; ð2Þ where i j; 0: Note that with this definition of local coordinate systems the angle H is defined between the global horizontal plane and the plane of the crystal. Analogically, angle N is the angle between the global axis OX and the local axis O 0 X 0 , thus defining the plane of the steepest descent of the crystal from the horizontal plane. Also note that if we transform the coordinates of any vector in a certain coordinate system using a appropriate rotation matrix, this procedure can mathematically be interpreted in two ways: • the coordinate system remains unchanged and the vector rotates around its center, • the vector remains constant (eg relative to the global coordinate system), and the local coordinate system changes-but in the opposite direction. In this paper second interpretation of rotation matrix has been used.
We are going to write the vector W ! j in the local coordinate system 0X 0 Y 0 Z 0 related to the central crystal. To do this, we first multiply the vector W ! j by an angle N about O 00 Y 00 axis. As we have already established, after this transformation, the coordinate system 0 0 X 0 Y 0 Z 0 is virtually rotated by an angle H j . Due to the fact that the local axis O 00 X 00 is denoted as a direction of the steepest descent on the interacting crystal, the O 00 Y 00 axis is always horizontally in the global coordinate system. If the angle H j is defined as the angle of inclination of the crystal relative to the level, then after rotation by angle H j the coordinate system O 00 X 00 Y 00 Z 00 remains horizontally in relation to the global coordinate system. Angle between axis OX and axis O 00 X 00 will be denoted as N j , hence, to write the vector W ! j in a global coordinate system we multiply it by the rotation matrix R x ðN i Þ. Analogous procedure is used to transform this vector from the global system into a local coordinate system related to the central crystal. We apply here the transposed matrices of rotation and we obtain the vector W ! j finally written in the local coordinate system related to the central crystal: Analogous reasoning is performed for the vector W ! d connecting the centers of both crystals and we obtain it written in O 0 X 0 Y 0 Z 0 coordinate system: Wherein the vector W ! d depends on the position coordinates of the interacting crystal j: where 2d, the characteristic size of the grid of crystals (Fig. 1); a, b, d, integers that determine the position of the interacting crystal in the space. Finally, we can write the vector connecting the two infinitely small pieces of crystals as the sum of the two vectors: By substituting Eqs. (1)-(6) to (7) we get: where This vector will be used to calculate moments from the electrostatic forces acting on the central crystal.

Equilibrium equations of a single crystal
in local axis OX 0 Y 0 Z 0 Our discussion will start from the equations of motion for the rotation of a single crystal of the individual scalar components: where j, the axis of rotation of the crystal in O 0 X 0 Y 0 Z 0 the local coordinate system; M j , the j-th component of the total moment M from the sum of electrostatic forces acting on the central crystal; e j , the j-th component of the angular acceleration e about local axis of rotation; crystals moments of inertia about local axis of rotation; R, the radius of the disk (Fig. 1); h, the crystals height (Fig. 1); q ice density.
Note that Eq. (11) represents only the instantaneous values, as with the rotation, the local coordinate system and crystals' orientation change, as a result, the values and directions of all moments and angular accelerations components change.
The moment M is defined as the sum of the moments from the electrostatic forces acting on the censtral crystal from all interacting crystals. Due to the assumption 2 and 6 it will be given by a general vector equation: where N, the number of interacting crystals; x j Þ, the force of electrostatic interaction between two infinitesimal parts of central and inter acting crystal (Fig. 1). This force is given by the Coulomb formula: where dq, the electrostatic load of infinitesimal parts of central and interacting crystals; e 0 , the vacuum permittivity; e c , the relative permittivity. Due to the assumption 5 we will consider only two of the three equations of motion defining a crystal rotation-around the local axes O 0 X 0 and O 0 Y 0 . Using the properties of the vector product and substituting the formulas (11)-(13) we finally obtain: 2.4 Numerical derivation of discrete equations of motion The direct, analytical calculation of quadruple integrals (14) and (15) is extremely difficult. Therefore, in order to obtain the solution of these integrals and derive the discrete equations of rotation of the crystal, we first apply the expansion of the functions M x and M y to the Taylor series. We will consider Eq. (14). The Moment M y described by Eq. (15) will be considered analogously, for this reason, we will first focus on the first Eq. (14) describing the moment M x , then, in an analogous way, we will consider Eq. (15) and M y moment.
Moment M x is a function of 4 independent variables: where j, the number of the considered interacting crystal; N, the number of interacting crystals considered.
Given the assumption 8, we can assume that the difference between the angles are always small enough: Let the space be four-dimensional u ¼ ðH 0 ; H j ; N 0 ; N j Þ and we assume a certain point in this space u 0 ¼ ðH 0 ; H j ; N 0 ; N j Þ We use Taylor expansion of the function M xj around point in variables H j , N j up to first order: where R 2 ðuÞ is negligible. Then we use the finite difference method and we expand the derivatives found in the formula (18).
By substituting (19) and (20) to (18) and treating constant parameters H 0 , N 0 as variables we obtain: where Let's note that assuming some finite but small values of parameters dH j and dN j we can numerically calculate the value of integrals (14) for a fixed pair of angles (H 0 , N 0 ). In this way, we can numerically determine functions M 0 Þfor each interacting crystal j separately. Figure 2 shows graphs of the functions (22)-(23) obtained for the first two crystals interacting with centers at points with absolute coordinates presented in Table 1.
The functions obtained are cyclical with a period of 2p for both variables. It is also worth noting that for the crystals interacting on opposite sides of the central crystal the functions obtained are identical, which means: for every x, y, z and for all 3 functions: Property (25) is used to derive continuous equations of motion. 72 function coefficients (8 interacting crystals, 3 directions x, y and z, and 3 types of coefficients 0, H, N) were found.

Numerical derivation of continous equations of motion
In Chapter 2.4, we obtain equation for moment (21) with functional coefficients, but it is defined  Fig. 2 Examples of graphs of the coefficients from the Eq. (21). a M 0 independently for each interacting crystal j. In this chapter, using the finite difference method and the property (25) we will construct a continuous equation defined for the entire cloud of crystals. In this chapter, using the finite difference method and the property (25) we will construct a continuous equation defined for the entire cloud of crystals. Let us first derive the formulas for the differential derivatives for the three-dimensional FDM grid we have. The assumed distribution of FDM grid points is identical with the distribution of crystals in space (see Fig. 1 and Table 1). We assume function F to be any continuous and differentiable to the second power function in the space R 3 : F x; y; z ð Þ. We use Taylor expansion of the function F x; y; z ð Þaround point X 0 ¼ x 0 ; y 0 ; z 0 ð Þ¼ 0; 0; 0 ð Þin 3 variables x, y, z up to second order: where parameters a, b, d have the values according to the Table 1.
Let us note that if we substitute the data from the Table 1 for formula (26), then some combinations of sums and differences of functions defined for each point give us simple derivatives on the right side of the equation: where in the symbols are defined as it is shown in Table 2.
In this way, we have derived formulas for partial derivatives of the 2nd order for our 3D FDM grid. Let's now define functions that are combinations of functional coefficients M 0 Þ and H j angles analogous to (27)-(28) formulas: In an analogous way, we could define the formulas for variable y and the moment (11) and the angle N. Let us also note that after taking into account properties (25) there are equalities: and similarly we obtain: By summing the Eq. (21) for all the 8 interacting crystals, and using Eqs. (31)-(32) and finite difference formulas (27)-(28) applied to the definition (29)-(30) we finally get: Table 2 Description of used symbols Note that the functional coefficients M H found in formula (33) can be found numerically as defined by (29)-(30).
Note that the functions shown in Fig. 3 have much more symmetry and antismetry axes, than the functions in Fig. 2. These functions can also be well approximated by trigonometric functions. As a base approximating functions assume the function of the It is important to note that after adopting such a basis of approximation we obtain with great accuracy (relative error less then 0.01%) integral coefficients before some expressions and a significant number of them are equal to zero.
The approximation procedure is performed for all coefficients occurring in Eq. (33) and the corresponding equation for M y moment. Thus, there is the possibility of writting the Eq. (33) in the form of analytical formulas. The dependence of coefficient (34) on the distance between crystals was determined on the basis of a series of numerical experiments.
Let us now consider the connection between the individual components of the moment (11) and the time derivatives of the unknown angles H and N. Due to the fact that the axis O 0 Y 0 is oriented always horizontally in the global coordinate system then the rotation around this axis is rotated by an angle H 0 . Therefore, considering the situation of the Eq. (11) we have: In a similar way we can say that since the O 0 X 0 axis is always the axis of steepest descent in the global coordinate system rotation around this axis by a small angle dH changes the direction of the axis O 0 X 0 ie rotation of angle dN equals to value of dH. Therefore, we can write an equation for the moment M x in analogy to the previous: Fig. 3 Examples of functional coefficients of formula (33). a M H Finally, we get the equations in the form: where a; b, the model constants that value depend only on the amount of considered interacting The Eqs. (29), (30) of motion are second order relative to time and space and they have continuous, trigonometric variable coefficients. To get numerical constants in the equation taht do not depend on the amount of crystals used for the modeling the amount of crystals considered and the distance from the central crystal must tend to infinity.

Extension of the model into infinite number of interacting crystals
The coefficients in Eqs. (37), (38) depend on the number of interacting crystals considered. Previous considerations were carried out for 8 crystals of X ¼ AE1; AE1; AE1 ½ coordinates. It should be noted, however, that if we appropriately group the further interacting crystals, their influence will have an identical character, and will differ only by the proportionality coefficient. Table 3 lists the next 234 crystals together with the analysis of their influence on the a coefficient. In order to obtain the coordinates of all n crystals from a given group one should permutate the x, y, z coordinates.
An analogous list can be made for the b coefficient (Table 4).
In both cases we can see that the influence of successive groups of crystals decreases approximately in the fifth power of the distance, which results in a rapid stabilization of the results. Additionally, subsequent groups of crystals have a partially increasing, partially decreasing effect, which results in even faster stabilization of the sum of this sequence. As a result, after considering 234 crystals, we can determine the values of a and b coefficients with the accuracy of 0.5% as follows:

Example of the solution of model equations
Equations (37) and (38) with its coefficients are in general a system of two non-linear differential equations of the 2nd order in relation to time and space. Finding general solutions for such a system of equations is practically impossible, but we can find Columns x, y, z describe coordinates of crystals acting from a given group, n-number of crystals in the group, bcoefficient obtained for a given group, d-distance of crystals from a given group to the central crystal a Including all permutations some specific solutions which, despite their simplicity have a practical meaning, as we will show in the next subsection. To simplify the solution, we assume that the vibrations of all crystals take place in a single plane, i.e., Equations (37) and (38) are then reduced to one equation of the form: Equation (40) is a certain extension of the known sine-Gordon equation in 3D. It is worth comparing them with the equations considered in the works of other authors, e.g.: [11][12][13][14][15]. It should be noted that the equations analysed in the above works differ from Eq. (40) only by the following coefficient cos 4H ð ÞÀ0:2 ½ . In all cases, regardless of the adopted method of analytical or numerical solution and regardless of the number of spatial dimensions taken into account, the authors state the existence of solutions of the sine-Gordon equation in the form of kink or breather solitons. Suggesting the conclusions from [11][12][13][14][15] we will look for the solution of the Eq. (40) in the form of breather and we will limit ourselves to one spatial dimension (3D, 2D and 1D solutions have the same character) and we will adapt the following initial and boundary conditions. As a result, the numerical problem that we solve can be presented as: Hðx; tÞj x¼AE1 ¼ 0; Hðx; tÞj t¼0 Let us assume an example of probable, real values of dimensions of ice crystals occurring in high clouds as [16]: Remembering about (39) and the above values we can adapt to the calculation example: We will perform numerical calculations using the Maple software and the finite difference method, assuming a spatial step Dx ¼ 0:2m and a time step Dt ¼ 0:2s. We should note that the solution is a quasibreather that disappears over a relatively long period of time. In the centre we have practically cyclic x, y, z columns describe coordinates of crystals acting from a given group, n-number of crystals in the group, bcoefficient obtained for a given group, d-distance of crystals from a given group to the central crystal a Including all permutations vibrations with slowly decreasing amplitude, which generate a certain wave propagating and disappearing with increasing distance from the centre. This situation is shown in the following space-time graphs for the time intervals on Fig. 4. The solution obtained depends to a great extent on the adoption of specific physical values related to the size of crystals, their charge and cloud density. Only for some specific combinations of these values can we obtain a stable or quasi-stable solution. If the numerical factor on the left side of Eq.  Too high a coefficient quickly leads to an unstable solution.

Comparison of the obtained solutions with observations and other theoretical results
We can check the obtained sample solution in two ways: first, by comparing it with other theoretical solutions of equations of a similar nature available in literature and second, by comparing it with observations and relations of witnesses of real optical atmospheric phenomena. The literature presents commonly available solutions of the sine-Gordon equation, both numerical and analytical [11][12][13][14][15]. The solutions presented in these works are completely stable, theoretically not fading away in a long period of time. Adding coefficient causes spatial propagation of disturbances. Indeed, if in our model we would analyze the influence of b factor on the shape of the solution, we would notice that at b ! 0 we are really receiving classical literature solutions of the breather type. A comparison of the results obtained with the observations of witnesses of the atmospheric optical phenomenon called the ''miracle of the sun'' leads to much more interesting conclusions. In the films available on the Internet and analysed in the paper [4] we deal with quasi-periodic changes of apparent brightness of the sun. However, the recorded vibration sequences of the phenomenon are relatively short, the phenomenon occurs most frequently in short sequences after 10-60 s, during which the period of vibration ranges from fractions to several seconds. This fact is consistent with the obtained numerical solution, where the vibration period was about 4 s and the vibrations significantly decreased with time. Since when their amplitude drops below the value of sun height above the horizon, according to work [4] the observer will stop seeing the phenomenon, we can say that the numerically obtained time of disappearance of the phenomenon is probable. Amplitude of crystal vibrations, according to the work [4] must be connected with the height of the sun above the horizon. The maximum amplitude of vibrations leading to a quasi stable solution is 0.3-0.4 radian, i.e. not more than approx. 22 . This fact is also consistent with the observations that are made in the vast majority of evenings at sunset when the sun is relatively low above the horizon. Finally, let us also note that from the theoretical point of view, this phenomenon requires very specific ranges of physical input parameters (the size of crystals, their charges and cloud density), for which we have a quasi stable solution. This leads to the obvious conclusion that the occurrence of this phenomenon in nature is relatively rare and requires very exceptional physical-geographical conditions in order to occur in reality. This conclusion is also consistent with observations, because in contrast to other optical phenomena occurring in the atmosphere (rainbows, halo, arches) the phenomenon of ''miracle of the sun'' occurs very rarely and only in certain places in the world.

Summary
As a result of the analysis and modeling of cloud vibrations of electrically charged ice crystals we have built a model of dynamic performance of such a cloud described by a system of partial differential equations of the 2nd order in relation to time and space (37)-(38). As a result of taking into account a large number of interacting crystals, coefficients occurring in this equation were found numerically. Then an exemplary solution for a one-dimensional case was shown and compared with literature solutions and analyses by observations of a real optical phenomenon called the ''miracle of the sun''. Therefore, the presented model can be successfully used for further analysis of this phenomenon, study of the detailed influence of individual cloud parameters on the possibility of its occurrence and its characteristics.