format: text

text:
# Non-symplectic symmetries and bi-Hamiltonian structures of the rational Harmonic Oscillator

## 1 Non-symplectic symmetries

It is well known that there is a close relation  between integrability and the existence of alternatives structures (see e.g.  for a recent paper) and also that integrable systems are systems endowed with a great number of symmetries. The purpose of this letter is to analyze, in the particular case of the $`n=2`$ harmonic oscillator, how these additional structures arise from the existence of dynamical symmetries of non-symplectic (non-canonical) type.

Let $`(M,\omega _0,H)`$ be a Hamiltonian system and $`\mathrm{\Gamma }_H`$ the associated Hamiltonian vector field, defined by $`i(\mathrm{\Gamma }_H)\omega _0=dH`$. A (infinitesimal) dynamical symmetry of this system is a vector field $`YX(M)`$ such that $`[Y,\mathrm{\Gamma }_H]=0`$. When $`Y`$ is a dynamical but non-symplectic symmetry of the system, then we have that (i) the dynamical vector field $`\mathrm{\Gamma }_H`$ is bi-Hamiltonian, and (ii) the function $`Y(H)`$ is the new Hamiltonian, and therefore it is a constant of motion.

A sketch of the proof - of this statement is as follows: The vector field $`Y`$ does not preserve $`\omega _0`$ and, as it is a non-canonical transformation, it determines a new $`2`$-form $`\omega _Y=_Y\omega _0`$ ($`_Y`$ denotes de Lie derivative with respect to $`Y`$). As $`Y`$ is a symmetry, $`[Y,\mathrm{\Gamma }_H]=0`$, then $`_Yi_{\mathrm{\Gamma }_H}=i_{\mathrm{\Gamma }_H}_Y`$, and, consequently,

$$i_{\mathrm{\Gamma }_H}\omega _Y=i_{\mathrm{\Gamma }_H}_Y\omega _0=_Yi_{\mathrm{\Gamma }_H}\omega _0=_Y(dH)=d(YH).$$

Therefore, the 2-form $`\omega _Y`$ is admissible for the dynamical vector field $`\mathrm{\Gamma }_H`$, i.e. $`_{\mathrm{\Gamma }_H}\omega _Y=0`$, which is weakly bi-Hamiltonian with respect to the original symplectic $`2`$-form $`\omega _0`$ and the new structure $`\omega _Y`$. Of course the particular form of $`\omega _Y`$ depends on $`Y`$ and, in some cases, it can be just a constant multiple of $`\omega _0`$ (trivial bi-Hamiltonian system). In some other cases $`\omega _Y`$ may be a degenerate $`2`$-form with a nontrivial kernel. In any case, the vector field $`\mathrm{\Gamma }_H`$ is a dynamical system solution of the following two equations

$$i(\mathrm{\Gamma }_H)\omega _0=dH,\mathrm{and}i(\mathrm{\Gamma }_H)\omega _Y=d[Y(H)]$$

Therefore the function $`H_Y=Y(H)`$, that must be a constant of motion, can be considered as a new Hamiltonian for $`\mathrm{\Gamma }_H`$.

## 2 Bi-Hamiltonian structures of the rational Harmonic Oscillator

The two-dimensional harmonic oscillator

$$H=\frac{1}{2}(p_x^2+p_y^2)+\frac{1}{2}(\lambda _1^2x^2+\lambda _2^2y^2)$$
(1)
has the two one-degree of freedom energies, $`I_1=E_x`$ and $`I_2=E_y`$, as fundamental constants of motion. The superintegrability of the rational case, $`\lambda _1=m\lambda _0`$, $`\lambda _2=n\lambda _0`$, with $`m,n`$, can be proved by making use of a complex formalism . Let $`K_x`$, $`K_y`$, be the following two functions $`K_x=p_x+\mathrm{i}m\lambda _0x`$ and $`K_y=p_y+\mathrm{i}n\lambda _0y`$; then the Hamiltonian $`H`$ and the canonical symplectic form $`\omega _0`$ become

$$H=\frac{1}{2}\left(K_xK_x^{}+K_yK_y^{}\right),$$

and

$$\omega _0=\frac{\mathrm{i}}{2m\lambda _0}dK_xdK_x^{}+\frac{\mathrm{i}}{2n\lambda _0}dK_ydK_y^{}.$$

We have

$$\{K_x,K_x^{}\}=2\mathrm{i}m\lambda _0,\{K_y,K_y^{}\}=2\mathrm{i}n\lambda _0,$$

and therefore, the evolution equations are

$$\frac{d}{dt}K_x=\mathrm{i}m\lambda _0K_x,\frac{d}{dt}K_y^{}=\mathrm{i}n\lambda _0K_y^{}.$$

Hence, the complex function $`J`$ defined as

$$J=K_x^n(K_y^{})^m$$
(2)
is a constant of motion that determines two different real first integrals, $`I_3=Im(J)`$ and $`I_4=Re(J)`$, which are polynomials in the momenta of degree $`m+n1`$ and $`m+n`$, respectively. As an example, for the Isotropic case, $`\lambda _1=\lambda _2=\lambda _0`$, we obtain

$`Re(J)`$ $`=`$ $`p_xp_y+\lambda _{0}^{}{}_{}{}^{2}xy,`$ (3)
$`Im(J)`$ $`=`$ $`\lambda _0(xp_yyp_x).`$
$`Im(J)`$ is just the angular momentum, and $`Re(J)`$ is the non-diagonal component of the Fradkin tensor . For the first non-isotropic case, $`\lambda _1=\lambda _0`$, $`\lambda _2=2\lambda _0`$, we arrive to

$`Re(J)`$ $`=`$ $`p_x^2p_y+\lambda _{0}^{}{}_{}{}^{2}(4yp_xxp_y)x,`$ (4)
$`Im(J)`$ $`=`$ $`(xp_yyp_x)p_x+\lambda _{0}^{}{}_{}{}^{2}x^2y.`$
This complex procedure provides not just the fundamental constant $`I_3`$, but the pair $`(I_3,I_4)`$; although the ‘partner’ function $`I_4`$ is not independent (is a function of $`I_1`$, $`I_2`$, $`I_3`$), we will see that it plays an important rôle, since it is closely concerned with the bi-Hamiltonian formalism. In fact, we will take the complex function $`J`$ as our starting point for the search of symmetries, but $`J`$ means not only one but two functions, $`I_3`$ and $`I_4`$.

The Noether theorem in the Hamiltonian formalism states that all constants of motion arise from canonical symmetries of the Hamiltonian function. Moreover, in differential geometric terms, the infinitesimal symmetries are simply those corresponding to the Hamiltonian vector fields, with respect to the canonical structure $`\omega _0`$, defined by the constants of motion. In this particular case, the above complex function $`J`$ given by (2) arises from a symmetry of (1) represented by the complex vector field $`X_J`$ defined by

$$i(X_J)\omega _0=dJ,X_J(H)=0.$$
(5)
In the following, and for easy of notation, we will suppose $`\lambda _0=1`$.

###### Proposition 1

The complex vector field $`X_J`$, defined by (5) as the canonical infinitesimal symmetry associated to $`J`$, can be written as a linear combination of two dynamical but non-symplectic symmetries of $`\mathrm{\Gamma }_H`$.

Proof: Let us denote by $`Y_{xm}`$ and $`Y_{yn}`$ the Hamiltonian vector fields of $`K_x`$ and $`K_y`$

$$i(Y_{xm})\omega _0=dK_x,i(Y_{yn})\omega _0=dK_y,$$

with coordinate expressions

$$Y_{xm}=\frac{}{x}\mathrm{i}m\frac{}{p_x},Y_{yn}=\frac{}{y}\mathrm{i}n\frac{}{p_y}.$$

Notice that, as $`H=I_1+I_2`$ with $`|K_x|^2=2I_1`$ and $`|K_y|^2=2I_2`$, we have

$$\mathrm{\Gamma }_H=Re(K_x^{}Y_{xm}+K_yY_{yn}^{}).$$

Then, the complex vector field $`X_J`$, canonical infinitesimal symmetry of the Harmonic oscillator, can be written as the following linear combination

$$X_J=nY+mY^{},$$

where the $`Y`$, $`Y^{}`$, are given by

$$Y=(K_x^{(n1)}K_y^m)Y_{xm},Y^{}=(K_x^nK_y^{(m1)})Y_{yn}^{}.$$

The important point is that these two vector fields, $`Y`$ and $`Y^{}`$, are neither locally-Hamiltonian with respect to $`\omega _0`$

$$_Y\omega _00,_Y^{}\omega _00,$$

nor infinitesimal symmetries of the Hamiltonian

$$_YH0,_Y^{}H0.$$

Concerning the Lie bracket of $`Y`$ with the dynamical vector field $`\mathrm{\Gamma }_H`$, it is given by

$$[Y,\mathrm{\Gamma }_H]=(K_x^{n1}K_y^m)[Y_{xm},\mathrm{\Gamma }_H]\mathrm{\Gamma }_H(K_x^{n1}K_y^m)Y_{xm}$$

but as

$$[Y_{xm},\mathrm{\Gamma }_H]=\mathrm{i}Y_{xm},[Y_{yn}^{},\mathrm{\Gamma }_H]=\mathrm{i}Y_{yn}^{},$$

and

$`\mathrm{\Gamma }_H(K_x^{n1}K_y^m)`$ $`=`$ $`(n1)(\mathrm{i}m)(K_x^{n1}K_y^m)+m(\mathrm{i}n)(K_x^{n1}K_y^m)`$
$`=`$ $`\mathrm{i}m(K_x^{n1}K_y^m),`$
we arrive to

$$[Y,\mathrm{\Gamma }_H]=0.$$

Thus $`Y`$ is a dynamical but non-symplectic (non-canonical) symmetry of $`\mathrm{\Gamma }_H`$. It can be proved, in a similar way, that this property is also true for $`Y^{}`$. Notice that, in the language of 1-forms, this property arises from the fact that $`dJ`$ splits as a sum of two non-closed 1-forms that, nevertheless, remain invariant under $`\mathrm{\Gamma }_H`$, that is, $`dJ=n\varphi _1+m\varphi _2`$, $`d\varphi _r0`$, $`_{\mathrm{\Gamma }_H}(\varphi _r)=0`$, $`r=1,2`$.

Two new structures can be obtained from $`\omega _0`$ by Lie derivation with respect to $`Y`$ and $`Y^{}`$. If we denote by $`\omega _Y`$ and $`\omega _Y^{}`$ these two new 2-forms, $`\omega _Y=_Y\omega _0`$ and $`\omega _Y^{}=_Y^{}\omega _0`$, then we obtain

$$\omega _Y=m(K_x^{(n1)}K_y^{(m1)})dK_xdK_y^{},\omega _Y^{}=n(K_x^{(n1)}K_y^{(m1)})dK_xdK_y^{}.$$

In the following we will denote by $`\mathrm{\Omega }`$ the complex 2-form defined as

$$\mathrm{\Omega }=dK_xdK_y^{}=\mathrm{\Omega }_1+\mathrm{i}\mathrm{\Omega }_2$$

where the two real 2-forms, $`\mathrm{\Omega }_1=Re(\mathrm{\Omega })`$ and $`\mathrm{\Omega }_2=Im(\mathrm{\Omega })`$, take the form

$$\mathrm{\Omega }_1=mndxdy+dp_xdp_y,\mathrm{\Omega }_2=mdxdp_y+ndydp_x.$$

Notice that $`\omega _Y`$ and $`\omega _Y^{}`$ satisfy the relation $`n\omega _Y+m\omega _Y^{}=0`$. Actually, this can be considered as a consequence of the fact that $`X_J`$ is locally Hamiltonian with respect to the canonical form $`\omega _0`$.

###### Proposition 2

The dynamical vector field $`\mathrm{\Gamma }_H`$ of the rational Harmonic Oscillator is a bi-Hamiltonian system with respect to $`(\omega _0,\omega _Y)`$.

Proof: Notice that

$$i(\mathrm{\Gamma }_H)\omega _Y=m(K_x^{(n1)}K_y^{(m1)})i(\mathrm{\Gamma }_H)\mathrm{\Omega },$$

and as

$$i(\mathrm{\Gamma }_H)\mathrm{\Omega }=\mathrm{\Gamma }_H(K_x)dK_y^{}\mathrm{\Gamma }_H(K_y^{})dK_x=\mathrm{i}mK_xdK_y^{}+\mathrm{i}nK_y^{}dK_x,$$

we obtain that

$$i(\mathrm{\Gamma }_H)\omega _Y=\mathrm{i}md(K_x^nK_y^m).$$

Thus, $`\mathrm{\Gamma }_H`$ is Hamiltonian vector field with respect to $`\omega _Y`$ with $`K_x^nK_y^m`$ as Hamiltonian function. Moreover, we can also compute the action of $`Y`$ on $`H`$; a direct calculation gives

$$H_YY(H)=\mathrm{i}m(K_x^nK_y^m).$$

To conclude, we have found that the integral of motion $`J`$ determines the following bi-Hamiltonian system

$$i(\mathrm{\Gamma }_H)\omega _0=dH,i(\mathrm{\Gamma }_H)\omega _Y=dH_Y.$$

Remark first that $`\mathrm{\Gamma }_H`$ is bi-Hamiltonian with respect to two different structures: the canonical symplectic form $`\omega _0`$ and another one, $`\omega _Y`$, which is complex. If we write $`\omega _Y=\omega _4+\mathrm{i}\omega _3`$, then $`\mathrm{\Gamma }_H`$ can be considered as a bi-Hamiltonian system with respect to the following three real forms $`(\omega _0,\omega _3,\omega _4)`$ (i.e. it is a three-Hamiltonian system). The $`\omega _0`$-Hamilton equation determined by $`J`$,

$$i(X_J)\omega _0=dJ,$$

is also complex; thus it determines two real Hamiltonian equations

$$i(X_4)\omega _0=dI_4,i(X_3)\omega _0=dI_3,$$

with $`X_4`$, $`X_3`$, given by $`X_J=X_4+\mathrm{i}X_3`$.

As a second remark, the complex 2-form $`\mathrm{\Omega }=dK_xdK_y^{}`$ is well defined but it is not symplectic. In fact, it can be proved that $`\mathrm{\Omega }_1=Re(\mathrm{\Omega })`$ and $`\mathrm{\Omega }_2=Im(\mathrm{\Omega })`$ satisfy

$$\mathrm{\Omega }_1\mathrm{\Omega }_1=\mathrm{\Omega }_2\mathrm{\Omega }_2=mn(dxdydp_xdp_y),\mathrm{and}\mathrm{\Omega }_1\mathrm{\Omega }_2=0,$$

so we obtain

$$\mathrm{\Omega }\mathrm{\Omega }=(\mathrm{\Omega }_1\mathrm{\Omega }_1\mathrm{\Omega }_2\mathrm{\Omega }_2)+2\mathrm{i}\mathrm{\Omega }_1\mathrm{\Omega }_2=0.$$

Thus, the degenerate character of $`\mathrm{\Omega }`$ is directly related with its complex nature. Moreover, the kernel of $`\mathrm{\Omega }`$ is the distribution generated by $`Y_{xm}`$ and $`Y_{yn}^{}`$,

$$Ker\mathrm{\Omega }=\{fY_{xm}+gY_{yn}^{}f,g:^2\times ^2\},$$

therefore it satisfies

$$[Ker\mathrm{\Omega },\mathrm{\Gamma }_H]Ker\mathrm{\Omega }.$$

Finally, the 2-form $`\omega _Y`$ is also degenerate. We obtain, $`Ker\omega _Y=Ker\mathrm{\Omega }`$, because of the relation between $`\mathrm{\Omega }`$ and $`\omega _Y`$. However $`\omega _3`$ and $`\omega _4`$, defined as $`\omega _Y=\omega _4+\mathrm{i}\omega _3`$, are symplectic real forms. Moreover, the form $`\omega _0+\mathrm{\Omega }`$ is symplectic because of $`\{K_x,K_y^{}\}=0`$.

## 3 Recursion operators

The bi-Hamiltonian structure $`(\omega _0,\omega _Y)`$ defines a complex recursion operator $`R_Y`$ by

$$\omega _Y(X,Y)=\omega _0(R_YX,Y),X,YX(M),$$

or, equivalently, $`R_Y=\widehat{\omega }_0^1\widehat{\omega }_Y`$. Since it is complex, it can be written as $`R_Y=R_4+\mathrm{i}R_3`$, so that $`R_4`$ and $`R_3`$ satisfy the relations

$$\omega _3(X,Y)=\omega _0(R_3X,Y),\mathrm{and}\omega _4(X,Y)=\omega _0(R_4X,Y).$$

Thus, we have that $`R_3`$ and $`R_4`$ are given by $`R_3=\widehat{\omega }_0^1\widehat{\omega }_3`$ and $`R_4=\widehat{\omega }_0^1\widehat{\omega }_4`$.

The important point is that the complex 2-form $`\mathrm{\Omega }=dK_xdK_y^{}`$ can be decomposed as $`\mathrm{\Omega }=\mathrm{\Omega }_1+\mathrm{i}\mathrm{\Omega }_2`$, where both 2-forms, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, are symplectic. Hence, we have, in addition to $`R_3`$ and $`R_4`$, two other recursion operators $`R_1`$ and $`R_2`$ associated with the bi-Hamiltonian structures provided by $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, respectively.

###### Proposition 3

The tensor fields $`R_1`$ and $`R_2`$ are invertible operators which anticommute and satisfy $`R_2^2=R_1^2`$.

Proof: As $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ are symplectic forms, the operators $`R_1`$ and $`R_2`$ are invertible. Their coordinate expressions are given by

$`R_1`$ $`=`$ $`{\displaystyle \frac{}{y}}dp_x{\displaystyle \frac{}{x}}dp_y+mn\left({\displaystyle \frac{}{p_x}}dy{\displaystyle \frac{}{p_y}}dx\right),`$ (7)
$`R_2`$ $`=`$ $`m\left({\displaystyle \frac{}{y}}dx+{\displaystyle \frac{}{p_x}}dp_y\right)+n\left({\displaystyle \frac{}{x}}dy+{\displaystyle \frac{}{p_y}}dp_x\right).`$ (8)
Therefore,

$$R_1^2=mnId,R_2^2=mnId.$$

Moreover we have

$$R_2R_1=n\left(\frac{}{x}dp_xm^2\frac{}{p_x}dx\right)m\left(\frac{}{y}dp_yn^2\frac{}{p_y}dy\right),$$

and $`R_1R_2=R_2R_1`$; therefore $`R_2R_1+R_1R_2=0`$.

We recall that the relation between $`\omega _Y`$ and $`\mathrm{\Omega }`$ is $`\omega _Y=\omega _4+\mathrm{i}\omega _3=m\mathrm{lK}\mathrm{\Omega }`$, with the complex function $`\mathrm{lK}`$ given by $`\mathrm{lK}=K_x^{(n1)}K_y^{(m1)}=\mathrm{lK}_r+\mathrm{i}\mathrm{lK}_i`$. Thus, the above two tensor fields, $`R_3`$ and $`R_4`$, are given by

$`R_4`$ $`=`$ $`m\left(\mathrm{lK}_rR_1\mathrm{lK}_iR_2\right),`$ (9)
$`R_3`$ $`=`$ $`m\left(\mathrm{lK}_rR_2+\mathrm{lK}_iR_1\right),`$ (10)
and, making use of the preceding proposition, we arrive to

$$R_4^2=R_3^2=rId,\mathrm{with}r=m^2(mn)^2|\mathrm{lK}|^2,$$

and

$$R_4R_3=m^2|\mathrm{lK}|^2R_1R_2,R_3R_4=m^2|\mathrm{lK}|^2R_2R_1,$$

where the modulus of $`\mathrm{lK}`$ is function of the two first integrals, $`I_1`$ and $`I_2`$,

$$|\mathrm{lK}|^2=\mathrm{lK}_r^2+\mathrm{lK}_i^2=(2I_1)^{(n1)}(2I_2)^{(m1)}.$$

Thus, the two tensor fields, $`R_3`$ and $`R_4`$, anticommute as well.

###### Proposition 4

The complex operator $`R_Y=R_4+\mathrm{i}R_3`$ is such that $`Image(R_Y)=KerR_Y`$.

Proof: Let $`X`$ be a generic vector field on the phase space

$$X=a\frac{}{x}+b\frac{}{y}+c\frac{}{p_x}+d\frac{}{p_y},$$

with $`a,b,c,d`$, arbitrary functions. Then, since $`R_Y=R_4+\mathrm{i}R_3`$ is given by $`R_Y=m\mathrm{lK}(R_1+\mathrm{i}R_2)`$, the subspace Image of $`R_Y`$ is given by

$$Image(R_Y)=\left\{ZX(^2\times ^2)\right|Z=m\mathrm{lK}\left(R_1(X)+\mathrm{i}R_2(X)\right)\}.$$

We obtain

$$R_1(X)+\mathrm{i}R_2(X)=(d\mathrm{i}nb)Y_{xm}+(c+\mathrm{i}ma)Y_{yn}^{}.$$

Thus, $`Image(R_Y)`$ is made up of linear combinations of $`Y_{xm}`$ and $`Y_{yn}^{}`$ with arbitrary complex functions as coefficients. But, since $`KerR_Y=Ker\omega _Y`$ and $`Ker\omega _Y`$ coincides with $`Ker\mathrm{\Omega }`$, which is also spanned by $`Y_{xm}`$ and $`Y_{yn}^{}`$, we arrive to $`Image(R_Y)=KerR_Y`$. Consequently $`R_Y^2=R_YR_Y=0`$.

Given a bi-Hamiltonian system on a manifold $`M`$, $`i(\mathrm{\Gamma })\omega _0=dH_0`$ and $`i(\mathrm{\Gamma })\omega _1=dH_1`$, the point is that the tensor field $`R`$, that was just defined by the relation between $`\omega _1`$ and $`\omega _0`$, induces a sequence of structures. Starting with the basic Hamiltonian system $`(\omega _0,\mathrm{\Gamma }_0=\mathrm{\Gamma },dH_0)`$ we can construct a sequence of 2-forms $`\omega _k`$, of vector fields $`\mathrm{\Gamma }_k`$, and of 1-forms $`\alpha _k`$, $`k=1,2,\mathrm{}`$, defined by $`\widehat{\omega }_k=\widehat{\omega }_0R^k`$, $`\mathrm{\Gamma }_k=R^k(\mathrm{\Gamma }_0)`$, and $`\alpha _k=R^k(dH_0)`$. Then it follows that

$`i(\mathrm{\Gamma }_0)\omega _1=i(\mathrm{\Gamma }_1)\omega _0=dH_1,`$ (11)
$`i(\mathrm{\Gamma }_0)\omega _2=i(\mathrm{\Gamma }_1)\omega _1=i(\mathrm{\Gamma }_2)\omega _0=\alpha _2,`$
where

$`\widehat{\omega }_1=\widehat{\omega }_0R,`$ $`\widehat{\omega }_2=\widehat{\omega }_1R,`$ (12)
$`\mathrm{\Gamma }_1=R(\mathrm{\Gamma }_0),`$ $`\mathrm{\Gamma }_2=R(\mathrm{\Gamma }_1),`$ (13)
$`dH_1=R^{}(dH_0),`$ $`\alpha _2=R^{}(dH_1),`$
The 1-form $`\alpha _2`$ is not necessarily exact, but if there is $`H_2`$ such that $`\alpha _2=dH_2`$, then the vector field $`\mathrm{\Gamma }_1`$ is a bi-Hamiltonian system as well. An interesting case is when $`\alpha _2`$ is not exact but there exist a nonvanishing function $`F_2`$ and another function $`H_2`$ such that $`\alpha _2=F_2dH_2`$. Then $`F_2^1`$ is an integrating factor for $`\alpha _2`$, and the vector field $`\mathrm{\Gamma }_1`$ is said to be quasi-bi-Hamiltonian ,.

Coming back to the rational harmonic oscillator as a bi-Hamiltonian system, $`i(\mathrm{\Gamma }_H)\omega _0=dH_0`$ and $`i(\mathrm{\Gamma }_H)\omega _Y=dH_Y`$, the situation is as follows:

(i) The action of $`R_Y`$ is such that $`\mathrm{\Gamma }_H\mathrm{\Gamma }_0`$ becomes $`\mathrm{\Gamma }_1=R_Y(\mathrm{\Gamma }_H)=\mathrm{i}mX_J`$,

(ii) $`dH_0`$ transforms into $`dH_Y=R_Y^{}(dH_0)=\mathrm{i}dJ`$, and

(iii) $`\omega _0`$ becomes $`\omega _Y`$ such that $`\widehat{\omega }_Y=\widehat{\omega }_0R_Y`$.

We have proved that $`R_Y^2=0`$ because of Proposition 4; therefore, $`\mathrm{\Gamma }_1`$ transforms into the new field, $`\mathrm{\Gamma }_2=R_Y(\mathrm{\Gamma }_1)=R_Y^2(\mathrm{\Gamma }_H)=0`$, while $`dH_Y`$ transforms into $`\alpha _2=R_Y^{}(dH_Y)=R_Y^2(dH_0)=0`$. Hence, it follows that the equation

$$i(\mathrm{\Gamma }_H)\omega _2=i(\mathrm{\Gamma }_1)\omega _1=i(\mathrm{\Gamma }_2)\omega _0=\alpha _2$$

becomes

$$i(\mathrm{\Gamma }_1)\omega _1=0.$$

Notice that this last equation corresponds to the property $`i(X_J)\omega _Y=0`$.

The Harmonic Oscillator can be considered as a complex and weakly bi-Hamiltonian system, or alternatively, as endowed with two different real bi-Hamiltonian structures

$$i(\mathrm{\Gamma }_H)\omega _0=dH_0,i(\mathrm{\Gamma }_H)\omega _4=mdI_3,i(\mathrm{\Gamma }_H)\omega _3=mdI_4.$$

One real structure gives rise to $`\mathrm{\Gamma }_{13}`$ defined by $`\mathrm{\Gamma }_{13}=R_3(\mathrm{\Gamma }_H)`$, and the other one to $`\mathrm{\Gamma }_{14}=R_4(\mathrm{\Gamma }_H)`$. They are such that

$$i(\mathrm{\Gamma }_{14})\omega _0=i(\mathrm{\Gamma }_H)\omega _4=mdI_3,i(\mathrm{\Gamma }_{13})\omega _0=i(\mathrm{\Gamma }_H)\omega _3=mdI_4.$$

Moreover taking into account that $`R_4^2=rId`$, we arrive to

$`\mathrm{\Gamma }_{24}`$ $`=`$ $`R_4(\mathrm{\Gamma }_{14})=R_4^2(\mathrm{\Gamma }_H)=r\mathrm{\Gamma }_H,`$ (14)
$`\alpha _{24}`$ $`=`$ $`R_4^{}(mdI_3)=R_4^2(dH_0)=rdH_0,`$ (15)
$`\widehat{\omega }_{24}`$ $`=`$ $`\widehat{\omega }_4R_4=\widehat{\omega }_0R_4^2=r\widehat{\omega }_0,`$
and similar results for $`R_3`$.

Now, making use of all these relations, we can prove the following final proposition concerning the properties of the vector fields $`X_3`$ and $`X_4`$.

###### Proposition 5

Let $`X_3`$ and $`X_4`$ denote the two infinitesimal canonical symmetries generating the two constants of motion $`I_3`$ and $`I_4`$. Then $`X_3`$ and $`X_4`$ are quasi-bi-Hamiltonian systems. Moreover, $`\omega _4(X_3,\mathrm{\Gamma }_H)=\omega _3(X_4,\mathrm{\Gamma }_H)=0`$.

Proof: The rational Harmonic Oscillator is endowed with the two constants of motion $`I_3`$ and $`I_4`$ which means, via the Hamiltonian Noether theorem, the existence of two symmetries. They are geometrically represented by two vector fields, $`X_3`$ and $`X_4`$, that can be uniquely determined as solutions of the following two equations

$$i(X_3)\omega _0=dI_3,i(X_4)\omega _0=dI_4.$$

Then we have

$$\mathrm{\Gamma }_{13}=R_3(\mathrm{\Gamma }_H)=mX_4,\mathrm{\Gamma }_{14}=R_4(\mathrm{\Gamma }_H)=mX_3.$$

Hence, if we denote by $`f_{34}`$ the function $`f_{34}=m(mn)^2|\mathrm{lK}|^2,`$ we arrive to

$$i(X_3)\omega _0=dI_3,i(X_3)\omega _4=f_{34}dH_0,$$

and

$$i(X_4)\omega _0=dI_4,i(X_4)\omega _3=f_{34}dH_0.$$

So, both $`X_3`$ and $`X_4`$ are quasi-bi-Hamiltonian systems.

A direct consequence of this property is that the dynamical vector field $`\mathrm{\Gamma }_H`$ is orthogonal to $`X_3`$ with respect to the symplectic structure $`\omega _4`$,

$$i(X_3)i(\mathrm{\Gamma }_H)\omega _0=0,\mathrm{and}i(X_3)i(\mathrm{\Gamma }_H)\omega _4=0.$$

Similarly, we obtain

$$i(X_4)i(\mathrm{\Gamma }_H)\omega _0=0,\mathrm{and}i(X_4)i(\mathrm{\Gamma }_H)\omega _3=0.$$

Finally, $`X_3`$ and $`X_4`$ are orthogonal vector fields with respect to both structures, $`\omega _3`$ and $`\omega _4`$:

$$i(X_3)i(X_4)\omega _3=0,\mathrm{and}i(X_3)i(X_4)\omega _4=0.$$

## Acknowledgements.

This work is dedicated to Prof. Domingo González, from the University of Zaragoza, on the occasion of his retirement. Support of Spanish DGI projects, BFM-2000-1066-C03-01 and FPA-2000-1252, is acknowledged.


text:
# Resonance widening in spherical GW detectors: model descriptions of the dissipation processes

## 1 Introduction

Spherical Gravitational Wave (GW) detectors will almost certainly be, in one of its variants (solid, hollow or dual), the next generation of acoustic antennae, due to their multimode-multifrequency capabilities , as well as their potentially enhanced sensitivity relative to their currently operating cylindrical counterparts . Conviction that this is going to be the case has encouraged a remarkable research effort within the GW community, and a variety of important topics have been addressed, ranging from theoretical to practical aspects of the problem. Several countries worldwide are currently developing projects to build and operate spherical detectors .

Some of the salient properties of a spherical body as a detector of GWs can be accurately established by means of an idealised model, where dissipative effects can be safely overlooked . Such effects are however very important when it comes to determining the ultimate sensitivity of the device, due to their fundamental bearing on the noise characteristics of the antenna via the fluctuation–dissipation theorem, more specifically its spectral density.

In a solid elastic body there are many oscillation modes in addition to the fundamental one chosen for GW detection, which therefore add noise to the frequency band of operation. It has been suggested  that normal mode analysis of this contribution to the noise level can be either inconvenient or too expensive in terms of computational cost if heat is e.g. generated by a laser spot on the oscillating solid, such as happens in an interferometric GW detector or in a dual sphere . Experimental measurements clearly show that this is the case .

The naïve approach, and actually the one most often taken, to account for dissipative effects is to add an ad hoc term proportional to the velocity in the solid’s eigenfrequency equation. This results in a damped oscillatory behaviour of the sort

$$e^{\omega t/Q}\mathrm{sin}\omega t$$
(1)
where the quality factor Q is assumed constant (i.e., time independent), and accounts for the linewidth of the mode with frequency $`\omega `$. It is well known  that this quality factor is different for different oscillation modes, so the set of $`Q`$’s is very large, and the procedure appears highly artificial. On the other hand, there has been criticism of the concept of dissipation forces proportional to velocity in terms of actual noise spectrum in an elastic solid , and suggestion that a better description may rather derive from the use of complex frequencies.

In the specialised literature on the subject (see  and references therein), viscoelastic effects are often described by means of so called constituent equations. These are extensions of Hooke’s law relating the stress and strain tensors in the solid. Much like in a simple one-dimensional spring, internal friction forces can be considered proportional to the instantaneous velocity of the oscillating mass, hence constituent equations usually contain time derivatives of those tensors, and depend on a small number of viscosity parameters, to be added to the elastic Lamé coefficients $`\lambda `$ and $`\mu `$ —see below. There is however no unique way in which this idea can be carried over from a one-dimensional system into a constituent equations set for a three-dimensional solid, and so different alternatives result in different models for the purpose.

In this paper we propose to study and discuss the results of applying to a spherical GW detector the equations of various such phenomenological models, in view to determine how quality factors change from mode to mode in each case. This, it is hoped, will generate some insight into the nature of the viscous processes which take place in a specific spherical GW detector, and help to assess on the basis of spectral measurements which particular class of viscoelastic solid a given material belongs to. In turn, better understanding of material’s macroscopic properties should also contribute relevant information to the currently important issue of spherical GW detector design. In section 2 we present the general equations, and in the subsequent sections we successively consider the Kelvin–Voigt, Maxwell, Standard Linear and Genaralised Mechanical models. As we shall see, the sphere’s vibration eigenmodes always group into the usual families of toroidal (purely torsional) and spheroidal modes, but different quality factor dependences arise in different models. In section 6 we match our theoretical results to a number of experimental measurements of quality factors performed on various spheres, and discuss which is the most appropriate model to fit such data. In the final section we present a summary of conclusions, and two appendices are added to clarify a few mathematical technicalities.

## 2 The general equations

GWs bathing the Earth are known to be extremely weak —see e.g.  for a review— so that the classical equations of linear Elasticity are very good to describe the GW induced motions of a spherical antenna in the expected frequency range, roughly 102–104 Hz. These equations are

$$\rho \frac{^2s_i}{t^2}\frac{\sigma _{ij}}{x_j}=f_i$$
(2)
where $`𝐬(𝐱,t)`$ is the field of displacements in the elastic body, and $`𝐟(𝐱,t)`$ is the GW induced density of forces acting on the solid —see  for full technical details. $`\sigma _{ij}`$ is the the stress tensor, and is related to the strain tensor

$$s_{ij}\frac{1}{2}\left(\frac{s_i}{x_j}+\frac{s_j}{x_i}\right)$$
(3)
through a set of constituent equations. In the case of a non-dissipative (ideal) solid these are simply the expression of Hooke’s law. In a dissipative one, such equations include time derivatives of both $`s_{ij}`$ and $`\sigma _{ij}`$ to account for internal friction effects. Constituent equations are of the following general type:

$$L(s_{ij},\dot{s}_{ij},\ddot{s}_{ij},\mathrm{};\sigma _{ij},\dot{\sigma }_{ij},\ddot{\sigma }_{ij},\mathrm{})=0$$
(4)
In this paper we shall limit ourselves to linear constituent equations, an excellent approach for a GW detector, as already stressed. In the simplest instance, only first order derivatives will appear in (4), and we shall consider this first. Then we shall also devote some attention to more complicted models.

The equations of motion (2) must be supplemented with suitable boundary conditions. We shall prescribe the usual ones

$$\sigma _{ij}n_j=0\mathrm{at}r=R$$
(5)
where n is a unit outward-pointing vector, expressing that the surface of the sphere is free from any tensions and/or tractions.

## 3 Kelvin–Voigt model

This model assumes that the solid is homogeneous and isotropic, and is characterised by the following constituent equations :

$$\sigma _{ij}=\left(\lambda +\lambda ^{}\frac{}{t}\right)s_{kk}\delta _{ij}+2\left(\mu +\mu ^{}\frac{}{t}\right)s_{ij}$$
(6)
The constants $`\lambda `$ and $`\mu `$ are the usual Lamé coefficients describing the purely elastic behaviour of the body , while the positive coefficients $`\lambda ^{}`$ and $`\mu ^{}`$ parametrise its viscous properties111 These coefficients are actually analogous to those which describe the viscosity of fluids in Hydrodynamics: shear viscosity ($`\mu ^{}`$) and bulk viscosity ($`2\mu ^{}+3\lambda ^{}`$)., which are proportional to the change rate of the strain tensor, $`_ts_{ij}`$.

If equations (6) are replaced into (2) we obtain the equations of motion for $`𝐬(𝐱,t)`$:

$$\rho \frac{^2𝐬}{t^2}=\left(\mu +\mu ^{}\frac{}{t}\right)^2𝐬\left[(\lambda +\mu )+(\lambda ^{}+\mu ^{})\frac{}{t}\right](𝐬)+𝐟(𝐱,t)$$
(7)
The solution to this system of coupled equations can be expressed in terms of a Green function, the construction of which requires explicit knowledge of its eigenmode solutions. As is well known, the latter correspond to the free oscillations of the solid —no density of external forces in the rhs of (7). As usual, we attempt to find such eigen-solutions in the factorised form

$$𝐬(𝐱,t)=T(t)𝐬(𝐱)$$
(8)
which results in the following (dots on symbols stand for time derivatives):

$`\rho \ddot{T}(t)𝐬(𝐱)`$ $`=`$ $`\left[\mu T(t)+\mu ^{}\dot{T}(t)\right]^2𝐬(𝐱)`$ (9)
$`+`$ $`\left[(\lambda +\mu )T(t)+(\lambda ^{}+\mu ^{})\dot{T}(t)\right]\left(𝐬(𝐱)\right)`$
Next we use the well known decomposition of a three dimensional vector field into its irrotational a divergence free components :

$$𝐬(𝐱)=𝐬_l(𝐱)+𝐬_t(𝐱),𝐬_t(𝐱)=\times 𝐬_l(𝐱)=0$$
(10)
By the methods described in A, it can be seen that the equations satisfied by $`𝐬_t(𝐱)`$, $`𝐬_l(𝐱)`$, and $`T(t)`$ are (11lnyaaaubabibkbmbncbcg)–(11lnyaaaubabibkbmbncbcj), i.e.,

$`^2𝐬_t+𝒦^2𝐬_t=0`$ (11a)
$`\mu T+\mu ^{}\dot{T}+𝒦^2\ddot{T}=0`$ (11b)
and

$`^2𝐬_l+𝒬^2𝐬_l=0`$ (11la)
$`\mu T+\mu ^{}\dot{T}+𝒬^2\ddot{T}=0`$ (11lb)
Since $`T(t)`$ must fulfil both equations (11b) and (11lb), the separation constants $`𝒦^2`$ and $`𝒬^2`$ are not independent. The binding relationship is established after it is realised that the solution to those equations is of the form

$$T(t)=e^{\gamma t}$$
(11lm)
where $`\gamma `$ is in general a complex quantity. We readily obtain

$`𝒬^2`$ $`=`$ $`\rho \gamma ^2\left[\lambda +2\mu +\gamma (\lambda ^{}+2\mu ^{})\right]^1`$ (11lna)
$`𝒦^2`$ $`=`$ $`\rho \gamma ^2\left[\mu +\gamma \mu ^{}\right]^1`$ (11lnb)
The values $`\gamma `$ can possibly take on are determined by the boundary conditions, equations (5). The reader is referred to B for a detailed description of the eigenvalue algebra of this problem. Just as in the non-dissipative case, there are seen to be two families of eigenmodes: toroidal (purely torsional) and spheroidal. Viscous effetcs are however small in practice, as inferred from the narrow linewidth of the measured resonances. This means that the following inequalities hold in any cases of interest to us:

$$\frac{\mu ^{}}{\mu },\frac{\lambda ^{}}{\mu }\frac{1}{\omega }$$
(11lno)
where $`\omega `$ is the frequency of the mode considered. We shall use these inequalities to estimate the roots of the eigenvalue equation (11lnyaaaubabibkbmbncbcy) (see B) perturbatively from the non-dissipative ones, already known —see  for full details. Clearly thus, our procedure will be valid for the lower frequency modes. We proceed sequentially for the two families of eigenmodes.

### 3.1 Toroidal modes

These correspond to the solutions to (see (11lnyaaaubabibkbmbncbcy))

$$\beta _1(𝒦R)=0$$
(11lnp)
which is formally identical to the toroidal eigenvalue equation for a non-dissipative solid sphere . If we call $`k_{nl}^T`$ the toroidal wave numbers of the latter, we have

$$𝒦_{nl}^T=k_{nl}^T=\sqrt{\frac{\rho }{\mu }}\omega _{nl}^T$$
(11lnq)
and hence, by equation (11lnb),

$$\gamma _{nl}^T=(\omega _{nl}^T)^2\frac{\mu ^{}}{2\mu }+i\omega _{nl}^T\sqrt{1\left(\frac{\omega _{nl}^T\mu ^{}}{2\mu }\right)}i\omega \omega ^2\frac{\mu ^{}}{2\mu }$$
(11lnr)
where the last approximation depends on the validity of the assumption (11lno). Expression (11lnr) nicely shows how this Kelvin-Voigt model predicts exponentially damped eigenmode oscillations. If we recall that such damping is expediently described in terms of a quality factor Q —see (1) above— then we discover that the prediction of the model is that

$$Q_{nl}^T=\frac{2\mu }{\mu ^{}}\left(\omega _{nl}^T\right)^1$$
(11lns)
i.e., the quality factor for toroidal modes is inversely proportional to the frequency of the mode. The amplitudes of these modes have the form

$$𝐬_{KV}^T(𝐱,t)=𝐬_E^T(𝐱,t)e^{\omega t/Q},Q=\frac{2\mu }{\mu ^{}\omega }$$
(11lnt)
where the subindex $`KV`$ stands for ‘Kelvin-Voigt’, while $`E`$ refers to the standard frictionless case, whose amplitudes are those given e.g. in reference .

### 3.2 Spheroidal modes

A second alternative to find a non-trivial solution of the linear system (11lnyaaaubabibkbmbncbcs) is to impose the condition

$$\beta _4(𝒬R,\frac{\lambda +\gamma \lambda ^{}}{\mu +\gamma \mu ^{}})\beta _3(𝒦R)l(l+1)\beta _1(𝒬R)\beta _1(𝒦R)=0$$
(11lnu)
This is characteristic of the spheroidal eigenmodes. By virtue of equations (11lna) and (11lnb), this relationship can be translated into a condition to be fufilled by $`\gamma `$, and which depends on the ratios $`\lambda /\mu `$, $`\lambda ^{}/\lambda `$ and $`\mu ^{}/\mu `$, as well as on the multipole index $`l`$. In this case, as we are not dealing with an eigenvalue problem of a selfadjoint operator, complex solutions to equation (11lnu) are allowed —indeed, expected. An exact solution of that equation implies a separation of its real and imaginary parts, followed by numerical calculations which determine the angular frequency and quality factor of the quasinormal mode at hand. We are interested in materials with long decay times, so we shall set up a perturbative solution to equation (11lnu), using

$$ϵ\frac{\mu ^{}\omega }{\mu }$$
(11lnv)
as the small perturbative parameter. In other words, we assume that the approximation (11lno) holds. Here, $`\omega `$ stands for a generic spheroidal eigenfrequency of the non-dissipative solid. Obviously, the unperturbed solution, i.e., that corresponding to $`ϵ`$ = 0, is the elastic solid’s solution, already discussed in reference . We thus introduce the perturbative expansion

$$\gamma =\gamma _0+\gamma _1ϵ+O(ϵ^2),\gamma _0=i\omega ,$$
(11lnw)
Using equations (11lna), (11lnb) and (11lnw), we obtain perturbative expansions for the parameters $`𝒦`$ and $`𝒬`$ which can be written as

$$𝒦=k_0+k_1ϵ+O(ϵ^2),𝒬=q_0+q_1ϵ+O(ϵ^2),$$
(11lnx)
where

$`k_0=k=\omega \sqrt{{\displaystyle \frac{\rho }{\mu }}},k_1=i\sqrt{{\displaystyle \frac{\rho }{\mu }}}\left(\gamma _1+{\displaystyle \frac{\omega }{2}}\right)`$ (11lnya)
$`q_0=q=\omega \sqrt{{\displaystyle \frac{\rho }{\lambda +2\mu }}},q_1=i\sqrt{{\displaystyle \frac{\rho }{\lambda +2\mu }}}\left(\gamma _1+{\displaystyle \frac{h^{}+2}{h+2}}{\displaystyle \frac{\omega }{2}}\right)`$ (11lnyb)
In the above equations, $`k`$ and $`q`$ are the parameters appearing in the elastic sphere’s case, and we have introduced the dimensionless ratios

$$h\frac{\lambda }{\mu },h^{}\frac{\lambda ^{}}{\mu ^{}}$$
(11lnyz)
which are both zero order quantities. We can now perform the perturbative expansion of the eigenvalue equation. In order to ease the resulting expressions, let us introduce the following notation:

$`l(l+1)\beta _1(𝒦R)`$ $`=`$ $`l(l+1)\beta _1(kR)+l(l+1)\beta _1^{}(kR)k_1Rϵ`$ (11lnyaaa)
$``$ $`B_0+B_1k_1Rϵ`$
$`\beta _1(𝒬R)`$ $`=`$ $`\beta _1(qR)+\beta _1^{}(qR)q_1Rϵ`$ (11lnyaab)
$``$ $`C_0+C_1q_1Rϵ`$
$`\beta _3(𝒦R)`$ $`=`$ $`\beta _3(kR)+\beta _3^{}(kR)k_1Rϵ`$ (11lnyaac)
$``$ $`D_0+D_1k_1Rϵ`$
$`\beta _4(𝒬R,{\displaystyle \frac{\lambda +\gamma \lambda ^{}}{\mu +\gamma \mu ^{}}})`$ $`=`$ $`\beta _4(qR,h)+\left[\beta _4^{}(qR,\lambda /\mu )q_1R{\displaystyle \frac{i}{2}}(h^{}h)j_l(qR)\right]ϵ`$ (11lnyaad)
$``$ $`A_0+(A_1q_1RiA_1^{})ϵ`$
where a prime on a $`\beta `$ function denotes differentiation with respect to its first argument. We note that the functions $`A_0,\mathrm{},D_1`$ (defined in every second line of the above set of equations) are real. With this notation, the zero-th order form of equation (11lnu) is

$$A_0D_0C_0B_0=0$$
(11lnyaaab)
which is simply the condition that $`\omega `$ be a spheroidal eigenvalue of the purely elastic case. On the other hand, the first order expansion of (11lnu) yields

$$(A_0D_1C_0B_1)k_1R+(A_1D_0C_1B_0)q_1R=iA_1^{}D_0,$$
(11lnyaaac)
whence, using the form of $`k_1`$ and $`q_1`$, the value of $`\gamma _1`$ ensues. It can be written as

$$\gamma _1=\frac{\omega }{2}[f(kR,h,h^{})]^1$$
(11lnyaaad)
where we define the dimensionless function $`f`$ by222 The case in which $`f`$ takes its simplest form is that of monopole modes. We know (see ) that when $`l`$ = 0 equation (11lnu) is no longer valid, and must be replaced by
$$\beta _4(𝒬R,\frac{\lambda +\gamma \lambda ^{}}{\mu +\gamma \mu ^{}})=0$$
(11lnyaaae)  Using now the expansion (11lnyaad) and the fact that $`A_0=0`$ for the uperturbed monopole eigenfrequencies, we obtain for the first order correction to $`𝒬`$,
$$q_1=iA_1^{}A_1^1$$
(11lnyaaaf)  and therefore, using the notation of equation (11lnyaaad) and the relationships (11lnya)-(11lnyb),
$$f(qR,h,h^{})=\left[\frac{h^{}+2}{h+2}\frac{2A_1^{}}{qRA_1}\right]^1.$$
(11lnyaaag)

$`f(kR,h,h^{})`$
$`{\displaystyle \frac{A_0D_1C_0B_1+(A_1D_0C_1B_0)(h+2)^{1/2}}{2A_1^{}D_0(kR)^1A_0D_1+C_0B_1(A_1D_0C_1B_0)(h^{}+2)(h+2)^{3/2}}}`$ (11lnyaaah)
We note that the first order correction obtained for $`\gamma `$ is real. Therefore, to this order of approximation, the frequencies of vibration remain unaltered, and are the same as those obtained for the elastic solid. Moreover, $`k_1`$ and $`q_1`$ happen to be purely imaginary. Therefore, the modulus of the radial functions appearing in the spatial part of spheroidal quasinormal modes of vibration will also be the same as those of the ideally elastic solid, for the corrections to $`k`$ and $`q`$ will just introduce, to first order, a complex phase factor.

Summing up, we have shown that while the spheroidal normal modes of vibration of an elastic solid are given by an expression of the form

$$𝐬_E^P(𝐱,t)=e^{i\omega _{nl}^Pt}\left[A_nl(r)Y_{lm}(\theta ,\phi )𝐧B_{nl}(r)i𝐧\times 𝐋Y_{lm}(\theta ,\phi )\right]$$
(11lnyaaai)
the spheroidal quasinormal modes $`𝐬_{KV}^P`$ of a Kelvin–Voigt solid are obtained from the normal modes of the elastic solid according to the following

$`𝐬_{KV}^P(𝐱,t)`$ $`=`$ $`e^{i\omega _{nl}^Pt\omega _{nl}^Pt/Q}[e^{i\chi _1(r)}A_{nl}(r)Y_{lm}(\theta ,\phi )𝐧`$ (11lnyaaaj)
$``$ $`e^{i\chi _2(r)}B_{nl}(r)i𝐧\times 𝐋Y_{lm}(\theta ,\phi )]`$
the qualitity factor being given by

$$Q_{nl}^P=\frac{2\mu }{\mu ^{}\omega _{nl}^P}f(k_{nl}^PR,h,h^{})$$
(11lnyaaak)
where, it is recalled, $`f`$ is given by (11lnyaaah) as a function of the mode and the coefficients characterising the viscoelastic properties of the solid. The real phases $`\chi _{1,2}(r)`$ can be computed from equations (11lnyaaaubabibkbmbncbcp)-(11lnyaaaubabibkbmbncbcr), (11lnya), and (11lnyb). Nevertheless, the explicit (and cumbersome) form of these phases is largely irrelevant, and whe shall omit its explicit form here . They merely introduce a position dependent shift in the phase of the vibrations which is of order $`ϵ`$, therefore not likely to give rise to measurable effects. More interesting and physically relevant is the behaviour of the function $`f`$ giving the precise dependence of the quality factor on frequency. First of all, it is easily seen that, for the special case $`h`$ = $`h^{}`$, $`f`$ is equal to 1, and thus the quality factor is proportional to $`\omega ^1`$, as was the case with toroidal modes. But when the aforementioned equality does not hold, numerical calculations are needed. Figure 1 shows that this function has a rather irregular dependence on frequency. We have represented $`f`$ for the first 18 eigenvalues of the spheroidal spectrum, where the inequality (11lno) safely holds. It is clear from equation (11lnyaaah) that holding fixed $`h`$ and $`kR`$ leaves us with a linear function of $`h^{}`$, whose slope varies from root to root.

## 4 Maxwell model

In this section we shall consider constituent equations given by the so called Maxwell model. Like the Kelvin-Voigt, these equations only involve first order time derivatives. As we shall see, quite different predictions will be obtained for the quality factor dependence on the mode frequencies.

The Maxwell model is also isotropic and homogenous, and is characterised by the following constituent equations :

$$_t\sigma _{ij}+\alpha \sigma _{kk}\delta _{ij}+\beta \sigma _{ij}=\frac{}{t}\left(\lambda s_{kk}\delta _{ij}+2\mu s_{ij}\right)$$
(11lnyaaal)
Here, the constants $`\lambda `$ and $`\mu `$ are again the Lamé coefficients describing the elastic behaviour of the body, while the constants $`\alpha `$ and $`\beta `$ parametrise the effects due to internal friction. To construct factorised solutions, we must also factorise both stress and strain333 Due to the equations of motion (2), if $`𝐬(𝐱,t)`$ is assumed to be separable in the fashion (8) then the strain tensor is separable too, and due to the constituent equation it is easily seen that the only possible time dependence is of the form $`\mathrm{exp}(\gamma t)`$.:

$$\sigma _{ij}(𝐱,t)=e^{\gamma t}\sigma _{ij}(𝐱),s_i(𝐱,t)=e^{\gamma t}s_i(𝐱),$$
(11lnyaaam)
The constituent equation is thus written, after separation of variables and contraction of its free indices, as

$$(\gamma +3\alpha +\beta )\sigma _{jj}(𝐱)=\gamma (2\mu +3\lambda )s_{jj}(𝐱)$$
(11lnyaaan)
and hence we have the following relationship between the spatial parts of stress and strain:

$$\left(1+\frac{\beta }{\gamma }\right)\sigma _{ij}(𝐱)=\left(\lambda \alpha \frac{2\mu +3\lambda }{\gamma +3\alpha +\beta }\right)s_{kk}(𝐱)\delta _{ij}+2\mu s_{ij}(𝐱)$$
(11lnyaaao)
Like in the Kelvin-Voigt model, we shall be mainly interested in the case of small internal friction, so we shall assume

$$\frac{\beta }{\gamma },\frac{\alpha }{\gamma }1$$
(11lnyaaap)
whence the following constituent relationship results:

$$\sigma _{ij}(𝐱)=\lambda \left(1\frac{\delta }{\gamma }\right)s_{kk}(𝐱)\delta _{ij}+2\mu \left(1\frac{\beta }{\gamma }\right)s_{ij}(𝐱)$$
(11lnyaaaq)
where we have introduced a new constant, $`\delta `$, defined by

$$\delta \frac{2\mu +3\lambda }{\lambda }\alpha +\beta ,$$
(11lnyaaar)
so that we can take as the parameters characterizing the Maxwell solid the set consisting of the Lamé coefficients $`\lambda `$, $`\mu `$, and the parameters $`\beta `$, $`\delta `$ which describe internal friction.

Let us compare equation (11lnyaaaq) with that of the Kelvin–Voigt model —cf. (6)— once the separation of variables has been performed:

$$\sigma _{ij}(𝐱)=(\lambda +\gamma \lambda ^{})s_{kk}(𝐱)\delta _{ij}+2(\mu +\gamma \mu ^{})s_{ij}(𝐱)$$
(11lnyaaas)
Comparing equations (11lnyaaaq) and (11lnyaaas), we observe that the solution of the Maxwell model can be carried out, as regards the spatial part of $`𝐬`$, following the same method used in the previous section for the Kelvin–Voigt model. In fact, we can directly take the expressions there derived, and make the substitutions

$$\mu ^{}\mu \beta \gamma ^2,\lambda ^{}\lambda \delta \gamma ^2$$
(11lnyaaat)
which transform equation (11lnyaaas) into (11lnyaaaq). Thus the form of the solutions and boundary conditions for a Maxwell viscoelastic sphere are those of B with constants $`𝒦`$ and $`𝒬`$ now given by the following functions of the parameter $`\gamma `$:

$`𝒬`$ $`=`$ $`i\sqrt{{\displaystyle \frac{\rho }{\lambda +2\mu }}}\left[\gamma +{\displaystyle \frac{1}{h+2}}\left({\displaystyle \frac{h}{2}}\delta +\beta \right)\right]`$ (11lnyaaaua)
$`𝒦`$ $`=`$ $`i\sqrt{{\displaystyle \frac{\rho }{\mu }}}\left(\gamma +{\displaystyle \frac{\beta }{2}}\right)`$ (11lnyaaaub)
where the approximation (11lnyaaap) has been taken into account, and $`h`$$``$$`\lambda /\mu `$.

Thus the two families of quasinormal modes of vibration are also present in this model, and we describe them in the following subsections.

### 4.1 Toroidal modes

As already discussed, the allowed values for $`\gamma `$ are again those making the linear system (11lnyaaaubabibkbmbncbcs) compatible, and there are two alternative ways to accomplish this. The first possibility yields purely tangential ($`C_t=C_l=0`$) vibrations satisfying once more the condition

$$\beta _1(𝒦R)=0𝒦=\sqrt{\frac{\rho }{\mu }}\omega _{nl}^T$$
(11lnyaaauav)
$`\omega _{nl}^T`$ being a toroidal eigenfrequency of the elastic sphere. Using the relationship between $`\gamma `$ and $`𝒦`$ for a Maxwell solid given by equations (11lnyaaaua)-(11lnyaaaub), we obtain the allowed values $`\gamma _{nl}^T`$ as

$$\gamma _{nl}^T=i\omega _{nl}^T\frac{\beta }{2}$$
(11lnyaaauaw)
Again, toroidal quasinormal modes have two fundamental properties: they have the same set of eigenfrequencies as the elastic sphere (to first order in the parameters describing internal friction, $`\beta `$ in this case), and also exactly the same spatial part (for all values of the viscosity parameters). The only difference between Kelvin–Voigt and Maxwell solids as regards toroidal modes appears in the dependence of the quality factor on $`\omega `$: as equation (11lnyaaauaw) shows, the quality factor in a linear Maxwell solid increases linearly with frequency. We can express all these properties by means of the following formulæ:

$$𝐬_M^T(𝐱,t)=𝐬_E^T(𝐱,t)e^{\omega _{nl}^Tt/Q},Q_{nl}^T=\frac{2\omega _{nl}^T}{\beta }$$
(11lnyaaauax)
relating Maxwell quasi-normal modes of vibration, $`𝐬_M^T(𝐱,t)`$, to elastic normal modes, $`𝐬_E^T(𝐱,t)`$, for the toroidal family.

### 4.2 Spheroidal modes

In order to handle the spheroidal family, we shall resort again to the perturbative expansions already used in the Kelvin–Voigt case, and also in the toroidal family, just described. The Maxwell model trivially reduces to the perfect elastic case when $`\beta `$ = $`\delta `$ = 0, hence we can take as the perturbative parameter

$$ϵ=\frac{\beta }{\omega }$$
(11lnyaaauay)
where $`\omega `$ is the elastic eigenfrequency to which $`\gamma `$ approaches when both $`\beta `$ and $`\delta `$ approach zero. Perturbative expansions in the fashion of section 3.2 can now be introduced:

$$\gamma =i\omega +\gamma _1ϵ,𝒦=k+k_1ϵ,𝒬=q+q_1ϵ$$
(11lnyaaauaz)
where the first order corrections $`k_1`$ and $`q_1`$ are given by equations (11lnyaaaua)-(11lnyaaaub) as functions of $`\gamma _1`$:

$`k_1`$ $`=`$ $`i\sqrt{{\displaystyle \frac{\rho }{\mu }}}\left(\gamma _1+{\displaystyle \frac{\omega }{2}}\right)`$ (11lnyaaaubaa)
$`q_1`$ $`=`$ $`i\sqrt{{\displaystyle \frac{\rho }{\lambda +2\mu }}}\left(\gamma _1+{\displaystyle \frac{h^{}+2}{h+2}}{\displaystyle \frac{\omega }{2}}\right)`$ (11lnyaaaubab)
The zero–order ratio $`h^{}`$ is now given by

$$h^{}=\frac{\alpha }{\beta }h$$
(11lnyaaaubabb)
With this definition, together with that of the perturbative parameter, the expressions at hand are formally identical to those of the Kelvin–Voigt model, and therefore the solutions to the Maxwell model share all their properties with their Kelvin–Voigt counterparts; the exception is the dependence of the quality factor on frequency: the product $`\gamma _1ϵ`$, which gives the exponential decay, is now independent of $`\omega `$.

Summing up, spheroidal quasinormal modes of the Maxwell solid, $`𝐬_M^P(𝐱,t)`$, are related to spheroidal normal modes of a perfectly elastic sphere by the equations

$$𝐬_M^P(𝐱,t)=𝐬_E^P(𝐱,t)e^{\omega _{nl}^Pt/Q+\chi (r)},Q_{nl}^P=\frac{2\omega _{nl}^P}{\beta }f(k_{nl}^PR,h,h^{})$$
(11lnyaaaubabc)
where the function $`f(kR,h,h^{})`$ is again given by (11lnyaaah).

As we see, the only difference between the behaviour of Maxwell and Kelvin–Voigt viscoelastic solids, when the internal friction effects can be considered small, appears in the dependence of $`Q`$ on frequency. We must however stress that, under other conditions (e.g. static load), both models show larger divergences in their physical properties .

## 5 Other Models

In this section we review other models which have been proposed to address the dynamics of a viscoelastic solid. They are generalisations of those in the two previous sections. We shall however not attempt to find complete solutions to all of them, as it eventually becomes too cumbersome. We shall however discuss in this section some of their most relevant traits.

### 5.1 The Standard Linear Model

The Standard Linear Model (SLM) for a viscoelastic solid is a generalised combination of the Kelvin–Voigt and Maxwell models. The constituent equations take here the form:

$$\sigma _{ij}+\frac{}{t}(\alpha \sigma _{kk}\delta _{ij}+2\beta \sigma _{ij})=\left(\lambda +\alpha ^{}\frac{}{t}\right)s_{kk}\delta _{ij}+2\left(\mu +\beta ^{}\frac{}{t}\right)s_{ij}$$
(11lnyaaaubabd)
where the effects of internal friction are described in this case with the aid of four constant parameters: $`\alpha `$, $`\beta `$, $`\alpha ^{}`$, and $`\beta ^{}`$. When looking for factorised solutions, the equations of motion and the above relationship force a time dependence of the form $`e^{\gamma t}`$ for both stress and strain. When such a dependence is introduced in equation (11lnyaaaubabd), we obtain the following relationship between the spatial parts of the stress and strain tensors:

$`(1+2\gamma \beta )\sigma _{ij}(𝐱)`$ $`=`$ $`\left[\lambda +\alpha ^{}\gamma \alpha \gamma {\displaystyle \frac{3\lambda +2\mu +\gamma (3\alpha ^{}+2\beta )^{}}{1+\gamma (3\alpha +2\beta )}}\right]s_{kk}(𝐱)\delta _{ij}`$ (11lnyaaaubabe)
$`+`$ $`2(\mu +\beta ^{}\gamma )s_{ij}(𝐱)`$
The case of small internal friction is treated by first order approximation in the quantities parametrising viscous processes, i.e.,

$$\alpha ,\beta ,\alpha ^{},\beta ^{}\gamma ^1$$
(11lnyaaaubabf)
When such approximation is made, the above equation reduces to

$$\sigma _{ij}(𝐱)=(\lambda +\lambda ^{}\gamma )s_{kk}(𝐱)\delta _{ij}+2(\mu +\mu ^{}\gamma )s_{ij}(𝐱)$$
(11lnyaaaubabg)
where we have introduced two new constants given by

$$\mu ^{}\beta ^{}2\beta \mu ,\lambda ^{}\alpha ^{}2\lambda \beta (3\lambda +2\mu )\alpha $$
(11lnyaaaubabh)
Therefore when equation (11lnyaaaubabf) holds, the SLM reduces to a Kelvin–Voigt model, i.e., for small internal friction, both models have the same set of quasinormal modes of vibration, which are characterized by the constants $`\lambda `$, $`\mu `$, $`\lambda ^{}`$ and $`\mu ^{}`$, the latter being given, for the Standard Linear solid, by equations (11lnyaaaubabh).

### 5.2 Generalised mechanical models

The models analyzed so far are the simplest ones obtained by three dimensional generalisations of mechanical viscoelastic models composed of linear springs. They give rise to differential constituent relations, with time derivatives up to the first order. Considering more involved generalisations yields differential relations involving higher order time derivatives of strain and stress —see e.g. . Thus, quite independently of any reference to the underlying mechanical model, we can consider general differential relations between stress and strain including any number of time derivatives. To ease the formulation of such differential constituent equations for the case of isotropic and homogenous solids, we shall introduce the trace-free parts of the strain and stress tensors, $`s_{ij}^{}`$, $`\sigma _{ij}^{}`$ (usually termed deviatoric components in the literature on viscoelasticy ), and their traces, $`s`$ and $`\sigma `$ (dilational components), defined by

$`s_{ij}^{}`$ $`=`$ $`s_{ij}{\displaystyle \frac{1}{3}}s\delta _{ij},ss_{kk}`$ (11lnyaaaubabia)
$`\sigma _{ij}^{}`$ $`=`$ $`s_{ij}{\displaystyle \frac{1}{3}}\sigma \delta _{ij},\sigma \sigma _{kk}.`$ (11lnyaaaubabib)
In terms of the above quantities, the linear Hooke law for an elastic solid takes the form

$$\sigma =(3\lambda +2\mu )s\sigma _{ij}^{}=2\mu s_{ij}^{}$$
(11lnyaaaubabibj)
while the constituent equation of an SLM is written

$`\left[1+(3\alpha +2\beta ){\displaystyle \frac{}{t}}\right]\sigma =\left[(3\lambda +2\mu )+(3\alpha ^{}+2\beta ^{}){\displaystyle \frac{}{t}}\right]s`$ (11lnyaaaubabibka)
$`\left(1+2\beta {\displaystyle \frac{}{t}}\right)\sigma _{ij}^{}=\left(2\mu +2\beta ^{}{\displaystyle \frac{}{t}}\right)s_{ij}^{}`$ (11lnyaaaubabibkb)
This equation can now be generalised to include higher order time derivatives. We can thus consider viscoelastic models whose constituent equation is given by

$$R(_t)\sigma =S(_t)s,R^{}(_t)\sigma _{ij}^{}=S^{}(_t)s_{ij}^{}$$
(11lnyaaaubabibkbl)
where we have introduced formal polynomials

$`R(x){\displaystyle \underset{l=0}{\overset{N1}{}}}r_lx^l,R^{}(x){\displaystyle \underset{l=0}{\overset{N1}{}}}r_l^{}x^l`$ (11lnyaaaubabibkbma)
$`S(x){\displaystyle \underset{l=0}{\overset{N1}{}}}s_lx^l,S^{}(x){\displaystyle \underset{l=0}{\overset{N1}{}}}s_l^{}x^l`$ (11lnyaaaubabibkbmb)
so that a general differential model is given for each set of 4$`N`$ real constants $`r_l`$, $`r_l^{}`$, $`s_l`$, and $`s_l^{}`$ characterising the solid. Some of these constants may vanish. The Kelvin-Voigt, Maxwell and SL models considered above are of course special cases within this general class. Several procedures have been proposed in the literature to solve the general equations —see  for a review and further reference—, which can be applied to the solid viscoelastic sphere problem. We now sketch how they work in this case of our interest.

Let $`\stackrel{~}{\sigma }_{ij}(𝐱,\mathrm{\Omega })`$ and $`\stackrel{~}{s}_{ij}(𝐱,\mathrm{\Omega })`$ be the Fourier transforms of the stress and strain tensors, and $`\stackrel{~}{s}_i(𝐱,\mathrm{\Omega })`$ that of the displacement vector field:

$`\stackrel{~}{\sigma }_{ij}(𝐱,\mathrm{\Omega })`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\sigma _{ij}(𝐱,t)e^{i\mathrm{\Omega }t}𝑑t`$ (11lnyaaaubabibkbmbna)
$`\stackrel{~}{s}_i(𝐱,\mathrm{\Omega })`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}s_i(𝐱,t)e^{i\mathrm{\Omega }t}𝑑t`$ (11lnyaaaubabibkbmbnb)
$`\stackrel{~}{s}_{ij}(𝐱,\mathrm{\Omega })`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}s_{ij}(𝐱,t)e^{i\mathrm{\Omega }t}𝑑t`$ (11lnyaaaubabibkbmbnc)
In terms of these, Fourier transforms the constituent equations (11lnyaaaubabibkbl) read:

$$\stackrel{~}{\sigma }=\frac{R(i\mathrm{\Omega })}{S(i\mathrm{\Omega })}\stackrel{~}{s},\stackrel{~}{\sigma }_{ij}^{}=\frac{R^{}(i\mathrm{\Omega })}{S^{}(i\mathrm{\Omega })}\stackrel{~}{s}_{ij}$$
(11lnyaaaubabibkbmbnbo)
while the equations of motion are

$$\mathrm{\Omega }^2\stackrel{~}{𝐬}=\frac{1}{3}\left(\frac{R(i\mathrm{\Omega })}{S(i\mathrm{\Omega })}\frac{R^{}(i\mathrm{\Omega })}{S^{}(i\mathrm{\Omega })}\right)(\stackrel{~}{𝐬})+\frac{R^{}(i\mathrm{\Omega })}{2S^{}(i\mathrm{\Omega })}^2\stackrel{~}{𝐬}$$
(11lnyaaaubabibkbmbnbp)
Comparing the above equations with the corresponding ones for normal modes of vibration of elastic solids, and the constituent relation (11lnyaaaubabibkbmbnbo) with 11lnyaaaubabibj, we note that the problem of finding solutions to the equation of motion of a general viscoelastic differential model reduces to that of finding the normal modes of vibration of an elastic solid having complex Lamé coefficients given by

$$\stackrel{~}{\lambda }(\mathrm{\Omega })=\frac{1}{3}\left(\frac{R(i\mathrm{\Omega })}{S(i\mathrm{\Omega })}\frac{R^{}(i\mathrm{\Omega })}{S^{}(i\mathrm{\Omega })}\right),\stackrel{~}{\mu }(\mathrm{\Omega })=\frac{1}{2}\frac{R^{}(i\mathrm{\Omega })}{S^{}(i\mathrm{\Omega })}$$
(11lnyaaaubabibkbmbnbq)
where the allowed values of $`\mathrm{\Omega }`$ are obtained as the solutions to the elastic solid’s eigenfrequency equation when the above complex coefficients are used instead of the real, constant Lamé coefficients $`\lambda `$, $`\mu `$. Generally, $`\mathrm{\Omega }`$ will have complex values, thus giving rise to damped system oscillations. After solving for $`\mathrm{\Omega }`$, the spatial part of the solutions is obtained from that of the normal modes by simply substituting the old, real–valued constants $`\omega `$, $`\lambda `$ and $`\mu `$ by the new complex values $`\mathrm{\Omega }`$, $`\stackrel{~}{\lambda }`$ and $`\stackrel{~}{\mu }`$. This method for solving the viscoelasticity is often termed in the literature on the subject the Correspondence Principle , and as a matter of fact our previous derivations of the form of the quasinormal modes for Kelvin-Voigt, Maxwell and SL models can be seen to be special cases of its application. The method is applicable to any boundary value problem whose elastic counterpart is solvable. The case of small internal friction (i.e., first order approximation in the coefficients of the polynomials (11lnyaaaubabibkbma)-(11lnyaaaubabibkbmb) has been considered by Graffi  for one dimensional wave propagation.

The three dimensional spherical case is also solvable, as we know. The toroidal modes are relatively straightforward to obtain from their elastic counterparts due to the simple form of their eigenvalue equation, while the spheroidal ones demand more complex algebra, which becomes increasingly cumbersome as the order $`N`$ of the model increases. We shall present here the general solution for the toroidal modes for any differential viscoelastic model, whereby we shall obtain the dependence of their $`Q`$ on frequency. This will also be the approximate dependence for the spheroidal modes, if friction effects are small, as was the case with the first order models analysed so far. A complete solution for the latter modes can also be systematically found, but will be omitted due to its scarcely useful algebraic complexity .

#### 5.2.1 Toroidal modes

As discussed above, the boundary equation for the toroidal modes in a general viscoelastic model is obtained from the eigenvalue equation of the elastic model:

$$\beta _1(kR)=0,k=\sqrt{\frac{\rho }{\mu }}\omega $$
(11lnyaaaubabibkbmbnbr)
Upon substitution of $`\mu `$ by $`\stackrel{~}{\mu }`$, we obtain

$$\beta _1(𝒦R)=0,𝒦=\sqrt{\frac{2\rho S^{}(i\mathrm{\Omega })}{R^{}(i\mathrm{\Omega })}}\mathrm{\Omega }$$
(11lnyaaaubabibkbmbnbs)
We know that the only solutions to the eigenvalue equation 11lnyaaaubabibkbmbnbr are the real eigenfrequencies of the elastic sphere $`\omega _{nl}^T`$, and therefore the allowed values for $`\mathrm{\Omega }`$ are given by the implicit relationship

$$\sqrt{\frac{2\mu S^{}(i\mathrm{\Omega })}{R^{}(i\mathrm{\Omega })}}\mathrm{\Omega }=\omega _{nl}^T$$
(11lnyaaaubabibkbmbnbt)
Let us now write the polynomials $`S^{}`$ and $`R^{}`$ in the form

$$S^{}(x)=1+ϵ\underset{l=1}{\overset{N1}{}}s_l^{}x^l,R^{}(x)=2\mu \left(1+ϵ\underset{l=1}{\overset{N1}{}}r_l^{}x^l\right)$$
(11lnyaaaubabibkbmbnbu)
as a suitable one to imply that internal friction effects are small, letting

$$ϵ1$$
(11lnyaaaubabibkbmbnbv)
The quantities $`r_l^{}\omega ^l`$ and $`s_l^{}\omega ^l`$ are thus zero order in $`ϵ`$ and dimensionless, $`\omega `$ being a toroidal eigenvalue of the elastic case. We then introduce an expansion for $`\mathrm{\Omega }`$ in the small parameter $`ϵ`$, whose zeroth order term corresponds to a given toroidal eigenfrequency $`\omega `$ of the elastic solid:

$$\mathrm{\Omega }=\omega +ϵ\mathrm{\Omega }_1$$
(11lnyaaaubabibkbmbnbw)
Under the above conditions, we have

$$\frac{2S^{}(i\mathrm{\Omega })}{\mu R^{}(i\mathrm{\Omega })}=1+2iϵ\underset{l=0}{\overset{N1}{}}t_l\omega ^l,t_l=i^{l1}(s_l^{}r_l^{})/2$$
(11lnyaaaubabibkbmbnbx)
and the value of $`\mathrm{\Omega }_1`$ follows when introducing the above expansion into equation (11lnyaaaubabibkbmbnbt), yielding

$$\mathrm{\Omega }_1=\underset{l=1}{\overset{N1}{}}t_l\omega ^{l+1}$$
(11lnyaaaubabibkbmbnby)
In terms of $`\mathrm{\Omega }_1`$, the quality factor reads

$$Q=\frac{\omega }{ϵ}\frac{1}{\text{Im}[\mathrm{\Omega }_1]}$$
(11lnyaaaubabibkbmbnbz)
where Im\[$``$\] denotes the imaginary part of its argument. Thus we observe that, as regards toroidal modes, using a general differential model gives us a polynomial in $`\omega `$ for $`1/Q`$, with no independent term, so that constant $`Q`$ is not allowed by these models. The polynomial only contains odd powers of the unperturbed frequency $`\omega `$. In general, whenever $`t_l`$$``$ 0 for even $`l`$, the real part of $`\mathrm{\Omega }_1`$ will not vanish, and the angular frequency of the periodic component of the quasinormal modes shall undergo first order corrections. Hence, in order to preserve the elastic spectrum to first order, our model must satisfy the conditions

$$t_l=0\text{(}l\text{ even)}$$
(11lnyaaaubabibkbmbnca)
Provided the preceeding equation holds, the corrections to $`𝒦`$ will be purely imaginary, and therfore the modulus of the spatial part of the modes will remain unaltered, the only effect of viscosity being the addition of a point dependent phase in the fashion of equation (11lnyaaaubabc)444 This correction will only appear provided that $`N`$$``$ 3. This is why it was absent in the toroidal families of the previously discussed models..

The calculation for the spheroidal quasinormal modes can be performed along the same lines but, as we have seen, the algebra is considerably more involved already in the simplest models. It does naturally become more cumbersome as the order of the model increases, so we omit a detailed discussion of its technicalities here.

## 6 Experimental data

We now present a confrontation of the above theoretical models with available experimental measurements, some of them published , and others unpublished as of this date .

It is very important, in order to rightly assess the suitability of the theoretical models, to keep in mind that measurements of $`Q`$ are intrinsically difficult. The most important source of problems is the fact that almost any experimental environment conditions (suspensions, readout, electronics, etc.) significantly couple to the primary oscillator (the sphere itself), and almost invariably result in a *degradation* of the measured $`Q`$. The way and intensity in which this happens is however utterly unpredictable. As a rule, theoretical $`Q`$ values should thus be expected *larger* than measured ones.

Then, the mode frequencies slightly split up as an unavoidable consequence of the necessity to *suspend* the sample in the laboratory —a fact not considered in the above models. This, at times, generates a forest of frequency peaks when different multipole modes happen to be close to one another. For example, within a mere 0.8 % of their nominal value, we find the three fundamental frequencies: $`\omega _{22}^P`$, $`\omega _{14}^P`$, and $`\omega _{14}^T`$ which, if split up, result in as many as 23 Fourier peaks in a very reduced frequency interval. On the other hand, however, the $`Q`$’s of the various multiplet components often group around clearly distinct values, close to those predicted by theory. In such cases, therefore, the viscoelastic model is helpful to tell the different multiplet members from one another, an otherwise very difficult task.

Rather than characterising a given viscoelastic model by the two parameters $`\lambda ^{}`$ and $`\mu ^{}`$, or $`\alpha `$ and $`\beta `$ —see sections 3 and 4—, it is expedient for the purposes of this section to use a different parameter pair. This will be chosen as the $`Q`$ of the lowest frequency toroidal mode, $`Q_0Q_{12}^T`$, and the viscoelastic ratio $`h^{}`$ of sections 3.2 and 4.2. The following expressions thus hold:

Kelvin–Voigt: $`\{\begin{array}{cc}Q_{nl}^T=\left(\frac{\omega _{12}^T}{\omega _{nl}^T}\right)Q_0\hfill & \\ Q_{nl}^P=\left(\frac{\omega _{12}^T}{\omega _{nl}^P}\right)f(k_{nl}^PR,h,h^{})Q_0\hfill & \end{array}`$ (11lnyaaaubabibkbmbncbc)
Maxwell: $`\{\begin{array}{cc}Q_{nl}^T=\left(\frac{\omega _{nl}^T}{\omega _{12}^T}\right)Q_0\hfill & \\ Q_{nl}^P=\left(\frac{\omega _{nl}^P}{\omega _{12}^T}\right)f(k_{nl}^PR,h,h^{})Q_0\hfill & \end{array}`$ (11lnyaaaubabibkbmbncbf)
obviously a rewrite of the corresponding ones in the respective sections.

In figure 2 we plot the measured versus theoretically calculated values of $`Q`$ for two sample spheres of Al 5056 and a 7 % alloy of CuAl —see more technical data in the caption to the figure. Agreement is quite good, when a Maxwell model is used, except for the higher frequency mode, reported in the case of CuAl. It appears that a better fit can be accomplished when a combination of both Maxwell and Kelvin–Voigt models is used, such as we see in figure 3, where $`Q`$’s of a sample sphere of steel are considered. The latter plot contains more experimental points than theoretical, and we see here a clear example of frequency splittings due to symmetry breaking caused by suspension, normally a diametral or semi-diametral bore drilled across the sphere .

## 7 Conclusions

In this paper we have addressed the problem of whether it is possible to systematically characterise the linewidths or, equivalently, the $`Q`$’s of the oscillation eigenmodes of a given spherical GW detector. To this end we have considered various phenomenological models, selected from the specialised literature on the subject, and solved the equations of motion in the cases of our interest. Different models are seen to predict different frequency dependences of the quality factors for the lower modes, which are the ones we have paid attention to, and the ones relevant for GW detection purposes. For example, in a Kelvin–Voigt solid the $`Q`$ of a given mode appears to be inversely proportional to its frequency, while in a Maxwell solid it is directly proportional to it.

As we have seen, however, the behaviour of a given material is generally not well described by a *single* such model but rather by a suitable combination of e.g. Maxwell and Kelvin–Voigt models, depending on the range of frequencies considered. This is not too surprising, as these models are *phenomenological*, i.e., they are not based on a detailed analysis of the *physics* of the various dissipation mechanisms which are responsible for the loss processes , but on more or less plausible (though not unique) generalisations of the “friction is proportional to velocity” principle, which so accurately holds in a simple unidimensional harmonic oscillator.

In spite of these limitations, and also in spite of the experimental difficulties inherent in the determination of $`Q`$, it appears that our treatment of the problem gives results which can be considered quite good, in view of the above limitations. A significant bonus of the present analysis is the possibility it offers to identify nearby frequencies in the sphere spectrum by their distinct value of $`Q`$, as predicted by theory. This is useful in GW detection science since, for example, the *second* spheroidal quadrupole frequency, $`\omega _{22}^P`$, is very important due to large cross section for GW energy absorption in that mode  yet it is very close (in frequency) to other modes with different $`Q`$’s.

Improvements on the results presented in this paper are certainly possible. But we feel they would have to be based on a different methodology, with more emphasis *ab initio* on Solid State lore, and/or on Materials Science.

We express our thanks to J.A. Ortega, who contributed most of the analytic work reported in this paper. Three of us (JAL, GF and AdW) acknowledg hospitality in the INFN Frascati Laboratories, and support from the TARI programme, contract number HPRI-CT-1999-00088.

## Appendix A

A well-behaved three-dimensional vector field $`𝐬(𝐱)`$ can be expressed as the sum of an irrotational, $`𝐬_l(𝐱)`$, and a divergence free, $`𝐬_t(𝐱)`$, vector fields, respectively called the longitudinal and transverse components of $`𝐬(𝐱)`$ :

$$𝐬(𝐱)=𝐬_l(𝐱)+𝐬_t(𝐱),𝐬_t=\times 𝐬_l=0$$
(11lnyaaaubabibkbmbncbcc)
We now replace this decomposition into equation (9) to find

$`\rho \ddot{T}(t)(𝐬_t+𝐬_l)`$ $`=`$ $`\left[\mu T(t)+\mu ^{}\dot{T}(t)\right]^2(𝐬_t+𝐬_l)`$ (11lnyaaaubabibkbmbncbcd)
$`+`$ $`\left[(\lambda +\mu )T(t)+(\lambda ^{}+\mu ^{})\dot{T}(t)\right]\left(𝐬_l\right)`$
Taking the rotational of this equation,

$$\times \left[\rho \ddot{T}𝐬_t(\mu T+\mu ^{}\dot{T})^2𝐬_t\right]=0$$
(11lnyaaaubabibkbmbncbce)
The vector between square brackets is thus divergence–free and irrotational, so it vanishes. We have therefore:

$$^2𝐬_t=\left\{\frac{\rho \ddot{T}}{\mu T+\mu ^{}\dot{T}}\right\}𝐬_t$$
(11lnyaaaubabibkbmbncbcf)
Since the left hand side of the above equation does not depend on time, the term between braces in the right hand side must equal a (complex) constant, say $`𝒦^2`$. Thus,

$`^2𝐬_t+𝒦^2𝐬_t=0`$ (11lnyaaaubabibkbmbncbcg)
$`\mu T+\mu ^{}\dot{T}+𝒦^2\ddot{T}=0`$ (11lnyaaaubabibkbmbncbch)
An analogous procedure, after taking the divergence of equation (11lnyaaaubabibkbmbncbcd), gives us the corresponding formulæ for the longitudinal part:

$`^2𝐬_l+𝒬^2𝐬_l=0`$ (11lnyaaaubabibkbmbncbci)
$`(\lambda +2\mu )T+(\lambda ^{}+2\mu ^{})\dot{T}+𝒬^2\ddot{T}=0`$ (11lnyaaaubabibkbmbncbcj)
where $`𝒬^2`$ stands for another complex separation constant.

## Appendix B

We describe in this appendix the algebraic operations which lead to the solution to the eigenvalue problem in a viscoelastic sphere. Equations (2) ought to be solved, subject to the boundary conditions (5). The latter can be cast in explicit vector form:

$$(\lambda +\gamma \lambda ^{})\left[𝐬(𝐱)\right]𝐧+2(\mu +\gamma \mu ^{})(𝐧)𝐬(𝐱)+2(\mu +\gamma \mu ^{})𝐧\times \left[\times 𝐬(𝐱)\right]=0$$
(11lnyaaaubabibkbmbncbck)
The irrotational and divergence free components, $`𝐬_t(𝐱)`$ and $`𝐬_l(𝐱)`$, can be expressed by means of auxiliary functions $`\varphi (𝐱)`$ and $`\psi (𝐱)`$:

$`𝐬_l(𝐱)`$ $`=`$ $`𝒬^1C_0\varphi (𝐱)`$ (11lnyaaaubabibkbmbncbcl)
$`𝐬_t(𝐱)`$ $`=`$ $`i𝒦^1C_1\times 𝐋\psi (𝐱)+iC_2𝐋\psi (𝐱)`$ (11lnyaaaubabibkbmbncbcm)
where $`C_0`$, $`C_1`$ and $`C_2`$ are (so far) undetermined integration constants, and L$``$$`i`$x$`\times `$$``$ is the “angular momentum” operator. Upon substitution of (11lnyaaaubabibkbmbncbcl)-(11lnyaaaubabibkbmbncbcm) into (2) it is readily seen that the functions $`\varphi `$ and $`\psi `$ are themselves also solutions to corresponding Helmholtz equations,

$$^2\varphi (𝐱)+𝒬^2\varphi (𝐱)=0,\mathrm{and}^2\psi (𝐱)+𝒦^2\psi (𝐱)=0$$
(11lnyaaaubabibkbmbncbcn)
They have therefore the general form, using spherical coordinates ($`r`$,$`\theta `$,$`\phi `$) for the vector x,

$$\varphi (𝐱)=j_l(𝒬r)Y_{lm}(\theta ,\phi ),\mathrm{and}\psi (𝐱)=j_l(𝒦r)Y_{lm}(\theta ,\phi )$$
(11lnyaaaubabibkbmbncbco)
where $`j_l`$ are spherical Bessel functions of the first kind and $`Y_{lm}`$ are spherical harmonics. The solutions (11lnyaaaubabibkbmbncbco) are those possessing regularity properties in the whole interior and boundary of the solid. We thus have

$`\varphi `$ $`=`$ $`{\displaystyle \frac{dj_l(𝒬r)}{dr}}Y_{lm}(\theta ,\phi )𝐧{\displaystyle \frac{j_l(𝒬r)}{r}}i𝐧\times 𝐋Y_{lm}(\theta ,\phi )`$ (11lnyaaaubabibkbmbncbcp)
$`\times 𝐋\psi `$ $`=`$ $`l(l+1){\displaystyle \frac{j_l(𝒦r)}{r}}Y_{lm}(\theta ,\phi )𝐧`$ (11lnyaaaubabibkbmbncbcq)
$`+`$ $`\left[{\displaystyle \frac{j_l(𝒦r)}{r}}+{\displaystyle \frac{d}{dr}}j_l(𝒦r)\right]i𝐧\times 𝐋Y_{lm}(\theta ,\phi )`$
$`𝐋\psi `$ $`=`$ $`j_l(𝒦r)i𝐋Y_{lm}(\theta ,\phi )`$ (11lnyaaaubabibkbmbncbcr)
These expressions ought to be substituted now into (11lnyaaaubabibkbmbncbcl)-(11lnyaaaubabibkbmbncbcm), and then into (11lnyaaaubabibkbmbncbck) —recall that s = $`𝐬_t`$ \+ $`𝐬_l`$. It is found thet these are equivalent to the following homogeneous linear system

$$\left(\begin{array}{ccc}\beta _4(𝒬R,\frac{\lambda +\gamma \lambda ^{}}{\mu +\gamma \mu ^{}})& l(l+1)\frac{𝒦}{𝒬}\beta _1(𝒦R)& 0\\ \beta _1(𝒬R)& \frac{𝒦}{𝒬}\beta _3(𝒦R)& 0\\ 0& 0& \frac{𝒦}{𝒬}𝒦R\beta _1(𝒦R)\end{array}\right)\left(\begin{array}{c}C_0\\ C_1\\ C_2\end{array}\right)=0$$
(11lnyaaaubabibkbmbncbcs)
where

$`\beta _0(z)`$ $``$ $`{\displaystyle \frac{j_l(z)}{z^2}}`$ (11lnyaaaubabibkbmbncbct)
$`\beta _1(z)`$ $``$ $`{\displaystyle \frac{d}{dz}}\left[{\displaystyle \frac{j_l(z)}{z}}\right]`$ (11lnyaaaubabibkbmbncbcu)
$`\beta _2(z)`$ $``$ $`{\displaystyle \frac{d^2}{dz^2}}\left[j_l(z)\right]`$ (11lnyaaaubabibkbmbncbcv)
$`\beta _3(z)`$ $``$ $`\frac{1}{2}\beta _2(z)+\left\{\frac{l\left(l+1\right)}{2}1\right\}\beta _0(z)`$ (11lnyaaaubabibkbmbncbcw)
$`\beta _4(z,A)`$ $``$ $`\beta _2(z){\displaystyle \frac{A}{2}}j_l(z)`$ (11lnyaaaubabibkbmbncbcx)
The system (11lnyaaaubabibkbmbncbcs) is to be satisfied by the constants $`C_0`$, $`C_1`$ and $`C_2`$, but has no meaningful solution unless the system matrix is singular, i.e., if its determinant vanishes. Therefore

$$\beta _1(𝒦R)det\left(\begin{array}{cc}\beta _4(𝒬R,\frac{\lambda +\gamma \lambda ^{}}{\mu +\gamma \mu ^{}})& l(l+1)\frac{𝒦}{𝒬}\beta _1(𝒦R)\\ \beta _1(𝒬R)& \frac{𝒦}{𝒬}\beta _3(𝒦R)\end{array}\right)=0$$
(11lnyaaaubabibkbmbncbcy)
This is an equation for the parameter $`\gamma `$, on which $`𝒦`$ and $`𝒬`$ depend through equations (11lna) and (11lnb). Clearly, there are two families of solutions, or eigenmodes, to (11lnyaaaubabibkbmbncbcy) associated to the vanishing of either of the two factors in its lhs, i.e., $`\beta _1(𝒦R)`$ or the determinant of the displayed 2$`\times `$2 matrix. They are called toroidal and spheroidal solutions, respectively.

## References


text:
# Formation of magnetic moments in the cuprate superconductor Hg0.8Cu0.2Ba2Ca2Cu3O8+δ below 𝑇_𝑐 seen by NQR

## 1 Introduction

A possible co-existence of antiferromagnetism and superconductivity remains one of the most interesting problem in high-$`T_c`$ cuprates. Despite that in several theoretical scenarios of high-$`T_c`$ superconductivity in cuprates the Cooper-pairing arises due to an exchange of antiferromagnetic spin fluctuations present in the paramagnetic phase , it is commonly believed that bulk superconductivity and the antiferromagnetism do not coexist in the cuprates at the same doping concentration and temperature. However, recently possible formation of antiferromagnetism below the superconducting transition temperature was found by several experimental techniques in underdoped YBa2Cu3Ox and La2-xSrxCuO4 . The relatively small values of the observed magnetic moments  ($`0.01\mu _B`$, where $`\mu _B`$ is the Bohr magneton) have indicated an orbital rather than an electronic origin of the observed antiferromagnetism. Interestingly, some of these findings were interpreted in terms of the time-reversal symmetry breaking state or $`id`$-density wave state  that may occur in the underdoped cuprates below the pseudogap formation temperature, $`T^{}`$. Most recently, the time-reversal symmetry breaking below $`T^{}`$ has been indeed observed in underdoped Bi2Sr2CaCu2O8+δ by means of the angle-resolved-photoemission spectroscopy (ARPES) .

Nuclear magnetic resonance (NMR) and nuclear quadrupole resonance (NQR) can be one of the most powerful tools to study the possible occurrence of the intrinsic antiferromagnetism in high-$`T_c`$ cuprate superconductors. Very recently the maximum in the transverse relaxation rate of Cu(2) nuclei observed in YBa2Cu3O7 below $`T_c`$ at T = 35 K was interpreted in terms of critical fluctuations of the $`id`$-density wave state . A similar conclusion was drawn from the analysis of the spin-freezing processes in underdoped Y2-xCaxBa2Cu3O6 and La2-xSrxCuO4. On the other hand, the NQR study provides a more direct prove for the formation of orbital magnetism since it is performed in zero magnetic field. Thus, the internal magnetic moments if they are present will result in the NQR line splitting. So far such a splitting was not observed in high-$`T_c`$ cuprates. Instead, one mainly observes an anomalous broadened Cu NQR spectra at low doping . Associating the broad Cu NQR spectra with doping inhomogeneities is leading to a few contradictions. The spectral width should depend on the preparation conditions and should have strongly decreased with improved preparation techniques developed in the past decade. Especially for the case of Hg-cuprates, that possess extremely broad NQR spectra, a dependency of the width of the NQR spectra on preparation conditions has never been reported. The second and more relevant reason against a doped inhomogeneity as the origin of the broad Cu NQR spectra is that it will spread the superconducting transition. An inhomogeneously distributed electric field gradient (EFG) is directly associated to an inhomogeneously distributed charge density in the CuO layers. The charge density again is associated to the superconducting transition temperature. Thus, a distribution of spheres with distinct EFGs inside a crystal has to be associated to a distribution of $`T_c`$’s. If one estimates the distributions of $`T_c`$’s of HgBaCuO from the phenomenological model by Gippius et al. , the superconducting transition should be broadened by several tenth of K, what is clearly not the case.

The most interesting reason for additional NQR lines is the removing of the Kramers degeneracy of the Cu quadrupole levels due to a formation of magnetic fields at the Cu sites. Such magnetic fields may arise from an antiferromagnetic state formatting spontaneously below $`T_c`$.

Additional lines may also occur from impurities and Cu outside the CuO layers. However, lines from impurities should be visible at any temperature. Furthermore, the relaxation times of Cu inside the CuO layers are extremely short, in comparison with metallic Cu, and driven by antiferromagnetic fluctuations (AF)  inside the CuO layers, strongly decreasing with distance. If a line from an impurity becomes visible at a certain temperature it can be easily identified by their distinct relaxation times with respect to Cu inside CuO layers.

In this work using NQR technique we present the first experimental evidence for the formation of the internal magnetic moments in the optimally doped three-CuO2-layer Hg0.8Cu0.2Ba2Ca2Cu3O8+δ (Hg-1223) and Hg0.8Cu0.2Ba2Ca2Cu3O8+δ$`F`$ high-temperature cuprate superconductors (HTSC) below $`T_c`$ ($`T_c`$ = 134 K). The synthesis of the samples is described in Ref. . An impurity content of less than 5 % within the sample was found from neutron powder diffraction and x-ray powder diffraction . The real part of the ac-susceptibility is shown in Fig 1. Both samples possess different ways of preparation and types of dopands. Thus, any common properties have to attributed to intrinsic properties of the superconducting state. In particular, we show that 63Cu NQR-lines split below $`T_c`$ due to a Zeeman splitting originating from the internal magnetic fields within the CuO2-layers.

Our measurements were carried out on a home-built NQR spectrometer. The spectra have been recorded by point-by-point technique with a Hahn spin-echo pulse sequence with pulse lengths of 4 $`\mu `$s and 8 $`\mu `$s. The recovery of the longitudinal magnetization was measured by the method of saturation recovery. Various numbers of saturating pulses (from 3 to 101) have been used to study the dependency of the spin-lattice relaxation on the initial conditions. The decay of the transverse magnetization was measured by varying the pulse separation ($`\tau `$) of the spin-echo sequence. In order to detect lines with strongly different $`T_2`$ times the spectra were scanned with first increasing pulse separation time $`\tau _2`$ of the Hahn echo sequence. To detect lines with strongly different $`T_1`$ times the point to be measured was saturating by a saturating comp and after a time $`\tau _1`$ measured by the Hahn spin-echo sequence.

## 2 Zero field NQR

Nuclei with a spin $`I1`$ possess a quadrupole moment and interacting with an EFG. The EFG is described by the traceless tensor $`V_{ij}=\frac{\delta ^2V}{\delta x_i\delta x_j}`$. Since the EFG tensor is a traceless tensor, only the values $`V_{zz}`$ and $`V_{xx}V_{yy}`$ are needed to characterize the EFG tensor. The EFG tensor is commonly described by the values $`eq`$ and $`\eta `$ defined by the equations $`eq=V_{zz}`$ and $`\eta =\frac{V_{xx}V_{yy}}{V_{zz}}`$. In high-$`T_c`$ superconductors the elongated CuO octahedrons leading to the case of axial symmetry of the EFG at the Cu site, giving $`\eta =0`$.

With the restriction to the case of axial symmetry of the EFG tensor, the interaction energy is defined by the well known Hamiltonian

$$H=\frac{e^2qQ}{4I(2I1)}(3I_z^2I^2)$$
(1)
with eQ the quadrupole moment of the nucleus. The Eigenvalues of this unperturbed Hamiltonian are

$$E=\frac{e^2qQ}{4I(2I1)}(3m^2I(I+1).$$
(2)
In general, there is a twofold degeneracy of the $`\pm m`$ eigenvalues corresponding to the symmetry of the nucleus charge distribution and of the EFG. For half integer spins a deviation from the EFGs rotational symmetry, $`\eta 0`$, will not remove the degeneracy corresponding to Kramers theorem that states that the degeneracy of half integer spins can only completely be removed by a magnetic field.

For the case of nuclei like Cu with spin $`I=3/2`$ the set of levels form a virtual two level system due to the Kramers degeneracy and only one transition will be observed. Note, that for this particular case, an observed Cu NQR line splitting can only be due to mesoscopic charge segregation and/or the appearance of a magnetic field at the nucleus site.

### 2.1 Zeeman perturbed NQR

The twofold degeneracy of the quadrupole level is corresponding to the fact, that turning the nuclei end for end will leave the electrostatic energy unchanged. Applying a magnetic field at the nuclei will lower the symmetry of the system and lifting the degeneracy. In general the situation is described by adding a Zeeman term $`H_z=\gamma _nB_0hI_z`$ to the unperturbed Hamiltonian (1)

$$H=\frac{e^2qQ}{4I(2I1)}(3I_z^2I^2)+\gamma _nB_0hI_z$$
(3)
and treating the problem in first order perturbation theory. The solutions can be divided into three cases: The magnetic field $`B_0`$ oriented in parallel, second perpendicular to the symmetry axis of the EFG, and third the magnetic field oriented with angles $`0^{},90^{}`$ with respect to the symmetry axis of the EFG. In the cases $`B_0`$ parallel to the symmetry axis, the mixing of the $`|\pm 1/2`$ level vanishes and only two lines are observed. In the case $`B_0`$ perpendicular to the symmetry axis the splitting of the $`|\pm 3/2`$ vanishes in first order and again two lines are observed. For case three the quadrupole levels split into four levels with a mixing of the $`|\pm 1/2`$ levels and four lines are observed. The three different cases are corresponding to different orientations of the magnetic field direction with respect to the symmetry axis of the EFG and can be distinguished by investigating the spin-lattice relaxation behavior under various initial conditions.

### 2.2  Spin-lattice Relaxation

In the simplest case, of a two level system, or a virtual two level system (degenerated system), the concept of spin temperature can be applied to the problem and the establishing of equilibrium magnetization $`m(\mathrm{})`$ is described by the single exponential function $`m(t)=m(\mathrm{})(1\mathrm{exp}(t/\tau ))`$. Thus, a line splitting in coincidence with a single exponential spin-lattice recovery function characterizes a line splitting due to a mesoscopic charge segregation and the absence of a magnetic field. The Kramers degeneracy is left unchanged and every line in the spectrum corresponds to a transition between degenerated $`m=\pm 3/2`$ levels. In case of unequally spaced levels (applied magnetic field) the system will come into internal equilibrium and into equilibrium with the lattice with distinct time constants. Furthermore, the recovery is not exponential and depend on the initial conditions i.e. the saturating time. In general, the recovery of the central transition of a four level system into equilibrium is described by a linear combination of three exponential functions

$`{\displaystyle \frac{M(\mathrm{})M(t)}{M(\mathrm{})}}={\displaystyle \underset{i=1}{\overset{3}{}}}a_i\mathrm{exp}(W_it)`$ (4)
whereas the $`a_i`$ depend on the relaxation mechanism and saturating time. In order to extract the correct recovery law one has to know the dominant relaxation mechanism. In some particular cases the recovery law for spin-lattice relaxation can be given and are well known in magnetic resonance. For quadrupole splitting and magnetic relaxation for example a $`I=3/2`$ nuclei like Cu the recovery of the central transition after short saturation time is described by the recovery law

$`{\displaystyle \frac{M(\mathrm{})M(t)}{M(\mathrm{})}}`$ $`=`$ $`0.1\mathrm{exp}(2Wt)`$
$`+`$ $`0.9\mathrm{exp}(12wt)`$
and

$`{\displaystyle \frac{M(\mathrm{})M(t)}{M(\mathrm{})}}`$ $`=`$ $`0.4\mathrm{exp}(2Wt)`$
$`+`$ $`0.6\mathrm{exp}(12wt)`$
for long saturating time . If magnetic and quadrupole relaxation are allowed no analytic solution can be given. The situation in NQR is slightly different since $`|\pm 1/2`$ $``$ $`|\pm 3/2`$, are excited, however, the general facts hold true that the spin-lattice relaxation is described by a multi exponential recovery law and that the component of the recovery

law with the longer relaxation time will increase with increasing saturating time. This characteristic behavior can be utilized to distinguish not only between a magnetic and electric line splitting. In case of a magnetic line splitting it can be utilized to distinguish between the various orientations of the magnetic to the symmetry axis of the EFG tensor.

In case of the orientation of the magnetic field perpendicular to EFG’s symmetry axis the $`|\pm 3/2`$ levels coincidence and the system consist of three levels. The relaxation is now described by two transitions between this three levels and thus, in general, described by a two exponential recovery law. In the two other cases the recovery law is, in general, described by a three exponential law.

In most cases the time constants of the recovery law are close to each other and the three distinct exponential functions can hardly be resolved. In such a case one can not distinguished a three exponential recovery law from a two exponential recovery law by a fitting procedure. However, three and two exponential recovery laws can be distinguish by their answers to various saturating times. In case of a two exponential recovery law varying the saturating time will only affect the coefficients of the individual exponential functions and not the time constants. If, however, a three exponential recovery trace is fitted with a two exponential recovery law, varying the saturating time will obtain different time constants for various saturating times.

Thus we safely conclude, that we can distinguish a line splitting due to a mesoscopic charge segregation and due to a removal of the Kramers degeneracy, by recording the spin-lattice relaxation curve under different initial conditions. Furthermore, in case of a removal of the Kramers degeneracy the orientation of the magnetic with respect to the symmetry axis of the EFG can be deduced from the spin-lattice relaxation curve. Different initial conditions may established e. g. by using the method of saturation recovery and applying different numbers of saturation pulses. Under the condition of mesoscopic charge segregation, one expects a single exponential recovery independent of the number of saturating pulses. However, a removal of the Kramers degeneracy is leading to a multi exponential recovery curve depending on the initial conditions.

## 3 Results

### 3.1 Cu NQR Spectra

In Fig. 2 we show the evolution of 63/65Cu spectra with temperatures from T = 140 K to T = 4.2 K for the originated sample. The corresponding spectra for the fluorinated sample are

shown in Fig. 3. For all temperatures the spectra lie between 13 MHz and 19 MHz. Above $`T_c`$ (T = 140 K) one clearly sees two pairs of lines corresponding to the resonance transition between $`|\pm 3/2>|\pm 1/2>`$ quadrupole levels of the copper nuclear spins belonging to the inner and outer CuO2-layers. The spectra above $`T_c`$ are fitted with equation (3.1) according to the various natural abundance and quadrupole moments of the 63/65Cu isotopes.

$`Y`$ $`=`$ $`y_0+A\mathrm{exp}{\displaystyle \frac{(xx_0)^2}{2W^2}}`$
$`+`$ $`2,3A\mathrm{exp}{\displaystyle \frac{(x1,082x_0)^2}{2(1,082W)^2}}.`$
In accordance with previous NMR studies the pair of lines with higher frequency and larger intensity corresponds to the two outer CuO2-layers while the other pair of lines to the inner one . Most importantly, the 63/65Cu NQR of both sample exhibit the same evolution with temperature.

Below $`T_c`$ the copper NQR spectra are strongly broadened and can not be fitted with two pairs of lines. Instead, a good agreement can be achieved by using six pairs of lines according to equation (3.1). It is important to note that the spectra are fitted with line pairs of fixed relative frequency spacings and relative intensities. This restricts the degree of

freedom of the fits despite their large number of free parameters. If one, for example, fits the lower end of the spectrum and subsequently the upper end the resulting curve has to fit the middle part of the spectrum. Under this condition the resulting fits are reliable if the measured spectra exhibits pronounced peaks.

The occurrence of additional lines within the spectrum is not an evidence for a line splitting. These lines may arise from impurities, Cu outside the CuO layers or Cu from distinct doped phases. However, additional lines from impurities should be clearly visible above $`T_c`$ within the spectrum and should strongly differ from their relaxation times. The same can be stated for lines arising from distinct doped phases.

### 3.2 Spin-lattice relaxation

The spin-lattice relaxation time of Cu in superconducting high-$`T_c`$s is extremely short and driven by antiferromagnetic (AF) fluctuation within the CuO layers . The spin-lattice relaxation times of Cu and Hg, in HgBaCuO for example, differ by two orders of magnitudes . From the ratio $`\gamma (Cu)/\gamma (Hg)`$ of their gyromagnetic ratios one would expect a ratio of the spin-lattice relaxation times of the order of two. The difference is explained by the AF fluctuations within the CuO layers strongly decreasing with distance. As demonstrated by Figs. 4 and 5 the relaxation times vary only smoothly within the spectra what can be interpreted as a proof for Cu inside CuO layers as the origin of the additional lines. Furthermore, these lines have to belong to the same superconducting phase as proofed by the sharp superconducting transition (Fig. 1). From the doping dependency of the Cu NQR frequency in the one layer superconductor HgBaCaCuO one can estimate the difference in $`T_c`$ of phases possessing a Cu NQR frequency difference of 1 MHz to be several tenth of K what is clearly not seen in Fig. 1. So far it is plausible that the additional are due to a line splitting of the two NQR lines visible at 140 K. In order to distinguish between the possible origins of the observed splitting we present in Fig. 6 the 63Cu spin-lattice and spin-spin relaxation curves measured at $`\nu =`$ 16.6 MHz. Various numbers of saturating pulses were applied according to the dependency of the relaxation process illustrated in Sec. 2.2.

As seen in Fig. 6 the spin-lattice relaxation follows a two exponential behavior strongly depending on the number of saturating pulses. The deduced relaxation times are shown in Tab. 1. It is remarkable that the deduced relaxation times depend on the initial conditions. As discussed in section2.2 this can be understood under the assumption of three exponential relaxation behavior fitted with a two exponential recovery law. This assumption is supported by the Cu spin-lattice relaxation in the fluorinated sample which is shown in the lower picture of Fig. (6). Increasing the number of saturating pulses reveals a third relaxation function. Exactly such a behavior is expected from the relaxation theory and equation (2.2) and (2.2). The coefficient of the component with the longer relaxation time increases with increasing saturating time. Furthermore, the transverse relaxation is purely single exponential. Thus, a formation of the internal magnetic fields as an origin of the observed splitting can be concluded.

## 4 Deduced magnetic fields

A magnetic field present at the copper nuclei will split the two line pairs observed at 140 K into four pairs of lines. However, the spectra can not be fitted with four pairs of lines within the whole temperature range. Instead, a reasonable agreement could be achieved by using six pairs of lines. This can be understood by the following assumption. If both outer CuO2-layers would not

be equal below $`T_c`$, the NQR spectrum will consist of three pairs of lines. In this case both outer CuO2-layers differ with respect to their EFG on the copper nuclei. If these three NQR line pairs split due to the formation of magnetic moments in the

CuO2-layers, the spectrum will consist of six pairs of lines. This is indeed observed. One may argue that an electrical inequality is against the crystal symmetry because the inner CuO layer is a symmetry layer. However, the symmetry is immediately broken when a magnetic field appears in the inner layer with a direction parallel to the c-axis. So far the rules of symmetry are not violated. As discussed above the deduced three exponential spin-lattice relaxation is as an indicator for exactly that field direction. It can be stated that the interpretation of the data is consistent.

Associating the line pairs to the CuO2-layers has been done by comparing the line intensities at 4.2 K under the condition that the sum of the intensities of the outer CuO2-layers lines have to be twice of the the intensity of the inner CuO2-layer. Subsequently, the magnetic field has been deduced from the equation $`H=\mathrm{\Delta }\nu /2\gamma `$, with $`\gamma `$ as the gyromagnetic ratio of the 63Cu nucleus and $`\mathrm{\Delta }\nu `$ as the frequency spacing of the 63Cu lines. The results are shown in Fig. 7. The values of the internal magnetic fields vary from the inner to the outer CuO2-layer and the maximum field is of order of 700 Gauss for the originated sample and 1500 gauss for the fluorinated sample .

Now we would like to discuss the possible origin of these internal magnetic fields. We can exclude the presence of the water as a source of the weak antiferromagnetism in the CuO2-plane as was suggested recently, since our measurements were performed on the fresh samples. Furthermore, the values of the deduced magnetic fields in our case is much larger than those measured in the water-inserted YBaCuO compounds. Similarly, the formation of the time-reversal symmetry breaking state  can be excluded since it would not give a magnetic moment on the copper site.

Note that one of the most probable scenarios for the microscopical origin of the observed magnetic fields could be the formation of $`id`$-density wave  or the short-range antiferromagnetic spin correlations remnant from the parent antiferromagnetic compounds and suppressed by the doped carriers. In the first case, one expects that the $`id`$-density wave(DDW) state with the d$`_{x^2y^2}`$-wave symmetry of the order parameter ($`\eta _𝐤=i\frac{\eta _0}{2}(\mathrm{cos}k_x\mathrm{cos}k_y)`$) is formed around the pseudogap formation temperature $`T^{}`$ which for optimally-doped cuprates is close to Tc. Below $`T^{}`$ the staggered orbital currents are formed due to the imaginary character of the order parameter . Furthermore, these currents induce the formation of the internal magnetic moments perpendicular to the CuO2-layer and their magnitudes are increasing upon decreasing temperature. These magnetic moments also would result in the splitting of the NQR lines due to the magnetic fields they produce. In this scenario the value of the internal magnetic fields, DDW gap, and pseudogap formation temperature will be ultimately related to each other. In particular, using the relation obtained earlier between the orbital currents and the DDW order parameter one gets $`j_0=\frac{4et}{\mathrm{}J}\eta _0`$ where $`t=250`$meV is a hopping integral between nearest neighbors in the CuO2-plane and $`J=120`$meV is a corresponding superexchange interaction. One can easily estimate the values of the internal magnetic fields as $`H_{int}=j(T=0)/cr`$ ($`r2\AA `$). For example, taking the value of $`\eta _0=30meV`$ one obtains $`H_{int}`$ of the order of 100G. This value qualitatively agrees with those we obtain from the splitting of NQR lines at T$`=4.2`$ K for the inner CuO2-layer. Note, the corresponding internal magnetic moments are of the order $`l_c0.03\mu _B`$ that agrees with neutron scattering experiments. On the other hand, despite some consistency we notice several problems. For example, the magnitude of the magnetic fields varies and enhances for the outer CuO2-layers in comparison to the inner CuO2-layer. This seems to contradict the $`id`$-density wave scenario, since it is expected that the doping concentration of the holes in the inner and the outer CuO2-layers are different and the latter are more doped. Then, the magnitude of the magnetic moments should be smaller there. This is in contrast to our observation. Therefore, further theoretical and experimental studies are needed in order to identify the microscopical origin of the internal magnetic fields.

## 5 Summary

In summary, we perform copper NQR measurements in the optimally-doped three-CuO2-layer Hg-1223 high-temperature cuprate superconductor. We show that the 63Cu NQR-lines in all three CuO2-layers broaden below $`T_c`$ due to a Zeeman splitting originating from the internal magnetic fields (magnetic moments) within the CuO2-layers. Analyzing the NQR spectra we obtain the values of the internal magnetic fields within the CuO2-layer. For the outer CuO2-layers its values vary from $`400`$G and $`1400`$G, respectively, while for the inner CuO2-layer its value is approximately $`200`$ G. Note, the formation of $`id`$-density wave or short-range antiferromagnetic spin correlations remnant from the parent antiferromagnetic compounds and suppressed by the doped carriers could be possible scenarios as a microscopic origin of the observed fields.

D.M. and I.E. are supported by the INTAS Grant No. 01-0654.