Stochastic Dissipative Euler's equations for a free body
J. A. de la Torre J. Sánchez-Rodríguez Pep Español





 Dept.   Física Fundamental, Universidad Nacional
  de Educación a Distancia, Madrid, Spain








  Intrinsic  thermal fluctuations  within a  real solid  challenge the
  rigid body assumption  that is central to Euler's  equations for the
  motion of a  free body.  Recently, we have  introduced a dissipative
  and stochastic  version of Euler's equations  in a thermodynamically
  consistent way (European Journal of Mechanics - A/Solids 103, 105184
  (2024)).  This framework describes the evolution of both orientation
  and shape  of a  free body, incorporating  internal thermal
  fluctuations and  their concomitant dissipative mechanisms.   In the
  present  work,  we  demonstrate  that, in  the  absence  of  angular
  momentum, the theory predicts that  principal axis unit vectors of a
  body undergo an anisotropic Brownian motion on the unit sphere, with
  the anisotropy arising  from the body's varying  moments of inertia.
  The resulting equilibrium time correlation function of the principal
  eigenvectors decays  exponentially.  This theoretical  prediction is
  confirmed in  molecular dynamics  simulations of small  bodies.  The
  comparison of theory and equilibrium MD simulations allow us to measure
  the orientational diffusion tensor.  We then use this information in
  the  Stochastic   Dissipative  Euler's  Equations,  to   describe  a
  non-equilibrium  situation of  a body  spinning around  the unstable
  intermediate axis.  The agreement  between theory and simulations is
  excellent, offering a validation of the theoretical framework.

  
Stochastic Dissipative Euler's equations for a free body
    J.A. de la Torre, J. Sánchez-Rodríguez, Pep Español
    April 29, 2024
========================================================



§ INTRODUCTION


The  motion of  a rigid  body  in free  space is  governed by  Euler's
equations  which  assume that  the  distance  between the  constituent
particles                                                           is
constant<cit.>.
This idealization overlooks the inherent  elasticity of solids as well
as the thermal  motion that all particles within a  body experience at
finite  temperatures.   While Euler's  equations  have  been used  for
centuries  to  describe  rotational   dynamics,  they  fall  short  in
explaining  certain  dissipative  phenomena observed  in  nature  like
precession relaxation  or tidal locking. Under  precession relaxation,
free spinning bodies  always end up rotating around the  major axis of
inertia,  explaining  why  98%  of asteroids  are  spinning  in  pure
rotation   <cit.>.     Tidal   locking
describes how a gravitationally bounded  body ends up sincronizing its
orbital and  rotational periods, as in  the case of our  Moon.  From a
conceptual point of view, Euler's  equations can be obtained under the
rigid body  idealization from the relationship  between two orthogonal
reference systems in  motion, the observation that  the inertia tensor
diagonalizes  in   the  principal  axis  reference   system,  and  the
conservation   laws,   in   particular  that   of   angular   momentum
<cit.>.  However,  until recently,  the link
between the Hamiltonian motion of  the particles constituting the body
and Euler's  equations for the motion  of the body was  lacking.  This
missing  link  has been  addressed  in  Ref.  <cit.>.   To
“derive” Euler's  equations from  Hamilton's equations one  needs to
resort to  the powerful  machinery of  the Theory  of Coarse-Graining,
also known  as Non-Equilibrium Statistical Mechanics  and, some times,
as      the      Mori-Zwanzig      projector      operator      method
<cit.>.
In           this           theory,          following           Gibbs
<cit.>  one   introduces  the
notions  of   microstates  and   macrostates.   The   macrostates,  or
coarse-grained  (CG) variables,  define  the  coarse-grained level  of
description of the system.  The theory makes the fundamental modelling
assumption that the selected CG  variables evolve in two distinct time
scales,  one  characterized by  small  and  fast contribution  due  to
collisions/vibrations which  is modelled  as white noise,  and another
slow  contribution  due  to  the cummulative  effect  of  these  rapid
contributions.  In this way, the theory describes the evolution of the
CG      variables     as      a      diffusive     Markov      process
<cit.>.  One possible  level of description for
the  motion  of a  quasi-rigid  body  would  consider  the body  as  a
continuum,  and  use  the  displacement  and  velocity  fields  as  CG
variables, to produce a viscoelastic field description (including free
boundary conditions) of the free  body.  In the present work, however,
we consider  a coarser  level of description,  closer to  the original
Euler's description,  which takes the  gyration tensor as  the primary
variable.   The gyration  tensor  is closely  related  to the  inertia
tensor and its eigenvectors (defining the principal axis) describe the
orientation of the body, while the eigenvalues (referred to as central
moments) describe the overall shape of the body.


In Ref.  <cit.>, we have constructed from first principles
the stochastic differential equations that govern the evolution of the
orientation  and  shape  of  a  quasi-rigid  body.   These  equations,
referred  to  as  Stochastic  Dissipative  Euler's  Equations  (SDEE),
generalize Euler's equations in  order to include thermal fluctuations
and  its  associated  dissipation.    The  reversible  part  of  these
equations is  given by  the usual Euler's  equations, while  the newly
derived  dissipative part  of  the dynamics  contains two  dissipative
mechanisms: an  orientational diffusion of  the principal axis,  and a
dilational friction  damping the oscillations in  the central moments.
In addition,  the SDE describe  the random motion of  the eigenvectors
and eigenvalues due to thermal fluctuations.

For     macroscopic     bodies,    thermal     fluctuations     become
negligible. However,  dissipation cannot be neglected  in general.  In
Ref. <cit.>,  we have  studied numerically  the ordinary
differential   equations   that   result   from   neglecting   thermal
fluctuations  in the  SDEE.  By  switching  on and  off the  different
dissipative mechanisms  involved, we  gave support  to the  claim that
precession relaxation  is due  to orientational diffusion  rather than
dilational  friction.  In  fact,  Euler's equations  predict that  the
motion of a  body spinning around the intermediate axis  of inertia is
unstable
<cit.>,
a result known  as the tennis racket theorem or  the intermediate axis
theorem. As  a consequence,  a free  body experiences  the Dzhanibekov
effect, a striking  phenomena in which the  intermediate axis performs
180^0 degrees flips in a  periodic way. When dissipation is included
in the description, the Dzhanibekov effect dissapears as a consequence
of precession  relaxation.  In  the process,  “organized” rotational
kinetic energy  is transformed into “disorganized”  internal thermal
energy while  increasing the entropy  of the system,  and consequently
heating the  body.  The resulting  minimum kinetic energy  occurs when
the body  is spinning  around the  axis of  larger moment  of inertia.

In  the  present  paper,  we  turn our  attention  to  the  Stochastic
Dissipative  Euler's  equations   including  thermal  fluctuations  as
formulated in  <cit.>.  The objective  of the paper  is to
compare  the  predictions  of  this  theory  with  Molecular  Dynamics
simulations of a  free body composed of bonded  interacting atoms, see
Fig. <ref>. This  constitutes a stringent test  of the theory
and  consitutes a  necessary  validation step.   The SDEE  encompasses
several  parameters  that  must  be  determined  to  enable  numerical
solutions  and generate  predictions.   The list  of these  parameters
include the  average central moments giving  the size and geometry of the  body, an
elastic  matrix defined  in terms  of  the covariance  of the  central
moments, and  two sets  of dissipative coefficients:  an orientational
diffusion matrix and a dilational friction matrix. The values of these
parameters are  obtained from  equilibrium MD  simulations of
the body at zero angular momentum.  The test phase of the procedure is
to   predict   non-equilibrium   results   about   precession
relaxation using these previously determined set of parameters.












In this paper,  we further elucidate a  theoretical insight concerning
the SDEE derived in <cit.>.   We demonstrate that when the
angular momentum  of the body  vanishes, the principal axis  follow an
anisotropic Brownian motion  on the unit sphere.  In this  way, a body
of nanoscopic dimensions  with zero angular momentum  will explore all
possible  orientations <cit.>.   Loosely speaking  the
body “spins” without  angular momentum. The key  observation is that
the  movement of  the principal  axes aligns  with the  mathematically
precise      concept      of       spherical      Brownian      motion
<cit.>.  The  anisotropic behavior is  due to
the distinct moments of inertia, which introduce a directional bias in
the Brownian  motion.  This  bias reflects  the principle  that larger
moments  of inertia  result  in slower  rotational  motion around  the
corresponding  axis.  We  stress  that  the  physical  origin  of  the
rotational Brownian  motion of  a free body  is the  intrinsic thermal
fluctuations of the particles that  constitute it.  This is physically
different from the usual Brownian  Rotor, where a particle immersed in
a  fluid experiences  Brownian motion  due to  the bombardment  of the
surrounding                                                  molecules
<cit.>.
    


As the motion of the principal axis is a spherical Brownian motion, we
predict that  the time-correlation  function of the  principal vectors
decays in a matrix exponential  form.  We compare this prediction with
the  MD   simulation  results,   finding  excellent   agreement.   The
exponential matrix contains  the orientational diffusion coefficients,
and  the  fitting  against  the  MD simulation  results  allow  us  to
determine these parameters.


The  paper is  organized  as follows.   In  Sec.  <ref>  we
summarize  the main  result of  Ref. <cit.>  which is  the
SDEE. In  Sec <ref> we  consider the SDEE when  the angular
momentum  of the  system is  zero,  and show  that the  motion of  the
principal  axis correspond  to  an anisotropic  Brownian motion.   The
details of the mathematical description of anisotropic Brownian motion
are given  in Appendix <ref>.  In  Sec.  <ref>
we run equilibrium MD simulations of a small body, and measure all the
required  parameters  in  the  SDEE.  In  Sec.   <ref>,  we
compare  the numerical  predictions of  the mesoscopic  SDEE with  the
results of  non-equilibrium MD  simulations in which  the body  is set
into  motion through  an “angular  kick” that  suddenly sets  a body
initially  at rest  into  rotation along  the  intermediate axis.   We
compare  the  evolution of  the  rotation  kinetic energy,  reflecting
precession relaxation, in the  mesoscopic and microscopic simulations,
obtaining excellent agreement.  Finally, in Sec. <ref>
we present our conclusions.


§ THE MOTION OF A QUASI-RIGID FREE BODY


In  this   section,  we  summarize   the  theory  presented   in  Ref.
<cit.> for  the motion  of a  quasi-rigid free  body. This
allow us to set up the notation used in the paper. The body, which is
composed  of particles  bonded with  a potential,  is described  at a
coarse-grained level  with the center  of mass position  R̂
and  gyration tensor  Ĝ,  closely related  to the  inertia
tensor Î.  These CG variables capture how the particles of
the body distribute in space, the first giving the “location” of the
body and the  latter giving a sense of its  “shape and orientation”.
These  variables  are the  ones  used  to  describe  a rigid  body  in
Classical Mechanics under  the rigid constraint assumption,  and it is
natural  to us  these phase  functions as  CG variables.  The gyration
tensor is defined as

    Ĝ   ≡1/4∑_im_i( r_i-R̂)( r_i-R̂)^T

where the superscript T denotes the transpose matrix.
The prefactor  1/4 in  the definition  (<ref>) allowed  us in
<cit.>    to    interpret   directly    the    eigenvalues
M̂_α  of  Ĝ as  “dilational  masses”,  or
inertia to dilations.  The inertia tensor Î is defined as

    Î   ≡∑^N_im_i
    [ r_i-R̂]_×^T
    [ r_i-R̂]_×

where [  u]_× is  the cross-product  matrix formed  from a
vector  u.   The  action of  the cross  product  matrix on  an
arbitrary vector  v  gives the cross  product of  both vectors
[ u]_× v= u× v.

The  tensors
Ĝ and Î are both symmetric, positive definite,
and  they  commute  with  each  other.   As  a  result,  they  can  be
simultaneously  diagonalized   in  the  same  reference   system.  The
principal  axis system,  denoted as  S^0,  is defined  as the
reference system with its origin at  the center of mass, in which both
tensors   can  be   diagonalized.   Let   ê_α   with
α=1,2,3  be  the  basis  vectors  of  the  inertial  laboratory
reference system S and ê^0_α be the basis of
the  the non-inertial  principal axis  reference system  S^0.
The components of the rotation matrix  of S^0 with respect to
S are defined as

    R̂_αβ   =ê^T_βê^0_α

This shows that ê^0_α  are the  rows of the  rotation matrix
R.
In S^0 the gyration and inertia tensor take the form
 
    R̂ĜR̂^T   = 𝔾̂
    R̂ÎR̂^T   = 𝕀̂

where 𝔾̂ is  a diagonal matrix whose  elements are the
central  moments  M̂_1,M̂_2,M̂_3,  which  we  write
compactly           as           a           list           M̂=(M̂_1,M̂_2,M̂_3). Also,  𝕀̂ is
a   diagonal  matrix   whose  elements   are  the   principal  moments
Î_1,Î_2,Î_3.  In this paper, diagonal matrices are
represented with  voided fonts,  as in  𝔸.  As  the inertia
tensor (<ref>) can be expressed in terms of the gyration
tensor (<ref>) in a linear  way, the eigenvalues of the inertia
tensor are given in terms of the eigenvalues of the gyration tensor as

    Î_α   =  4(M̂_1+M̂_2+M̂_3-M̂_α)

Finally,  observe  that   according  to  (<ref>),  the
vectors ê^0_α are  the unit eigenvectors of the
gyration and inertia tensors.

The  rotation matrix  can be  expressed  in terms  of the  exponential
matrix

    R̂=e^-[Λ̂]_×

where  Λ̂ are  the  attitude parameters,  or
orientation    for    short.      The    angular    velocity    vector
ω is defined as <cit.>

    d/dt R   ≡- R[ω]_×

The time derivative of the orientation and the angular velocity are related through <cit.>

    dΛ/dt   = Bω

where B is the Attitude Kinematic Operator <cit.> given
by

    B   =1+p [ n]_×+q[ n]_×[ n]_×

where        n=Λ/Λ        with
Λ=|Λ|  and  p,q   are  the  following
functions of the modulus of the orientation Λ

    p   =-Λ/2,
          q=1-Λ/2sinΛ/(1-cosΛ)


From (<ref>) and (<ref>), the gyration tensor can be written as

    Ĝ   =e^[Λ̂]_×𝔾̂ e^-[Λ̂]_×

which can be  understood as a change of variables
from the six independent components  of the symmetric gyration tensor,
and         the         six         degrees         of         freedom
Λ̂,M̂.   We    have   chosen   in
Ref.  <cit.> the  latter as  the primary  CG variables  to
describe the motion of a quasi-rigid free body.

The  CG  variables  are  phase functions,  denoted  with  circumflexed
symbols. As  the microscopic  state of the  body evolves  according to
Hamilton's        equations,         the        phase        functions
Λ̂,M̂ also  evolve  in time.   The
theory of  CG offers a modelling  of the evolution of  the CG
variables  in terms  of  a  diffusive Markov  process.    The corresponding Ito
SDE   for  the   orientation   Λ   is  given   by
<cit.>

    dΛ   =  B(Ω
        -Γ^_Λ B^-T
        (Ω× S))dt
        + k_BT F^ th dt+dΛ̃

where the  spin velocity Ω is obtained  from the
conserved  angular  momentum S  of  the body  as

    Ω   = I^-1 S

The spin velocity  Ω is a  dynamic quantity which
is   different  from   the   angular  velocity   ω
introduced in (<ref>) which is  a purely kinematic quantity. For
a rigid  body both quantities  coincide but for  a real body  they are
different    <cit.>.    The    dissipative   matrix    in
(<ref>) is given by

    Γ^_Λ   =     B
        e^[Λ]_× D_0
         e^-[Λ]_× B^T

where D_0 is  the orientational  diffusion matrix

The component α of the thermal drift vector F^ th
is given  in terms of  the Kinematic  Operator as (see  Eq. (135) of
<cit.>)

    F^ th_α   =    (∂ B_α'α/∂Λ_β)
         D^α'β'_0 B_β'β

The noise term in (<ref>) has the form

    dΛ̃   =(2k_BT)^1/2 B
        e^[Λ]_× D_0^1/2 dW̃

where dW̃  is a vector of independent increments  of the Wiener
process satisfying the mnemotechnical Ito rule

    dW̃dW̃^T=1dt

and the square root matrix D_0^1/2 satisfies

    D_0^1/2( D_0^1/2)^T   = D_0

Now,  all the  symbols appearing  in  the SDE  (<ref>) for  the
orientation are defined.


The SDEs for the central moments obtained in <cit.> are

    d M   =Πdt
                  
    
      dΠ   = Kdt-Γνdt
                         +   dΠ̃

Here Π is the  dilational momentum, defined by the first equation
as the time derivative of  the central moments. The α component
of the dilational force is defined as

    K_α   = M_α( 1/2ν_α^2
        +2(    Ω_p^TΩ_p -Ω_pα^2 ).
        
          . +k_BT/2 M_α
    - [Σ^-1]_αβ( M_β- M_β^ rest) )

Here,  repeated   indices  are   summed  over   except  if   they  are
underlined. The need of  breaking this Einstein's summation convention
arises from  the fact that the  central moments do not  transform as a
vector.           The          dilational         velocity          is
ν_α=Π_α/  M_α and the
spin     velocity    in     the     principal     axis    frame     is
Ω_p=e^-[Λ]_×Ω.
The  elasticity matrix  Σ is  proportional to  the
covariance of central moments fluctuations, given by

    Σ   =1/k_BT  (M̂- M^ rest)(M̂- M^ rest)^T^ E

where ⋯^  E is  an average with  the rest
microcanonical ensemble  <cit.>.  M^  rest is
the  average  of the  central  moments  with the  rest  microcanonical
ensemble.   We refer  to the  contribution to  K
quadratic in ν as the convective term, whose physical
meaning  has been  discussed in  <cit.>. The  contribution
quadratic in  Ω is referred to  as the centrifugal
term,   and   the   last   term  involving   the   elasticity   matrix
Σ  as   the  elastic   term.   Observe   that  the
dilational momentum  equation in (<ref>) involves  a dissipative
force   -Γ·ν  that   is  to   be
interpreted as a dilational friction, where Γ is a
dilational friction matrix.


The   temperature   of   the   body   T   appearing   in   the   SDE
(<ref>) and (<ref>) depends in general on the thermal
energy E of the body. For the solid bodies considered in this work,
typically

    T=   E/C

where the  heat capacity is  given by the Dulong-Petit  law C=3Nk_B,
where N  is the number  of atoms  of the body.   The thermal  energy is
defined as

    E   =E-K^ rot-K^ dil

where E  is the total conserved  energy of the body,  the rotational
kinetic energy is

    K^ rot   =1/2 S I^-1 S

and the dilational kinetic energy is

    K^ dil   =1/2Π𝔾^-1Π


The   SDE  (<ref>)   for  the   orientation  and   the  SDE
(<ref>) for central moments are coupled through different terms.
For   example,  the   spin  velocity   Ω  in   the
orientation  equation (<ref>)  contains the  inertia tensor
I that  depends on  the  central moments.   Observe that,  in
principle,       the        orientational       diffusion       tensor
D_0= D_0( M, E) and it
depends on  the instantaneous value  of the central moments  M
and the  thermal energy  E  <cit.>. On  the other
hand the equation of the  central moments (<ref>) depends on the
orientation  through  the  centrifugal  term  quadratic  in  the  spin
velocity.





§ THE SDE AT ZERO ANGULAR MOMENTUM


In this section,  we show that when the angular  momentum of the body
vanishes,  the  theory  in  Ref <cit.>  outlined  in  Sec.
<ref> characterizes the evolution  of the principal axis as
an anisotropic Brownian motion on the unit sphere.  We offer a detailed description of
anisotropic   Brownian  motion   on   the  sphere   in  the   Appendix
<ref>.



 §.§ The evolution of orientation

When the  angular momentum  vanishes, S=0, (<ref>)
with (<ref>), reduces to the following Ito SDE

    dΛ   =k_BT F^ thdt +      B^T C dW̃

where 

    C   ≡  (2k_BT)^1/2 D_0^1/2

We wish  to obtain the  SDE governing  the evolution of  the principal
axis or, equivalently, of the rotation matrix.  This requires deriving
an                             equation                            for
R=e^-[Λ]_×  through
Ito Calculus.  Calculations vastly  simplify by using the Stratonovich
SDE corresponding  to (<ref>) and use  of ordinary calculus.
In   Appendix  <ref>   we   show  that   the  Stratonovich   SDE
corresponding to the Ito SDE (<ref>) is

    dΛ   =      B^T C∘ dW̃

From  (<ref>) and  by  using the  chain  rule in  ordinary
calculus, we  show in Appendix  <ref> that the SDE  governing the
rotation matrix is

    d R   =     -[ C∘ d W_t]_× R

Therefore, the column vectors  of the rotation matrix satisfy
the Stratonovich SDE

    d c_α   = -[ C∘ d W_t]_× c_α

or in terms of the vector product

    d c_α    = c_α×( C∘ d W_t)

This  equation  is  identical   to  (<ref>)  in  the  Appendix
<ref> that describes the  anisotropic Brownian motion of
a particle moving on the surface of a unit sphere.  Therefore, the SDE
(<ref>) for the orientation predicts that the unit principal
vectors describe an  anisotropic Brownian motion on the  sphere. It is
straightforward  to show  that  the SDEs  (<ref>) conserve  the
scalar  products  c_α  c_β  (see  Appendix
<ref>).   Mantaining  these  conservation  laws  with  a
numerical integrator  with finite  time step requires  special methods
<cit.>. Instead  of using  (<ref>), in  the present
paper   we   update    the   orientation   Λ   in
(<ref>)  with  a  predictor-corrector squeme.  This  ensures
automatically that the eigenvalues remain unitary at all times.
  


















 §.§ The evolution of the shape

Let  us move  now to  the form  of the  SDE (<ref>)  for central
moments  when  the  angular  momentum vanishes.   In  this  case,  the
dilational  force  has  no  centrifugal  component  and  (<ref>)
reduces to

    d M_α   =Π_α dt
        
    
      dΠ_α   = M_α( 1/2ν_α^2
        +k_BT/2 M_α
        - [Σ^-1]_αβ( M_β- M_β^ rest) )dt
        
       -Γ_αβν_β dt
        +   dΠ̃_α

Observe that when the body has zero angular momentum, the evolution of
the  central   moments  is  uncoupled   from  the  evolution   of  the
orientation.

A good approximation to the  set of equations (<ref>) is given
by the following linearized set of equations

    d M_α   =Π_α dt
        
    
      dΠ_α   ≃ - M^ eq_α[Σ^-1]_αβ( M_β- M_β^ rest)dt
        
       -Γ'_αβΠ_β dt
        +   dΠ̃_α

where      we      have      approximated     in      some      places
M≃  M^ eq,  have neglected the  convective term
quadratic in  the dilational velocity ν  and the small term
k_BT/     M_α.      Finally,     we     have     redefined
Γ'_αβ
=Γ_αβ/            M^
  eq_β. With these approximations, the evolution of
the central  moments is an  Ornstein-Ulenbeck (OU) process,  for which
the equilibrium time-correlation of  central moments can be explicitly
computed.  The  equilibrium time correlation functions  oscillate with
frequencies determined by  the elasticity matrix Σ
which  decay in  a time  scale determined  by the  dilational friction
matrix Γ'.



§ EQUILIBRIUM MD SIMULATIONS AT S=0


In  this  section, we  consider  equilibrium MD simulations of a small body directed  to
measure    the   different    parameters    that    enter   the    SDE
(<ref>), (<ref>). These parameters  will then be used
in non-equilibrium simulations in order to validate the theory.

A parallepiped  crystal made of 6×5×3  atoms that interact
with  a  Lennard-Jones potential  (LJ),  non-linear  bonds, and  angle
contributions is  simulated with LAMMPS  (Large-scale Atomic/Molecular
Massively  Parallel  Simulator)<cit.>.    We  take  ϵ,
σ and  m (the mass of  one atom) as fundamental  LJ units, so
that τ=√(mσ^2/ϵ) is our unit of time.  Details of
the numerical simulations are given in the Supplemental Material.  The
initial microstate is selected by first carrying out an NVT simulation
to equilibrate  the crystal  to a prescribed  temperature T  using a
Nosé-Hoover thermostat.  In order to average different simulations, we
uniformly choose  30 equilibrated microstates  and we subject  each of
them to  a transformation leading to  a fixed value E  of the energy
and  to a  zero  value of  the angular  momentum  vector.  Then,  each
microscopic configuration evolves in an  NVE ensemble until the system
is fully equilibrated  at the prescribed temperature.  Once  we have a
typical equilibrium state at energy  E and zero angular momentum, we
consider two  types of  simulations, equilibrium  and non-equilibrium.
In  the equilibrium  simulations we  continue in  the NVE  ensemble to
enter a  production phase in  which we measure  different observables.
In  non-equilibrium  simulations,  we  produce an  angular  kick  that
transforms  the   velocities  in   the  initial   typical  equilibrium
microstate  in such  a  way that  the  system begins  to  rotate at  a
particular initial angular velocity  around the intermediate principal
axis.


 We set
the total  energy for  the NVE  simulations to  E=2342ϵ. This
value  is a  typical  one  observed in  the  NVT equilibration  phase.
Through the equipartion theorem <cit.>

    ∑_im_i/2 v_i^2^ E_ rest    =3(N-2)k_BT/2

this corresponds to a temperature of T=8.7 ϵ/k_B. We have checked
that, to  a very  good approximation  the temperature  scales linearly
with the  total energy, E=C T  with the heat capacity  following the
Dulong and Petit law C=3Nk_B.

The gyration tensor Ĝ(t) of the body as a function of time is
computed  from (<ref>).   The eigenvalues  provide the  central
moments   M̂_α(t)    and   the   unit   eigenvectors
ê^0_α(t) give the direction of the principal axis.

The  measured equilibrium  averages  of the  central  moments at  rest
M^  rest and the elasticity  matrix Σ
of the central moments (<ref>) are

    M^  rest   =(91.2,62.5,21.0)  m  σ^2
      
    Σ   =(
                           [ 2.635 0.003 0.001; 0.003 1.273 0.001; 0.001 0.001 0.162 ]) τ^4 ϵ

We have checked that the  elasticity matrix is practically independent
on the temperature in the range T=4-8 in LJ units.  Observe that the
elasticity  matrix  is very  approximately  diagonal.  This entails  a
simplification   in    the   dynamics    of   the    central   moments
(<ref>).  By  assuming  that the  dilational  friction  matrix
Γ' is  also diagonal,  the evolution  equations of
the central  moments (<ref>)  take the form  (repeated indices
are not summed over here)

    d M_α   =Π_α dt
        
    
      dΠ_α   ≃ -ω_α^2( M_α- M_α^ rest)dt
    -Γ_α'Π_αdt
        +   dΠ̃_α

where
ω_α^2=                                          M^
  eq_α[Σ^-1]_αα    are    the
frequencies of oscillation. Using  the values (<ref>), the
theoretical frequencies are

    ω_α→(5.873, 6.998, 11.363) τ^-1

    The     equilibrium
time-correlation   for   central   moments,  as   predicted   by   the
Ornstein-Ulhenbeck process is

    M_α(t)  M_α^ eq   =   M_α M_α^ eq
        e^-Γ_α' tcos(ω_α t)

In Fig. <ref> we show  the equilibrium time correlation functions
for  the  central  moments.   To  these curves,  we  have  fitted  the
expressions     (<ref>)     with    the     fitting     parameters

    ω_α   →(5.870, 7.008, 11.357)τ^-1
    Γ'
         =(
        [ 0.0295      0      0;      0 0.0393      0;      0      0 0.0930 ]) τ^-1

The  agreement  of  the  fitted values  for  ω_α  and  the
predicted values (<ref>) is very good, and shows that
the modelling  of the  central moments  with a  simple O-U  process is
quite accurate and allows us to measure the dilational friction matrix
Γ'.

Let us move now to the dynamics  of the unit eigenvectors. As shown in
Fig <ref>, the unit eigenvectors trace a random trajectory as
a  consequence   of  thermal   fluctuations.   We  claim   that  these
MD trajectories  can be modelled  as realizations  of a
Brownian motion on  the unit sphere. To confirm that  this is the case
and  that,  consequently,  the  theory presented  describes  well  the
observed microscopic dynamics of  the orientation  of  the body,  we compute  the
equilibrium time correlation matrix of the eigenvectors

    E(t)   ≡ê^0(t)ê^0T^  eq

As  shown in  the Appendix  <ref>, a  crucial analytical
prediction  of the  SDE (<ref>)  is that  the equilibrium  time
correlation matrix of the eigenvectors is

    E(t)
         =1/3e^- At

where  e^-  At is  the  exponential  matrix, and  the  matrix
A is given by

    A   ≡1/2( Tr[ D]1 - D)

Here, the matrix D is related to the noise amplitude through

    D   ≡ C C^T(<ref>)=
               2 k_BT D_0

In terms of D_0 the matrix A takes the form

    A   
    =k_BT
      ( Tr[ D_0]1 - D_0)

while  D_0 in terms of A is

    D_0
         =1/k_BT[1/2 Tr[ A]1- A]

as can be easily deduced by taking the trace of (<ref>).

The  measured   off-diagonal  elements   of  the   correlation  matrix
E(t)  are vanishingly  small.   A  diagonal time  correlation
matrix  E(t)=𝔼(t)  indicates   that  the  matrix  in
(<ref>) is also diagonal  A=𝔸. The autocorrelation
function that resides in the diagonal  of E(t) is predicted to
decay as a simple exponential function

    𝔼_1(t)   =1/3e^-𝔸_1t
    𝔼_2(t)   =1/3e^-𝔸_2t
    𝔼_3(t)   =1/3e^-𝔸_3t

We  plot in  Fig  <ref> the  autocorrelation  function of  the
principal unit  vectors which  decay exponentially, in  full agreement
with  the  theoretical  prediction.    From  the  fitting  values  for
𝔸_α  we  extract  the orientational  diffusion  matrix
D_0,  as follows. Because A=𝔸 is
diagonal (<ref>) implies

    D_0
         =1/k_BT[1/2 Tr[𝔸]1-𝔸]

which is the diagonal matrix

    D_0
         =(
        [ 𝔸_2+𝔸_3-𝔸_1/2k_BT                 0                 0;                 0 𝔸_1+𝔸_3-𝔸_2/2k_BT                 0;                 0                 0 𝔸_1+𝔸_2-𝔸_3/2k_BT ])

From the  fitted values  of 𝔸_α in  Fig. <ref>
the matrix takes the value

    D_0
         =(
        [ 2.2338      0      0;      0 0.8587      0;      0      0 1.9129 ]) × 10^-6 (τϵ)^-1

In  summary, the  exponential  decay of  the  autocorrelation of  unit
eigenvalues allow us to measure the orientational diffusion matrix.


§ NON-EQUILIBRIUM SIMULATIONS WITH S≠ 0


In  this  section,  we   will  compare  results  from  non-equilibrium
microscopic MD simulations with non-equilibrium mesoscopic simulations
of  the SDEs  (<ref>),(<ref>) for  a spinning  body.
The  parameters  to  be  used  in  the  SDEs  have  been  obtained  in
(<ref>),(<ref>),(<ref>).  Because the
size  of the  system  is  small (the  crystal  has  90 atoms)  thermal
fluctuations are rather large and the  signal to noise ratio is small.
Therefore, in  order to compare  the noisy signals in  both simulation
methods, and  to be able  to validate the  theory, it is  necessary to
perform some  sort of averaging.  On  one hand, the average  in the MD
simulations will be  over initial microstates compatible  with a given
(non-equilibrium) macrostate, i.e., with the same energy, angular momentum,
and central  moments.  The first  two observables (energy  and angular
momentum) are conserved  by construction of the  algorithm.  To obtain
microscopic configurations with the same value for the central moments
M, we carry out a long simulation in the NVE ensemble, compute
the  central moments  along the  entire simulation,  and select  those
microstates  which have  similar values  of M.  This  is done
through  a Kmeans  algorithm  that cluster  the  central moments  into
groups of equal variance, minimizing the within-cluster sum-of-squares
distance to the centroid of each group.  <cit.> Following
this  procedure,  we select  a  cluster  of 50  different  microscopic
configurations  with equal  energy and  angular momentum,  and central
moments differing  by less  than 1%.  Note that we  do not  group the
microstates   with   similar   values  of   the   dilational   momenta
Π, which fluctuates  around zero and seems  to have no
big effect on the results.

On the other hand, the average in  the SDE will be done over different
realizations  of  the  stochastic  process starting  from  an  initial
non-equilibrium macrostate  identical to  the one  ocurring in  the MD
simulations.



 §.§ MD simulations





The  non-equilibrium  MD  simulations  are performed  as  follows.   A
microstate  z  is selected  from  the  equilibrium simulations  with
prescribed values  of E,  S, M.  Then,  a rotation  of the
crystal is performed to reorient it  with the intermediate axis in the
y axis,  and an “angular  kick” is  imparted, as described  in the
Supplemental Material.   This angular  kick transforms  the microstate
z       to      z'       with       an      angular       velocity
Ω=(0,Ω,0) and the  same original total energy
E.  We choose a value of  Ω=0.75 τ^-1. If we run now the
MD simulation starting from the  new microstate z', the crystal will
rotate  around   the  intermediate  axis  with   an  angular  momentum
S   =      I·Ω≡
  S^eq.  After the applied  angular kick, the centrifugal
force  produces  an  inital  expansion  of  the  crystal  that  starts
oscillating, producing  an oscillatory  motion of the  central moments
which  is   superimposed  to   their  thermal   motion  as   shown  in
Fig. <ref>.  Eventually these  oscillations damp out and the
dilational  kinetic energy  K^ dil(t)  defined in  (<ref>)
goes   down  to   its  thermal   noise   level,  as   shown  in   Fig.
<ref>.    After  the   initial   angular   kick  around   the
intermediate  axis, the  system rotates  and displays  the Dzhanibekov
effect in  very much the same  way as in the  fully deterministic case
<cit.>. This  effect is  appreciable from  the flip-flop
evolution   of    the   intermediate    axis   as   shown    in   Fig.
<ref>.  At  long times,  the crystal ends  up spinning
around  the  major  axis,  which aligns  with  the  conserved  angular
momentum vector.  To keep track of the evolution of the system towards
the equilibrium  final state, we  choose as observable  the rotational
kinetic energy K̂^ rot  defined in (<ref>).  A single
realization  of   the  MD   simulation  is  shown   in  gray   in  Fig
<ref>, where it is apparent  that the noise level does
not  allow for  a clear  distinction  between signal  and noise,  thus
justifying the averaging  over 50 different simulations  with the same
E, S and similar central moments.




 §.§ SDE simulations

We carried out 50 simulations of the SDEs with the same initial values
of  the orientation,  central moments  and dilational  momenta of  the
corresponding     non-equilibrium    MD     simulation,    using     a
predictor-corrector     explicit      scheme<cit.>.      In
Fig. <ref> we  compare the  average of  the rotational
kinetic energy K^rot  obtained from the SDE  (red) and MD
(green) simulations.  The rotational kinetic energy starts and ends at
the same values in both SDE and  MD, as a consequence of the identical
initial  macrostate selected.   The  stochastic precession  relaxation
produces a  spinning about  the major axis.   The initial  damping and
long-term decay  time of  the MD  signal is well  captured by  our SDE
model. The  very good agreement of  the MD and SDE  results provides a
further validation of the proposed theory.



§ CONCLUSIONS



In this paper we have compared the predictions  of the
Stochastic           Dissipative           Euler's           Equations
(<ref>),(<ref>)   with  the   results  of   Molecular
Dynamics simulation.   Fitting parameters at equilibrium  allows us to
predict  the  non-equilibrium  precession   relaxation  of  the  body,
reflected  in   the  evolution  of  the   average  rotational  kinetic
energy. In addition,  we confirm through the exponential  decay of the
autocorrelation function  of the principal unit  eigenvectors that the
motion of  these vectors  can be modelled  as an  anisotropic Brownian
motion, where  the stochasticity is  intrinsic, rather than due  to an
external bath (gas or liquid).   The agreement between the theoretical
model and the  simulations is excellent, thus  confirming the validity
of the former.

We have observed that orientational  diffusion strongly depends on the
value of the rest central moments.   Changing the size and geometry of
the  body changes  strongly  the value  of  the coefficients.   Larger
bodies display  very small orientation fluctuations.  In addition, the
correlation time  of the principal  axis also depends strongly  on the
body's  dimensions,  increasing  with  size.   As  a  consequence,  MD
simulations are readily  unfeasible to access the long  time scales in
which a  non-spinning body  will change  its overall  orientation.  In
order to estimate the orientational diffusion coefficients from MD, we
can only  consider small bodies.  We are currently attempting  to find
scaling relationship between geometry  and the orientational diffusion
coefficients. This requires extremely  long simulation times, and will
be presented  elsewhere.  On  the other hand,  using too  small bodies
(like molecules) may break one of the implicit assumptions made in the
Theory of Coarse-Graining used to derive the SDEE, which is the mixing
assumption  of  the  Hamiltonian  flow ensuring  the  existence  of  a
well-defined  equilibrium  state.    Not  all  interaction  potentials
between the particles of the  body ensure this property.  For example,
using almost  linear harmonic  springs does not  give, in  general, an
ergodic  system, as  the  famous  Fermi-Pasta-Ulam-Tsingou problem  of
thermalization                                                  showed
<cit.>.  We  have dealt with
this problem  by ensuring sufficiently non-linear  interactions, based
on  the Lennard-Jones  potential.  The  model atomic  potential energy
considered in the  present paper does not correspond  necessarily to a
real material.  Future  work should focus on  realistic potentials for
diamond and  silica, for example.  Further  comparisons between theory
and  MD simulations  are strongly  constrained  to a  small window  of
system sizes: smaller sizes may not be ergodic, larger body sizes lead
to  extremely long  correlation times.   It is  possible to  play with
temperature,  as   we  have  observed  that   orientational  diffusion
increases with the temperature of the  body.  In the present paper, in
order to have observable effects that allow for a quantitative precise
estimation  of d_α  we have  considered very  high temperatures
(k_BT≃9ϵ) that  may not be entirely  realistic.  Because
the model selected does not correspond necessarily to a real material,
we are not  yet able to provide accurate estimates  of the time scales
involved  in the  precession  relaxation rate  for  real materials  of
realistic sizes.

The present theory describes  orientational diffusion due to intrinsic
thermal fluctuations.  The phenomena is very different from rotational
diffusion due  to the  interaction with an  ambient gas  that produces
additional friction. This latter phenomenon  should be modelled with a
Brownian          rotor         with          viscous         friction
<cit.>  or with
kinetic  theory  <cit.>.   For   a  dense  gas,  it  is
anticipated  that the  friction with  the gas  will dominate  over the
intrinsic orientational  diffusion.  It is  an open question  at which
gas densities the intrinsic rotational diffusion, as discussed in this
paper, will be surpassed by  the effects of gas friction.  Interaction
with radiation  as it is encountered  in optolevitodynamic experiments
is  also a  source  of stochasticity  due to  the  discrete nature  of
photons. The  comparison of the  magnitude of these  different effects
clearly deserves further study.  The  significance of the present work
lies in its  validation of the theoretical model by  comparing it with
Molecular  Dynamics   (MD)  simulations   in  situations   where  such
comparisons are feasible.


















§ ACKNOWLEDGMENTS

We  thank  Mark  Thachuck  for  his  useful  comments  regarding  this
manuscript.  This  research has  been  supported  through MCIN  grants
PDC2021-121441-C22  and  PID2020-117080RB-C54.  JSR  acknowledges  the
support of  the Ministerio de  Ciencia, Innovación y  Universidades of
Spain  under  a  Margarita  Salas  contract  funded  by  the  European
Union-NextGenerationEU.   We acknowledge  the computational  resources
and  assistance  provided  by  the   Centro  de  Computación  de  Alto
Rendimiento CCAR-UNED.



§ APPENDIX: BROWNIAN MOTION ON A SPHERE


In  this  appendix, we  summarize  some  mathematical results  on  the
Brownian  motion of  a particle  on the  surface of  a sphere  and its
anisotropic generalization.



 §.§ Ito vs Stratonovich

We  first  recall the  connection  between  the Ito  and  Stratonovich
interpretation of an SDE.  The following Ito SDE

    d x= A( x)dt+Θ( x) d W_t

with A( x)  the drift, Θ( x) the
noise  amplitude matrix,  and  d W_t  a  vector of  independent
increments  of  the  Wiener  process,  corresponds  to  the  following
Stratonovich SDE <cit.>

    d x=[ A( x)- V( x)]dt+Θ( x)∘ d W_t

where the stochastic drift is 

    V_i
         =1/2Θ_kj∂/∂ x_kΘ_ij

and  repeated  indices are  summed  over.   Both SDE  (<ref>),
(<ref>) correspond to the following Fokker-Planck Equation (FPE)

    ∂_tP( x,t)
         =-∂/∂ x AP( x,t)
    +1/2∂/∂ x∂/∂ xΘΘ^T P( x,t)





 §.§ Anisotropic Brownian motion on the sphere

It is well-known <cit.>  that the evolution of
a particle with position r∈ℝ^3 on the surface of a
sphere due to Brownian motion is described by the following Stochastic
Differential Equation (SDE)

    d r   =[ r]_×∘ d W_t

with the corresponding SDE in the Ito interpretation

    d r   =- r+[ r]_× d W_t

where  d  W_t∈ℝ^3  is  a  a  vector  of  independent
increments  of  the  Wiener  process and  [  r]_×  is  the
cross-product matrix formed from the  vector r.  

A natural generalization of  the SDE (<ref>) for describing
 non-isotropic Brownian motion on the sphere is

    d r   =[ r]_× C∘ d W_t

where C∈ℝ^3×3 is a constant symmetric  matrix.
We write this equation in the form (<ref>)

    d r=Θ( r)∘ d W_t

where

    Θ( r)=[ r]_× C

and C is a constant matrix.  The
corresponding Ito SDE  is given by

    d r   = V dt + Θ d W_t

where V is given in (<ref>). Noting that the matrix Θ has components

    Θ_ij=ϵ_ilm r_l C_mj

then

    V_i( x)
         =1/2Θ_kj∂/∂ r_kΘ_ij
       =1/2ϵ_kl'm' r_l' C_m'j∂/∂ r_kϵ_ilm r_l C_mj
     
     
     
     
     
     
     
     
     
     
     
        
        =-1/2( Tr[ C C^T] r_i
        - C_ij C_mj r_m)

which means

    V   =  -1/2( Tr[ D]1 - D) r,

where we have defined the symmetric positive definite diffusion tensor 

    D≡ C C^T

and the  superscript T
denotes the transpose  matrix. The Ito SDE  (<ref>) is then

    d r   =- A r+
                [ r]_× C d W_t

where the matrix A∈ℝ^3×3 is

    A≡1/2( Tr[ D]1 - D)

Here 1 is the unit matrix.  The matrix A is symmetric
and positive  definite. This  can be  easily seen  from the  fact that
A, D commute and, hence diagonalize in the same basis. If
the eigenvalues of D are (d_1,d_2,d_3), then the eigenvalues
of      A      are       easily      seen      to      be
(d_2+d_3/2,d_2+d_3/2,d_1+d_2/2). Because
d_α>0, all  the eigenvalues  of A are  positive.  When
C=1=  D,   then   A=1,   and
(<ref>)  recovers  the  SDE (<ref>)  for  isotropic
Brownian   motion.



  

 §.§ The FPE for anisotropic Brownian motion

  To get  the FPE that corresponds  to the Ito SDE  (<ref>) we
  translate the connection between the FPE (<ref>) and the Ito SDE
  (<ref>), this is

    ∂_tP( r,t)
         =∂/∂ r[1/2( Tr[ D]1 - D) r]
        P( r,t)
        +1/2∂/∂ r∂/∂ r[ r]_× C C^T[ r]^T_× P( r,t)
      
      
      
      
      
      
        
       =∂/∂ r[1/2( Tr[ D]1 - D) r]
        P( r,t)
        +1/2∂/∂ r
       [ r]_× D [ r]^T_×∂/∂ rP( r,t)
        +1/2∂/∂ r(
        ∂/∂ r[ r]_× D [ r]^T_×)P( r,t)


Let us compute the term within rounded parenthesis

    (
      ∂/∂ r_j[[ r]_× D [ r]^T_×]_ij)
         =-
        ∂/∂ r_jϵ_iki' D_i'j'ϵ_j'k'j∂/∂ r_j r_k r_k'
        =-ϵ_iki'ϵ_j'k'j D_i'j'(δ_kj r_k'+δ_k'j r_k)
        
        =-ϵ_iji'ϵ_j'k'j D_i'j' r_k'-ϵ_iki'ϵ_j'jj_=0 D_i'j' r_k
      
       =ϵ_ii'jϵ_j'k'j D_i'j' r_k'
        =[δ_ij'δ_i'k'-δ_ik'δ_i'j']
         D_i'j' r_k'
        =δ_ij'δ_i'k' D_i'j' r_k'-δ_ik'δ_i'j' D_i'j' r_k'
        = D_i'i r_i'- Tr[ D] r_i


The last term in the rhs of (<ref>) then cancels the first term
and the FPE (<ref>) for anisotropic Brownian motion is simply

    ∂_tP( r,t)
         =
        1/2∂/∂ r
       [ r]^T_× D [ r]_×∂/∂ rP( r,t)

The diffusion matrix of this  FPE is clearly positive (semi) definite,
because, for any arbitrary vector v we have

    v^T    [ r]^T_× D [ r]_× v   =   v^T    [ r]^T_× C C^T [ r]_× v
           =( C^T [ r]_× v)^2≥0




 §.§ The equilibrium distribution for spherical Brownian motion

Let us  find the equilibrium  distribution of the  FPE (<ref>).
Observe that  any probability of  the form P(  r)=ϕ(r), with
r=|  r|, is  a  stationary solution  of  the FPE  (<ref>)
because

    ∂/∂ rϕ(r)=ϕ'(r) r/r

At   the  same   time,   any   average  of   a   function  F(r)   is
time-independent. This is shown as follows

    d/dt∫ d rF(r)P( r,t)
        
       =∫ d rF(r)
        1/2∂/∂ r
        [ r]^T_× D [ r]_×∂/∂ rP( r,t)  
        
       =- 1/2∫ d r∂ F(r)/∂ r
        [ r]^T_× D [ r]_×∂/∂ rP( r,t)  
        =0    
       =- 1/2∫ d rF'(r)
         r/r
        [ r]^T_×_=0 D [ r]_×∂/∂ rP( r,t)  
        =0

Now  consider  the probability  that  the  particle has  a  particular
modulus r. By definition it is given by the pushforward

    P(r,t)=∫ d rδ(r-| r|)P( r,t)

This  quantity  is, in  fact,  time-independent.  It takes  the  value
P(r,0)  that  it has  at  the  initial  time  and, also,  the  value
P(r,∞)=P(r,0) that has at equilibrium. Therefore, we have

    P(r,0)   =∫ d rδ(r-| r|)P^ eq( r)
              
       =ϕ(r)∫ d rδ(r-| r|)=ϕ(r)4π r^2

This means that the equilibrium solution must be given by

    P^ eq( r)=P(r,0)/4π r^2

Observe  that the  stationary solution  of  the FPE  depends on  the
initial distribution. For example, if the initial condition is  one in which
we are  certain that  the particle is  at a particular point  r_0  on the
surface  of  the  unit  sphere,  i.e.  |  r_0|=1,  the  initial
condition   is   P(   r,0)=δ(   r-   r_0).   The
corresponding  stationary  solution  of  the  FPE  with  this  initial
condition is

    P^ eq( r)=δ(| r|-1)/4π




 §.§ Correlation function

The equilibrium correlation  function of  the position of  a particle  perfoming a
Brownian motion on  the unit sphere can be  analytically computed, as follows. The
FPE (<ref>) can be written as

    ∂_tP( r,t)
         = L
        P( r,t)

where the Fokker-Planck operator is given by

    L   =1/2∂/∂ r
        [ r]^T_× D [ r]_×∂/∂ r

The stationary correlation function can be expressed as <cit.>

    r_α(t) r_β =
      ∫ d rP^ eq( r) r_β e^ Lt r_α

where the exponential of the Fokker-Planck operator is defined in terms of its series expansion

    e^ Lt=∑_n=0^∞t^n/n! L^n

We wish to compute the action  of the exponential operator on the unit
vector,  e^  Lt  r.   The second  term  in  the  series
expansion is L r that we write  explicitly in component
form as

    L r_γ   =-1/2∂/∂ r_μϵ_μαμ' r_α D_μ'ν'ϵ_ν'βν r_β∂/∂ r_ν r_γ
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
       
       =1/2 D_μ'ν'ϵ_μαμ'ϵ_μν'γ r_α
       =1/2 D_μ'ν'(δ_αν'δ_μ'γ-δ_αγδ_μ'ν')
         r_α
       =1/2 D_γα r_α
        -1/2 Tr[ D]  r_γ

Therefore, 

    L r=- A r

where the  matrix A  has been  defined in (<ref>).   It is
then obvious that

    e^ Lt r=e^- At r

where the matrix exponential is introduced through the series expansion.

The correlation matrix (<ref>) is now

    r_α(t) r_β   =
        ∫ d rP^ eq( r) r_α[e^- At]_ββ' r_β'
       =
        [e^- At]_ββ'∫ d rP^ eq( r) r_α r_β'
       (<ref>)=[e^- At]_ββ'1/3δ_αβ'

and we  conclude that the  equilibrium time correlation matrix  of the
position of the Brownian particle on the unit sphere is given by

    r(t) r^T=1/3e^- At




 §.§ Representation in terms of a rotation matrix

We consider now a time-dependent rotation matrix R

    R   
        =([ c_1, c_2, c_3 ])

where  c_α  are   column  vectors.

Because R is  an orthogonal matrix,  its columns
c_α form  an orthonormal basis set. We propose  now a SDE
for the rotation matrix in such  a way that the basis vectors describe
a Brownian motion on the unit sphere.  In order to reach this goal, we
introduce   the  angular   velocity  ω_0   in  the
principal axis frame through

    d/dt R   ≡-[ω_0]_× R

to be compared with (<ref>).  
Both angular velocities (<ref>) and (<ref>) are related through

    ω_0= Rω

We may use  (<ref>) as an inspiration  to construct an
SDE in the space of rotation  matrices.  To achieve this objective, we
transform   the  angular   velocity  ω_0   into  a
stochastic process, this is

    ω_0    → C∘d W_t/dt

and postulate the  following Stratonovich  SDE for
the rotation matrix

    d R   =-[ C∘ d W_t]_× R

In  terms of  the
columns c_α, we may write (<ref>) in the form

    d c_1   =-[ C∘ d W_t]_× c_1
                  
    
      d c_2   =-[ C∘ d W_t]_× c_2
                  
    
      d c_3   =-[ C∘ d W_t]_× c_3

where in each equation we have  exactly the same value of the independent increment of the Wiener process d W_t. Therefore, the
columns  of the  rotation  matrix experience  an anisotropic  Brownian
motion   on   the  surface   of   the   unit   sphere  of   the   form
(<ref>). We conclude that  the postulated SDE (<ref>)
for the rotation matrix is an alternative representation of anisotropic
Brownian motion.

Observe  that the  SDE (<ref>)  has a  number of  conserved
quantities. By scalarly multiplying with c_α we have

    c^T_α d c_β   = c^T_α( c_β× C∘ d W_t)
        =( c_α× c_β )^T  C∘ d W_t

and this implies

    d( c_α^T c_β)
     =     c^T_α d c_β+   c^T_β d c_α
       =
      ( c_α× c_β)^T   C∘ d W_t
        +( c_β× c_α)^T  C∘ d W_t
     =0

This means that both the modulus  of the vectors and the angle between
them    are    conserved    by     the    dynamics.    If    initially
c^T_α    c_β=δ_αβ,   these
orthonormality conditions are mantained at all times, as expected.

For  future  reference, we  also  construct  the  Ito version  of  the
Stratonovich  SDE  <ref>).   By  analogy  with  (<ref>),
(<ref>), the  Stratonovich SDE (<ref>) is  equivalent to
the following Ito SDE

    d c_α   =- A c_α
        -[ C d W]_× c_α

where A is given  by (<ref>).  Equation (<ref>) is
equivalent to the following Ito SDE for the rotation matrix

    d R   =- A R
        -[ C d W]_× R

This concludes our review of the anisotropic Brownian motion of a particle on the unit sphere. 






§ APPENDIX: DERIVATION OF (<REF>) 

We start from the Stratonovich SDE (<ref>) for the orientation
  
    dΛ   =k_BT F^ thdt +      B^T C dW̃

where 

    C   ≡  (2k_BT)^1/2 D_0^1/2

This is of the form of the Ito SDE (<ref>) with

    A   =k_BT F^ th
    Θ   = B^T(2k_BT)^1/2 D_0^1/2

The corresponidng Stratonovich SDE (<ref>) is now 

    dΛ   =(k_BT F^ th- V)dt +      B^T C∘ dW̃

where the stochastic drift V is given  by (<ref>), which in the present case is

    V_i   =k_BT[ B R^T D_0^1/2]_kj∂/∂Λ_k[ B R^T D_0^1/2]_ij
       =k_BT [ B R^T]_kk'[ D_0]_k'i'∂/∂Λ_k[ B R^T]_ii'

The rotation matrix and the kinematic operator satisfy

    B R^T= B^T

as can  be shown  from the  results in Sec.   K.1 of  the Supplemental
Material  of   Ref.   <cit.>.    By  using   the  property
(<ref>), (<ref>) becomes

    V_i
       =k_BT [ D_0]_k'i' B_k'k∂/∂Λ_k B_i'i(<ref>)=k_BT F^ th_i

Therefore, (<ref>) becomes

    dΛ   =      B^T C∘ dW̃

which is (<ref>), as we wished to demonstrate.


§ APPENDIX: DERIVATION OF (<REF>) 


By using the chain rule of ordinary calculus with the Stratonovich SDE
(<ref>)  we wish  to  obtain a  Stratonovich  SDE for  the
rotation matrix. Indeed

    d R   =∂ R/∂Λ dΛ(<ref>)=∂ R/∂Λ B^T C∘ dW̃

From the  results in Sec.   K.1 of  the Supplemental Material  of Ref.
<cit.>, the derivative of the rotation matrix with respect
to the orientation is given by
  
    ∂ R_μν/∂Λ_α   =- B^-1_αβϵ_μβδ R_δν

Using this derivative in (<ref>) leads to 

    d R_μν   =  ∂ R_μν/∂Λ_α B_αα' R^T_α'β[ C dW̃_t]_β
       =-   B^-1_αβϵ_μβδ R_δν B_αα' R^T_α'β'[ C dW̃_t]_β'
       =-ϵ_μβδ R_δν B^-1_αβ B_αα' R^T_α'β'_[ B^-T· B· R^T]_ββ'[ C dW̃_t]_β'

The property (<ref>) implies 

    B^-T B R^T=1

and, therefore,

    d R_μν   =-  ϵ_μβδ R_δν B^-1_αβ B_αα' R^T_α'β'_δ_ββ'[ C dW̃_t]_β'
       =      -ϵ_μαβ R_βν[ C d W_t]_α

which is the desired result (<ref>).


36
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Goldstein(1983)]goldsteinClassicalMechanics1983
authorH. Goldstein,
  titleClassical Mechanics
  (publisherAddison-Wesley, Massachusetts, year1983).

[V. I. Arnold(1989)]v.i.arnoldMathematicalMethodsClassical1989
authorK. V. V. I. Arnold,
  A. Weinstein, titleMathematical Methods Of
  Classical Mechanics, Graduate Texts in Mathematics
  (publisherSpringer, year1989), edition2nd
  ed.

[Lamy and Burns(1972)]Lamy1972
authorP. L. Lamy and
  authorJ. A. Burns,
  journalAmerican Journal of Physics volume40,
  pages441 (year1972).

[Warner et al.(2009)Warner, Harris, and
  Pravec]Warner2009
authorB. D. Warner,
  authorA. W. Harris,
  and authorP. Pravec,
  journalIcarus volume202,
  pages134 (year2009).

[Breiter et al.(2012)Breiter, Rozek, and
  Vokrouhlický]Breiter2012
authorS. Breiter,
  authorA. Rozek, and
  authorD. Vokrouhlický,
  journalMonthly Notices of the Royal Astronomical Society
  volume427, pages755 (year2012).

[Descamps(2008)]Descamps2008
authorS. X. Descamps,
  journalQuaderns d'Historia de l Enginyeria
  volume9, pages270 (year2008).

[Gautschi(2008)]Gautschi2008
authorW. Gautschi,
  journalSIAM review volume50,
  pages3 (year2008).

[Español et al.(2024)Español,
  Thachuk, and De La Torre]espanol2024
authorP. Español,
  authorM. Thachuk, and
  authorJ. De La Torre,
  journalEuropean Journal of Mechanics - A/Solids
  volume103, pages105184
  (year2024).

[Delong et al.(2014)Delong, Sun, Grittith, Vanden-Eijnden, and Donev]Delong2014
  authorS. Delong,
  authorY. Sun,
  authorB.E. Griffith,
  authorE. Vanden-Eijnden, and
  authorA. Donev,
  journalPhysical Review E
  volume90, pages063312
  (year2014).

[Einstein(1905)]Einstein1905
authorA. Einstein,
  journalAnn. Phys. (Leipzig) volume19,
  pages549 (year1905).

[Green(1952)]Green1952
authorM. Green, journalJ.
  Chem. Phys. volume20, pages1281
  (year1952).

[Zwanzig(1961)]Zwanzig1961
authorR. Zwanzig,
  journalPhysical Review volume124,
  pages983 (year1961).

[Grabert(1982)]Grabert1982
authorH. Grabert,
  titleProjection Operator Techniques in Nonequilibrium
  Statistical Mechanics (publisherSpringer,
  year1982).

[Öttinger(2005)]Ottinger2005
authorH. C. Öttinger,
  titleBeyond Equilibrium Thermodynamics
  (publisherJ. Wiley & Sons, year2005).

[Gibbs(1960)]gibbsElementaryPrinciplesStatistical1960
authorJ. W. Gibbs,
  titleElementary Principles in Statistical
  Mechanics (publisherYale Univ. Press, 1902. Dover,,
  addressNew York,, year1960).

[De La Torre and Español(2024)]delatorre2024
authorJ. De La Torre and
  authorP. Español,
  journalEuropean Journal of Mechanics - A/Solids p.
  pages105298 (year2024).

[Poinsot and Whitley(2022)]poinsot2022outlines
authorL. Poinsot and
  authorC. Whitley,
  titleOutlines of a New Theory of Rotatory Motion
  (publisherCreative Media Partners, LLC, year2022).

[Landau and
  Lifshitz(1960)]landauMechanicsThirdEdition1960
authorL. D. Landau and
  authorE. M. Lifshitz,
  titleMechanics (Third Edition)
  (publisherPergamon Press, year1960).

[Ashbaugh et al.(1991)Ashbaugh, Chicone,
  and Cushman]Ashbaugh1991
authorM. S. Ashbaugh,
  authorC. C. Chicone,
  and authorR. H.
  Cushman, journalJournal of Dynamics and
  Differential Equations volume3, pages67
  (year1991).

[Saporta Katz and Efrati(2019)]SaportaKatz2019
authorO. Saporta Katz
  and authorE. Efrati,
  journalPhysical Review Letters volume122,
  pages024102 (year2019).

[Price and Williams(1983)]price1983
authorG. C. Price and
  authorD. Williams, in
  booktitleSéminaire de Probabilités XVII
  1981/82, edited by
  editorJ. Azéma and
  editorM. Yor
  (publisherSpringer Berlin Heidelberg, addressBerlin,
  Heidelberg, year1983), vol. volume986, pp.
  pages194–197.

[Van Den Berg and Lewis(1985)]vandenberg1985
authorM. Van Den Berg
  and authorJ. T. Lewis,
  journalBulletin of the London Mathematical Society
  volume17, pages144 (year1985).

[Perrin(1934)]perrin1934
authorF. Perrin,
  journalJournal de Physique et le Radium
  volume5, pages497 (year1934).

[Furry(1957)]furry1957
authorW. H. Furry,
  journalPhysical Review volume107,
  pages7 (year1957).

[Favro(1960)]favro1960
authorL. D. Favro,
  journalPhysical Review volume119,
  pages53 (year1960).

[Hubbard(1972)]hubbard1972
authorP. S. Hubbard,
  journalPhysical Review A volume6,
  pages2421 (year1972).

[Höfling and Straube(2024)]hofling2024
authorF. Höfling and
  authorA. V. Straube,
  titleLangevin equations and a geometric integration scheme
  for the overdamped limit of homogeneous rotational Brownian motion
  (year2024), arXiv 2403.04501.

[Díaz(2019)]diaz2019
authorE. O. Díaz,
  title3D Motion of Rigid Bodies A Foundation for
  Robot Dynamics Analysis (publisherSpringer,
  year2019).

[Thompson et al.(2022)Thompson, Aktulga,
  Berger, Bolintineanu, Brown, Crozier, in 't Veld, Kohlmeyer, Moore, Nguyen
  et al.]LAMMPS
authorA. P. Thompson,
  authorH. M. Aktulga,
  authorR. Berger,
  authorD. S. Bolintineanu,
  authorW. M. Brown,
  authorP. S. Crozier,
  authorP. J. in 't Veld,
  authorA. Kohlmeyer,
  authorS. G. Moore,
  authorT. D. Nguyen,
  et al., journalComp. Phys. Comm.
  volume271, pages108171
  (year2022).

[Pedregosa et al.(2011)Pedregosa,
  Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss,
  Dubourg et al.]scikit-learn
authorF. Pedregosa,
  authorG. Varoquaux,
  authorA. Gramfort,
  authorV. Michel,
  authorB. Thirion,
  authorO. Grisel,
  authorM. Blondel,
  authorP. Prettenhofer,
  authorR. Weiss,
  authorV. Dubourg,
  et al., journalJournal of Machine Learning Research
  volume12, pages2825 (year2011).

[Fermi et al.()Fermi, Pasta, and
  Ulam]fermi
authorE. Fermi,
  authorJ. Pasta, and
  authorS. Ulam,
  titleStudies of nonlinear problems. I. Los Alamos
  report LA-1940 (1955), published later in Collected Papers of
  Enrico Fermi, E. Segré (Ed.) (University of Chicago
  Press, Chicago 1965); also in Non- linear Wave Motion, Newell
  A. C. Ed., Lecture in Applied Mathematics 15 (AMS,
  Providence, Rhode Island, 1974); also in The Many-Body Problem,
  Mattis C. C. Ed. (World Scientific, Singapore, 1993).

[Berman and Izrailev(2005)]Berman2005
authorG. P. Berman and
  authorF. M. Izrailev,
  journalChaos: An Interdisciplinary Journal of Nonlinear Science
  volume15, pages15104 (year2005).

[Dauxois et al.(2005)Dauxois, Peyrard,
  and Ruffo]Dauxois2005
authorT. Dauxois,
  authorM. Peyrard, and
  authorS. Ruffo,
  journalEuropean journal of physics volume26,
  pagesS3 (year2005).

[Dauxois(2008)]Dauxois2008a
authorT. Dauxois,
  journalPhysics Today volume61,
  pages55 (year2008).

[Martinetz et al.(2018)Martinetz,
  Hornberger, and Stickler]martinetz2018a
authorL. Martinetz,
  authorK. Hornberger,
  and authorB. A.
  Stickler, journalPhysical Review E
  volume97, pages052112
  (year2018).

[Gardiner(1983)]Gardiner1983
authorC. W. Gardiner,
  titleHandbook of Stochastic Methods
  (publisherSpringer-Verlag, Berlin, year1983).

[Risken(1984)]riskenFokkerPlanckEquation1984
authorH. Risken,
  titleThe Fokker-Planck Equation
  (publisherSpringer Verlag, Berlin, year1984).



