# Modulation of LISA free-fall orbits due to the Earth-Moon system

Massimo Cerdonio1, Fabrizio De Marchi2777Corresponding author:
<EMAIL_ADDRESS>Roberto De Pietri3, Philippe Jetzer4, Francesco
Marzari1, Giulio Mazzolo5, Antonello Ortolan6 and Mauro Sereno7,8 1 Department
of Physics, University of Padova and INFN Padova, via Marzolo 8, I-35131
Padova, Italy  2 Department of Physics, University of Trento and INFN Trento,
I-38100 Povo (Trento), Italy  3 Department of Physics, University of Parma and
INFN Parma I-43100 Parma, Italy  4 Institute of Theoretical Physics,
University of Zürich, Winterhurerstrasse 190, 8057 Zürich, Switzerland  5 Max
Planck Institut für Gravitationsphysik, Callinstrasse 38, 30167 Hannover,
Germany  6 INFN Laboratori Nazionali di Legnaro, Viale dell’Università 35020
Legnaro (Padova), Italy  7 Dipartimento di Fisica, Politecnico di Torino,
Corso Duca degli Abruzzi 24, 10129 Torino, Italia  8 INFN, Sezione di Torino,
Via Pietro Giuria 1, 10125, Torino, Italia

###### Abstract

We calculate the effect of the Earth-Moon (EM) system on the free-fall motion
of LISA test masses. We show that the periodic gravitational pulling of the EM
system induces a resonance with fundamental frequency 1 $yr^{-1}$ and a series
of periodic perturbations with frequencies equal to integer harmonics of the
synodic month ($\simeq 3.92\times 10^{-7}\ Hz$). We then evaluate the effects
of these perturbations (up to the 6th harmonics) on the relative motions
between each test masses couple, finding that they range between $3\ mm$ and
$10\ pm$ for the 2nd and 6th harmonic, respectively.

If we take the LISA sensitivity curve, as extrapolated down to $10^{-6}\ Hz$
in [1], we obtain that a few harmonics of the EM system can be detected in the
Doppler data collected by the LISA space mission. This suggests that the EM
system gravitational near field could provide an additional crosscheck to the
calibration of LISA, as extended to such low frequencies.

###### pacs:

04.80.nn , 95.10.Eg

##  1 Introduction

LISA (Laser Interferometer Space Antenna) is a ten years long NASA-ESA space
mission to detect gravitational waves in the frequency range $10^{-4}-10^{-1}\
Hz$ [2]. It consists of three spacecrafts whose mutual distances are about
$L=5\times 10^{6}\ km$. The LISA constellation will orbit around the Sun
following the same path of the Earth, $\phi_{0}=20^{\circ}$ behind [3].

The ideal configuration for LISA performances should be a rigid equilater
triangle [4]; however, the shape of the LISA constellation is subject to
significant variations because of the gravitational interaction due to Sun and
planets [5].

Much smaller perturbations can be induced by the presence of interplanetary
dust [6] or dark matter in the solar system [7]. The perturbations due to each
celestial body can be treated, at a first approximation, independently. The
gravitational effects are quite different in intensity (orders of magnitude)
and/or in behavior (stationary or time-dependent); the frequencies involved
are, in general, not commensurable and so resonance effects are not observed.
The only significative exception is the Earth perturbation [5], [8] and [9],
which gives a resonance because of the 1:1 commensurability between the Earth
and LISA orbits.

However, such a resonance will not have enough time to grow during the 10
years of LISA mission. In this paper we focus on the perturbations induced by
the EM system on LISA at the frequency of the synodic month and its harmonics,
which are much higher than $1\ yr^{-1}$ and so, in a first order
approximation, perturbative effects can be treated independently. We extend
the approach of [5] to include the time dependent perturbation of the Moon.
The plan of the paper is as follows. In Section 2 we shortly describe the
perturbative approach used to study the EM effects. Section 3 is devoted to
illustrate the EM system and the approximations we used to model its
gravitational near field. In Section 4 we calculate the modulations of the
distance between two LISA test masses due to EM system. In Section 5 we
estimate the perturbations induced by Venus and Jupiter, compared to those due
to the EM system calculated in the previous Sections. In Section 6 conclusions
are drawn and the future research potential of method for further studies and
application to LISA are given.

##  2 Perturbative dynamics of the LISA test masses

In order to simplify the notation, we use the Astronomical System of Units,
for length (AU), mass (M⊙) and time (days). However, the quantities that
affect the relative motion of the LISA spacecrafts will be reconverted in the
SI units.

By means of a F77 code based on the inverse 15th-order Runge-Kutta method
[10], we calculate the modulus of the force gradient between 2 LISA
spacecrafts due to the main Solar System bodies, i.e. Sun, Venus, Earth (more
precisely the EM system) and Jupiter (Figure 1). The effects are different by
orders of magnitude in amplitude and show also different frequencies.

In fact, the Earth and Jupiter tidal effects are $10^{4}$ and $10^{5}$ times
smaller than the Sun contribution, respectively; after 3 years, Sun, Earth and
Jupiter cause arms length changes (the so-called ”arm breathing”) of $\simeq
1.2\times 10^{5}\,km$, 4.8$\times 10^{4}\ km$ and 4$\times 10^{3}\ km$,
respectively [4], [5].

The Venus contribution oscillates in time by two orders of magnitude and it is
comparable with Earth effects for a short time interval every 584$\ d$ [11].
However, as we will show in Section 5, this perturbation is negligible with
respect to the Earth one.

Figure 1:  Gravity gradients between two spacecraft due to EM system (solid
line), Venus (dotted line), Jupiter (short dashed line) and Sun (long dashed
lines).

We therefore conclude that the celestial bodies affect the LISA arms length by
at most a few percent. As a consequence, at any time the distance of each
satellite from the LISA barycenter is $\rho_{0}=L/\sqrt{3}\simeq 1.9\times
10^{-2}\ AU$, and the distance between the Earth and the LISA barycenter is
$r_{g}=2\,\sin(\phi_{0}/2)\simeq 0.347\ AU$, within a relative fluctuation of
few percent.

###  2.1 Hill-Clohessy-Wiltshire (HCW) reference frame

To study the EM system effect, we made use of the Hill-Clohessy-Wiltshire
(HCW) reference frame $\\{x,y,z\\}$ [12] defined as follows (see Figure 2):

  1. 1.

the origin $O^{\prime}$ describes a circular orbit on the ecliptic plane at
$1\ AU$ from the Solar System Barycenter $O$ [13];

  2. 2.

the $xy$ plane coincides with the ecliptic;

  3. 3.

the $z$-axis is perpendicular to the ecliptic plane and parallel to the Solar
System total angular momentum111The Solar System total angular momentum is
perpendicular to a plane slightly inclined with respect to the ecliptic one.
For our scopes, such inclination is negligible.;

  4. 4.

the $x$-axis is tangent to the orbit and is antiparallel to the origin
$O^{\prime}$ velocity vector;

  5. 5.

the $y$-axis is directed radially outward.

Figure 2:  Hill-Clohessy-Wiltshire frame: $R_{0}$ is the radius of the orbit,
$\rho$ and $R$ are the the distances of $P$ from the origin of the rotating
and fixed frames, respectively. $\\{X,Y,Z\\}$ is the associated reference
frame (see details in Section 3). The EM system is also represented. The
figure is not to scale.

In this frame the coordinates of the EM system barycenter are
$(x_{g},y_{g},z_{g})=R_{0}(-\sin\phi_{0},\cos\phi_{0}-1,0)$. Of course, the
LISA spacecrafts can be considered as three proof masses. Their equations of
motion in the HCW frame read

$\eqalign{\ddot{x}-2\omega_{0}\dot{y}-\omega_{0}^{2}x=-\frac{\mu}{R^{3}}x\cr\ddot{y}+2\omega_{0}\dot{x}-\omega_{0}^{2}(y+R_{0})=-\frac{\mu}{R^{3}}(R_{0}+y)\cr\ddot{z}=-\frac{\mu}{R^{3}}z\
,}$ (1)

where $\mu=GM_{\odot}=\omega_{0}^{2}R_{0}^{3}$,
$\omega_{0}=2\pi/365.257\,d^{-1}$, $R_{0}=1\,AU$, and
$R=\sqrt{x^{2}+(y+R_{0})^{2}+z^{2}}$. Being the breathing length $\Delta L\ll
L$, each satellite is located at any time at a distance
$\rho(t)=\sqrt{x(t)^{2}+y(t)^{2}+z(t)^{2}}\simeq\rho_{0}$. Since $\rho_{0}\ll
R_{0}$, we expand the acceleration due to the Sun in terms of
$x/R_{0},y/R_{0},z/R_{0}$. Retaining the first term of the series, we have the
so-called Hill-Clohessy-Wiltshire equations [12] of relative motion

$\eqalign{\ddot{x}-2\omega_{0}\dot{y}=0\cr\ddot{y}+2\omega_{0}\dot{x}-3\omega_{0}^{2}y=0\cr\ddot{z}+\omega_{0}^{2}z=0\
,}$ (2)

with general solutions [14], [15]

$\eqalign{x(t)=x_{0}+2\
\frac{\dot{y}_{0}}{\omega_{0}}-3\,\left(\frac{\dot{x}_{0}}{\omega_{0}}-2y_{0}\right)t\,-2\frac{\dot{y}_{0}}{\omega_{0}}\cos\omega_{0}t+2\left(2\frac{\dot{x}_{0}}{\omega_{0}}-3y_{0}\right)\sin\omega_{0}t\cr
y(t)=2\left(2y_{0}-\frac{\dot{x}_{0}}{\omega_{0}}\right)+\left(2\frac{\dot{x}_{0}}{\omega_{0}}-3y_{0}\right)\cos\omega_{0}t+\frac{\dot{y}_{0}}{\omega_{0}}\sin\omega_{0}t\cr
z(t)=z_{0}\cos\omega_{0}t+\frac{\dot{z}_{0}}{\omega_{0}}\sin\omega_{0}t\ ,}$
(3)

where $x_{0},y_{0},z_{0}$ and $\dot{x}_{0},\dot{y}_{0},\dot{z}_{0}$ are the
initial positions and velocities respectively.

Since $x(t)$ contains a term proportional to $t$, after some time the
assumption $\rho\ll R_{0}$ is no more valid and the above approximation breaks
down. However, the divergent term can be cancelled by choosing
$\dot{x}_{0}=2\omega_{0}y_{0}$. For the LISA case, the constraints of rigid
and bounded relative motions lead to the solutions [4]

$\eqalign{x_{k}(t)=-\rho_{0}\sin\left[\omega_{0}t+\sigma_{k}\right]\cr
y_{k}(t)=-\frac{1}{2}\rho_{0}\cos\left[\omega_{0}t+\sigma_{k}\right]\cr
z_{k}(t)=-\frac{\sqrt{3}}{2}\,\rho_{0}\cos\left[\omega_{0}t+\sigma_{k}\right]\
,}$ (4)

where $\sigma_{k}=(1-k)\ 2\pi/3$ and $\,k=1,2,3$ is a label which enumerates
the LISA spacecrafts.

###  2.2 Rescaling and expansion of the HCW equations

For our calculation we rewrite (2) by means of the coordinate transformations
$\hat{x}=x/\rho_{0},\hat{y}=y/\rho_{0},\hat{z}=z/\rho_{0}$ and
$\hat{t}=\omega_{0}t$.

$\eqalign{{\hat{x}}^{\prime\prime}-2{\hat{y}}^{\prime}=0\cr{\hat{y}}^{\prime\prime}+2{\hat{x}}^{\prime}-3{\hat{y}}=0\cr{\hat{z}}^{\prime\prime}+{\hat{z}}=0}$
(5)

with the notation $\ {}^{\prime}=d/d\hat{t}$. We will refer to (5) as the HCW1
equations. We will show that the higher order terms in the expansions are
negligible for calculating the perturbation of the LISA rigid and bounded
orbits.

It is worth noticing that the right-hand side of (5) is zero only at the first
order in the force expansion. In general, the right-hand side of (5) is a
polynomial of degree $n$, where $n$ is the order of the expansion [16]. In our
case, the order of magnitude of the neglected terms is
$\varepsilon_{HCW}=\rho_{0}/R_{0}\simeq 1.9\times 10^{-2}$.

##  3 Gravitational near field of the EM system

In our model, the EM system is constituted by 2 point masses: $m_{1}\simeq
3.0\times 10^{-6}M_{\odot}$ (Earth) and $m_{2}\simeq 3.7\times
10^{-8}M_{\odot}$ (Moon) located at a constant distance $l\simeq 2.57\times
10^{-3}\ AU$ and describing a circular orbit with angular velocity
$\omega_{M}=2\pi/P_{M}$, where $P_{M}\simeq 29.53\ d$ is the synodic month,
around their common barycenter. We assume that the barycenter of the EM system
makes a circular orbit with radius $1\ AU$ around the Sun, i.e. we neglect the
eccentricity $\simeq 0.0167$ of the Earth orbit around the Sun.

In addition, we disregard the eccentricity of the Moon orbit around the Earth
($\simeq 0.054$), its inclination to the ecliptic plane ($\simeq
5.14^{\circ}$), the motion of the perigee of the Moon ($\simeq 8.85\ yr$) and
the precession of the Moon orbit plane ($\simeq 18.03\ yr$) [17], [18].

We consider a non rotating reference frame $\\{X,Y,Z\\}$ centered on the EM
system barycenter, with the $X$ axis along the line joining Earth and Moon at
$t=0$, and the $Z$ axis perpendicular to the ecliptic plane (see Figure 2).
The gravitational potential due to the Earth and the Moon at a point $(X,Y,Z)$
is given by

$\displaystyle
U=-\frac{Gm_{1}}{\sqrt{\left[X-\frac{m_{2}}{m_{1}+m_{2}}l\cos(\omega_{M}t)\right]^{2}+\left[Y-\frac{m_{2}}{m_{1}+m_{2}}l\sin(\omega_{M}t)\right]^{2}+Z^{2}}}+$
$\displaystyle-\frac{G{m_{2}}}{\sqrt{\left[X+\frac{m_{1}}{m_{1}+m_{2}}l\cos(\omega_{M}t)\right]^{2}+\left[Y+\frac{m_{1}}{m_{1}+m_{2}}l\sin(\omega_{M}t)\right]^{2}+Z^{2}}}\
.$

###  3.1 Multipole expansion of the EM gravitational potential

We are interested in the effect of the EM system on the LISA constellation,
located at a distance $r_{g}$. The size $l$ of the EM system is small
relatively to $r_{g}$, $l/r_{g}=7\times 10^{-3}$. We therefore expand the
total potential in series of $\varepsilon_{M}=l/r_{g}$.

$U(X,Y,Z,t)\equiv\sum_{n=0}^{\infty}\frac{1}{n!}\ \varepsilon_{EM}^{n}\
U_{n}(X,Y,Z,t)\ ,$ (6)

where $U_{n}$ are the well known multipole terms.

At the zeroth order we have the monopole term

$U_{0}(X,Y,Z)=-\frac{G({m_{1}}+{m_{2}})}{\sqrt{X^{2}+Y^{2}+Z^{2}}}\ .$ (7)

The first order term, the dipole term, is equal to zero, due to the
conservation of the linear momentum. The second order term, the quadrupole
term, is

$U_{2}(X,Y,Z,t)=\frac{m_{1}m_{2}}{m_{1}+m_{2}}\ r_{g}^{2}\
\left[\frac{X^{2}+Y^{2}+Z^{2}-3\left[X\cos\omega_{M}t+Y\sin\omega_{M}t\right]^{2}}{(X^{2}+Y^{2}+Z^{2})^{5/2}}\right]$
(8)

and so on.

Each $U_{n}$ term contains sinusoidal terms as $\sin k\omega_{M}t\ $ and $\cos
k\omega_{M}t$, with $k=0,2,\dots,n$ for even $n$ and $k=1,3,\dots,n$ for odd
$n$.

Defining $\hat{\omega}_{M}=\omega_{M}/\omega_{0}\simeq 12.3687$ and operating
the substitution $X=x+x_{g},\ Y=y+y_{g},\ Z=z$ in all $U_{n}$, we obtain the
EM potential in the HCW frame. The corresponding force per unit mass is

$\mathbf{F}_{n}=-\frac{1}{n!}\ \varepsilon_{EM}^{n}\,\nabla U_{n}(x,y,z,t)\ .$

The resulting functions $\mathbf{F}_{n}(x,y,z,t)$ are still too complex to be
treated analytically, but it is worth noticing that the values of the
coordinates range within $-\rho_{0}$ and $+\rho_{0}$, while $r_{g}$ is about
20 times larger ($\rho_{0}/r_{g}=\varepsilon_{L}\simeq 5.5\times 10^{-2}$). We
can therefore expand $\mathbf{F}_{n}$ in terms of $x/r_{g},y/r_{g},z/r_{g}$
around the origin of the HCW frame

$\mathbf{F}_{n}(x,y,z,t)=\sum_{m=0}^{\infty}\frac{\varepsilon_{L}^{m}}{m!}\
\left[\left(x\frac{\partial}{\partial\xi}+y\frac{\partial}{\partial\eta}+z\frac{\partial}{\partial\zeta}\right)^{m}\mathbf{F}_{n}(\xi,\eta,\zeta,t)\right]_{\xi=0,\eta=0,\zeta=0}\
.$ (9)

The above formula can be expressed in a more useful way (after the rescaling)

$\mathbf{F}_{n}(\hat{x},\hat{y},\hat{z},\hat{t})=\sum_{i=0}^{\infty}\frac{1}{i\
!}\sum_{j=0}^{\infty}\frac{1}{j\ !}\ \sum_{k=0}^{\infty}\ \frac{1}{k\ !}\
\varepsilon_{L}^{i+j+k}\ \mathbf{a}_{n,ijk}(\hat{t})\
\hat{x}^{i}\hat{y}^{j}\hat{z}^{k}\ ,$ (10)

where the $\mathbf{a}_{n,ijk}(\hat{t})$ are

$\mathbf{a}_{n,ijk}(\hat{t})=\sum_{m}\left[\mathchoice{\hbox{\boldmath$\displaystyle\alpha$}}{\hbox{\boldmath$\textstyle\alpha$}}{\hbox{\boldmath$\scriptstyle\alpha$}}{\hbox{\boldmath$\scriptscriptstyle\alpha$}}_{n,ijk,m}\sin
m\hat{\omega}_{M}\hat{t}+\mathchoice{\hbox{\boldmath$\displaystyle\beta$}}{\hbox{\boldmath$\textstyle\beta$}}{\hbox{\boldmath$\scriptstyle\beta$}}{\hbox{\boldmath$\scriptscriptstyle\beta$}}_{n,ijk,m}\cos
m\hat{\omega}_{M}\hat{t}\,\right]$ (11)

with
$\mathchoice{\hbox{\boldmath$\displaystyle\alpha$}}{\hbox{\boldmath$\textstyle\alpha$}}{\hbox{\boldmath$\scriptstyle\alpha$}}{\hbox{\boldmath$\scriptscriptstyle\alpha$}}_{n,ijk,m},\mathchoice{\hbox{\boldmath$\displaystyle\beta$}}{\hbox{\boldmath$\textstyle\beta$}}{\hbox{\boldmath$\scriptstyle\beta$}}{\hbox{\boldmath$\scriptscriptstyle\beta$}}_{n,ijk,m}$
numerical coefficients and $m=0,2,\dots,n$ for even $n$ and $m=1,3,\dots,n$
for odd $n$. Note that for even $n$, (11) contains a constant term (i.e.
$\mathchoice{\hbox{\boldmath$\displaystyle\beta$}}{\hbox{\boldmath$\textstyle\beta$}}{\hbox{\boldmath$\scriptstyle\beta$}}{\hbox{\boldmath$\scriptscriptstyle\beta$}}_{n,ijk,0}$)
plus sinusoids, while for odd $n$ it presents only sinusoids
($\mathchoice{\hbox{\boldmath$\displaystyle\beta$}}{\hbox{\boldmath$\textstyle\beta$}}{\hbox{\boldmath$\scriptstyle\beta$}}{\hbox{\boldmath$\scriptscriptstyle\beta$}}_{n,ijk,1}$
is multiplied by $\cos\,\hat{\omega}_{M}\hat{t}$).

Let us estimate the intensity of the acceleration due to the EM system : the
most important contribution is given by the monopole term ($n=0$), and we
assume its value at the origin of the coordinates as an indicator of its
intensity. This term is

$\varepsilon_{0}=\frac{G(m_{1}+m_{2})}{\omega_{0}^{2}\,r_{g}^{2}\,\rho_{0}}\simeq
1.31\times 10^{-3}\simeq 7\times 10^{-2}\varepsilon_{HCW}$

therefore the EM system influence is about the 7% of the contribution of the
Sun (Section 2.2).

At the first order in the coordinates the monopole term (normalized to
$\varepsilon_{0}$) is

$\displaystyle\eqalign{f_{x}=\frac{x_{g}}{r_{g}}-\hat{x}\
\frac{\rho_{0}(r_{g}^{2}-3x_{g}^{2})}{r_{g}^{3}}+\hat{y}\
\frac{3\rho_{0}x_{g}y_{g}}{r_{g}^{3}}\\\ f_{y}=\frac{y_{g}}{r_{g}}+\hat{x}\
\frac{3\rho_{0}x_{g}y_{g}}{r_{g}^{3}}-\hat{y}\
\frac{\rho_{0}(r_{g}^{2}-3y_{g}^{2})}{r_{g}^{3}}\\\ f_{z}=-\hat{z}\
\frac{\rho_{0}}{r_{g}}\ ,\\\ }$ (15)

where $\hat{x},\hat{y},\hat{z}$ are the scaled coordinates.

The most important multipole term is the quadrupole ($n=2$) which is mostly
constituted by the zeroth order term of its expansion in spatial coordinates,
or equivalently, the force at the origin. This term is periodic with period
$P_{M}/2=\pi/\omega_{M}$, and its mean value is

$\varepsilon_{2}=\frac{3}{4}\ \frac{G\,l^{2}}{\omega_{0}^{2}\ r_{g}^{4}\
\rho_{0}}\ \frac{m_{1}m_{2}}{m_{1}+m_{2}}\simeq 6.48\times 10^{-10}\ \simeq\
5\times 10^{-7}\ \varepsilon_{0}\ .$

The monopole force has an identical time-independent polynomial structure as
the expansion in coordinates of the Sun force. The difference is that this
latter starts with second order terms [16], while the EM monopole presents
also linear terms and a constant term. Moreover, the monopole force is quite
smaller than the Sun force.

The multipole force with even $n$ also contains time-independent terms but the
largest ones (corresponding to $n=2$) have order of magnitude
$\simeq\varepsilon_{2}\ll\varepsilon_{HCW}$.

We can separate the total force in a more useful way:

  1. 1.

A component $\mathbf{F_{0}}(\hat{x},\hat{y},\hat{z})$ independent of time due
to the terms of the Sun force ($\propto\varepsilon_{HCW}$) [16], plus the EM
system monopole ($\propto\varepsilon_{0}$, (10) with $n=0$) plus terms
contained in the $n$-even multipole expansion ($\propto\varepsilon_{2}$, (11)
with $m=0$);

  2. 2.

A component $\mathbf{F_{1}}(\hat{t})$ that depends only on time due to the EM
system multipole terms ($\propto\varepsilon_{2}$, (10) with $i,j,k=0$ and
$n\geq 2$). Physically this is the force at the origin of the HCW system;

  3. 3.

A component
$\mathbf{F_{2}}(\hat{x},\hat{y},\hat{z},\hat{t})\propto\varepsilon_{2}$ that
depends on both time and coordinates (all other cases: (10) with $n\geq 2$ and
$m\neq 0$).

##  4 Perturbation of the LISA orbits due to EM system

If we indicate with $\mathbf{r}_{0i}(t)$ the unperturbed trajectory of the
$i^{th}$ LISA test mass ($i=1,2,3$) and with $\mathbf{r}_{1i}(t)$ its small
perturbation, the difference between two perturbed trajectories is simply

$\Delta\mathbf{r}_{ij}(t)\equiv\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)=\Delta\mathbf{r}_{0ij}(t)+\Delta\mathbf{r}_{1ij}(t)\
,$

where $\Delta\mathbf{r}_{0ij}(t)=\mathbf{r}_{0i}(t)-\mathbf{r}_{0j}(t)$,
$\Delta\mathbf{r}_{1ij}(t)=\mathbf{r}_{1i}(t)-\mathbf{r}_{1j}(t)$.

The perturbation of the relative displacement between the pair $i,j$ of LISA
test masses can be written as

$\Delta
L_{ij}(t)\simeq\frac{\Delta\mathbf{r}_{0ij}\cdot\Delta\mathbf{r}_{1ij}}{L_{0ij}}\
,$ (16)

where $L_{0ij}=|\Delta\mathbf{r}_{0ij}|$ is the distance between $i$ and $j$
test masses in the unperturbed case, and
$|\mathbf{r}_{1i}|\ll|\mathbf{r}_{0i}|$. We also define the perturbation to
the differential distances between each pair of LISA as

$\delta L_{ijk}(t)\equiv\Delta L_{ij}(t)-\Delta L_{jk}(t)$ (17)

which represents the variation of $L_{ij}(t)-L_{jk}(t)$ due to a small
perturbation and relates directly to the LISA sensitivity curve.

###  4.1 Effects of the EM monopole

At the first order, under the effect of the Sun, the LISA motion is described
by (4) and the EM monopole perturbation is contained in the time-independent
force per unit mass $\mathbf{F_{0}}$.

We write the motion of the $k$-spacecraft under the effect of $\mathbf{F_{0}}$
in the following form:
$\hat{\mathbf{r}}_{k}=\hat{\mathbf{r}}_{0k}+\varepsilon_{HCW}\hat{\mathbf{r}}_{1k}+\varepsilon_{0}\hat{\mathbf{r}}_{2k}$,
where $\hat{\mathbf{r}}_{0k}$ is the unperturbed motion ((4), rescaled), and
$\hat{\mathbf{r}}_{1k}$ and $\hat{\mathbf{r}}_{2k}$ are the perturbations due
to the Sun force terms [16] and the EM monopole, respectively. Being
$\varepsilon_{HCW}^{2}\ll\varepsilon_{0}$, to calculate
$\hat{\mathbf{r}}_{2k}$ it is not necessary to know $\hat{\mathbf{r}}_{1k}$,
[8].

The equations for
$\hat{\mathbf{r}}_{2k}=(\hat{x}_{2k},\hat{y}_{2k},\hat{z}_{2k})$ are [5], [8]
and [9]

$\eqalign{{\hat{x}_{2k}}^{\prime\prime}-2{\hat{y}_{2k}}^{\prime}=f_{x}(\hat{x}_{0k},\hat{y}_{0k},\hat{z}_{0k})\cr{\hat{y}_{2k}}^{\prime\prime}+2{\hat{x}_{2k}}^{\prime}-3{\hat{y}_{2k}}=f_{y}(\hat{x}_{0k},\hat{y}_{0k},\hat{z}_{0k})\cr{\hat{z}_{2k}}^{\prime\prime}+{\hat{z}_{2k}}=f_{z}(\hat{x}_{0k},\hat{y}_{0k},\hat{z}_{0k})\
,}$ (18)

where ($f_{x},f_{y},f_{z}$) are evaluated along the trajectory
$\hat{\mathbf{r}}_{0k}$, using (15).

The solution can be written as

$\eqalign{\hat{x}_{2k}=c_{1,k}+c_{2,k}\hat{t}+c_{3,k}\hat{t}^{2}+(c_{4,k}+c_{5,k}\hat{t})\sin\hat{t}+(c_{6,k}+c_{7,k}\hat{t})\cos\hat{t}\cr\hat{y}_{2k}=c_{8,k}+c_{9,k}\hat{t}+(c_{10,k}+c_{11,k}\hat{t})\sin\hat{t}+(c_{12,k}+c_{13,k}\hat{t})\cos\hat{t}\cr\hat{z}_{2k}=(c_{14,k}+c_{15,k}\hat{t})\sin\hat{t}+(c_{16,k}+c_{17,k}\hat{t})\cos\hat{t}\
,}$ (19)

where $c_{1,k}\dots c_{17,k}$ are constants that depend on the initial
conditions and on the geometric parameters ($\rho_{0},x_{g},y_{g}$).

The solution contains terms $\propto\hat{t}$ and $\propto\hat{t}^{2}$,
moreover there are also mixed perturbations as $\hat{t}\sin\hat{t}$, i.e.
perturbations increases with time.

In particular, the coefficient of $\hat{t}^{2}$,
$c_{3,k}=-3\,x_{g}/(2\,r_{g})$ is positive number and the same for all
spacecrafts. This means that the entire constellation is ”pushed away” by the
EM system ($r_{g}$ increases with time). The variation of the LISA arms length
can be calculated using (4), (19) and (16) and the result is represented in
Figure 3 (right panel).

The indefinite increasing of the perturbation (19) is not physical because the
perturbative regime would not be valid anymore, after few years. This is a
direct consequence of the force linearization: the terms proportional to
$\hat{t}\sin\hat{t}$ and $\hat{t}\cos\hat{t}$ are first-order terms of the
real, bounded, solution for small $\hat{t}$. Using the F77 code [10], we found
that the perturbative pulls lead to a complete dismembering of the
constellation and a successive recombination will occur after several tens of
thousand year. In this scenario, the distance of each spacecraft from the HCW
frame origin ranges from zero to $2\ AU$. Such motion is not a solution of the
HCW equations, which are valid only if $\rho\ll R_{0}$. An all-time valid
solution to our perturbative problem can be obtained applying the Lindstedt-
Poincaré method [19]. However, as we are interested in the LISA motion during
a few complete orbits, the Lindstedt-Poincaré method is not necessary.

The increase of $r_{g}$ is shown in Figure 3 (left panel), where its time
evolution is represented, during the hypothetical first 10 years of the
mission. The perturbation $\Delta L_{ij}$ to the relative motion between the
pair $i,j$ of LISA test masses due to the monopole perturbation, calculated
both analytically (dashed line) and via numerical integration [10] (solid
line), is plotted on the right panel of Figure 3. It is worth noticing the
good agreement during the 10 years of the LISA mission.

Figure 3: _L_ eft panel: Evolution of $r_{g}$ during 10 years: real (solid
line) and simplified case (Earth describes a circular orbit around the Sun).
The periodic component in the solid line is due to the eccentricity of the
Earth orbit, while the trend is due to the EM system influence. _R_ ight
panel: perturbation of the LISA arms length due to the EM system monopole
effect (rescaled): comparison between the numerical and the analytical
calculation (solid and dashed lines, respectively).

###  4.2 Effects of the EM multipoles

We now search for a perturbative solution to the HCW1 equations in presence of
$n\geq 2$ multipole terms. The intensity of this force is of the order of
$\varepsilon_{2}$. We have already shown that the multipole force is composed
by a periodic and a polynomial component independent of time. The polynomial
component is not important, as it can be added to the EM system monopole
expansion and solved. The solution has the same structure as (19), with
different coefficients (rescaled of a factor
$\varepsilon_{2}/\varepsilon_{0}$), and the motion described in the previous
section is therefore a very good approximation of the EM system polynomial
component influence.

####  4.2.1 $\mathbf{F_{1}}(t)$: periodic solutions in $n\hat{\omega}_{M}t$.

At the zeroth order ($i,j,k=0$) the multipole force does not depend on the
coordinates. Therefore, for each order $n$ we have periodic terms in
$n\hat{\omega}_{M}\hat{t}$. These terms represent the ”pure” oscillations of a
test mass due to the EM system that are not involved with $\omega_{0}$
harmonics. The solutions relative to these frequencies are equal for each test
mass, being independent of its position (we can interpret this as a common
motion). Since we are interested in the relative motion of the LISA
satellites, we know a priori that these terms are subtracted when one measures
the distance between two satellites.

The equations to be solved are

$\eqalign{\hat{x}^{\prime\prime}-2\hat{y}^{\prime}=\sum_{n}\left[a_{xn}\sin
n\hat{\omega}_{M}\hat{t}+b_{xn}\cos
n\hat{\omega}_{M}\hat{t}\,\right]\cr\hat{y}^{\prime\prime}+2\hat{x}^{\prime}-3\hat{y}=\sum_{n}\left[a_{yn}\sin
n\hat{\omega}_{M}\hat{t}+b_{yn}\cos
n\hat{\omega}_{M}\hat{t}\,\right]\cr\hat{z}^{\prime\prime}+\hat{z}=\sum_{n}\left[a_{zn}\sin
n\hat{\omega}_{M}\hat{t}+b_{zn}\cos n\hat{\omega}_{M}\hat{t}\,\right]\ ,}$
(20)

where $a_{xn,yn,zn},b_{xn,yn,zn}$ are constants. Being $n\hat{\omega}_{M}\neq
1$ for each $n$, particular solutions can be written as

$\eqalign{\hat{x}=\sum_{n}a^{\prime}_{xn}\sin
n\hat{\omega}_{M}\hat{t}+b^{\prime}_{xn}\cos
n\hat{\omega}_{M}\hat{t}\cr\hat{y}=\sum_{n}a^{\prime}_{yn}\sin
n\hat{\omega}_{M}\hat{t}+b^{\prime}_{yn}\cos
n\hat{\omega}_{M}\hat{t}\cr\hat{z}=\sum_{n}a^{\prime}_{zn}\sin
n\hat{\omega}_{M}\hat{t}+b^{\prime}_{zn}\cos n\hat{\omega}_{M}\hat{t}}$ (21)

and the corresponding coefficients are

$\displaystyle\eqalign{a^{\prime}_{xn}=-\frac{2b_{yn}n\hat{\omega}_{M}+a_{xn}(3+n\hat{\omega}_{M}^{2})}{n\hat{\omega}_{M}^{2}(n\hat{\omega}_{M}^{2}-1)}&b^{\prime}_{xn}=-\frac{-2a_{yn}n\hat{\omega}_{M}+b_{xn}(3+n\hat{\omega}_{M}^{2})}{n\hat{\omega}_{M}^{2}(n\hat{\omega}_{M}^{2}-1)}\\\
a^{\prime}_{yn}=\frac{2b_{xn}-a_{yn}n\hat{\omega}_{M}}{n\hat{\omega}_{M}(n\hat{\omega}_{M}^{2}-1)}&b^{\prime}_{yn}=\frac{-2a_{xn}-b_{yn}n\hat{\omega}_{M}}{n\hat{\omega}_{M}(n\hat{\omega}_{M}^{2}-1)}\\\
a^{\prime}_{zn}=-\frac{{a_{zn}}}{n\hat{\omega}_{M}^{2}-1}&b^{\prime}_{zn}=-\frac{{b_{zn}}}{n\hat{\omega}_{M}^{2}-1}\
.}$ (25)

Inserting the numerical values, it comes out that the most important
contribution is due to the $2\,\hat{\omega}_{M}$ frequency and it corresponds
to an amplitude of about $1~{}cm$ for $x,y$ coordinates, while the
coefficients $a^{\prime}_{zn},b^{\prime}_{zn}$ are all equal to zero.

In Table 1 we report the coefficients $a^{\prime}_{xn,yn},b^{\prime}_{xn,yn}$
of the fluctuations in meters (i.e. multiplied by $\rho_{0}=2.89\times 10^{9}\
m$), relatively to the first six harmonics of the fundamental frequency
$\hat{\omega}_{M}$.

Table 1:  Rescaled coefficients
$\rho_{0}\,a^{\prime}_{xn,yn},\rho_{0}\,b^{\prime}_{xn,yn}$ of the sinusoidal
terms with frequency $n\hat{\omega}_{M}$ with $1\leq n\leq 6$.

$n$ | Frequency | $\rho_{0}\,a^{\prime}_{xn}$ | $\rho_{0}\,b^{\prime}_{xn}$ | $\rho_{0}\,a^{\prime}_{yn}$ | $\rho_{0}\,b^{\prime}_{yn}$
---|---|---|---|---|---
| [Hz] | [m] | [m] | [m] | [m]
1 | 3.920$\times 10^{-7}$ | -4.5$\times 10^{-5}$ | -1.8$\times 10^{-4}$ | 6.5$\times 10^{-5}$ | -4.4$\times 10^{-5}$
2 | 7.840$\times 10^{-7}$ | 4.4$\times 10^{-3}$ | 8.5$\times 10^{-3}$ | -5.7$\times 10^{-3}$ | 3.9$\times 10^{-3}$
3 | 1.176$\times 10^{-6}$ | -2.1$\times 10^{-5}$ | -2.7$\times 10^{-5}$ | 1.9$\times 10^{-5}$ | -1.8$\times 10^{-5}$
4 | 1.568$\times 10^{-6}$ | 1.1$\times 10^{-7}$ | 1.0$\times 10^{-7}$ | -7.6$\times 10^{-8}$ | 9.8$\times 10^{-8}$
5 | 1.960$\times 10^{-6}$ | -6.4$\times 10^{-10}$ | -3.9$\times 10^{-10}$ | 3.0$\times 10^{-10}$ | -5.7$\times 10^{-10}$
6 | 2.352$\times 10^{-6}$ | 3.8$\times 10^{-12}$ | 1.5$\times 10^{-12}$ | -1.1$\times 10^{-12}$ | 3.4$\times 10^{-12}$

####  4.2.2 $\mathbf{F_{2}}(x,y,z,t)$: periodic solutions in
$(n\hat{\omega}_{M}\pm m)\,\hat{t}$.

In Section 4.2.1 we showed that the $\mathbf{F}_{1}$ term corresponds to the
particular solutions of (20) and that $\mathbf{F}_{1}$ does not affect $\Delta
L_{ij}$ being independent of coordinates.

The motion associated with the coordinate–dependent term $\mathbf{F}_{2}$,
will be different between each pair of test masses, and so the relative
displacements $\Delta L_{ij}$ will be different from zero.

The most direct approach to solve the equation of motion is to write HCW1
equations with $\mathbf{F}_{2}$ evaluated along the unperturbed trajectories
given by (4). Thus, we have obtained only the amplitudes relative to
frequencies $(n\hat{\omega}_{M}\pm 1)$, which represent the main effect of
$\mathbf{F}_{2}$. In order to have the complete spectrum
$(n\hat{\omega}_{M}\pm m)$ we should consider the solutions of HCW with the
complete expansion of the Sun force per unit mass (Section 2.2).

The solution can be written as sum of sinusoidal terms with frequencies
$(n\hat{\omega}_{M}\pm 1)$. In addition, using (16), a similar relation can be
written also for the perturbation $\Delta L_{ij}$ that in SI units reads

$\Delta
L_{ij}(t)=\rho_{0}\sum_{n}\sum_{m=-1}^{1}\left[a_{nm,ij}\sin\omega_{nm}t+b_{nm,ij}\cos\omega_{nm}t\,\right]\
,$ (26)

where $\omega_{nm}=n\omega_{M}+2m\omega_{0}$, with $n=1,2,\dots,+\infty$ and
$m=-1,0,+1$.

The quantity $\delta L_{ijk}$, which can be directly compared with the LISA
sensitivity curve, reads

$\delta
L_{ijk}(t)=\rho_{0}\sum_{n}\sum_{m=-1}^{1}\left[a_{nm,ijk}\sin\omega_{nm}t+b_{nm,ijk}\cos\omega_{nm}t\,\right]\
,$ (27)

where $a_{nm,ijk}=a_{nm,ij}-a_{nm,jk}$ and $b_{nm,ijk}=b_{nm,ij}-b_{nm,jk}$.

The coefficients of (26) and (27) are reported in Table 2 and Table 3,
respectively, relatively to a certain number of frequencies.

In Figure 4 we plot the $\delta L_{123}$ amplitudes (filled circles)
superimposed to two LISA sensitivity curves (straight lines) corresponding to
integration times of 12 days (the upper one), period below which there should
not be disturbances or if present should be removable and 1 year (the lower
one), respectively. The straight lines were obtained by extrapolating the LISA
sensitivity curve down to $10^{-6}\ Hz$, as discussed in [1].

Each amplitude is subject to time variations, because i) the Earth and the
Moon orbits are not circular; ii) the orbital plane of the Moon is slightly
inclined; and iii) the LISA constellation is not rigid.

We estimated uncertainties for $x_{g},y_{g},l$ and $\rho_{0}$ by taking into
account all these effects; in particular, we assumed the error relative to
$\rho_{0}$ equal to the amplitude of the Sun induced breathing, i.e. $\simeq
2$% of $\rho_{0}$ [4].

The resulting relative uncertainties of the $\delta L_{ijk}$ amplitudes are
$\sim 30$%.

Figure 4:  The filled circles give the amplitudes of the differential relative
displacement $\delta L_{123}$ between the LISA test masses $1,2$ and $2,3$
induced by multipoles of the EM system; relative errors are also plotted. The
two straight lines represent the LISA sensitivity curve for 12 days (upper
straight line) and 1 year (lower straight line) of integration time. Table 2:
Coefficients relative to $\Delta L_{ij}(t)$. Amplitudes and frequencies are in
SI units and the overall accuracy on amplitudes is $\simeq 30\%$.

($n,m$) | Frequency | $a_{nm,12}$ | $b_{nm,12}$ | $a_{nm,23}$ | $b_{nm,23}$ | $a_{nm,31}$ | $b_{nm,31}$
---|---|---|---|---|---|---|---
| [Hz] | [m] | [m] | [m] | [m] | [m] | [m]
(1,-1) | 3.286$\times 10^{-7}$ | 5.0$\times 10^{-5}$ | -8.0$\times 10^{-6}$ | -1.8$\times 10^{-5}$ | 4.8$\times 10^{-5}$ | -3.2$\times 10^{-5}$ | -4.0$\times 10^{-5}$
(1, 0) | 3.920$\times 10^{-7}$ | 6.8$\times 10^{-6}$ | 1.6$\times 10^{-5}$ | 6.8$\times 10^{-6}$ | 1.6$\times 10^{-5}$ | 6.8$\times 10^{-6}$ | 1.6$\times 10^{-5}$
(1, 1) | 4.553$\times 10^{-7}$ | 1.9$\times 10^{-5}$ | -1.4$\times 10^{-5}$ | 2.7$\times 10^{-6}$ | 2.4$\times 10^{-5}$ | -2.2$\times 10^{-5}$ | -9.5$\times 10^{-6}$
(2,-1) | 7.206$\times 10^{-7}$ | -2.1$\times 10^{-3}$ | -2.1$\times 10^{-4}$ | 1.2$\times 10^{-3}$ | -1.7$\times 10^{-3}$ | 8.8$\times 10^{-4}$ | 1.9$\times 10^{-3}$
(2, 0) | 7.839$\times 10^{-7}$ | -6.4$\times 10^{-4}$ | -7.9$\times 10^{-4}$ | -6.4$\times 10^{-4}$ | -7.9$\times 10^{-4}$ | -6.4$\times 10^{-4}$ | -7.9$\times 10^{-4}$
(2, 1) | 8.473$\times 10^{-7}$ | -3.9$\times 10^{-4}$ | 5.0$\times 10^{-4}$ | -2.3$\times 10^{-4}$ | -5.9$\times 10^{-4}$ | 6.2$\times 10^{-4}$ | 9.1$\times 10^{-5}$
(3,-1) | 1.113$\times 10^{-6}$ | 8.9$\times 10^{-6}$ | 2.6$\times 10^{-6}$ | -6.7$\times 10^{-6}$ | 6.4$\times 10^{-6}$ | -2.2$\times 10^{-6}$ | -9.0$\times 10^{-6}$
(3, 0) | 1.176$\times 10^{-6}$ | 3.6$\times 10^{-6}$ | 3.0$\times 10^{-6}$ | 3.6$\times 10^{-6}$ | 3.0$\times 10^{-6}$ | 3.6$\times 10^{-6}$ | 3.0$\times 10^{-6}$
(3, 1) | 1.239$\times 10^{-6}$ | 1.1$\times 10^{-6}$ | -2.2$\times 10^{-6}$ | 1.4$\times 10^{-6}$ | 2.0$\times 10^{-6}$ | -2.5$\times 10^{-6}$ | 1.6$\times 10^{-7}$
(4,-1) | 1.505$\times 10^{-6}$ | -4.4$\times 10^{-8}$ | -2.2$\times 10^{-8}$ | 4.2$\times 10^{-8}$ | -2.7$\times 10^{-8}$ | 2.6$\times 10^{-9}$ | 5.0$\times 10^{-8}$
(4, 0) | 1.568$\times 10^{-6}$ | -2.3$\times 10^{-8}$ | -1.3$\times 10^{-8}$ | -2.3$\times 10^{-8}$ | -1.3$\times 10^{-8}$ | -2.3$\times 10^{-8}$ | -1.3$\times 10^{-8}$
(4, 1) | 1.631$\times 10^{-6}$ | -3.0$\times 10^{-9}$ | 1.2$\times 10^{-8}$ | -8.5$\times 10^{-9}$ | -8.4$\times 10^{-9}$ | 1.2$\times 10^{-8}$ | -3.2$\times 10^{-9}$
(5,-1) | 1.896$\times 10^{-6}$ | 2.3$\times 10^{-10}$ | 1.8$\times 10^{-10}$ | -2.7$\times 10^{-10}$ | 1.1$\times 10^{-10}$ | 3.7$\times 10^{-11}$ | -2.9$\times 10^{-10}$
(5, 0) | 1.960$\times 10^{-6}$ | 1.5$\times 10^{-10}$ | 5.4$\times 10^{-11}$ | 1.5$\times 10^{-10}$ | 5.4$\times 10^{-11}$ | 1.5$\times 10^{-10}$ | 5.4$\times 10^{-11}$
(5, 1) | 2.023$\times 10^{-6}$ | 3.9$\times 10^{-12}$ | -6.5$\times 10^{-11}$ | 5.4$\times 10^{-11}$ | 3.6$\times 10^{-11}$ | -5.8$\times 10^{-11}$ | 2.9$\times 10^{-11}$
(6,-1) | 2.288$\times 10^{-6}$ | -1.2$\times 10^{-12}$ | -1.3$\times 10^{-12}$ | 1.7$\times 10^{-12}$ | -3.8$\times 10^{-13}$ | -5.4$\times 10^{-13}$ | 1.7$\times 10^{-12}$
(6, 0) | 2.320$\times 10^{-6}$ | -9.9$\times 10^{-13}$ | -1.7$\times 10^{-13}$ | -9.9$\times 10^{-13}$ | -1.7$\times 10^{-13}$ | -9.9$\times 10^{-13}$ | -1.7$\times 10^{-13}$
(6, 1) | 2.415$\times 10^{-6}$ | 5.0$\times 10^{-14}$ | 3.7$\times 10^{-13}$ | -3.5$\times 10^{-13}$ | -1.4$\times 10^{-13}$ | 3.0$\times 10^{-13}$ | -2.3$\times 10^{-13}$

Table 3:  Coefficients relative to $\delta L_{ijk}(t)$. Amplitudes and
frequencies are in SI units and the overall accuracy on amplitudes is $\simeq
30\%$.

($n,m$) | Frequency | $a_{nm,123}$ | $b_{nm,123}$ | $a_{nm,231}$ | $b_{nm,231}$ | $a_{nm,312}$ | $b_{nm,312}$
---|---|---|---|---|---|---|---
| [Hz] | [m ] | [m ] | [m ] | [m ] | [m ] | [m ]
(1,-1) | 3.286$\times 10^{-7}$ | 6.9 $\times 10^{-5}$ | -5.6 $\times 10^{-5}$ | 1.4 $\times 10^{-5}$ | 8.7 $\times 10^{-5}$ | -8.3 $\times 10^{-5}$ | -3.2 $\times 10^{-5}$
(1, 0) | 3.920$\times 10^{-7}$ | 0 | 0 | 0 | 0 | 0 | 0
(1, 1) | 4.553$\times 10^{-7}$ | 1.6 $\times 10^{-5}$ | -3.8$\times 10^{-5}$ | 2.5 $\times 10^{-5}$ | 3.3$\times 10^{-5}$ | -4.1 $\times 10^{-5}$ | 4.7$\times 10^{-6}$
(2,-1) | 7.206$\times 10^{-7}$ | -3.4 $\times 10^{-3}$ | 1.5 $\times 10^{-3}$ | 3.6 $\times 10^{-4}$ | -3.7 $\times 10^{-3}$ | 3.0$\times 10^{-3}$ | 2.1 $\times 10^{-3}$
(2, 0) | 7.839$\times 10^{-7}$ | 0 | 0 | 0 | 0 | 0 | 0
(2, 1) | 8.473$\times 10^{-7}$ | -1.6$\times 10^{-4}$ | 1.1 $\times 10^{-3}$ | -8.6 $\times 10^{-4}$ | -6.8 $\times 10^{-4}$ | 1.0$\times 10^{-4}$ | -4.0$\times 10^{-4}$
(3,-1) | 1.113$\times 10^{-6}$ | 1.6 $\times 10^{-5}$ | -3.8$\times 10^{-6}$ | -4.6$\times 10^{-6}$ | 1.5 $\times 10^{-5}$ | -1.1 $\times 10^{-5}$ | -1.2 $\times 10^{-5}$
(3, 0) | 1.176$\times 10^{-6}$ | 0 | 0 | 0 | 0 | 0 | 0
(3, 1) | 1.239$\times 10^{-6}$ | -2.8$\times 10^{-7}$ | -4.3$\times 10^{-6}$ | 3.8$\times 10^{-6}$ | 1.9$\times 10^{-6}$ | -3.6$\times 10^{-6}$ | 2.4$\times 10^{-6}$
(4,-1) | 1.505$\times 10^{-6}$ | -8.6$\times 10^{-8}$ | 4.5$\times 10^{-9}$ | 3.9$\times 10^{-8}$ | -7.7$\times 10^{-8}$ | 4.7$\times 10^{-8}$ | 7.2$\times 10^{-8}$
(4, 0) | 1.568$\times 10^{-6}$ | 0 | 0 | 0 | 0 | 0 | 0
(4, 1) | 1.631$\times 10^{-6}$ | 5.5$\times 10^{-9}$ | 2.0$\times 10^{-8}$ | -2.0$\times 10^{-8}$ | -5.3$\times 10^{-9}$ | 1.5$\times 10^{-8}$ | -1.5$\times 10^{-8}$
(5,-1) | 1.896$\times 10^{-6}$ | 5.0$\times 10^{-10}$ | 6.5$\times 10^{-11}$ | -3.0$\times 10^{-10}$ | 4.0$\times 10^{-10}$ | -1.9$\times 10^{-10}$ | -4.6$\times 10^{-10}$
(5, 0) | 1.960$\times 10^{-6}$ | 0 | 0 | 0 | 0 | 0 | 0
(5, 1) | 2.023$\times 10^{-6}$ | -5.0$\times 10^{-11}$ | -1.0$\times 10^{-10}$ | 1.1$\times 10^{-10}$ | 6.8$\times 10^{-12}$ | -6.2$\times 10^{-11}$ | 9.4$\times 10^{-11}$
(6,-1) | 2.288$\times 10^{-6}$ | -3.0$\times 10^{-12}$ | -9.4$\times 10^{-13}$ | 2.3$\times 10^{-12}$ | -2.1$\times 10^{-12}$ | 6.6$\times 10^{-13}$ | 3.0$\times 10^{-12}$
(6, 0) | 2.352$\times 10^{-6}$ | 0 | 0 | 0 | 0 | 0 | 0
(6, 1) | 2.415$\times 10^{-6}$ | 4.0$\times 10^{-13}$ | 5.2$\times 10^{-13}$ | -6.5$\times 10^{-13}$ | 8.6$\times 10^{-14}$ | 2.5$\times 10^{-13}$ | -6.0$\times 10^{-13}$

##  5 Effect of Venus and Jupiter on LISA motion

To make a comparison with the effect on the LISA motion due to EM system
effect, we have numerically evaluated the effect of Venus, Jupiter and the EM
system (its monopole contribution) using a F77 code implementing the algorithm
described in [10]. In practice, we independently computed the LISA motion
under the effect of Sun+Venus, Sun+Jupiter and Sun+EM system and we
subtracted, to each of them, the unperturbed motion due to the Sun. In Figure
5 we plot the $\delta L_{123}$ perturbations (in km) once the modulation due
to the Sun has been subtracted. We have found that the monopole contribution
of the EM system (solid line) is much larger than that of Venus (dotted line)
and Jupiter (dashed line). Figure 5 also shows the onset of the resonant
effect of the EM system monopole after 2-3 years.

Figure 5: $\delta L_{123}$ perturbation (in km) due to EM system (solid line),
Venus (dotted line) and Jupiter (dashed line). $\delta L_{123}$ perturbation
(in km) due to monopole contribution of the EM system (solid line), Venus
(dotted line) and Jupiter (dashed line). See Section 5.

##  6 Conclusions

We calculated the relative motion of LISA test masses due to the effect of the
EM system monopole and we found $\Delta L_{ij}\simeq 3\times 10^{5}\ km$ after
a period of 10 years. We also found that the Jupiter and Venus influences are
at least 10 times smaller than the EM system one.

The perturbations of the differential relative motion of LISA test masses
$\delta L_{ijk}$ are in the $10^{-6}$ to $10^{-7}\ Hz$ decade. This is a very
low frequency range in which LISA residual acceleration noise may be much
larger than extrapolated on grounds of known effects in [1]. However, it might
not be completely hopeless to get an interesting sensitivity also at such low
frequencies, provided that one will have to face the problem of i) loss of
signal coherence over a time scale of one month; and ii) low frequency range
calibration.

As discussed in Pollack [20] it is possible to extract a signal from the LISA
data, even in presence of disturbances. These latter arise due to
environmental effects, such as cosmic rays induced by solar flares, and the
telecommunication antenna which periodically has to be rotated. Pollack showed
how these disturbances can be identified and subsequently removed from the
data even at low frequencies. Assuming, for instance, that a disturbance
appears every 19 days (see Table 6 in [20]) the resulting error in the signal
frequency of $3\times 10^{-6}\,Hz$ is only $\sim 1.3\times 10^{-9}\,Hz$. Thus
it should not be a problem to extend to such low frequencies the calibration
from the verification binaries, and, by using their signals, ensure the
continuity of data over time spans of many weeks. Still the signal from the EM
system, as understood at the level of accuracy given in the present paper, can
give a relevant additional crosscheck to such an extension of the calibration.
It thus may help in improving our knowledge of the LISA acceleration noise at
very low frequencies and contribute to extend to such low frequencies the
capabilities of LISA.

We are indebted to Peter Bender for a critical reading of the initial version
of the manuscript, together with helpful suggestions. We thank Oliver Jennrich
and Gerard Gómez, for useful discussions. Mauro Sereno was supported in the
early stages of this work by the Swiss National Science Foundation.

## References

## References

  * [1] Bender P L 2003 Class. Quantum Grav.20 301-310
  * [2] LISA Pre-Phase A Report, LISA Project internal report number Max Planck Institut für Quantenoptik 233 (July 1998)
  * [3] LISA: System and Technology Study Report ESA document ESA-SCI 2000 11, July 2000
  * [4] Dhurandhar S V, Rajesh Nayak K, Koshti S, Vinet K 2005 Class. Quantum Grav.22 481-487
  * [5] Dhurandhar S V, Vinet J Y; Rajesh Nayak K 2008 Class. Quantum Grav.25 245002
  * [6] Cerdonio M, De Marchi F, De Pietri R, Jetzer P, Marzari F, Mazzolo G, Ortolan A and Sereno M 2010 (arXiv: gr-qc/1002.0489v1)
  * [7] Cerdonio M, De Pietri R, Jetzer P, Sereno M 2009 Class. Quantum Grav.26 094022
  * [8] Povoleri A, Kemble S 2006 in Laser Interferometer Space Antenna, AIP Conf. Proc. 873, (Amer. Inst. of Physics, Melville, NY) pp 702-706
  * [9] Dhurandhar S V 2009 J. Phys. Conf. Ser. 154 012047
  * [10] Everhart E 1985, ”An efficient integrator that uses Gauss?-Radau spacings”, in Dynamics of Comets: Their Origin and Evolution, A. Carusi and G. B. Valsecchi, Eds., Dordrecht, Reidel (1985), pp. 185-?202.
  * [11] Dixon R T 1971 Dynamic astronomy (Prentice-Hall, New Jersey)
  * [12] Clohessy W H and Wiltshire R S 1960 J. Aerospace Sci. 653-58
  * [13] Shirley J H & Fairbridge R W 1997 Encyclopedia of planetary sciences, (Springer Us)
  * [14] Bocaletti D and Pucacco G 2001 Theory of orbits (Springer, Berlin)
  * [15] Bakulin P I, Kononovich E V, Moroz V I 1984 Course in General Astronomy
  * [16] Nayak K R, Koshti S, Dhurandhar S V and Vinet J-Y 2006 Class. Quantum Grav.23 1763-78
  * [17] Abhyankar K D 1999 Astrophysics of the solar system (Sangam Books Ltd)
  * [18] Murray C D & Dermott S F 1999 Solar System Dynamics (Cambridge University Press p 184)
  * [19] Gómez G, Marcote M AAS/AIAA Astrodynamics Specialists Conference Aug, 7-11, 2005 Paper AAS 05-359
  * [20] Pollack S E, 2004 Class. Quantum Grav. 21 3419-3432


# Quarkonium spin alignment in a vortical medium

Paulo Henrique De Moura1,Kayman J. Gonçalves1, Giorgio Torrieri1,2
$\phantom{A}{}^{1}$ Universidade Estadual de Campinas - Instituto de Fisica
”Gleb Wataghin”
Rua Sergio Buarque de Holanda, 777
CEP 13083-859 - Campinas SP
$\phantom{A}{}^{2}$ Institute of Physics, Jan Kochanowski University, Ul.
Uniwersytecka 7, 25-406 Kielce, Poland.

###### Abstract

We use a potential model to investigate the phenomenology of quarkonium in a
thermal rotating medium, where vorticity and spin density are not necessarily
in equilibrium. We find that the quarkonium spin density matrix, as well as
the binding energy and melting temperature, are sensitive to both the
vorticity and the lack of equilibrium between vorticity and spin. This means
that quarkonium spin alignment is a sensitive probe for vorticity and spin
within the hydrodynamic phase. Information unequivocally pointing to spin-
orbit non-equilibrium dynamics can be obtained from a combined study of
quarkonium relative abundance and spin alignment, as well as experimentally
obtainable off-diagonal density matrix elements.

##  I Introduction

Quarkonium as been used as a probe of thermodynamic properties of the Quark-
Gluon Plasma (QGP) since the seminal paper from MATSUI_SATZ . The heavy mass
scale means that it is a probe that can both be examined reliably by
theoretical calculations and susceptible to non-trivial in-medium effects,
serving as a “thermometer”, sensitive to the interplay of thermal fluctuations
and the QCD medium, and as a probe for non-equilibrium behavior review ;
classic1 ; classic2 ; classic3 .

The advent of the study of vorticity in heavy ion collisions bec  added a
potentially new arena where quarkonium could be used. In fact, we shall argue
that quarkonium provides unique opportunities for the phenomenology of
hydrodynamics with spin. First of all, quarkonium can be formed early in the
collision and can survive throghout the quark gluon plasma evolution. Unlike
polarized $\Lambda$s bec  and spin-aligned vector mesons alices ; star ; wang
; xia ; kayman ; oliva , they are potentially sensitive to the entire dynamics
of the fluid and not just to the freeze-out.

Furthermore, the long-lived quarkonium state is spin 1, having a 3$\times$3
density matrix with 8 degrees of freedom, 6 of which are accessible to spin
alignment measurements. Thus, if, as seems to be theoretically highly likely
uscausality ; uscas2 ; usdiss ; relax ; jeon ; gursoy ; hongo ; hongo2 , spin
and vorticity are not in equilibrium, this lack of equilibration will can be
imprinted on the density matrix’s measurable off-diagonal elements
$\rho_{0,\pm 1},\rho_{\pm 1,\mp 1}$ kayman  (No equivalent elements exist for
the “qubit” density matrix of a fermion). In fact, these elements have
recently been measured ALICE .

Last but not least, bottomonium and to a certain extent charmonium states can
be viewed as solutions to a Schrodinger equation with heavy Quark
wavefunctions moving around a QCD potential (including a weakly coupled and a
confined part) brambilla ; boundquark . It is therefore possible to
understand, both analytically tuchin  and via effective theory kim , effect
that vorticity will have to the properties of quarkonium.

In this work, we shall go in this direction, combining the insights developed
in kayman  with a potential model for quarkonium solved using standard methods
extended to rotating frames NU ; ARF ; Anan ; Anan2  and finite temperature
AHM ; WONG . In section II we shall assess currently available experimental
data, and point out what would be necessary to probe spin-orbit non-
equilibrium. Then in section III we calculate the quarkonium wavefunction
properties, namely binding energy and a semi-classical estimate of the melting
temperature, in a rotating frame. Finally, in section IV we calculate
quarkonium observables that could indicate a lack of equilibrium between spin
and vorticity.

##  II The Quarkonium spin density matrix

In kayman , we have argued that vector spin alignment contains crucial
information on the still-unknown spin hydrodynamic evolution in heavy ion
collisions, provided that not just the $\rho_{00}$ coefficient but the “off-
diagonal” coefficients are measured. The former, $\rho_{00}$ was measured in
alices ; star  and is sensitive to $\theta$, the angle w.r.t. the spin
alignment direction, which for heavy ions would be mainly the reaction plane,
reflecting the vorticity structure in analogy to the global polarization
measurement bec .

The other coefficients, related to density matrix elements $\rho_{0\pm 1,\pm
1\pm 1,\pm 1\mp 1}$ (see equation (2) of kayman ) would depend on a second
“reference” angle $\phi$, whose most logical definition is in terms of the
beam axis.

Figure 1:  A schematic representation of how the vorticity in heavy ion
collisions looks like at early and late times Blue dashed arrows represent
vorticity, cyan dot-dashed ones polarization, red solid ones flow. The mis-
alignment of spin in the right panel is due to the lack of equilibrium between
vorticity and polarization uscausality ; usdiss .

In such a set-up, the non-equilibrium between spin and vorticity would
manifest itself by the interplay of transverse vs. longitudinal polarization
bec . The heavy ion system actually lends hope that this interplay could occur
(Fig. 1 ): Transverse polarization is thought to be present from the
beginning, since it is the result of the angular momentum of the initial state
of an off-central collision. Longitudinal polarization forms from azimuthal
gradients, on a time-scale comparable to the formation of momentum anisotropy
causing elliptic flow. If polarization “relaxes to vorticity” uscausality ;
usdiss , one expects polarization to be aligned in the transverse direction,
but comparatively mis-aligned in the longitudinal one usdiss , which indeed
seems to be out of phase with some dynamical calculations bec . Hence, there
will be an angle in the $z-y$ plane between vorticity and polarization, whose
magnitude we do not yet know to estimate quantitatively but whose experimental
confirmation will add a phenomenological dimension to the recent theoretical
development of spin hydrodynamics uscausality ; uscas2 ; usdiss ; relax ; jeon
; gursoy ; hongo ; hongo2 .

The hope, therefore, is that off-diagonal (longitudinal-transverse)
coefficients are potentially important for they could signal deviations from
local equilibrium (included via a Cooper-Frye type formula zanna ) due to the
presence of two distinct axial currents, representing spin and vorticity,
evolving on different time-scales usdiss . In kayman  we illustrated this with
a coalescence type model.

Coalescence of only spin within a vortical background should not change the
coherence of the density matrix, since it is a unitary process and the
dynamics is symmetric around the vortical axis. But assuming vorticity and
pre-existing spin density are not in equilibrium and pointing in different
directions, this is no longer true kayman ; Vorticity is “classical”
background, interacting with the quantum spin state, so if coalescence happens
in a vortical background (i.e. if spin and vorticity are out of equilibrium)
one expects impurity of the density matrix. Mathematically, the loss of purity
is manifest in Eq. (23) of kayman  $P_{L}(w)$ representing the (unknown)
classical probability of a vortex $\omega$ giving an angular momentum $L$ to
the meson wavefunction. When this probability becomes uniform
($P_{L}(w)\rightarrow$ constant ) we recover a maximally impure
state111Through not quite the Cooper-frye ansatz of zanna . A $P_{L}=$
constant impure state can be regarded as a microcanonical density matrix
assuming the diquark quarkonium state is exact. The grand canonical matrix
inherent in the Cooper-Frye formula of zanna  and it’s vector extension would
arise if all values of $L$ up to $\infty$ were allowed due to angular momentum
fluctuations and a bath of degrees of freedom.. For vector bosons, this
impurity is manifest in the off-diagonal matrix components (see the discussion
between eqs 4 and 5 of kayman ).

While this data as yet does not exist for vector mesons, it does exist for
quarkonia $J/\psi$ and $\Upsilon$ states ALICE , since the ALICE Collaboration
measurement of the quarkonium polarization included the off-diagonal values of
the spin density matrix. Therefore, we can do direct connection between
polarization parameters $\lambda_{\theta}$, $\lambda_{\phi}$ and
$\lambda_{\theta\phi}$ and density matrix ALICE  and the coefficients used in
kayman  (table in eq. 1)

$\begin{array}[]{ccc}\mathrm{Variable}&\mathrm{Element}&\mathrm{coefficient}\times\frac{3}{4\pi}\\\
\rho_{00}&\rho_{00}&\cos^{2}\theta\\\
\frac{1-\rho_{00}}{2}&\frac{\rho_{11}+\rho_{-1-1}}{2}&\sin^{2}\theta\\\
r_{10}&Re[\rho_{-10}-\rho_{10}]&\sin(2\theta)\cos(\phi)\\\
\alpha_{10}&Im[-\rho_{-10}+\rho_{10}]&\sin(2\theta)\sin(\phi)\\\
r_{1,-1}&Re[\rho_{1,-1}]&\sin^{2}\theta\cos(2\phi)\\\
\alpha_{1,-1}&Im[\rho_{1,-1}]&\sin^{2}\theta\sin(2\phi)\end{array}$ (1)
$\rho_{00}=\frac{1+\lambda_{\theta}}{3+\lambda_{\theta}}\phantom{AA},\phantom{AA}r_{1,-1}=\frac{\lambda_{\phi}}{3+\lambda_{\theta}}\phantom{AA},\phantom{AA}r_{10}=\frac{\lambda_{\theta\phi}}{3+\lambda_{\theta}}$
(2)

This is possible by making the comparison between $J/\psi$ angular
distribution and the standard vector mesons angular distribution shown in Eq.
1 of kayman  (Eq. 2).

Thus, we can do the analyses presented in kayman  to relate
$\lambda_{\theta,\phi,\theta\phi}$ to the wave function coherence via the
parametrization in terms of Gell-Mann matrices. Choosing the $n_{3-8}$ basis
for this parametrization, we need to solve the following system of algebraic
equations derived in kayman  in terms of the frame relating the lab to the
spin direction (defined by angles $\theta_{r},\phi_{r}$)

$\displaystyle\frac{1}{12}\left(3\left(n_{8}-\sqrt{3}\;n_{3}\right)\cos\left(2\theta_{r}\right)-\sqrt{3}\;n_{3}+n_{8}+4\right)=\rho_{00}$
(3)
$\displaystyle\frac{\left(n_{8}-\sqrt{3}\;n_{3}\right)\sin\left(\theta_{r}\right)\cos\left(\theta_{r}\right)\cos\left(\phi_{r}\right)}{\sqrt{2}}=r_{10}$
(4)
$\displaystyle-\frac{\left(\sqrt{3}\;n_{3}+3n_{8}\right)\sin\left(\theta_{r}\right)\sin\left(\phi_{r}\right)}{3\sqrt{2}}=\alpha_{10}$
(5)
$\phi_{r}=-\frac{1}{2}\tan^{-1}\left(\frac{\alpha_{1,-1}}{r_{1,-1}}\right)$
(6)

Now, we will do the follow change variable $\tilde{n}=n_{8}-\sqrt{3}\;n_{3}$,
and knowing that variables are equal to zero $\alpha_{10}$ and
$\alpha_{1,-1}$. So, we can write this system equation the following form:

$\displaystyle\frac{1}{12}\left(3\tilde{n}\cos\left(2\theta_{r}\right)+\tilde{n}+4\right)=\rho_{00}$
(7)
$\displaystyle\frac{\tilde{n}\sin\left(\theta_{r}\right)\cos\left(\theta_{r}\right)\cos\left(\phi_{r}\right)}{\sqrt{2}}=r_{10}$
(8) $\displaystyle\phi_{r}=0$ (9)

Therefore, we have the following solution:

$\displaystyle\tilde{n}\left(\rho_{00},r_{10}\right)=-\frac{(1-3\rho_{00})^{2}+3\sqrt{(1-3\rho_{00})^{4}+4(1-3\rho_{00})^{2}r_{10}^{2}}}{6\rho_{00}-2}$
(11)
$\displaystyle\Theta_{1}\left(\rho_{00},r_{10}\right)=-\sqrt{\frac{2(1-3\rho_{00})^{2}-2\sqrt{(1-3\rho_{00})^{4}+4(1-3\rho_{00})^{2}r_{10}^{2}}+6r_{10}^{2}}{2(1-3\rho_{00})^{2}+9r_{10}^{2}}}$
(12)
$\displaystyle\Theta_{2}\left(\rho_{00},r_{10}\right)=\sqrt{\frac{(1-3\rho_{00})^{2}-\sqrt{(1-3\rho_{00})^{4}+4(1-3\rho_{00})^{2}r_{10}^{2}}+3r_{10}^{2}}{2(1-3\rho_{00})^{2}+9r_{10}^{2}}}$
(13)
$\displaystyle\Theta_{3}\left(\rho_{00},r_{10}\right)=\frac{\Theta_{2}\left(\rho_{00},r_{10}\right)\left(2(1-3\rho_{00})^{2}+2\sqrt{(1-3\rho_{00})^{4}+4(1-3\rho_{00})^{2}r_{10}^{2}}\right)}{2(3\rho_{00}-1)r_{10}}$
(14)
$\displaystyle\theta_{r}=\tan^{-1}\left(\Theta_{1}\left(\rho_{00},r_{10}\right),\Theta_{3}\left(\rho_{00},r_{10}\right)\right)$
(15)

So, using polarization parameters that were obtained from ALICE collaboration
ALICE  at different transverse momentum ($p_{T}$) ranges, we can determine
whether the density matrix represents a coherence state or not. To make it we
will use the equation 11 and reach the following figure 2.

Figure 2: The experimental results of $J/\psi$ and $\Upsilon$ polarization
measurements analyzed in terms of the Gell-Mann matrix representation of
kayman . These results were made using the helicity and Collins-Soper frame to
$J/\psi$ polarization ALICE . The other Collins-Soper values was not take in
this plot because they were above the Helicity $6<p_{T}<10~{}(GeV/C^{2})$
uncertainty band as were shown in the text.

Now, the coefficients in the frame Collins-Soper frame, given in the ranges
$4<p_{T}<6~{}(GeV/C^{2})$ and $6<p_{T}<10~{}(GeV/C^{2})$ respectively result
in $\tilde{n}=(0.09\pm 0.10,0.02\pm 0.07)$. Making the comparison with the
Helicity frame value $\tilde{n}=0.09\pm 0.11$ to $6<p_{T}<10~{}(GeV/C^{2})$.
Therefore, we can see that they are the same within error bar.

Looking at the figure 2, we can conclude that the density matrix from $J/\psi$
particle does not represent the pure state since none of the values for
$n_{3,8}$ obtained from the data intersects the black pink points, i.e. these
points represent the pure state in other words when the density matrix satisfy
$\rho^{2}=\rho$. This might indicate that statistical freeze-out advocated in
zanna ; classic1  is a good estimate of particle production in heavy-ion
collisions because the density matrix does not represent a coherent state as
argued in kayman . In the bottomonium $\Upsilon(1S)$ case, we can see in the
right panel of Fig. 2 that because of large uncertainty do not know whether
for these particles the density matrix represents a coherence state or not.

To fully assess the significance of the estimate above one must recall Fig. 4
of kayman  and the definition of $P_{L}(w)$ given in the introduction. We have
no idea what $P_{L}(w)$ is beyond the fact that it overall conserves momentum,
but it acts as a projector. One recovers a pure state when
$\rho_{ij}=\delta_{Ll}$, (there is a certainty of vorticity giving a certain
momentum) and a maximally mixed state when the momentum given by vorticity is
independent of $L$. So the measurement in Fig. 4 kayman  is directly connected
to how out of equilibrium vorticity and spin are, and how much vorticity vs
pre-existing spin influences the final spin of the vector meson. Linear
combination of the different $L$-values in Fig. 4 of kayman  are possible,
illustrating a probability of different spin configurations.

Note that these coefficients are given in terms of an angle $\theta$, which in
kayman  is related to $\theta_{r}$, the angle between the hadronization frame
and the lab frame. This angle of course depends on the detailed hydrodynamical
and spin-hydrodynamical evolution of the system, but it is obviously highly
dependent on the reaction plane angle $\Phi$. Note that since local
polarization necessarily averages to zero in symmetric collisions by parity,
so must the off-diagonal matrix elements (unless large event by event
fluctuations in transparency occur, something never seen in experimental
data). Because of this, a zero finding integrated over azimuthal anisotropy
can not count as conclusive evidence.

Considering the Harmonic behavior of the coefficients in Fig. 4 of kayman
w.r.t. $\theta_{r}$ (most coefficients average to zero for all angles),
therefore, it would be crucial to measure $\lambda_{\theta,\phi,\theta\phi}$
not as a function of $p_{T}$ as in ALICE  but as function of azimuthal
reaction plane angle (This can be measured directly at lower energies event1 ,
and via cumulants at the LHC event2 ). A modulated behavior would be a clear
signature for a non-trivial $P_{L}(w)$ which can then be harmonically
decomposed into $L=0,1,2$ components of Fig 4 of kayman  to obtain information
of the impact of spin vs vorticity in $J/\psi$ vs $\Upsilon$ hadronization. If
the dependence w.r.t. $\theta_{r}$ will be compatible with zero as it was for
$p_{T}$ in each $\Phi$ bin, this is good evidence for a statistical Cooper-
Frye freeze-out as in zanna ; classic1 . Schematically, these two alternatives
are illustrated in Fig. 3 .

The quantitative details of Fig. 3 would require a hydrodynamic simulation
with a hydrodynamic model where spin and vorticity are not in equilibrium,
given the interplay between spin, vorticity,radial flow and its anisotropies
in the quarkonia $p_{T}$ and rapidity distributions. We have used an
analytically solvable model subject to ongoing work blast1 ; blast2  to
provide an order of magnitude estimate of the maximum of the effect compatible
with Fig. 2. For reasons expanded on in the appendix this should be regarded
as just that, an order of magnitude estimate, although it is gratifying that
the possible off-diagonal terms are comparable in magnitude to
$\lambda_{\theta}$.

For this reason, in the rest of the paper we proceed to examine the
microscopic dynamics of the charmonium state using a potential model, with a
view of developing quantitative signatures of non-equilibrium between
charmonium spin and vorticity.

Figure 3:  A schematic illustration of how we expect the coefficients $\lambda_{\theta,\phi,\theta\phi}$ to evolve with reaction plane angle $\Phi-\Phi_{reaction}$ in the two scenarios. In ALICE , these coefficients are compatible with zero when integrated over this angle. A flat dependence (a) indicates a maximally incoherent Cooper-Frye/thermal type production. A harmonic dependence (b) might indicate a more non-trivial coalescence scenario involving vorticity. Note that both cases are consistent with zero when integrated over $\Phi$, as was done in ALICE . NB: The exact azimuthal shape will need more realistic modeling, but the amplitude can be qualitatively estimated in Table 1 using blast1 ; blast2  | $\lambda_{\theta}$ | $\lambda_{\phi}$ | $\lambda_{\theta\phi}$
---|---|---|---
$J/\psi$ | 0.021 | 0.023 | 0.012
$\Upsilon$ | 0.054 | 0.054 | 0.002
Table 1:  These are estimates of density matrix coefficients using Blast Wave
Model blast1 ; blast2  more details in Appendix Appendix B .

##  III The quarkonium state in rotating reference frames

###  III.1 The Schrödinger equation

We are interested in the angular momentum due to vortices couples with the
quark spin. So, we need to take an extra term,
$\boldsymbol{\omega}\cdot\mathbf{S}$ Anan ,

$\mathcal{H}=\frac{(\boldsymbol{p}-m\boldsymbol{\omega}\times\mathbf{r})^{2}}{2m}-\frac{m}{2}(\boldsymbol{\omega}\times\mathbf{r})^{2}-\boldsymbol{\omega}\cdot\mathbf{S}+V(r)$
(16)

Now, we can write the equation 16 for a two-body case. Thus we have the
following expression:

$\mathcal{H}=\sum_{i=1,2}\left(\frac{(\boldsymbol{p}_{i}-m_{i}\boldsymbol{\omega}_{i}\times\mathbf{r}_{i})^{2}}{2m_{i}}-\frac{m_{i}}{2}(\boldsymbol{\omega_{i}}\times\mathbf{r}_{i})^{2}-\boldsymbol{\omega_{i}}\cdot\mathbf{S_{i}}\right)+V(\lvert\mathbf{r_{1}}-\mathbf{r_{2}}\rvert)$
(17)

Using the relations:

$\left\\{\begin{array}[]{lcl}\mathbf{P}=\mathbf{p_{1}}+\mathbf{p_{2}},&&\mathbf{p}=\mu\left(\dfrac{\mathbf{p_{1}}}{m_{1}}-\dfrac{\mathbf{p_{2}}}{m_{2}}\right)\\\\[8.53581pt]
\mu=\dfrac{m_{1}\;m_{2}}{m_{1}+m_{2}},&&\\\\[8.53581pt]
\mathbf{r_{1}}=\mathbf{R}+\dfrac{m_{2}}{m_{1}+m_{2}}\mathbf{r},&&\mathbf{r_{2}}=\mathbf{R}-\dfrac{m_{1}}{m_{1}+m_{2}}\mathbf{r}\end{array}\right.$
(18)

As the two quarks that form the mesons are in the same vortical background, we
can suppose that
$\boldsymbol{\omega}_{1}=\boldsymbol{\omega}_{2}=\boldsymbol{\omega}$ and
$\mathbf{S}=\mathbf{S_{1}}+\mathbf{S_{2}}$, and then write:

$\mathcal{H}=\frac{P^{2}}{2M}+\frac{p^{2}}{2\mu}-\mathbf{P}\cdot(\boldsymbol{\omega}\times\mathbf{R})-\mathbf{p}\cdot(\boldsymbol{\omega}\times\mathbf{r})-\boldsymbol{\omega}\cdot(\mathbf{S}_{1}+\mathbf{S}_{2})+V(\lvert\mathbf{r}\rvert)$
(19)

Since $\mathbf{S}=\mathbf{S}_{1}+\mathbf{S}_{2}$ and we are interested just in
the reduced coordinate,

$\mathcal{H}=\frac{p^{2}}{2\mu}-\mathbf{p}\cdot(\boldsymbol{\omega}\times\mathbf{r})-\boldsymbol{\omega}\cdot\mathbf{S}+V(\lvert\mathbf{r}\rvert)$
(20)

We can rewrite this equation as

$i\frac{\partial}{\partial t}=-\frac{1}{2\mu}\boldsymbol{\nabla}^{2}-\omega
L_{z}-\omega S_{z}+V(\lvert\mathbf{r}\rvert)$ (21)

where we used $\mathbf{p}=-i\hbar\boldsymbol{\nabla}$,
$\mathcal{H}=i\frac{\partial}{\partial t}$,
$\boldsymbol{\omega}=\omega\mathbf{\hat{z}}$ and
$L_{z}=-i\frac{\partial}{\partial\phi}$.

In rotating frames, the contribution in Hamiltonian is only the product
between orbital angular momentum due to the vortices with the spin of meson.
Thus, in this non-inertial frame, the contribution is just to $L$ different
from zero. Then, we will write just the radial part of the equation 21.

$\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{dR(r)}{dr}\right)+\left[2\mu(E-V(r))-\frac{l(l+1)}{r^{2}}-2\mu\omega
L_{z}-2\mu\omega S_{z}\right]R(r)=0$ (22)

To go forward, we shall assume the rotation to be classical and related to the
hydrodynamic vorticity. Then we can define $\omega$ in terms of the a
conserved circulation

$C=\oint\mathbf{v}\cdot d\mathbf{l}$ (23)

and also assume the Cornell potential brambilla

$V(r)=b\;r-\frac{\alpha_{eff}}{r}$ (24)

so

$\frac{1}{2\mu
r^{2}}\frac{d}{dr}\left(r^{2}\frac{dR(r)}{dr}\right)+\left[(E-V(r))-\frac{l(l+1)}{2\mu
r^{2}}-\frac{\;m_{j}\;C}{\pi r^{2}}\right]R(r)=0$ (25)

Making $\chi(r)=rR(r)$ and change the variable $x=1/r$, then the equation 25,
we get:

$\frac{d^{2}\chi(x)}{dx^{2}}+\frac{2x}{x^{2}}\frac{d\chi(x)}{dx}+\frac{2\mu}{x^{4}}\left[E-\frac{b}{x}+\alpha_{eff}x-\frac{l(l+1)}{2\mu}x^{2}-\frac{m_{j}\;C}{2\pi}x^{2}\right]\chi(x)=0$
(26)

Now we will expand the variable $y$ ($y=x-\delta$) around zero with
$\delta=1/r_{0}$ where $r_{0}$ is the mean meson radius. Therefore we have:

$\frac{b}{x}\approx
b\left(\frac{3}{\delta}-\frac{3x}{\delta^{2}}+\frac{x^{2}}{\delta^{3}}\right)$
(27)

Thus we can write the equation 26 in following form:

$\frac{d^{2}\chi(x)}{dx^{2}}+\frac{2x}{x^{2}}\frac{d\chi(x)}{dx}+\frac{2\mu}{x^{4}}\left[E-\frac{3b}{\delta}+\left(\alpha_{eff}+\frac{3b}{\delta^{2}}\right)x-\left(\frac{l(l+1)}{2\mu}+\frac{m_{j}\;C}{2\pi}+\frac{b}{\delta^{3}}\right)x^{2}\right]\chi(x)=0$
(28)

The coefficients $H_{i}$ with $i=0,1,2$, are given by:

$H_{0}=-2\mu\left(E-\frac{3b}{\delta}\right)$ (29)

$H_{1}=2\mu\left(\alpha_{eff}+\frac{3b}{\delta^{2}}\right)$ (30)

$H_{2}=-2\mu\left(\frac{l(l+1)}{2\mu}+\frac{m_{j}\;C}{2\pi}+\frac{b}{\delta^{3}}\right)$
(31)

Thus, we can write:

$\frac{d^{2}\chi(x)}{dx^{2}}+\frac{2x}{x^{2}}\frac{d\chi(x)}{dx}+\frac{2\mu}{x^{4}}[-H_{0}+H_{1}\;x+H_{2}\;x^{2}]R(x)=0$
(32)

Now, we can compare the differential equation 32 with 57, we get:

$\sigma(x)=x^{2}\phantom{AA},\phantom{AA}\tilde{\tau}(x)=2x\phantom{AA},\phantom{AA}\tilde{\sigma}(x)=-H_{0}+H_{1}\;x+H_{2}\;x^{2}$
(33)

The function $\pi(s)$ from equation 61 is given by:

$\pi(x)=\pm\sqrt{(k-H_{2})x^{2}-H_{1}x+H_{0}}$ (34)

Choosing the negative solution the polynomial inside of square must have
discriminant equal to zero, so we get:

$k=\frac{1}{4H_{0}}\left(H^{2}_{1}+4H_{0}H_{2}\right)$ (35)

Then,

$\pi(x)=-\frac{1}{2\sqrt{H_{0}}}\left(H_{1}x-2H_{0}\right)$ (36)

and $\tau(x)$ is given by equation 64, so

$\tau(x)=2x-\frac{1}{\sqrt{H_{0}}}\left(H_{1}x-2H_{0}\right)$ (37)

Then equaling the equations 62 and 63 we have:

$\frac{H^{2}_{1}}{4H_{0}}+H_{2}-\frac{H_{1}}{2\sqrt{H_{0}}}=\frac{H_{1}}{\sqrt{H_{0}}}n-n(n+1)$
(38)

Solving the equation 38 to $H_{0}$, we get:

$\sqrt{H_{0}}=\frac{H_{1}}{(1+2n)\pm\sqrt{1-4H_{2}}}$ (39)

In this way, we can obtain the energy levels and wavefunctions as a function
of the rotation parameters via the Eigenvalue equation

$\sqrt{-2\mu\left(E-\frac{3b}{\delta}\right)}=\frac{2\mu\left(\alpha_{eff}+\frac{3b}{\delta^{2}}\right)}{(1+2n)+\sqrt{1+4l(l+1)+\frac{4\mu\;m_{j}\;C}{\pi}+\frac{8\mu
b}{\delta^{3}}}}$ (40)

###  III.2 Mass and vorticity

The previous equation gives the binding energy in the rotating frame

$E_{n,l,m}=\frac{3b}{\delta}-\frac{2\mu\left(\alpha_{eff}+\frac{3b}{\delta^{2}}\right)^{2}}{\left[(1+2n)+\sqrt{1+4l(l+1)+\frac{4\mu\;m_{j}\;C}{\pi}+\frac{8\mu
b}{\delta^{3}}}\right]^{2}}$ (41)

With he wave function determination is explicitly given in Appendix C using
the method outlined in Appendix A . We can obtain the quarkonium mass from the
equation 41:

$M=2m_{q}+E_{n,l,m}$ (42)
$M=2m_{q}+\frac{3b}{\delta}-\frac{2\mu\left(\alpha_{eff}+\frac{3b}{\delta^{2}}\right)^{2}}{\left[(1+2n)+\sqrt{1+4l(l+1)+\frac{4\mu\;m_{j}\;C}{\pi}+\frac{8\mu
b}{\delta^{3}}}\right]^{2}}$ (43)

Figure 4: Charmonium mass in relation to parameter $C$ for various energy and
momentum levels (notation follows atomic physics quantum numbers) using
$T=168$ MeV

Figure 5: Bottomonium mass in relation to parameter $C$ for various energy and
momentum levels. Notation follows atomic physics quantum numbers using $T=168$
MeV

Fig. 4 and Fig. 5 gives a quantitative value of the binding energy for
quarkonia polarized with and opposite the vorticity as a function of the
circulation. Unsurprisingly, one is the opposite of the other. This, however,
is not experimentally detectable since spin alignment measurements do not
distinguish between $m=\pm 1$ states. The mass difference in section 4 would
therefore appear as an impact parameter dependent widening of the quarkonium
state. However, since a similar widening occurs in any kind of in-medium
interaction, particularly in interactions leading to the melting of quarkonium
classic3 , a univocal proof of spin-orbit coupling can not be obtained by
measuring $E_{\pm 1}$ alone. As we will show in the next section section IV ,
however, off-diagonal matrix elements could be of help here.

###  III.3 Vorticity and melting

We can use the mass dependence calculated in the previous sub-section to
study, qualitatively, how the temperature for quarkonium melting (melting
temperature, $T_{melt}$) changes under a non-inertial frame. This calculation
can be considered to be a semi-classical estimate of the imaginary part of the
energy calculated in the previous subsection section III.2 .

Using a semi-classical analysis with $\left\langle p\right\rangle\sim 1/r$ and
$\left\langle p^{2}\right\rangle\sim 1/r^{2}$ we can write the energy by using

$E(r)=\frac{1}{2\mu r^{2}}-m_{j}\;\omega+V(r)$ (44)

For the $J/\psi$ melting, we shall use the potential that takes into
consideration the Debye screen where $\lambda_{D}$ is Debye screening length,
which is given by:

$V(r)=-\alpha_{eff}\frac{e^{-r/\lambda_{D}}}{r}$ (45)

Next, using the circulation theorem 23. Then,

$E(r)=\left[\frac{1}{2\mu}-\frac{m_{j}\;C}{2\pi}\right]\frac{1}{r^{2}}-\frac{\alpha_{eff}e^{-r/\lambda_{D}}}{r}$
(46)

The bound state is defined when the energy 46 has a minimum, so we can write:

$\frac{dE(r)}{dr}=-\left[\frac{1}{\mu}-\frac{m_{j}\;C}{\pi}\right]\frac{1}{r^{3}}+\alpha_{eff}e^{-r/\lambda_{D}}\left[1+\frac{r}{\lambda_{D}}\right]\frac{1}{r^{2}}=0$
(47)

, then

$r\left[1+\frac{r}{\lambda_{D}}\right]\alpha_{eff}e^{-r/\lambda_{D}}=\frac{1}{\mu}-\frac{m_{j}\;C}{\pi}$
(48)

Making the following variable change $\tilde{r}=r/\lambda_{D}$, we get:

$f(\tilde{r})=\tilde{r}\left(1+\tilde{r}\right)e^{-\tilde{r}}=\frac{1}{\alpha_{eff}\lambda_{D}}\left[\frac{1}{\mu}-\frac{m_{j}\;C}{\pi}\right]$
(49)

The maximum value of $f(\tilde{r})$ is $0.840$ at $\tilde{r}=1.92$ in this
limit we will have the non bound states. Thereby, we can write the inequality
for this limit:

$\frac{1}{\alpha_{eff}\lambda_{D}}\left[\frac{1}{\mu}-\frac{m_{j}\;C}{\pi}\right]>0.840$
(50)

At this point, we will utilize the Debye lenght $\lambda_{D}$ value from the
lowest-order perturbative QCD WONG  222The $\lambda_{D}$ value to non-inertial
frames is not same as the inertial case. Nevertheless, the Debye mass to an
inertial frame is an acceptable estimate because the non-inertial effect gives
just a second-order contribution.:

$\lambda_{D}=\sqrt{\frac{2}{9\pi\alpha_{eff}}}\frac{1}{T}$ (51)

We can obtain an estimate for the melting temperature $T_{melt}$ in non-
inertial frames using the equation 49 in the $f(\tilde{r})$ maximum value,
then:

$0.840\sqrt{\frac{2}{9\pi\alpha_{eff}}}\frac{1}{T_{melt}}=\frac{1}{\alpha_{eff}}\left[\frac{1}{\mu}-\frac{m_{j}\;C}{\pi}\right]$
(52)

So,

$T_{melt}=0.840\sqrt{\frac{2\alpha_{eff}}{9\pi}}\left[\frac{1}{\mu}-\frac{m_{j}\;C}{\pi}\right]$
(53)

We can see from figure 6 that the melting temperature $T_{melt}$ depends on
the spin polarization quantum number $m$. In particular for $m=1$ this
temperature increases considerably. For $m=0$ it predictably does not change
and for $m=-1$ it decreases but, and this is a fundamental point, it decreases
a lot less than the increase for $m=1$. Therefore, in a vortical medium
polarized quarkonia will be much more likely to survive, while un-polarized or
anti-polarized quarkonium’s survival probability does not change that much
w.r.t. quarkonia in a non-rotating medium.

This finding is qualitatively important since it shows that vorticity can link
quarkonium suppression and polarization via ”distillation”. In a vortical
medium, the melting probability of quarkonium states will depend strongly on
their polarization. Such a mechanism will result, analogously to the $\eta$
eta , to a strong and novel dependence of quarkonium abundance on centrality,
which could be investigated quantitatively by a hydrodynamic model.

One can be ”brave” and try to apply our potential model to the meson $\phi$,
which is formed of strange and anti-strange quarks. The strong-vorticity
dependent melting might be able to explain the strong spin alignment observed
in experiment, which seems incompatible with Cooper-Frye freezeout star . In
the picture described here, melting temperature for aligned $\phi$ rises a
lot, while melting temperature for anti-aligned and non-aligned $\phi$ stays
nearly the same. Thus, the large apparent spin alignment of $\phi$s comes from
a ”distillation” process where only polarized $\phi$s survive, and this
increases relative $\phi$ spin alignment. Of course, applying potential models
to $\phi$ is not justified theoretically and requires model-building, so such
a solution will need considerable quantitative and phenomenological
development.

Figure 6: An estimate of melting temperature $T_{melt}$ for $J/\psi$ (left)
and $\Upsilon$ (right) melting to different spin polarization $m$ values.

In the absence of a detailed hydrodynamic modeling of the correlation between
charmonium suppression and polarization, however, distinguishing the effects
outlined above from other variations of charmonium suppression with
centrality, namely the effect of centrality dependence of temperature, looks
complicated, with vorticity just adding an “event-by-event widening” to the
processes of melting and regeneration of quarkonia. The next section, however,
provides a direct experimentally measurable indication of spin-orbit non-
equilibrium in quarkonium polarization measurements

##  IV Density matrix elements and vorticity

Turning our attention back to the density matrix we can write this operator on
basis of energy in the following way:

$\hat{\rho}=e^{-\beta\hat{H}}\phantom{AA},\phantom{AA}\beta=\frac{1}{T}$ (54)

Where the eigenvalues is given by the equation 41. In this moment, we will
make a rotation to lab frame, so:

$\hat{\rho}^{r}=U(\theta_{r},\phi_{r})\;\hat{\rho}\;U^{-1}(\theta_{r},\phi_{r})$
(55)

Now, we can expand the equation 55 in following way:

$\rho^{r}_{m,m^{\prime}}=\sum_{m^{\prime\prime},m^{\prime\prime\prime}}e^{i(m^{\prime\prime\prime}-m)\phi_{r}}d^{j}_{m,m^{\prime\prime}}(\theta_{r})\rho_{m^{\prime\prime},m^{\prime\prime\prime}}\left[d^{j}_{m^{\prime\prime\prime},m^{\prime}}(\theta_{r})\right]^{-1}\phantom{AA},\phantom{AA}\rho_{m^{\prime\prime},m^{\prime\prime\prime}}=\frac{1}{Z}e^{\beta
E_{m^{\prime\prime\prime}}}\delta_{m^{\prime\prime},m^{\prime\prime\prime}}$
(56)

We can relate the density matrix coefficients from variable $C$ as both
density matrix coefficients and energy variation depend on the parameter $C$.
So we can relate these two values we obtain the figure 8. Then, we can relate
the energy $(\Delta E_{m,m^{\prime}}=E_{m}-E_{m^{\prime}})$ with off-diagonal
density matrix $Re[\rho_{-1,0}]$ and $Re[\rho_{-1,1}]$ using the values
$\phi=0$, $\theta_{r}=1.75\pm 0.10$ to Collins-Soper refer to transverse
momentum range $2<p_{T}<4$ (GeV).

Figure 7: The alignment factor $\rho_{00}$ in relation C for charmonium and
bottomonium using $T=168$ MeV.

Figure 8: The off-diagonal density matrix components in relation to energy variation $\Delta E_{m,m^{\prime}}$ using $T=168$ MeV. | $m\ (GeV)$ | $\alpha_{eff}$ | $b\ (GeV^{2})$ | $\delta\ (GeV)$ | $\mu\ (GeV)$ |  $\theta_{r}$ (rad)
---|---|---|---|---|---|---
$c\bar{c}$ | 1.209 | 1.244 | 0.2 | 0.231 | 0.6045 | $1.75\pm 0.10$
$b\bar{b}$ | 4.823 | 1.569 | 0.2 | 0.378 | 2.4115 | $1.58\pm 0.10$
Table 2:  Eigenvalue quarkonium parameters obtained from AHM  with
$\theta_{r}$ determined in system equation 9 to Collins-Soper frame to
$J/\psi$ and $\Upsilon$ with respectively ranges $4<p_{T}<6\ (GeV)$ and
$p_{T}<15\ (GeV)$.

As can be seen in Fig. 7 the obtained $\rho_{00}$ are compatible with the
estimates shown in section II

Figure 9: The off-diagonal density matrix components in relation to energy
variation $\Delta E_{m,m^{\prime}}$ using $T=168$ MeV.

We can obtain the off-diagonal density matrix components in relation to the
quarkonium energy using the parameters shown in Table 2. In figure 8 and 9, we
can note the relation between the alignment factor, $\rho_{00}$, with the
circulation parameter C. It is evident that an increase of C increments
$\rho_{00}$; however, this increase will depend on the type of meson. As we
can see the bottomonium alignment factor is larger than the charmonium when
compared under the same parameter C.

Unlike the figures in the previous sections, the $x,y$ axes of the plots in
this section are independently measurable. $\Delta E$ will be manifest event-
by-event as an invariant mass correlation. $\rho_{i,j}$ can be obtained wang
from $dN/d\Omega$ according to a harmonic analysis of Eq. (1). Barring
acceptance effects, there is no non-vortical dynamics capable of generating a
correlation between these two quantities. The most direct analysis to do is a
correlation between $\lambda_{\theta,\phi,\theta\phi}$ and the invariant mass
and width. A definite signal of the correlation can be used, using the
functions described in this section, to extract the angle between vorticity
and heavy quark spin polarization. This can then be used to constrain models
of spin hydrodynamics such as uscausality ; uscas2 ; usdiss ; relax ; jeon ;
gursoy ; hongo ; hongo2 . This section’s results, therefore, can be used as a
baseline for developing an experimental analysis capable of probing non-
equilibrium between spin and vorticity using quarkonium probes.

The most direct analysis to do is a correlation between
$\lambda_{\theta,\phi,\theta\phi}$ and the invariant mass and width. A
definite signal of the correlation can be used, using the functions described
in this section, to extract the angle between vorticity and heavy quark spin
polarization. This can then be used to constrain models of spin hydrodynamics
such as uscausality ; uscas2 ; usdiss ; relax ; jeon ; gursoy ; hongo ; hongo2
. This section’s results, therefore, can be used as a baseline for developing
an experimental analysis capable of probing non-equilibrium between spin and
vorticity using quarkonium probes.

##  V Conclusion

The results obtained in this work are certainly to be considered more like a
rough estimate than a quantitative analysis, since they use an early version
of the Cornell potential model. However we think that this model is good
enough to get a physical intuition of the problem of linking spin-vorticity
non-equilibrium to the quarkonium state in a rotating frame, and,
respectively, the quarkonium state to experimental observables. As section III
illustrates, our formalism in fact reproduces the reasonable expectation of
how the binding energy (mass) and melting probability (width) of the
quarkonium state respond to rotation.

Section IV therefore takes the model further and examines what happens when
rotation and spin are not aligned, and suggests experimental measurements. A
positive experimental observation of the observables described in the previous
section will provide evidence that quarkonium state polarization is not
thermally aligned, but in fact “remembers” a polarization state which is quite
distinct from the vorticity state.

As argued in section II and Fig. 1 , one in fact expects quarks to be aligned
transversely, which would also follow from initial state dynamics pqcd1 ;
pqcd2  and from the large transverse vorticity present in the beginning of the
collision bec . However, the vortex where the quarkonia will be could up could
well be aligned to the flow gradients that developed on the time-scale of a
hydrodynamic expansion, which will be in the longitudinal direction. The
combination of the two will result in a density matrix exhibiting large
corrections to the expectation value in zanna , with the magnitude of these
correcting correlated to the invariant mass of quarkonium.

A very interesting associated result is the strong dependence, non-monotonic
in both $T$ and $m$, of the dissociation temperature $T_{melt}$ of quarkonium
shown in section III.3 . It suggests a new ”distillation” mechanism for
quarkonium polarization, where aligned quarkonia have a stronger probability
to survive in medium while non-aligned and anti-aligned quarkonia simply melt.
Speculatively considering the $\phi$ a quarkonium state could explain the
strong alignment signal seen in star , though some phenomenological work is
needed to confirm such a model is viable.

So far, as seen in section II there is no direct experimental evidence for
such non-equilibrium. In the short term, an azimuthal modulation of
$\lambda_{\theta,\phi,\theta\phi}$ (Fig. 3 ) could provide a strong indication
that such non-equilibrium should be investigated. In the long term, an
experiment with sufficient resolution in invariant mass and capable of
reconstructing the spin alignment could be able to perform the analysis
advocated in section IV , providing more direct experimental evidence.

Quarkonium polarization, therefore, could well be a very promising observable
for experimentally probing spin hydrodynamics in the quark-gluon plasma, a
theoretically interesting and challenging topic that so far had little contact
with phenomenology.

GT thanks CNPQ bolsa de produtividade 306152/2020-7, bolsa FAPESP 2021/01700-2
and participation in tematic FAPESP, 2017/05685-2 and grant
BPN/ULM/2021/1/00039 from the Polish National Agency for Academic Exchange
K.J.G. is supported by CAPES doctoral fellowship 88887.464061/2019-00. This
work was supported in part by the Polish National Science Centre Grant Nos
2018/30/E/ST2/00432. We thank the hospitality of the Jagellonian university
when part of this work was performed.

##  Appendix A Nikiforov-Uvarov method

We will use the method developed by Nikiforov and Uvarov to solve a
differential equation where it is possible to reduce to hypergeometric
function. We can do it too to Schrödinger equation, which can write in the
following form:

$\psi^{{}^{\prime\prime}}(s)+\frac{\tilde{\tau}(s)}{\sigma(s)}\psi^{{}^{\prime}}(s)+\frac{\tilde{\sigma}(s)}{\sigma^{2}(s)}\psi(s)=0$
(57)

Using the variable separation technique,

$\psi(s)=\phi(s)y(s)$ (58)

Using the variable separation technique, we can see that the equation must
satisfy the following expression NU :

$\frac{\phi^{{}^{\prime}}(s)}{\phi(s)}=\frac{\pi(s)}{\sigma(s)}$ (59)

so we can simplify as

$\sigma(s)y^{{}^{\prime\prime}}(s)+\tau(s)y(s)+\lambda y(s)=0$ (60)

At this moment, we need to define the $\pi(s)$ function in the following form:

$\pi(s)=\frac{\sigma^{{}^{\prime}}(s)-\tilde{\tau}(s)}{2}\pm\sqrt{\left(\frac{\sigma^{{}^{\prime}}(s)-\tilde{\tau}(s)}{2}\right)^{2}-\tilde{\sigma}(s)+k\sigma(s)}$
(61)

and $\lambda$ parameter is given by:

$\lambda=k+\pi^{{}^{\prime}}(s)$ (62)

The Schrödinger equation eigenvalue obtained from this method is given by:

$\lambda_{n}=-n\tau^{{}^{\prime}}(s)-\frac{n(n-1)}{2}\sigma^{{}^{\prime\prime}}(s)$
(63)

where $\tau(s)$ is defined by:

$\tau(s)=\tilde{\tau}(s)+2\pi(s)$ (64)

The solution to hypergeometric equation is given by Rodrigues relation:

$y_{n}(s)=\frac{B_{n}}{\tilde{\rho}(s)}\frac{d^{n}}{ds^{n}}[\sigma^{n}\tilde{\rho}(s)]$
(65)

$B_{n}$ is the normalization coefficient and $\tilde{\rho}(s)$ must obey the
following relation:

$(\sigma(s)\tilde{\rho}(s))^{{}^{\prime}}=\tau(s)\tilde{\rho}(s)$ (66)

Then

$\frac{\tilde{\rho}^{{}^{\prime}}}{\tilde{\rho}}=\frac{\tau-\sigma^{{}^{\prime}}}{\sigma}$
(67)

We have the following solution:

$\tilde{\rho}(s)=\exp\left(\int\frac{\tau-\sigma^{{}^{\prime}}}{\sigma}\;ds\right)$
(68)

##  Appendix B Estimate of $\boldsymbol{\lambda_{\theta,\phi,\theta\phi}}$

The blast wave model, described in detail in blast1 , is to parametrize the
transverse plane flow profile of the fireball to obtain soft observables and
relate them to particle properties via a Cooper-Frye type formula zanna . The
fireball is assumed to be an ellipse in coordinate space

$x=r_{max}\sqrt{1-\epsilon}\cos\phi,\;\;\;\;\;\;\;y=r_{max}\sqrt{1+\epsilon}\sin\phi,$
(69)

where the $r_{max}$ and $\epsilon$ represents the parameters of model and
$\phi$ is the azimutal angle in the transverse plane. To obtain the thermal
vorticity, we need first define the fluid velocity that describes the
azimuthally asymmetric fluid flow, parametrized by a Hubble-type Ansatz

$u^{\mu}=\left(\frac{t}{N},\frac{x}{N}\sqrt{1+\tilde{\delta}},\frac{y}{N}\sqrt{1-\tilde{\delta}},\frac{z}{N}\right)$
(70)

Where the normalization factor is equal to
$N=\sqrt{\tau^{2}-(x^{2}-y^{2})\tilde{\delta}}$ and the proper time is
$\tau^{2}=t^{2}-x^{2}-y^{2}-z^{2}$ ensuring, $u^{\mu}u_{\mu}=1$. At this
point, we can obtain the thermal vorticity from the equation eq6 defined
above:

$\omega_{\mu\nu}=-\frac{1}{2T}(\partial_{\mu}u_{\nu}-\partial_{\nu}u_{\mu})-\frac{1}{2T^{2}}(u_{\mu}\partial_{\nu}T-u_{\nu}\partial_{\mu}T).$
(71)

The components are given explicitly in blast1 . Now, we will suppose that the
$\omega$ is defined in the following way:

$\omega=\sqrt{\omega^{2}_{x}+\omega^{2}_{y}+\omega^{2}_{z}}$ (72)

Where the spacial components to vorticity are given by
$\omega^{i}=\frac{1}{2}\epsilon^{ijk}\omega^{jk}$. Now, make the average on
the azimuthal angle, $\phi$, using the parameters to centrality $30–60\%$ in
midrapidity $\eta=0$ in equation 72. We have:

$\left<\omega\right>=\frac{1}{2\pi}\int^{2\pi}_{0}\omega(r_{max},\epsilon,\tilde{\delta})d\phi=24.08\;MeV$
(73)

We can make an estimative for this value for thermal vorticity 73 using
equation 23 in the following form: $C=2\pi\left<\omega\right>\;r_{0}^{2}$.
Now, using $r_{0}=1/\delta$, we have:

$C=\frac{2\pi\left<\omega\right>}{\delta^{2}}$ (74)

Using the numerical parameters reported in blast1 ; blast2  via a fit to
experimental data, we can obtain the vorticity.

To obtain an estimate of the maximum value of
$\lambda_{\theta,\phi,\theta\phi}$ we also need a numerical estimate of
$\theta_{r},\phi_{r}$ in Eq. 5. This can be obtained from inverting Fig. 2.
The final result, calculated by plugging in these values in equations 2, 74
and 56 is shown in Table 1.

This is a very rough estimate, since these parameters were fitted at a
different energy, and since only the second Fourier coefficient of the flow is
considered (in vorticity all coefficients mix non-trivially). The latter
reason means we cannot calculate the exact azimuthal dependence. However,
since $v_{n}$ changes slowly with energy we are confident this is a good
order-of-magnitude estimate of the effect experimentalists are looking for. If
azimuthal modulation is found, more realistic models will be needed to
describe it quantitatively.

##  Appendix C Quarkonium wave function determination

Now, we can obtain the wave function from the equation 68 and using the
equation 37 and $\sigma(x)=x^{2}$. From it, we get

$\tilde{\rho}(x)=C_{1}x^{-H_{1}/\sqrt{H_{0}}}e^{-2\sqrt{H_{0}}/x}$ (75)

We are able to obtain the value of function $\phi(x)$ from equation 36 and
using the values of equation 36. Then

$\frac{\phi^{{}^{\prime}}(x)}{\phi(x)}=-\frac{1}{2\sqrt{H_{0}}}\left(\frac{H_{1}}{x}-\frac{2H_{0}}{x^{2}}\right)$
(76)

Then

$\phi(x)=C_{2}\;x^{-H_{1}/2\sqrt{H_{0}}}e^{-\sqrt{H_{0}}/x}$ (77)

We can determine the $y_{n}(r)$ from equations 65, 75 and $\sigma(r)$. Then,

$y_{n}(r)=B_{n}x^{H_{1}/\sqrt{H_{0}}}e^{2\sqrt{H_{0}}/x}\frac{d^{n}}{dx^{n}}\left[x^{2n-H_{1}/\sqrt{H_{0}}}e^{-2\sqrt{H_{0}}/x}\right]$
(78)

From equation 77, 78 and putting in 58, we have that

$\psi(x)=B_{n}x^{H_{1}/2\sqrt{H_{0}}}e^{\sqrt{H_{0}}/x}\frac{d^{n}}{dx^{n}}\left[x^{2n-H_{1}/\sqrt{H_{0}}}e^{-2\sqrt{H_{0}}/x}\right]$
(79)

Thereby,changing the variable $x=1/r$ we can write the radial wave fuction
from equation 79 and relation $\psi(r)=rR(r)$. So:

$R_{n}(r)=B_{n}\;r^{-(H_{1}/2\sqrt{H_{0}})}e^{\sqrt{H_{0}}r}\left[-r^{2}\frac{d}{dr}\right]^{n}\left[x^{2n-H_{1}/\sqrt{H_{0}}}e^{-2\sqrt{H_{0}}/x}\right]$
(80)

Thus, we can write the wave fuction that following form:

$\Psi_{n,l,m}(r,\theta,\phi)=R_{n}(r)Y^{m}_{l}(\theta,\phi)$ (81)

Using the value of 80, we get:

$\Psi_{n,l,m}(r,\theta,\phi)=B_{n}\;r^{-(H_{1}/2\sqrt{H_{0}})}e^{\sqrt{H_{0}}r}\left[-r^{2}\frac{d}{dr}\right]^{n}\left[x^{2n-H_{1}/\sqrt{H_{0}}}e^{-2\sqrt{H_{0}}/x}\right]Y^{m}_{l}(\theta,\phi)$
(82)

## References

  * (1) Matsui, T. & Satz, H. $J/\psi$ suppression by quark-gluon plasma formation. Physics Letters B 178, 416 – 422 (1986)
  * (2) A. Andronic, F. Arleo, R. Arnaldi, A. Beraudo, E. Bruna, D. Caffarri, Z. C. del Valle, J. G. Contreras, T. Dahms and A. Dainese, et al. Eur. Phys. J. C 76 (2016) no.3, 107 doi:10.1140/epjc/s10052-015-3819-5 [arXiv:1506.03981 [nucl-ex]].
  * (3) Andronic, A., Braun-Munzinger, P., Köhler, M. K., Mazeliauskas, A., Redlich, K., Stachel, J., amp; Vislavicius, V. (2021). The multiple-charm hierarchy in the statistical hadronization model. Journal of High Energy Physics, 2021(7). https://doi.org/10.1007/jhep07(2021)035.
  * (4) R. L. Thews, M. Schroedter and J. Rafelski, Phys. Rev. C 63 (2001), 054905 doi:10.1103/PhysRevC.63.054905 [arXiv:hep-ph/0007323 [hep-ph]].
  * (5) A. Rothkopf, Phys. Rept. 858 (2020), 1-117 doi:10.1016/j.physrep.2020.02.006 [arXiv:1912.02253 [hep-ph]].
  * (6) F. Becattini and M. A. Lisa, Ann. Rev. Nucl. Part. Sci. 70 (2020), 395-423 doi:10.1146/annurev-nucl-021920-095245 [arXiv:2003.03640 [nucl-ex]].
  * (7) S. Acharya et al. [ALICE], [arXiv:1910.14408 [nucl-ex]].
  * (8) M. S. Abdallah et al. [STAR], Nature 614, no.7947, 244-248 (2023) doi:10.1038/s41586-022-05557-5 [arXiv:2204.02302 [hep-ph]].
  * (9) Z. T. Liang and X. N. Wang, Phys. Lett. B 629, 20-26 (2005) doi:10.1016/j.physletb.2005.09.060 [arXiv:nucl-th/0411101 [nucl-th]].
  * (10) X. L. Xia, H. Li, X. G. Huang and H. Zhong Huang, Phys. Lett. B 817 (2021), 136325 doi:10.1016/j.physletb.2021.136325 [arXiv:2010.01474 [nucl-th]].
  * (11) K. J. Gonçalves and G. Torrieri, Phys. Rev. C 105 (2022) no.3, 034913 doi:10.1103/PhysRevC.105.034913 [arXiv:2104.12941 [nucl-th]].
  * (12) X. L. Sheng, L. Oliva, Z. T. Liang, Q. Wang and X. N. Wang, [arXiv:2205.15689 [nucl-th]].
  * (13) D. Montenegro and G. Torrieri, Phys. Rev. D 100, no.5, 056011 (2019) doi:10.1103/PhysRevD.100.056011 [arXiv:1807.02796 [hep-th]].
  * (14) D. Montenegro and G. Torrieri, [arXiv:2004.10195 [hep-th]].
  * (15) D. Montenegro and G. Torrieri, Phys. Rev. D 107 (2023) no.7, 076010 doi:10.1103/PhysRevD.107.076010 [arXiv:2207.00537 [hep-th]].
  * (16) S. Bhadury, W. Florkowski, A. Jaiswal, A. Kumar and R. Ryblewski, Phys. Lett. B 814 (2021), 136096 doi:10.1016/j.physletb.2021.136096 [arXiv:2002.03937 [hep-ph]].
  * (17) S. Shi, C. Gale and S. Jeon, [arXiv:2008.08618 [nucl-th]].
  * (18) A. D. Gallegos, U. Gürsoy and A. Yarom, SciPost Phys. 11 (2021), 041 doi:10.21468/SciPostPhys.11.2.041 [arXiv:2101.04759 [hep-th]].
  * (19) M. Hongo, X. G. Huang, M. Kaminski, M. Stephanov and H. U. Yee, JHEP 11 (2021), 150 doi:10.1007/JHEP11(2021)150 [arXiv:2107.14231 [hep-th]].
  * (20) K. Hattori, M. Hongo, X. G. Huang, M. Matsuo and H. Taya, Phys. Lett. B 795 (2019), 100-106 doi:10.1016/j.physletb.2019.05.040 [arXiv:1901.06615 [hep-th]].
  * (21) Acharya S, Adamová D, Adler A, Adolfsson J, Aggarwal MM, Aglieri Rinella G, et al. First measurement of quarkonium polarization in nuclear collisions at the LHC. Physics Letters B. 2021;815:136146.
  * (22) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B 566 (2000), 275 doi:10.1016/S0550-3213(99)00693-8 [arXiv:hep-ph/9907240 [hep-ph]].
  * (23) Lucha, W. (1991). Bound states of Quarks. Physics Reports, 200(4), 127-240. https://doi.org/10.1016/0370-1573(91)90001-3
  * (24) M. Buzzegoli and K. Tuchin, Nucl. Phys. A 1030 (2023), 122577 doi:10.1016/j.nuclphysa.2022.122577 [arXiv:2209.03991 [nucl-th]].
  * (25) H. Kim, S. Cho and S. H. Lee, [arXiv:2212.14570 [hep-ph]].
  * (26) Yazdankish, E. (2020). Solving of the Schrodinger equation analytically with an approximated scheme of the woods-saxon potential by the systematical method of Nikiforov-Uvarov. International Journal of Modern Physics E, 29(06), 2050032. https://doi.org/10.1142/s0218301320500329
  * (27) Arfken, G. B. (2005). Mathematical Methods for Physicists (3rd ed.). Elsevier Academic Press.
  * (28) Anandan, J., Suzuki, J. (n.d.). Quantum Mechanics in a Rotating Frame. https://doi.org/https://arxiv.org/pdf/quant-ph/0305081.pdf
  * (29) Anandan, J. (1992). Comment on spin-rotation-gravity coupling. Physical Review Letters, 68(25), 3809-3810. https://doi.org/10.1103/physrevlett.68.3809
  * (30) Ahmadov, A. I., Abasova, K. H., amp; Orucova, M. S. (2021). Bound State Solution schrodinger equation for extended Cornell potential at finite temperature. Advances in High Energy Physics, 2021, 1–13. https://doi.org/10.1155/2021/1861946
  * (31) C. Y. Wong. Introduction to High-Energy Heavy-Ion Collisions, World Scientific Publisher, 1994\. https://doi.org/10.1142/1128
  * (32) F. Becattini, V. Chandra, L. Del Zanna and E. Grossi, Annals Phys. 338 (2013), 32-49 doi:10.1016/j.aop.2013.07.004 [arXiv:1303.3431 [nucl-th]].
  * (33) G. Torrieri, Phys. Rev. C 98, no.1, 014901 (2018) doi:10.1103/PhysRevC.98.014901 [arXiv:1802.09011 [nucl-th]].
  * (34) P. Faccioli, V. Knünz, C. Lourenco, J. Seixas and H. K. Wöhri, Phys. Lett. B 736 (2014), 98-109 doi:10.1016/j.physletb.2014.07.006 [arXiv:1403.3970 [hep-ph]].
  * (35) V. Cheung and R. Vogt, [arXiv:2203.10154 [hep-ph]].
  * (36) J. Adams, A. Ewigleben, S. Garrett, W. He, T. C. Huang, P. M. Jacobs, X. Ju, M. A. Lisa, M. Lomnitz and R. Pak, et al. Nucl. Instrum. Meth. A 968 (2020), 163970 doi:10.1016/j.nima.2020.163970 [arXiv:1912.05243 [physics.ins-det]]
  * (37) G. Giacalone, L. Yan, J. Noronha-Hostler and J. Y. Ollitrault, Phys. Rev. C 94 (2016) no.1, 014906 doi:10.1103/PhysRevC.94.014906 [arXiv:1605.08303 [nucl-th]].
  * (38) W. Florkowski, A. Kumar, R. Ryblewski and A. Mazeliauskas, Phys. Rev. C 100 (2019) no.5, 054907 doi:10.1103/PhysRevC.100.054907 [arXiv:1904.00002 [nucl-th]].
  * (39) K. J. Gonçalves and G. Torrieri, Baryon polarization and Spin alignment of vector mesons in a thermal model with spin-vorticity non-equilibrium, to be submitted