Stability and linear response

The core of the harmonic balance method is expressing the system's behaviour in terms of Fourier components or harmonics. For an $N$-coordinate system, we choose a set of $M_i$ harmonics to describe each coordinate $x_i$ :

$$x_i(t)=\sum _{j=1}^{M_i}{u}_{i,j}(T)\cos ({\omega }_{i,j}t)+v_{i,j}(T)\sin ({\omega }_{i,j}t),$$

This means the system is now described using a discrete set of variables $u_{i,j}$ and $v_{i,j}$. Constructing the vector

$$\mathbf{u}(T)=(u_{1,1}(T),v_{1,1}(T),\dots u_{N,M_N}(T),v_{N,M_N}(T)),$$

we may obtain the harmonic equations (see an example of this procedure)

$$\frac{d\mathbf{u}(T)}{dT}=\overline{\mathbf{F}}(\mathbf{u})$$

where $\overline{\mathbf{F}}(\mathbf{u})$ is a nonlinear function. A steady state ${\mathbf{u}}_0$ is defined by $\overline{\mathbf{F}}({\mathbf{u}}_0)=0$.

Stability

Let us assume that we found a steady state ${\mathbf{u}}_0$. When the system is in this state, it responds to small perturbations either by returning to ${\mathbf{u}}_0$ over some characteristic timescale (stable state) or by evolving away from ${\mathbf{u}}_0$ (unstable state). To analyze the stability of ${\mathbf{u}}_0$, we linearize the equations of motion around ${\mathbf{u}}_0$ for a small perturbation $\delta \mathbf{u}=\mathbf{u}-{\mathbf{u}}_0$ to obtain

$$\frac{d}{dT}[\delta \mathbf{u}(T)]=J({\mathbf{u}}_0)\delta \mathbf{u}(T),$$

where $J({\mathbf{u}}_0)={\mathrm{\nabla }}_{\mathbf{u}}\overline{\mathbf{F}}{|}_{\mathbf{u}={\mathbf{u}}_0}$ is the Jacobian matrix of the system evaluated at $\mathbf{u}={\mathbf{u}}_0$.

The linearised system is exactly solvable for $\delta \mathbf{u}(T)$ given an initial condition $\delta \mathbf{u}(T_0)$. The solution can be expanded in terms of the complex eigenvalues ${\lambda }_r$ and eigenvectors ${\mathbf{v}}_r$ of $J({\mathbf{u}}_0)$, namely

$$\delta \mathbf{u}(T)=\sum _r{c}_r{\mathbf{v}}_r{e}^{{\lambda }_{r}T}.$$

The dynamical behaviour near the steady states is thus governed by $e^{{\lambda }_{r}T}$: if $\mathrm{Re}({\lambda }_r)<0$ for all ${\lambda }_r$, the state ${\mathbf{u}}_0$ is stable. Conversely, if $\mathrm{Re}({\lambda }_r)>0$ for at least one ${\lambda }_r$, the state is unstable - perturbations such as noise or a small applied drive will force the system away from ${\mathbf{u}}_0$.

Linear response

The response of a stable steady state to an additional oscillatory force, caused by weak probes or noise, is often of interest. It can be calculated by solving for the perturbation $\delta \mathbf{u}(T)$ in the presence of an additional drive term.

$$\frac{d}{dT}[\delta \mathbf{u}(T)]=J({\mathbf{u}}_0)\delta \mathbf{u}(T)+\mathit{\xi }{e}^{i\mathrm{\Omega }T},$$

Suppose we have found an eigenvector of $J({\mathbf{u}}_0)$ such that $J(\mathbf{u})\mathbf{v}=\lambda \mathbf{v}$. To solve the linearised equations of motion, we insert $\delta \mathbf{u}(T)=A(\mathrm{\Omega })\mathbf{v}{e}^{i\mathrm{\Omega }T}$. Projecting each side onto $\mathbf{v}$ gives

$$A(\mathrm{\Omega })(i\mathrm{\Omega }-\lambda )=\mathit{\xi }\cdot \mathbf{v}⟹A(\mathrm{\Omega })=\frac{\mathit{\xi }\cdot \mathbf{v}}{-\text{Re}[\lambda ]+i(\mathrm{\Omega }-\text{Im}[\lambda ])}$$

We see that each eigenvalue $\lambda$ results in a linear response that is a Lorentzian centered at $\mathrm{\Omega }=\text{Im}[\lambda ]$. Effectively, the linear response matches that of a harmonic oscillator with resonance frequency $\text{Im}[\lambda ]$ and damping $\text{Re}[\lambda ]$.

Knowing the response of the harmonic variables $\mathbf{u}(T)$, what is the corresponding behaviour of the "natural" variables $x_i(t)$? To find this out, we insert the perturbation back into the harmonic ansatz. Since we require real variables, let us use $\delta \mathbf{u}(T)=A(\mathrm{\Omega })(\mathbf{v}{e}^{i\mathrm{\Omega }T}+{\mathbf{v}}^{\ast }{e}^{-i\mathrm{\Omega }T})$. Plugging this into

$$\delta x_i(t)=\sum _{j=1}^{M_i}\delta u_{i,j}(t)\cos ({\omega }_{i,j}t)+\delta v_{i,j}(t)\sin ({\omega }_{i,j}t)$$

and multiplying out the sines and cosines gives

$$\begin{array}{r}\delta x_i(t)=\sum _{j=1}^{M_i}\{(\text{Re}[\delta u_{i,j}]-\text{Im}[\delta v_{i,j}])\cos [({\omega }_{i,j}-\mathrm{\Omega })t]\\ +(\text{Im}[\delta u_{i,j}]+\text{Re}[\delta v_{i,j}])\sin [({\omega }_{i,j}-\mathrm{\Omega })t]\\ +(\text{Re}[\delta u_{i,j}]+\text{Im}[\delta v_{i,j}])\cos [({\omega }_{i,j}+\mathrm{\Omega })t]\\ +(-\text{Im}[\delta u_{i,j}]+\text{Re}[\delta v_{i,j}])\sin [({\omega }_{i,j}+\mathrm{\Omega })t]\}\end{array}$$

where $\delta u_{i,j}$ and $\delta v_{i,j}$ are the components of $\delta \mathbf{u}$ corresponding to the respective harmonics ${\omega }_{i,j}$.

We see that a motion of the harmonic variables at frequency $\mathrm{\Omega }$ appears as motion of $\delta x_i(t)$ at frequencies ${\omega }_{i,j}\pm \mathrm{\Omega }$.

To make sense of this, we normalize the vector $\delta \mathbf{u}$ and use normalised components $\delta {\hat{u}}_{i,j}$ and $\delta {\hat{v}}_{i,j}$. We also define the Lorentzian distribution

$$L(x{)}_{x_0,\gamma }=\frac{1}{{(x-x_0)}^2+{\gamma }^2}$$

We see that all components of $\delta x_i(t)$ are proportional to $L(\mathrm{\Omega }{)}_{\text{Im}[\lambda ],\text{Re}[\lambda ]}$. The first and last two summands are Lorentzians centered at $\pm \mathrm{\Omega }$ which oscillate at ${\omega }_{i,j}\pm \mathrm{\Omega }$, respectively. From this, we can extract the linear response function in Fourier space, $\chi (\tilde{\omega })$

$$\begin{array}{c}|\chi [\delta x_i](\tilde{\omega }){|}^2=\sum _{j=1}^{M_i}\{[{(\text{Re}[\delta {\hat{u}}_{i,j}]-\text{Im}[\delta {\hat{v}}_{i,j}])}^2+{(\text{Im}[\delta {\hat{u}}_{i,j}]+\text{Re}[\delta {\hat{v}}_{i,j}])}^2]L({\omega }_{i,j}-\tilde{\omega }{)}_{\text{Im}[\lambda ],\text{Re}[\lambda ]}\hfill \\ \hfill +[{(\text{Re}[\delta {\hat{u}}_{i,j}]+\text{Im}[\delta {\hat{v}}_{i,j}])}^2+{(\text{Re}[\delta {\hat{v}}_{i,j}]-\text{Im}[\delta {\hat{u}}_{i,j}])}^2]L(\tilde{\omega }-{\omega }_{i,j}{)}_{\text{Im}[\lambda ],\text{Re}[\lambda ]}\}\end{array}$$

Keeping in mind that $L(x{)}_{x_0,\gamma }=L(x+\mathrm{\Delta }{)}_{x_0+\mathrm{\Delta },\gamma }$ and the normalization $\delta {\hat{u}}_{i,j}^2+\delta {\hat{v}}_{i,j}^2=1$, we can rewrite this as

$$|\chi [\delta x_i](\tilde{\omega }){|}^2=\sum _{j=1}^{M_i}(1+{\alpha }_{i,j})L(\tilde{\omega }{)}_{{\omega }_{i,j}-\text{Im}[\lambda ],\text{Re}[\lambda ]}+(1-{\alpha }_{i,j})L(\tilde{\omega }{)}_{{\omega }_{i,j}+\text{Im}[\lambda ],\text{Re}[\lambda ]}$$

where

$${\alpha }_{i,j}=2(\text{Im}[\delta {\hat{u}}_{i,j}]\text{Re}[\delta {\hat{v}}_{i,j}]-\text{Re}[\delta {\hat{u}}_{i,j}]\text{Im}[\delta {\hat{v}}_{i,j}])$$

The above solution applies to every eigenvalue $\lambda$ of the Jacobian. It is now clear that the linear response function $\chi [\delta x_i](\tilde{\omega })$ contains for each eigenvalue ${\lambda }_r$ and harmonic ${\omega }_{i,j}$ :

• A Lorentzian centered at ${\omega }_{i,j}-\text{Im}[{\lambda }_r]$ with amplitude $1+{\alpha }_{i,j}^{(r)}$

• A Lorentzian centered at ${\omega }_{i,j}+\text{Im}[{\lambda }_r]$ with amplitude $1-{\alpha }_{i,j}^{(r)}$

Sidenote: As $J$ a real matrix, there is an eigenvalue ${\lambda }_r^{\ast }$ for each ${\lambda }_r$. The maximum number of peaks in the linear response is thus equal to the dimensionality of $\mathbf{u}(T)$.

The linear response of the system in the state ${\mathbf{u}}_0$ is thus fully specified by the complex eigenvalues and eigenvectors of $J({\mathbf{u}}_0)$. In HarmonicBalance.jl, the module LinearResponse creates a set of plottable Lorentzian objects to represent this.

Check out this example of the linear response module of HarmonicBalance.jl