text:
# Instabilities and disorder-driven first-order transition of the vortex lattice

## Abstract

Transport studies in a Corbino disk geometry suggest that the Bragg glass phase undergoes a first-order transition into a disordered solid. This transition shows a sharp reentrant behavior at low fields. In contrast, in the conventional strip configuration, the phase transition is obscured by the injection of the disordered vortices through the sample edges, which results in the commonly observed vortex instabilities and smearing of the peak effect in NbSe2 crystals. These features are found to be absent in the Corbino geometry, in which the circulating vortices do not cross the sample edges.

The nature of the disorder-driven solid-solid transition in the vortex matter in superconductors , and the associated instabilities , have recently attracted wide attention. A number of anomalous instability phenomena were reported, which include memory effects , frequency and bias dependence , low frequency noise , history dependent response , slow voltage oscillations , and negative dynamic creep . NbSe2 is a very convenient system for studying these phenomena since it displays a pronounced peak effect (PE) in the critical current $`I_c`$ below the upper critical field $`H_{c2}`$. The PE constitutes a transformation of the quasi-ordered Bragg glass phase  below the PE, to a highly disordered phase in the PE region . Although several sharp features have been observed at the PE  and various models for the PE suggested , there is currently no general consensus regarding the underlying nature of the disordered phase (DP) and of the corresponding order-disorder transition.

An important property of the disorder-driven transition is that the DP can be ‘supercooled’ to well below the transition by a field-cooling procedure, where it remains metastable in the absence of driving currents . In this Letter we address the question of how the metastable DP modifies the apparent transport behavior and what are the actual underlying vortex matter properties in the absence of the metastabilities. Our study is motivated by a recent model, according to which the metastable DP can be formed dynamically by an edge-contamination mechanism . In the presence of a driving current, the flowing vortices have to penetrate into the sample through the surface barriers at the edges . Since the barrier height is non-uniform due to material imperfections, the vortices are injected predominantly at the weakest points of the barrier, hence destroying the ordered phase (OP) and creating a metastable DP near the edges. Consequently, the common experimental procedures inadvertently cause a dynamic admixture of the OP with the metastable DP, thus preventing observation of unperturbed OP and of its behavior in the vicinity of the transition. We have carried out transport measurements in a Corbino disk geometry (inset to Fig. 1), in which the vortices circulate in the bulk without crossing the edges, and hence the injection of the DP could be prevented. These measurements are compared with the regular strip configuration in the same crystals. If the instability effects are a bulk phenomenon, no significant difference between the two geometries is expected, whereas if the sample edges do play a dominant role, qualitatively different behavior should be observed. This paper demonstrates that the two geometries yield strikingly different response and, significantly, the instability effects are eliminated in the Corbino geometry. Furthermore, by removing the complicating effects of the edges, the true bulk dynamics in the vicinity of the transition can be probed for the first time. The results suggest that the disorder-driven transition $`T_{DT}`$ is possibly a thermodynamic first-order phase transition. In addition, it is found that the $`T_{DT}`$ line displays a sharp reentrant behavior at low fields in a wide temperature range.

The results presented here were obtained on a Fe (200ppm) doped 2H-NbSe2 single crystal 2.2 $`\times `$ 1.5 $`\times `$ 0.04 mm3 with $`T_c`$=5.7 K, which display a significantly broader PE as compared to the pure crystals . Similar data were obtained on three additional samples including one pure crystal. The crystals were cut with a very fine wire-saw to ensure uniform edges. Ag$`/`$Au contacts were evaporated through a mask designed for a comparative study of the Corbino and strip geometries on the same crystal, using the same voltage contacts +V,-V (Fig. 1 inset). When measuring the Corbino, the current is applied to the +C,-C contacts, while for the strip configuration the +S,-S contacts are used. To prevent vortices from crossing the edges it is crucial to enforce a uniform radial current in the Corbino. The current uniformity is limited by the variations in the local contact resistance along the outer ring electrode. In order to improve the current uniformity, the ring electrode is divided into four quadrants, and 1/4 of the current is applied to each quadrant. Without such a careful current balance, some of the vortices are injected through the sample edges resulting in metastability effects similar to the strip configuration. The Corbino and the strip geometries were usually compared at the same applied current. In both geometries the current density is not uniform throughout the sample. A calibration above $`T_c`$ shows that the average current density between the +V,-V voltage contacts is higher by about 25% in the strip configuration, consistent with a geometrical calculation. Therefore, all voltage readings ($`V`$) of the strip were divided by about 1.25 in order to obtain the same voltage above $`T_c`$ in the two geometries. Similarly, the $`I_c`$ values of the strip were multiplied by the same factor. The ac measurements were limited to frequencies below 1 kHz to ensure that there were no effects of finite skin depth. The temperature reproducibility between different runs was about 30 mK.

Figure 1a shows the critical current vs. temperature measurements in the vicinity of the PE at 2 kOe. The striking difference between the two geometries is evident. While the dc $`I_c`$ in the strip geometry shows the usual smooth PE, the Corbino $`I_c`$ displays a very sharp drop at the maximum of the PE. Similarly, Fig. 1b shows the voltage at a constant current at 2.5 kOe, which displays a very sharp drop at $`T_{DT}`$ in the Corbino, for both dc and ac currents. Within our temperature reproducibility the response in the Corbino is practically frequency independent in the entire field and temperature range, as expected in the absence of the instability phenomena. The sharp transition in the Corbino geometry reflects a disorder driven transition $`T_{DT}`$ between two thermodynamically stable phases: an OP below $`T_{DT}`$, which is dominated by the elastic energy and is characterized by a low critical current $`J_c^{ord}`$, and a DP above $`T_{DT}`$, which is governed by the pinning energy and has a high $`J_c^{dis}`$.

We now analyze the transport behavior in the strip geometry. Figure 1 shows that the system has a conventional response with no geometry and no frequency dependence above $`T_{DT}`$. This is where only one stable phase, the DP, is present. In contrast, the region below $`T_{DT}`$ is where all the anomalous phenomena in the strip configuration are found . In this region, in the presence of a driving current, a metastable DP is injected through the sample edges instead of the equilibrium OP . The driven metastable DP has a finite lifetime, $`\tau _r`$, that defines a characteristic relaxation length $`L_r`$=$`v\tau _r`$, over which the metastable DP anneals into the OP ($`v`$ is the vortex drift velocity). As a result, the measured dc $`I_c`$ in the strip is a weighted superposition of $`J_c^{ord}`$ with the metastable $`J_c^{dis}`$, given by $`I_c=d_0^WJ_c(x)𝑑xdL_rJ_c^{dis}+d(WL_r)J_c^{ord}`$, where $`d`$ and $`W`$ are the sample thickness and width . For $`L_rW`$ we obtain $`I_c=dWJ_c^{ord}`$, and $`I_c=dWJ_c^{dis}`$ for $`L_rW`$. The relaxation length $`L_r`$ is very sensitive to the proximity to the transition. Near $`T_{DT}`$ the elastic and the pinning energies are comparable, and therefore the metastable DP has a long lifetime, and hence $`L_r`$ is large. Further below $`T_{DT}`$ the driven DP becomes progressively unstable, and therefore $`L_r`$ drops continuously at lower temperatures, resulting in a smooth decrease of $`I_c`$ in the strip . Since $`L_r`$ starts to decrease immediately below $`T_{DT}`$, the maximum in $`I_c`$ of the strip coincides with $`T_{DT}`$, as seen in Fig. 1a. The enhanced $`I_c`$ in the strip below $`T_{DT}`$ causes the dc voltage response to be significantly suppressed and shifted to lower temperatures compared to the Corbino, as seen in Fig. 1b. This voltage has a smoother temperature and field dependence due to the gradual evolution of the effective $`I_c`$ in the strip.

The edge contamination also results in a very pronounced frequency dependence below $`T_{DT}`$. An ac driving current limits the contamination by the DP to narrow regions, of width $`x_d^{ac}`$, near the edges of the strip, where the vortices penetrate and exit during the ac cycle . Since $`x_d^{ac}L_r`$, the $`I_c`$ measured by an ac current in the strip is significantly reduced relative to the dc measurement, as shown in Fig. 1a. Also, since $`x_d^{ac}`$ shrinks with frequency, the voltage increases with frequency and is shifted towards the response of the Corbino as shown in Fig. 1b. The combination of the observed differences in the ac and dc properties below $`T_{DT}`$ in the two geometries, as presented in Fig. 1, cannot be attributed to any known bulk vortex mechanism or inhomogeneities, and clearly demonstrates the key role of the sample edges.

The location of the disorder-driven transition line $`T_{DT}`$ or $`H_{DT}`$ on the $`HT`$ phase diagram, derived from the Corbino data, is shown in Fig. 2. Striking reentrant behavior of $`H_{DT}`$ is observed. At elevated fields, the elastic energy decreases with field resulting in the high-field $`H_{DT}`$ line when it becomes equal to the pinning energy. The recent theoretical studies of the disorder-driven transition have mainly focused on this high-field $`H_{DT}`$ transition . However, the elastic energy also decreases rapidly at low fields where the vortex interactions start to decrease exponentially, which leads to a mirror-like reentrant $`H_{DT}`$ line . Similar arguments involving thermal fluctuations lead to the well known prediction of reentrant melting . Although observations of a reentrant PE were originally interpreted in terms of reentrant melting , subsequent studies attribute this behavior to a reentrant disordering of the lattice . Our results demonstrate the first unambiguous reentrant phase diagram with very sharp and pronounced transition lines. In principle, by measuring the high-field $`H_{DT}`$ line and knowing all the microscopic parameters and the details of the disorder, one should be able to predict the location of the reentrant $`H_{DT}`$ line and compare it with the experiment. Such an extensive theoretical analysis is beyond the scope of this work and will be presented elsewhere.

Both the high-field and the reentrant $`H_{DT}`$ lines can be crossed in a single experimental run by a field sweep at a constant temperature, as presented in Fig. 3 for 4.3 K and 4.6 K. The transport data in the Corbino geometry show remarkably sharp resistive drops at both the low-field and high-field $`H_{DT}`$ points. Note the almost linear field dependence of the resistance in between the two transitions. Such a linear behavior is one of the hallmarks of the weakly pinned OP, as observed in clean NbSe2 crystals . The strip geometry shows markedly different behavior. Below the upper $`H_{DT}`$ in Fig. 3a the resistive onset is gradual and is shifted to a lower field of $`H_{on}^s`$ due to the injection of the metastable DP from the edges, similar to the behavior described in Fig. 1b. Surprisingly, we find that also the instability phenomena display a reentrant behavior, closely following the reentrant $`H_{DT}`$ line. In this case, however, the metastable DP is present at fields above the reentrant $`H_{DT}`$, instead of below the transition, as near the high-field $`H_{DT}`$. As a result, the reentrant $`H_{on}^s`$ of the strip resides above the reentrant $`H_{DT}`$. Also the ac voltage of the strip displays reentrant anomalous behavior, as shown by the dashed curve in Fig. 3a. Near the upper $`H_{on}^s`$ the ac data is shifted towards higher fields with respect to the dc strip response, whereas in the reentrant region a mirror-image like displacement to lower fields is obtained.

The proximity to the $`H_{DT}`$ transition is essential to the enhanced lifetime of the metastable DP. Accordingly, the relaxation length $`L_r`$ diverges upon approaching the high-field $`H_{DT}`$ from below and the reentrant $`H_{DT}`$ from above. In the Corbino geometry the entire area within the triangle formed by the $`H_{DT}`$ line in Fig. 2 represents a stable OP. In the strip case, however, in this triangle, a metastable DP dynamically coexists with the OP. In the belt area between $`H_{on}^s`$ and $`H_{DT}`$, $`L_r`$ is sufficiently large, such that the $`V_{dc}`$ of the strip in Fig. 3a is immeasurably small at the applied current of 20 mA. Within the triangular area of $`H_{on}^s`$, marked by a dashed line, significant vortex motion gradually builds up away from the transition. Yet even here $`L_r`$ remains finite since the full flux-flow vortex velocity of the Corbino is not attained. An interesting case is shown in Fig. 3b where the phase diagram is crossed vertically at 4.6 K, cutting through the tip of the $`H_{on}^s`$ triangle (see Fig. 2). In the strip geometry no voltage response is observed at 20 mA except a noisy behavior near 2 kOe. This is the characteristic noise associated with the described instability phenomena, as reported previously .

The inset to Fig. 2 shows $`I_c`$ at 4.2 K as a function of $`H`$ on crossing the reentrant $`H_{DT}`$. One may expect to see here a PE similar to the high-field peak effect. The high-field PE originates from the fact that $`I_c`$ increases upon crossing from OP to DP, but then drops gradually to zero near $`H_{c2}`$. In the reentrant case, however, this analogy is not complete, since $`I_c`$ does not go to zero upon approaching $`H_{c1}`$. Instead, the $`I_c`$ in the Corbino shows a sharp increase at $`H_{DT}`$, but then continues to grow as the field is decreased within the reentrant DP, since the diluted vortices remain strongly pinned individually . In contrast to the Corbino, the strip configuration shows once again a smooth behavior of $`I_c`$ at $`H_{DT}`$, like a mirror-image of Fig. 1. Thus, regardless of whether the order to disorder transition occurs upon increasing or decreasing the field, the injection of the metastable DP from the sample edges always occurs on the OP side of the transition in the vicinity of $`H_{DT}`$.

The present findings allow us to derive some conclusions regarding the thermodynamic nature of the $`H_{DT}`$ transition. Figures 1 and 3, as well as the inset to Fig. 2, show extremely sharp resistive transitions at $`H_{DT}`$ in the Corbino geometry, which are much sharper than the resistive transition at $`T_c`$. This observation is indicative of the first-order nature of the disorder-driven transition. Obviously, resistivity is not a thermodynamic probe, and further investigations are required. Yet historically, the first strong indications of a possible first-order melting transition came from similarly sharp resistive kinks in HTS crystals , which were confirmed thermodynamically only later . Our results thus imply that the disorder-driven destruction of the Bragg glass is possibly of first order, similar to the thermally driven destruction upon melting. Since both the disordered phase and the vortex liquid do not possess any long range order, it is plausible to expect that their transition into an ordered phase should be of first order, involving topological symmetry breaking in both cases. This means that the second peak transition in HTS , which is of the same nature as the PE in NbSe2, could be of first order as well, thus forming a unified first-order destruction line of the Bragg glass at all temperatures. Recent theoretical considerations seem to support this scenario .

In summary, we find that the ordered Bragg glass becomes unstable with respect to disorder at both high and low fields, resulting in a reentrant disorder-driven transition line. By using a Corbino geometry, and thus avoiding the contamination from the sample edges, this transition is found to be very sharp and apparently of first order. The vortex instability phenomena are caused by the injection of metastable disorder through the sample edges, and are therefore absent in the Corbino geometry. The instabilities are present in the strip geometry on the Bragg glass side of the transition along both the high-field and the low-field branches of the transition line.

This work was supported by the US-Israel Binational Science Foundation (BSF), by the Israel Science Foundation - Center of Excellence Program, and by the Alhadeff Research Award. EYA acknowledges support by the NSF.

FIGURE CAPTIONS

Fig. 1. (a) The critical current $`I_c`$ vs. temperature at 2 kOe measured by a dc current in the Corbino geometry ($`\mathrm{}`$), and by dc ($`\mathrm{}`$) and 172 Hz ac ($``$) current in the strip. The $`I_c`$ is defined at a voltage criterion of 0.5 $`\mu `$V. The dashed line is a schematic guide to the eye extrapolation of the $`I_c^{dis}`$. (b) Voltage vs. temperature at 2.5 kOe and 20 mA, using dc ($``$) and 765 Hz ac ($`\mathrm{}`$) current in the Corbino, and dc, 22, 172 and 765 Hz in the strip. Inset: the electrode configuration allowing measurements in both the Corbino and strip configurations by using +C,-C and +S,-S current contacts, respectively. The outer diameter of the Corbino electrode is 1.1 mm and the distance between the centers of the voltage contacts is 0.15 mm.

Fig. 2. $`HT`$ phase diagram showing the disorder-driven phase transition line $`H_{DT}`$ with reentrant behavior, as determined from the Corbino data. The ordered phase (OP) is present within the triangular region defined by the solid $`H_{DT}`$ line, and is surrounded by the DP both at high and low fields. The instability phenomena in the strip are most pronounced between $`H_{on}^s`$ and $`H_{DT}`$ lines. $`H_{c2}`$ is defined resistively at 10% of the normal state resistance. Inset: dc $`I_c`$ vs. field in the vicinity of the reentrant $`H_{DT}`$ line at 4.2 K in the strip ($``$) and Corbino ($``$) geometries.

Fig. 3. The voltage vs. field at 20 mA at $`T`$ = 4.3 K (a) and 4.6 K (b) measured with dc current in the Corbino ($``$) and strip ($``$), and with 772 Hz ac current in the strip (dashed line). $`H_{DT}`$ marks the position of the thermodynamic phase transition, and $`H_{on}^s`$ is the onset of an observable dc response in the strip geometry.


format: text

text:
# Untitled Document

A finite–dimensional representation of the quantum angular momentum operator

(Short title: Angular momentum in a finite linear space)

Rafael G. Campos

Escuela de Ciencias Físico–Matemáticas, Universidad Michoacana

58060 Morelia, Michoacán, México

e-mail: rcampos@zeus.ccu.umich.mx

L. O. Pimentel

Departamento de Física, Universidad Autónoma Metropolitana–Iztapalapa

Apdo. Postal 55-534, México DF, 09340 México

e-mail: lopr@xanum.uam.mx

Abstract

A useful finite–dimensional matrix representation of the derivative of periodic functions is obtained by using some elementary facts of trigonometric interpolation. This $`N\times N`$ matrix becomes a projection of the angular derivative into polynomial subspaces of finite dimension and it can be interpreted as a generator of discrete rotations associated to the $`z`$–component of the projection of the angular momentum operator in such subspaces, inheriting thus some properties of the continuum operator. The group associated to these discrete rotations is the cyclic group of order $`N`$.

Since the square of the quantum angular momentum $`L^2`$ is associated to a partial differential boundary value problem in the angular variables $`\theta `$ and $`\phi `$ whose solution is given in terms of the spherical harmonics, we can project such a differential equation to obtain an eigenvalue matrix problem of finite dimension by extending to several variables a projection technique for solving numerically two point boundary value problems and using the matrix representation of the angular derivative found before. The eigenvalues of the matrix representing $`L^2`$ are found to have the exact form $`n(n+1)`$, counting the degeneracy, and the eigenvectors are found to coincide exactly with the corresponding spherical harmonics evaluated at a certain set of points.

PACS: 03.65.Ca, 03.65.Lx, 02.60.Jh, 02.60.Lj

1. Introduction

A Galerkin–collocation–type method based in a $`N`$–dimensional matrix representation of $`\frac{d}{dx}`$ obtained through Lagrange’s interpolation has been used to solve one–dimensional boundary value problems (see Refs. –). This technique consists basically in the substitution of the variable $`x`$ and the derivative $`\frac{d}{dx}`$ by $`N\times N`$ matrices $`X`$ and $`D`$ in a certain form of the given differential equation. The diagonal matrix $`X`$ has the $`N`$ different values $`x_j`$, $`j=1,2,\mathrm{},N`$ along the main diagonal, and the matrix $`D`$ is given by

$$D=P\stackrel{~}{D}P^1,$$
$`(1)`$
where

$$\stackrel{~}{D}_{ij}=\{\begin{array}{cc}\underset{l=1}{\overset{N}{\stackrel{}{}}}\frac{1}{(x_ix_l)},\hfill & i=j\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \frac{1}{(x_ix_j)},\hfill & ij\text{,}\hfill \end{array}P_{ij}=p^{}(x_i)\delta _{ij},i,j=1,2,\mathrm{},N,$$
$`(2)`$
The symbol $`^{}`$ appearing in (2) means the sum over $`li`$ and the prime on $`p`$ means differentiation of the polynomial

$$p(x)=\underset{k=1}{\overset{N}{}}(xx_k).$$

The $`N`$ nodes $`x_j`$ can be chosen by imposing a condition where the coefficients of the differential equation and the boundary conditions play the main part. This condition is

$$\underset{k=1}{\overset{N}{\stackrel{}{}}}\frac{1}{(x_jx_k)}=\frac{\gamma ^{}(x_j)}{\gamma (x_j)},j=1,2,\mathrm{},N,$$
$`(3)`$
where $`\gamma (x)`$ is a function defined by the boundary conditions and the differential equation (see Refs. –).

More precisely, the projection scheme for the $`k`$–th derivative of a real function $`f`$ evaluated at different (but otherwise arbitrary) points $`x_i,i=1,2,\mathrm{},N`$ is given by

$$f^{(k)}(x_j)=\underset{l=1}{\overset{N}{}}D_{}^{k}{}_{jl}{}^{}f(x_l)+E_j,j=1,2,\mathrm{},N,$$
$`(4)`$
where $`E_j`$ is the $`j`$–th component of the residual vector depending on $`f`$, $`N`$ and $`k`$. Let us denote by $`\pi _{N1}`$ the space of polynomials of degree at most $`N1`$. Thus, if $`f\pi _{N1}`$, it is found that (4) is exact at the nodes, i.e., $`E_j=0`$, and therefore, any differential problem closed in $`\pi _{N1}`$ can be solved also in $`\text{R}^N`$ as a matrix problem yielding the same solution. An estimate of the error $`E_j`$ for other kind of functions is given in .

The formal application of this elementary technique to some multivariate problems is straightforward, as shown in the following section (see also ).

2. Multivariate case

To illustrate how this method should be applied to some boundary value problems in several variables we will take first two real variables, $`x`$ and $`y`$.

Let $`\stackrel{~}{\mathrm{\Pi }}=\pi _{M1}\pi _{N1}`$ be the space of bivariate tensor–product polynomials of degree $`N1`$ in the variable $`x`$ and $`M1`$ in $`y`$. Thus, if $`f\stackrel{~}{\mathrm{\Pi }}`$, $`f(x,y)`$ can be written as

$$f(x,y)=\underset{j=0}{\overset{N1}{}}\underset{k=0}{\overset{M1}{}}a_{jk}x^jy^k.$$
$`(5)`$
Now let us take two sets of different but otherwise arbitrary points $`\{x_1,x_2,\mathrm{},x_{N_1}\}`$ and $`\{y_1,y_2,\mathrm{},y_{N_2}\}`$, on the $`x`$ and $`y`$ axes respectively, and let $`D_x`$ and $`D_x`$ be the $`N\times N`$ and $`M\times M`$ matrix representations of $`\frac{d}{dx}`$ and $`\frac{d}{dy}`$ constructed according to (1) with such sets.

By deriving partially (5) $`n`$ times with respect to $`x`$ and $`m`$ times with respect to $`y`$, and evaluating the result at the cartesian nodes $`(x_j,y_k)`$ we obtain the vector $`\text{f}^{(n,m)}`$ of dimension $`\stackrel{~}{N}=NM`$ whose entries are given by

$$f_r^{(n,m)}=\frac{^{n+m}f(x,y)}{x^ny^m}|_{(x_j,y_k)}$$

where the indexes $`r`$, $`j`$ and $`k`$ are related through

$$r=j+(k1)N,j=1,2,\mathrm{},N,k=1,2,\mathrm{},M,$$
$`(6)`$
in such a form that $`r=1,2,\mathrm{},\stackrel{~}{N}`$. It is known that a bivariate interpolation on the grid $`(x_j,y_k)`$ is uniquely possible in $`\stackrel{~}{\mathrm{\Pi }}`$ (see ) so that $`f`$ can be taken as a bivariate and sufficiently differentiable real function other than a polynomial, and (5) as its corresponding (tensor) Taylor polynomial. By choosing $`j`$ (the $`x`$–index) to run faster than $`k`$ and using (4) according to the case, we can write down a matrix formula in $`\text{R}^{\stackrel{~}{N}}`$ for $`\text{f}^{(n,m)}`$ in terms of $`\text{f}=\text{f}^{(0,0)}`$, generalizing (4):

$$\text{f}^{(n,m)}=(D_{y}^{}{}_{}{}^{m}D_{x}^{}{}_{}{}^{n})\text{f}+\text{E}.$$
$`(7)`$
Here, the Kronecker product of matrices $`A=(a_{jk})`$ and $`B=(b_{jk})`$ of sizes $`N\times N`$ and $`M\times M`$, respectively, is defined by

$$AB=\left(\begin{array}{cccc}a_{11}B& a_{12}B& \mathrm{}& a_{1,N}B\\ a_{21}B& a_{22}B& \mathrm{}& a_{2,N}B\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ a_{N,1}B& a_{N,2}B& \mathrm{}& a_{N,N}B\end{array}\right).$$

The components of the residual vector E are zero if $`f\stackrel{~}{\mathrm{\Pi }}`$ and an expression for them is given in  for sufficiently differentiable functions.

Let $`1_x`$ and $`1_y`$ denote the identity matrices of dimension $`N`$ and $`M`$ respectively, and define the $`\stackrel{~}{N}\times \stackrel{~}{N}`$ matrices Dx and Dy by

$$\text{D}\text{x}=1_yD_x,\text{D}\text{y}=D_y1_x.$$
$`(8)`$
Due to the properties of the Kronecker product, these matrices commute:

$$\text{D}\text{x}\text{D}\text{y}=\text{D}\text{y}\text{D}\text{x}=D_yD_x.$$

More generally,

$$\text{D}\text{x}^n\text{D}\text{y}^m=(1_yD_x^n)(D_y^m1_x)=D_y^mD_x^n,$$

and therefore, (7) takes the form

$$\text{f}^{(n,m)}=\text{D}\text{x}^n\text{D}\text{y}^m\text{f}+\text{E},$$
$`(9)`$
indicating that the partial derivatives $`/x`$ and $`/y`$ take in $`\stackrel{~}{\mathrm{\Pi }}`$ the tensor–product forms given in (8).

By using the properties of the tensor product and defining $`\text{P}=P_yP_x`$ and

$$\stackrel{~}{\text{D}\text{x}}=1_y\stackrel{~}{D_x},\stackrel{~}{\text{D}\text{y}}=\stackrel{~}{D_y}1_x,$$

where $`P_x`$, $`P_y`$, $`\stackrel{~}{D_x}`$ and $`\stackrel{~}{D_y}`$ have the structure given in (2), it is possible to give the following alternate form of (9):

$$\text{f}^{(n,m)}=\text{P}\stackrel{~}{\text{D}\text{x}}^n\stackrel{~}{\text{D}\text{y}}^m\text{P}^1\text{f}+\text{E}.$$

On the other hand, the projection of the coefficient functions of the differential operator can be written as diagonal matrices since the the product of a function $`a(x,y)`$ by the partial derivatives of the unknown function evaluated at the nodes $`(x_j,y_k)`$ is

$$a_rf_r^{(n,m)}=a(x_j,y_k)f^{(n,m)}(x_j,y_k)$$

($`a(x_j,y_k)`$ should be well defined) and the indexes can be ordered according to (6), producing that in this scheme, the coefficient functions can be represented by $`\stackrel{~}{N}\times \stackrel{~}{N}`$ diagonal matrices whose non-zero elements are given by $`a_r=a(x_j,y_k)`$, where $`r`$, $`j`$ and $`k`$ are related by (6). Let us denote this generic coefficient matrix by A. Thus, the part of the differential operator consisting in the product $`a(x,y)f^{(m,n)}(x,y)`$ takes the matrix form

$$\text{A}\text{D}\text{x}^m\text{D}\text{y}^n\text{f}.$$

If $`a(x,y)`$ accept a Taylor expansion (this condition is too restrictive and it can be relaxed, but it is adequate for our illustrative purposes), we have that A can be defined by the same function $`a(x,y)`$ through

$$\text{A}=a(\text{X},\text{Y}),$$

where X and Y are the matrices given by

$$\text{X}=1_yX,\text{Y}=Y1_x,$$

where the $`N\times N`$ diagonal matrix $`X`$ has the set of points $`x_j`$, $`j=1,2\mathrm{},N`$ along the main diagonal whereas the $`M`$ points $`y_j`$ lie along the main diagonal of the $`M\times M`$ matrix $`Y`$. X and Y represent the variables $`x`$ and $`y`$, respectively.

The generalization to the case of $`q`$ variables $`x_1,x_2,\mathrm{},x_q`$ is straightforward. We will choose $`N_j`$ points on the $`x_j`$–axis, and the projection space as the tensor product of the subspaces of polynomials of degree at most $`N_j1`$ in $`x_j`$, i.e.,

$$\stackrel{~}{\mathrm{\Pi }}=\underset{j=1}{\overset{q}{}}\pi _{N_j1}.$$

The nodes will be ordered in such a way that a function $`f(x_1,\mathrm{},x_q)`$ evaluated at the nodes yields the vector f whose entries are

$$f_r=f(x_{j_1}^1,x_{j_2}^2,\mathrm{},x_{j_q}^q),j_k=1,2,\mathrm{},N_k,k=1,2,\mathrm{},q,$$

and the indexes $`r`$ and $`j_k`$ are related through

$$r=j_1+(j_21)N_1+(j_31)N_1N_2\mathrm{}+(j_q1)\stackrel{~}{N}/N_q,j_k=1,2,\mathrm{},N_k,k=1,2,\mathrm{},q,$$

and $`\stackrel{~}{N}=_{k=1}^qN_k`$. The index than runs faster is $`j_1`$, then $`j_2`$ and so on, yielding that $`r=1,2,\mathrm{},\stackrel{~}{N}`$.

The $`\stackrel{~}{N}\times \stackrel{~}{N}`$ matrix representation of $`/x^k`$, is now

$$\text{D}\text{k}=1_q\mathrm{}1_{k+1}D_k1_{k1}\mathrm{}1_1,$$

where $`1_j`$ is the $`N_j\times N_j`$ identity matrix and $`D_k`$ is a matrix of dimension $`N_k`$ having the structure given by (1). Similarly, the representation of the variable $`x^k`$ is

$$\text{X}\text{k}=1_q\mathrm{}1_{k+1}X^k1_{k1}\mathrm{}1_1.$$

In the next section we adapt this technique to other kind of subspaces and it will be applied to an important physical problem in Section 4.

3. Discrete rotations

In this section we show that a matrix representation of the angular derivative (the derivative of univariate periodic functions), can be related to a generator of discrete rotations (associated to a subgroup of the rotation group) and becomes a finite–dimensional matrix representation of the $`z`$–component of the angular momentum.

The notation and language used in this section is that of Quantum Physics, what can be seen as a digression in this numerical look–like paper, but we think that the reasons are obvious.

Let us begin by considering a complete set of quantum states in a space of finite dimension, $`|\phi _j`$, to be determined later. Here, $`\phi _j`$ indicates the $`j`$–th eigenvalue of the operator associated to the spatial observable $`\phi `$ (an angular variable). Now, we ask for the operator that produces a cyclic permutation, within a phase shift factor, of the complete set of states $`|\phi _1,\mathrm{},|\phi _N`$, that is, the operator $`\mathrm{\Delta }`$ that yields

$$\mathrm{\Delta }|\phi _j=\{\begin{array}{cc}e^{i\gamma _{j+1}}|\phi _{j+1},\hfill & jN\text{,}\hfill \\ \multicolumn{2}{c}{}\\ e^{i\gamma _1}|\phi _1,\hfill & j=N\text{.}\hfill \end{array}$$

Therefore, the representation of $`\mathrm{\Delta }`$ in the $`\phi `$–basis is the matrix of elements $`\phi _j|\mathrm{\Delta }|\phi _k=\mathrm{\Delta }_{jk}`$ given by

$$(\mathrm{\Delta }_{jk})=\left(\begin{array}{cccccc}0& 0& 0& \mathrm{}& 0& e^{i\gamma _1}\\ e^{i\gamma _2}& 0& 0& \mathrm{}& 0& 0\\ 0& e^{i\gamma _3}& 0& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& 0& 0\\ 0& 0& 0& \mathrm{}& e^{i\gamma _N}& 0\end{array}\right),$$
$`(10)`$
yielding that

$$\mathrm{\Delta }^N=e^{i\mathrm{\Gamma }}1,$$
$`(11)`$
where $`\mathrm{\Gamma }=_{j=1}^N\gamma _j`$, and 1 is the identity matrix of dimension $`N`$. Since the determinant of $`\mathrm{\Delta }`$ is

$$det\mathrm{\Delta }=(1)^{N+1}e^{i\mathrm{\Gamma }},$$

$`\mathrm{\Delta }`$ represents a proper rotation and generates a finite subgroup of the rotation group if

$$\mathrm{\Gamma }=\{\begin{array}{cc}2l\pi ,\hfill & N\text{ odd,}\hfill \\ \multicolumn{2}{c}{}\\ (2l+1)\pi ,\hfill & N\text{ even,}\hfill \end{array}$$

with $`l`$ integer. In such a case, (11) becomes

$$\mathrm{\Delta }^N=\{\begin{array}{cc}1,\hfill & \text{}\text{odd,}\hfill \\ \multicolumn{2}{c}{}\\ 1,\hfill & \text{}\text{even,}\hfill \end{array}$$
$`(12)`$
showing that the matrix set of powers of $`\mathrm{\Delta }`$ is a representation of the cyclic group of order $`N`$ or $`2N`$ for $`N`$ odd or even respectively (being a two–valued representation of our finite group of rotations in the latter case). Therefore, every eigenvalue of any power of $`\mathrm{\Delta }`$ is a rational root of unity. Let us make clear this point because it is related with the main result of this section. By fixing the phase shifts $`\gamma _j`$, we can write (10) in more tractable forms, but due to (12) it is necessary to consider the cases of $`N`$ odd or even separately. Thus, we can take the $`N\times N`$ basic circulant permutation matrix

$$\mathrm{\Delta }_{jk}=\{\begin{array}{cc}\delta _{jN},\hfill & \text{for }j=1\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \delta _{j,k+1},\hfill & \text{for }j2\text{,}\hfill \end{array}$$

for (10) in the odd case, whereas, for $`N`$ even,

$$\mathrm{\Delta }_{jk}=\{\begin{array}{cc}\delta _{jN},\hfill & \text{for }j=1\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \delta _{j,k+1},\hfill & \text{for }j2\text{.}\hfill \end{array}$$

It is not difficult to see that the characteristic polynomial of $`\mathrm{\Delta }`$ in the variable $`\lambda `$ is $`(\lambda )^N+1=0`$ for any $`N`$, and that the $`N`$–tuple of components $`e^{ika_j}/\sqrt{N}`$, $`k=1,2,\mathrm{},N`$, is the $`j`$–th normalized eigenvector of $`\mathrm{\Delta }`$ yielding the eigenvalue $`\lambda _j=e^{ia_j}`$, where

$$a_j=\frac{2\pi j}{N},jI_N,$$
$`(13)`$
and

$$I_N=\{\begin{array}{cc}0,\pm 1,\pm 2,\mathrm{},\pm n,\hfill & N=2n+1\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \pm 1/2,\pm 3/2,\mathrm{},\pm (2n1)/2,\hfill & N=2n\text{.}\hfill \end{array}$$

Being $`\mathrm{\Delta }`$ unitary, it defines an hermitian matrix $`A`$ through

$$\mathrm{\Delta }=e^{iA}$$
$`(14)`$
whose eigenvalues are given by (13) (note that the trace of $`A`$ vanishes) and its eigenvectors, denoted by $`|a_j`$, are those of $`\mathrm{\Delta }`$. Therefore, the unitary matrix diagonalizing simultaneously $`\mathrm{\Delta }`$ and $`A`$ is the one with entries

$$\phi _k|a_j=\frac{1}{\sqrt{N}}e^{ika_j},k,j=1,2,\mathrm{},N.$$
$`(15)`$
In the usual quantum case $`A`$ becomes proportional to the angle of rotation times $`L_z`$. To find a discrete analogue relation we need to calculate the elements of $`A`$ in the $`\phi `$–basis. This can be made through

$$A_{jk}=\phi _j|A|\phi _k=\underset{l=1}{\overset{N}{}}\phi _j|a_la_la_l|\phi _k=\frac{2\pi }{N^2}\underset{lI_N}{}le^{i(jk)a_l}.$$

Note that the diagonal elements vanish all of them. For $`N=2n+1`$ we have

$$A_{jk}=\frac{4\pi }{N^2}i\underset{l=1}{\overset{n}{}}l\mathrm{sin}\frac{2\pi l(jk)}{N},$$

whereas for $`N=2n`$,

$$A_{jk}=\frac{4\pi }{N^2}i\underset{l=1}{\overset{n}{}}(l1/2)\mathrm{sin}\frac{2\pi (l1/2)(jk)}{N}.$$

These sums can be calculated easily  and, for $`N`$ odd or even, they become

$$A_{jk}=\phi _j|A|\phi _k=\{\begin{array}{cc}0,\hfill & j=k\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \frac{i(1)^{j+k}\pi /N}{\mathrm{sin}\frac{\pi (jk)}{N}},\hfill & jk\text{.}\hfill \end{array}$$
$`(16)`$
We proceed now to establish a relation between $`A`$ and a $`N`$–dimensional matrix representation of the derivative of trigonometric polynomials. We choose the odd case $`N=2n+1`$ first.

3.1 Odd case

Let $`\tau _n`$ be the space of trigonometric polynomials of degree at most $`n`$. It is well known (see for example ) that any trigonometric polynomial $`f\tau _n`$ can be uniquely determined by its values at $`2n+1`$ arbitrary points $`\pi <x_1<x_2<\mathrm{}<x_{2n+1}\pi `$, (we change our notation for simplicity) through the formula

$$f(x)=\underset{k=1}{\overset{2n+1}{}}f(x_k)\frac{t_k(x)}{t_k(x_k)},$$

where the polynomials $`t_k\tau _n`$ are given by

$$t_k(x)=\underset{lk}{\overset{2n+1}{}}\mathrm{sin}\frac{(xx_l)}{2}.$$
$`(17)`$
To reach our goal we need to calculate $`\frac{df(x)}{dx}`$ at the nodes. The differentiation of (17) and algebraic manipulation of the result yields

$$t_k^{}(x_j)=\{\begin{array}{cc}t^{}(x_k)\underset{l=1}{\overset{2n+1}{\stackrel{}{}}}\mathrm{cot}\frac{(x_kx_l)}{2},\hfill & j=k\text{,}\hfill \\ \multicolumn{2}{c}{}\\ t^{}(x_j)\mathrm{csc}\frac{(x_jx_k)}{2},\hfill & jk\text{,}\hfill \end{array},$$

where

$$t(x)=t_k(x)\mathrm{sin}\frac{(xx_k)}{2}=\underset{l=1}{\overset{2n+1}{}}\mathrm{sin}\frac{(xx_l)}{2}.$$

Therefore, $`f^{}(x_j)`$ takes the form

$$f^{}(x_j)=\frac{1}{2}t^{}(x_j)\underset{k=1}{\overset{2n+1}{\stackrel{}{}}}\frac{f(x_k)}{\mathrm{sin}\frac{(x_jx_k)}{2}t^{}(x_k)}+\frac{1}{2}f(x_j)\underset{k=1}{\overset{2n+1}{\stackrel{}{}}}\mathrm{cot}\frac{(x_jx_k)}{2}.$$
$`(18)`$
Note that $`f^{}(x)`$ is again an element of $`\tau _n`$ and that this equation has the structure of (4) with the definitions

$$D=D_\phi =T\stackrel{~}{D}T^1,$$

$$\stackrel{~}{D}_{ij}=\{\begin{array}{cc}\underset{l=1}{\overset{N}{\stackrel{}{}}}\frac{1}{2}\mathrm{cot}\frac{(x_ix_l)}{2},\hfill & i=j\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \frac{1}{2}\mathrm{csc}\frac{(x_ix_j)}{2},\hfill & ij\text{,}\hfill \end{array}T_{ij}=t^{}(x_i)\delta _{ij},i,j=1,N.$$
$`(19)`$
Let choose the nodes such that

$$x_j=\pi +\frac{2\pi j}{N},j=1,2,\mathrm{},N.$$
$`(20)`$
Then, it is not difficult to see that regardless the parity of $`N`$, the product (17) and the sum given by $`\stackrel{~}{D}_{jj}`$, evaluated at (20), satisfy the equations

$$\frac{t_j(x_j)}{t_k(x_k)}=\frac{t^{}(x_j)}{t^{}(x_k)}=(1)^{j+k},$$

and

$$\stackrel{~}{D}_{jj}=\underset{k=1}{\overset{N}{\stackrel{}{}}}\mathrm{cot}\frac{(jk)\pi }{N}=0,j=1,2,\mathrm{},N.$$

Therefore, substituting these expressions in (18) and taking into account(16), we obtain for this choice of nodes and $`N=2n+1`$,

$$A=i\frac{2\pi }{N}D_\phi =\frac{2\pi }{N}L_z,$$
$`(21)`$
where we have defined the $`N`$-dimensional matrix $`L_z`$ as

$$L_z=iD_\phi $$
$`(22)`$
(we take $`\mathrm{}=1`$). Therefore, the equation (14) for the discrete rotation operator $`\mathrm{\Delta }`$ becomes

$$\mathrm{\Delta }=e^{{\displaystyle \frac{2\pi }{N}}D_\phi }=e^{iϵL_z},$$
$`(23)`$
where $`ϵ=2\pi /N`$. Like it comes out in the continuum case, the argument of the exponential is proportional to the angle of rotation (in this case $`ϵ`$) times the derivative $`D_\phi `$. Thus, the interpretation of $`L_z`$ given by (22) as generator of discrete rotations follows immediately. On account of (21), (13) and (15), $`L_z`$ has the $`2n+1`$ eigenvalues

$$\{n,n+1,\mathrm{},1,0,1,\mathrm{},n1,n\},$$

and the normalized eigenvectors $`|m`$, $`m=n,\mathrm{},n`$ whose components, in the $`\phi `$–basis are

$$\phi _k|m=\frac{1}{\sqrt{N}}e^{im\phi _k},$$

where $`\phi _k=kϵ`$, $`k=1,2,\mathrm{},N`$.

The fact that the quantum rule for the $`z`$–component of the angular momentum can be projected in a finite–dimensional space maintaining the form it has in the continuum case (a consequence of the use of $`D_\phi `$ as a projection of $`\frac{d}{d\phi }`$), reinforces the possibility of the construction of a finite–dimensional algebra for Quantum Mechanics as it appears in other problems –.

We consider now the case $`N=2n`$, but before we have an important remark. According to the uniqueness of the representation of trigonometric polynomials in terms of the Dirichlet kernel evaluated at differences of (20) (see for example ) we have that $`D_\phi `$ equals to $`(2n+1)/2`$ times the derivative of the Dirichlet kernel evaluated at $`x_jx_k`$. On the other hand, it should be pointed out the dependence of the explicit form of the matrix representation of $`\frac{d}{d\phi }`$ on the type of projection functions; this means that it is possible to construct other matrices representing the derivative $`\frac{d}{d\phi }`$, if we restrict the points to be in $`(0,\pi ]`$ and project on the cosine polynomials.

3.2. Even case

Essentially, the case of $`N=2n`$ differs from the odd case only by the space of functions where the projection takes place. Let us begin by considering a function of the form

$$g(x)=\mathrm{sin}(x/2)f(x),$$
$`(24)`$
where $`f\tau _{n1}`$. Let $`\pi <x_1^{}<x_2^{}<\mathrm{}<x_{2n1}^{}\pi `$ be $`N1`$ arbitrary points, different from zero. We can interpolate $`f`$ at these nodes to yield

$$g(x)=\underset{k=1}{\overset{N1}{}}f(x_k^{})\frac{\mathrm{sin}(x/2)t_k^{}(x)}{t_k^{}(x_k^{})}=\underset{k=1}{\overset{N1}{}}\left[\mathrm{sin}(x_k^{}/2)f(x_k^{})\right]\frac{\mathrm{sin}(x/2)t_k^{}(x)}{\mathrm{sin}(x_k^{}/2)t_k^{}(x_k^{})},$$

that is,

$$g(x)=\underset{k=1}{\overset{N1}{}}g(x_k^{})\frac{\mathrm{sin}(x/2)t_k^{}(x)}{\mathrm{sin}(x_k^{}/2)t_k^{}(x_k^{})}.$$
$`(25)`$
The functions $`t_k^{}(x)`$ are Gaussian polynomials, given by a product like (17) at the $`N1`$ nodes $`x_j^{}`$. Now, let us consider the set of points formed by zero and $`x_j^{}`$, $`j=1,2,\mathrm{},N1`$, and denote them by $`x_j`$, $`j=1,2,\mathrm{},2n`$. Let $`t_k(x)`$ be the corresponding basic interpolatory polynomials. Then, taking into account that $`g(0)=0`$, (25) becomes

$$g(x)=\underset{k=1}{\overset{2n}{}}g(x_k)\frac{t_k(x)}{t_k(x_k)},$$
$`(26)`$
and we have an interpolation formula for functions of the form (24) at the points $`x_j`$. Therefore, most of the reasonings made in the odd case can be applied to (26) to yield a formula similar to (18) for $`g(x)`$, so that, for $`2n`$ points given by any $`2n1`$ different points of $`(\pi ,\pi ]`$ and zero, the matrix

$$D=D_\phi =T\stackrel{~}{D}T^1,$$
$`(27)`$
with the definitions (19) is an exact representation of the derivative of functions given by (24). Moreover, taking into account that $`\mathrm{cos}x/2=\mathrm{sin}[(x+\pi )/2]`$ and that $`\pi `$ can be taken as one of our nodes and that this function vanishes just at $`x=\pi `$, we can manipulate the function $`h(x)=\mathrm{cos}(x/2)f(x)`$ in the same form as we did with (24) to conclude that the $`k`$-th power of (27), constructed with a set of $`2n`$ distinct points of $`(\pi ,\pi ]`$ that includes zero and $`\pi `$, is a $`2n`$–dimensional matrix projection of the $`k`$-th derivative on the subspace of functions of the type

$$e^{ix/2}f(x),$$

where $`f\tau _{n1}`$. Now, if we restrict these nodes to be equidistant as in (20), the same formulas (21)-(23) for the relations between $`A`$, $`D_\phi `$, $`L_z`$, and $`\mathrm{\Delta }`$ are yielded, but now for $`N=2n`$. This case makes up a two–valued or spin representation of the finite group of rotations.

Using the same definitions of the preceding case, we have that $`L_z`$ has the $`2n`$ eigenvalues

$$\{(2n1)/2,(2n3)/2,\mathrm{},1/2,1/2,\mathrm{},(2n3)/2,(2n1)/2\},$$

and the normalized eigenvectors $`|m`$, $`m=(2n1)/2,\mathrm{},(2n1)/2`$ are given again by

$$\phi _k|m=\frac{1}{\sqrt{N}}e^{im\phi _k},$$

where $`\phi _k=kϵ`$, $`k=1,2,\mathrm{},N`$.

To end this section let us remark that the validity of formulas (21)–(23) for both $`N`$ odd and even shows one discrete formulation of rotations in terms of an finite–dimensional matrix representation of the derivative of certain periodic functions that can be associated, according to quantum postulates, to one component of the angular momentum operator projected in $`\text{C}^N`$, giving thus, a finite subgroup of the rotation group that shares some of the properties of the full group.

4. The numerical eigenproblem of $`L^2`$.

Our purpose in this section is to obtain a finite–dimensional matrix representation of the angular momentum eigenproblem of a system described by three classical degrees of freedom by applying the results of sections 2 and 3.

As usual, we choose the spherical variables $`\theta `$ and $`\phi `$ to describe the problem. In order to apply the results of section 2, we need two matrices: $`D_\theta `$ of dimension $`N`$, and $`D_\phi `$ of dimension $`M`$, to represent $`d/d\theta `$ and $`d/d\phi `$ respectively. Since $`\pi <\phi \pi `$ and the functions to represent (the spherical harmonics $`Y_n^m(\theta ,\phi )`$) are trigonometric polynomials, we will use as matrix $`D_\phi `$, the one given by (19) with the change of notation $`x_j\phi _j`$, where the $`M`$ points $`\phi _j`$ are given by formula (20):

$$\phi _j=\pi +\frac{2\pi j}{M},j=1,2,\mathrm{},M.$$

Concerning the variable $`\theta `$, we have several alternatives to choose a $`N\times N`$ matrix for $`D_\theta `$. Among these, we present only two of them yielding exact results at the nodes. The first one follows the ideas given here and the second is taken from .

4.1. A matrix for $`L^2`$

Since $`0<\theta \pi `$, we can take any set of $`N`$ distinct points $`\theta _j`$ of $`(0,\pi )`$ to construct the matrix $`D_\theta `$ according to (19). Besides, the fact that $`Y_n^m(\theta ,\phi )`$ are defined for integral indexes excludes the use of the spin representations, i.e., the matrices $`D_\phi `$ and $`D_\theta `$ should be constructed with an odd number of points.

Now, according to the results of section 2, the differential eigenvalue problem for $`f_s(\theta ,\phi )=Y_n^m(\theta ,\phi )`$,

$$\frac{^2f_s}{^2\theta }+\mathrm{cot}\theta \frac{f_s}{\theta }+\frac{1}{\mathrm{sin}^2\theta }\frac{^2f_s}{^2\phi }=\lambda _sf_s=n(n+1)f_s,$$

takes the matrix form

$$\text{L}^2\text{f}_s^{}=\lambda _s^{}\text{f}_s^{},s=1,2\mathrm{},NM,$$
$`(28)`$
where $`\text{L}^2`$ is the $`NM\times NM`$ matrix given by

$$\text{L}^2=\left[\text{D}\text{θ}^2+\mathrm{cot}(\theta )\text{D}\text{θ}+\mathrm{sin}^2(\theta )\text{D}\text{φ}^2\right],$$
$`(29)`$
$`\text{f}_s^{}\text{C}^{NM}`$ and $`\lambda _s^{}`$ is in general a complex number. The matrices of (29) are given by

$$\theta =1_M\mathrm{\Theta },\text{D}\text{θ}=1_MD_\theta ,\text{D}\text{φ}=D_\phi 1_N,$$

where $`\mathrm{\Theta }`$ is a diagonal matrix with entries $`\mathrm{\Theta }_{jk}=\theta _j\delta _{jk},`$ and $`1_K`$ is the identity matrix of dimension $`K`$. Thus, (29) can be rewritten as

$$\text{L}^2=1_M[D_\theta ^2+\mathrm{cot}(\mathrm{\Theta })D_\theta ]+D_\phi ^2\mathrm{sin}^2(\mathrm{\Theta }).$$
$`(30)`$
Since $`Y_n^m(\theta ,\phi )`$ is the (tensor) product of the associated Legendre functions $`P_n^m(\theta )`$ and $`e^{im\theta }`$ and, on the other hand, we know (from the results of section 3) that the degree of the polynomial that can be differentiated through (19) yielding exact values at $`N`$ nodes ($`N`$ odd) is at most $`(N1)/2`$, formula (28) reproduces exactly the first

$$\underset{l=0}{\overset{(N1)/2}{}}(2l+1)=\left(\frac{N+1}{2}\right)^2$$

eigenvalues and the corresponding unnormalized functions $`Y_n^m(\theta ,\phi )`$ evaluated at the nodes $`(\theta _j,\phi _k)`$, provided $`MN`$. In this case $`\text{L}^2`$ will have necessarily $`(N+1)^2/4`$ real eigenvalues given by

$$\lambda _s^{}=n(n+1),$$

ordered according to

$$s=n^2+n+m+1,n=0,1,\mathrm{},(N1)/2,m=n,n+1,\mathrm{},n1,n.$$

The eigenvector $`\text{f}_s^{}`$ corresponding to $`\lambda _s^{}`$, has the components $`f_{rs}^{}`$ given by

$$f_{rs}^{}=c_{nm}P_n^m(\theta _j)e^{imkϵ},$$

where $`ϵ=2\pi /M`$, $`c_{nm}`$ is a normalization constant, $`\theta _j(0,\pi )`$, and the relation between $`r`$, $`j`$ and $`k`$, is given by (6).

It is possible to choose the $`\theta `$–nodes in such a way that $`\text{L}^2`$ becomes a positive semidefinite matrix (save to a similarity transformation). To see this, note that $`D_\phi ^2`$ is positive semidefinite whereas $`\mathrm{sin}^2(\mathrm{\Theta })`$ is positive definite since $`\theta _j(0,\pi )`$, $`j=1,2,\mathrm{},N`$. Therefore, according to (30) we only have to find the conditions which make the matrix

$$T_\theta ^1L_\theta ^2T_\theta =\stackrel{~}{D}_\theta ^2\mathrm{cot}(\mathrm{\Theta })\stackrel{~}{D}_\theta $$
$`(31)`$
positive semidefinite. To this end, let us separate the main diagonal of $`D_\theta `$ by writing

$$D_\theta =T_\theta D^{}T_\theta ^1+d,$$

where

$$D_{ij}^{}=\{\begin{array}{cc}0,\hfill & i=j\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \frac{1}{2}\mathrm{csc}\frac{(\theta _i\theta _j)}{2},\hfill & ij\text{,}\hfill \end{array}d_{ij}=\delta _{ij}\underset{l=1}{\overset{N}{\stackrel{}{}}}\frac{1}{2}\mathrm{cot}\frac{(\theta _i\theta _j)}{2}$$

\[cf. (19)\], and $`T_\theta `$ is given by (19). Thus,

$$L_\theta ^2=T_\theta [D_{}^{}{}_{}{}^{2}+D^{}d+dD^{}+d^2+\mathrm{cot}(\mathrm{\Theta })D^{}+\mathrm{cot}(\mathrm{\Theta })d]T_\theta ^1.$$
$`(32)`$
If we can find points $`\theta _j`$ such that $`d=\mathrm{cot}(\mathrm{\Theta })/2`$, i.e.,

$$\underset{l=1}{\overset{N}{\stackrel{}{}}}\mathrm{cot}\frac{(\theta _j\theta _l)}{2}=\mathrm{cot}(\theta _j)$$
$`(33)`$
\[cf. (3)\], equation (32) becomes

$$L_\theta ^2=T_\theta (D_{}^{}{}_{}{}^{2}+dD^{}D^{}d+d^2)T_\theta ^1=T_\theta (D^{}+d)(D^{}+d)T_\theta ^1=T_\theta \stackrel{~}{D_\theta ^t}\stackrel{~}{D_\theta }T_\theta ^1,$$

where $`\stackrel{~}{D_\theta ^t}`$ is the transpose of $`\stackrel{~}{D_\theta }`$. It remains to show that a solution of (33) always exists. Since $`\mathrm{cot}x`$ is the logarithmic derivative of $`\mathrm{sin}x`$, Eq. (33) is the condition for a critical point of the function of $`N`$ variables

$$U(z)=U(z_1,z_2,\mathrm{},z_N)=\underset{k=1}{\overset{N}{}}\mathrm{sin}(z_k)\underset{i>j}{\overset{N}{}}\mathrm{sin}\frac{(z_iz_j)}{2},$$

where $`0z_j\pi `$, $`j=1,2,\mathrm{},N`$. The existence (and uniqueness) of the solution can be proved along the same lines given in . It is worth to be noticed that the similarity transformation is not essential for an eigenproblem like (28) since the similarity matrices corresponding to $`\theta `$ and $`\phi `$ can be collected into a $`NM\times NM`$ diagonal matrix S to write (28) in the form

$$\text{L}_P^2\text{g}_s=\lambda _s^{}\text{g}_s,$$

where $`\text{g}_s=\text{S}^1\text{f}_s^{}`$ and $`\text{L}_P^2`$ is positive semidefinite. Thus, if we construct $`L_\theta ^2`$ with the set of nodes satisfying (33), $`\text{L}^2`$ is a positive semidefinite matrix (save a similarity transformation), a necessary property from the numerical and physical point of view in any projection scheme.

4.2. Other matrix for $`L^2`$

In , a matrix representation of the derivative for trigonometric polynomials of definite parity is given and, therefore, it can also be used to construct a representative of (31), yielding a matrix for $`L^2`$ with a higher degree of approximation than (29). The cost we have to afford for this is the lacking of simplicity: we can not use a single matrix for $`\frac{d}{d\theta }`$ to be substituted directly in (31) unless we accept to loose precision in the results (see ). According to this scheme, the representation of a differential operator is formed following a given rule where certain matrix D is involved. Such a matrix is constructed with $`N`$ arbitrary distinct points $`\theta _j(0,\pi )`$ through a formula similar to (19)

$$\text{D}=\text{D}_\theta =S\stackrel{~}{\text{D}}S^1,$$

$$\stackrel{~}{\text{D}}_{ij}=\{\begin{array}{cc}\underset{l=1}{\overset{N}{\stackrel{}{}}}\mathrm{cot}(\theta _i\theta _l),\hfill & i=j\text{,}\hfill \\ \multicolumn{2}{c}{}\\ \mathrm{cot}(\theta _i\theta _j),\hfill & ij\text{,}\hfill \end{array}S_{ij}=\delta _{ij}\underset{lj}{\overset{N}{}}\mathrm{sin}(\theta _j\theta _l).$$

The form that $`L_\theta ^2`$ adopts in this case is

$$L_\theta ^2=\text{D}_\theta ^2\mathrm{cot}(\mathrm{\Theta })\text{D}_\theta +NSOS^1,$$
$`(34)`$
where $`\mathrm{\Theta }`$ is again a diagonal matrix with entries $`\mathrm{\Theta }_{jk}=\theta _j\delta _{jk}`$ and $`O`$ is a projection matrix with ones everywhere, i.e., $`O_{jk}=1`$. According to , the degree of approximation of (34) is higher than that of (31). While the former yields exact results at $`N`$ nodes ($`N`$ odd or even) for trigonometric polynomials (of definite parity) of degree $`N`$, the latter produces exact results for polynomials of degree $`(N1)/2`$ ($`N`$ odd). Thus, the use of (34) in (30) gives the matrix

$$\text{L}^2=1_M[\text{D}_\theta ^2+\mathrm{cot}(\mathrm{\Theta })\text{D}_\theta NSOS^1]+D_\phi ^2\mathrm{sin}^2(\mathrm{\Theta })$$
$`(35)`$
with a higher degree of approximation: if $`N`$ is an odd integer and $`M=2N+1`$, then (35) produces the first

$$\underset{l=0}{\overset{N}{}}(2l+1)=(N+1)^2$$

exact eigenvalues (and eigenvectors) while (30) yields only the first $`(N+1)^2/4`$ for the same values of $`N`$ and $`M`$.

5. Final remarks

Summarizing, we have found a finite–dimensional representation of the square of quantum angular momentum with the following properties:

1. The coordinate representation of the $`z`$–component of the angular momentum operator is maintained in this discrete scheme and the projection of $`L_z`$ is the generator of a finite subgroup of the group of rotations.

2. The projection of $`L^2`$ is given in terms of finite–dimensional representations of the partial derivatives according to the well-known formula for $`L^2`$.

3. The spectrum of this projection contains the first eigenvalues of $`L^2`$ (counting the degeneracy) in such a form that the whole spectrum of $`L^2`$ can be reobtained when the number of nodes tends to infinity.

4. The eigenvectors of this projection, corresponding to the exact eigenvalues of $`L^2`$, can be converted into the exact eigenfunctions through an interpolation at the nodes. Again, this process yields the complete set of eigenfunctions of $`L^2`$ when the dimension of $`\text{L}^2`$ tends to infinity.

Finally, we note that these properties makes $`\text{L}^2`$ suitable for numerical applications to quantum problems as it will be shown in a subsequent work.

Acknowledgments

RGC wants to thank Prof. O. Obregón and Prof. L. Sabinin for very useful discussions. This research has been partially supported by Consejo Nacional de Ciencia y Tecnología (CONACYT).

REFERENCES

 F. Calogero, Lett. Nuovo Cimento 35, 273 (1982).

 F. Calogero and E. Franco, Nuovo Cimento B 89, 161 (1985).

 M. Bruschi, R.G. Campos and E. Pace, Nuovo Cimento B 105, 131 (1990).

 R.G. Campos, Bol. Soc. Mat. Mexicana, 3, 279 (1997)

 R.G. Campos y R. Muñoz B., Rev. Mex. Fís. 36, 1 (1990).

 F. Calogero, J. Math. Phys. 34, 4704 (1993)

 E.W. Cheney, Multivariate approximation theory: selected topics (SIAM Publications, Philadelphia, Pennsylvania, 1986)

 F. Calogero, Lett. Nuovo Cimento 38, 453 (1983).

 I.S. Gradshteyn and I.M. Ryzhik, Table of integrals, series, and products (Academic Press, London, 1994), 5th. Ed.

 A.F. Timan, Theory of Approximation of functions of a real variable (Dover Publications, Inc., New York, 1994)

 R.G. Campos, Rev. Mex. Fis, 29, 217 (1983).

 R.G. Campos, Rev. Mex. Fis., 32, 379 (1986).

 F. Calogero, Lett. Nuovo Cimento 39, 305 (1984).

 T. Popoviciu, Bull. Math. Soc. Roumaine Sci., 38, 73 (1936).


text:
# Phase-Induced (In)-Stability in Coupled Parametric Oscillators

## 1 Introduction

Parametric resonance is a phenomenon pervading several fields of science. It occurs when the modulation of a system parameter causes the system to become unstable. The literature on parametric resonance is enormous, so the best one can do is cite a range of different subjects where it may play a role. Examples include mechanical systems where such resonances were first identified , elementary particles , quantum dots , astrophysics , fluid mechanics , plasma physics , electronic networks , superconducting and laser devices , biomechanics , and even medicine . Connections with chaotic systems have been suggested recently . The simplest and perhaps most familiar example of parametric resonance occurs in a harmonic oscillator whose frequency varies periodically with time . For certain ranges of modulation parameters (frequency, amplitude) the oscillator is unstable while for others it is stable. Even for this seemingly simple system the stability boundary diagram is already quite rich and complex (see below).

A great deal of recent work has dealt with systems of coupled oscillators – again, the literature in this general area is enormous. However, very little attention has focused on systems of coupled parametric oscillators . Such coupled arrays are particularly intriguing because each single oscillator alone exhibits regions of stable or unstable behavior. Two questions arise naturally: 1) How does coupling modify the single-oscillator stability boundaries? 2) Are there collective parametric instabilities of the coupled system that are distinct from those experienced by single oscillators? In a recent paper Bena and Van den Broeck  addressed these questions for an infinite set of globally mean-field coupled harmonic oscillators subjected to time-periodic block pulses with quenched uniformly distributed random phases (“quenched” in this context means that the phase of frequency modulation of each oscillator is set at time $`t=0`$ and then remains unchanged). Within their mean field treatment they are able to deal with both of the questions posed above. In particular, they find wide ranges of parameter values that lead to collective instabilities.

Our interest in this paper is to identify generic coupling effects that are not specifically a consequence of the mean field analysis, and in this context to understand how a coupled system of (very) few oscillators with short-ranged interactions might carry in it the seeds of the infinite/infinitely cross-coupled array. In particular, we seek the seeds of the collective parametric instabilities. We do this by studying a model of two coupled oscillators. The control parameter in our study is the phase difference in the periodic modulation of the frequencies of the two oscillators.

In Section 2 we introduce the model and its partially analytic solution. Section 3 is a presentation of results as a large collection of stability boundary diagrams that convey the way in which stability boundaries shift, appear, and disappear with parameter changes. In Section 4 we analyze the ways in which the seeds of the collective parametric instabilities found in the mean field model already appear in the two-oscillator system. A brief summary of our conclusions is reiterated in Section 5.

## 2 The model

Our system consists of two linearly coupled parametric oscillators whose equations of motion are

$`\ddot{x}_1`$ $`=`$ $`\omega _0^2[1+\varphi _1(t)]x_1k(x_1x_2)\gamma \dot{x}_1`$
$`\ddot{x}_2`$ $`=`$ $`\omega _0^2[1+\varphi _2(t)]x_2k(x_2x_1)\gamma \dot{x}_2.`$ (1)
Here $`\omega _0`$ is the natural frequency of each (uncoupled) oscillator, $`k`$ is the coupling constant between them and $`\gamma `$ is the damping coefficient. The parametric modulations $`\varphi _1(t)`$ and $`\varphi _2(t)`$ are square waves with period $`T2\pi /\omega _p`$, amplitude $`A`$ and phase difference $`\theta `$:

$`\varphi _1(t)`$ $`=`$ $`A\text{sgn}\left(\mathrm{sin}(\omega _pt)\right)`$
$`\varphi _2(t)`$ $`=`$ $`A\text{sgn}\left(\mathrm{sin}(\omega _pt+\theta )\right).`$ (2)
With $`A=0`$ one is left with coupled ordinary damped harmonic oscillators whose total energy decays exponentially to zero. In the presence of the parametric modulation, however, energy is periodically pumped into the system, which may or may not lead to parametric resonance, i.e., to an infinite growth of the amplitudes of the oscillators. Our goal is to determine the boundaries between these two behaviors, which will be referred to as unstable and stable. As a point of reference we point the reader to the appendix where we review the results for a single parametric oscillator.

It is worth noting that a simple rescaling of time to dimensionless units, $`t^{}=t\omega _0`$, shows that Eqs. (2) are governed by the dimensionless parameter combinations $`r\omega _0/\omega _p`$, $`A`$, $`k/\omega _0^2`$, $`\gamma /\omega _0`$ and $`\theta `$. Moreover, the behavior of the system is invariant with respect to the replacement of $`\theta `$ by $`\theta ^{}2\pi \theta `$ (reflection around $`\pi `$) since this just amounts to an exchange of indices between the oscillators.

### 2.1 Floquet theory

The linearity of the equations allows one to make use of Floquet theory to solve the problem . Defining $`𝑿(x_1,\dot{x}_1,x_2,\dot{x}_2)^\mathrm{T}`$ where the superscript $`\mathrm{T}`$ denotes the transpose (i.e., $`𝑿`$ is a column matrix), one can rewrite Eqs. (2) in matrix form as $`\dot{𝑿}(t)=\widehat{D}(t)𝑿(t)`$ with a matrix $`\widehat{D}`$ satisfying $`\widehat{D}(t)=\widehat{D}(t+T)`$. Given a solution $`𝑿(t)`$ to Eqs. (2), the time periodicity of the parametric modulations thus implies that $`𝑿(t+T)`$ is also a solution. Therefore one can find a matrix $`\widehat{F}`$ such that $`𝑿(t+T)=\widehat{F}(T)𝑿(t)`$ and hence more generally through a repetition of this solution $`𝑿(t+nT)=\widehat{F}^n(T)𝑿(t)`$. The long time behavior of the system is thus clearly determined by the eigenvalues $`\{\lambda _i\}`$ of the Floquet operator $`\widehat{F}`$, which propagates the system in phase space for one period of the modulation. The eigenvectors $`𝒗_j`$ of $`\widehat{F}`$ satisfy $`|𝒗_j(t+nT)|=|\lambda _j|^n|𝒗_j(t)|`$, so in the limit $`n\mathrm{}`$ parametric resonance occurs if $`\mathrm{max}_j\{|\lambda _j|\}>1`$.

The model is appealing because for any number of oscillators (a single oscillator, or the coupled oscillators considered here, or even the mean field version of an infinitely cross-coupled infinite chain ) the piecewise constant parametric modulation leads to a piecewise linear system whose piecewise solution is known analytically. One can therefore construct the Floquet operator explicitly by simply multiplying together piecewise linear Floquet operators. This provides an immense reduction in computational effort by skipping what is usually the most time-consuming task in obtaining the regions of parametric resonance, namely, the numerical evaluation of the Floquet operator itself. For an isolated oscillator this procedure is well known  and reviewed in the appendix.

For the sake of clarity, let us first analyze the frictionless case $`\gamma =0`$. Each linear interval is characterized by one of the four possible states of the modulations $`(\varphi _1,\varphi _2)=(+A,+A)`$, $`(+A,A)`$, $`(A,+A)`$, or $`(A,A)`$. For each of these states, the solutions are of the form $`X_j(t)=A_je^{iPt}+B_je^{iMt}`$, where the eigenfrequencies $`P`$ and $`M`$ follow directly from the diagonalization of $`\widehat{D}`$:

$`P^2`$ $`=`$ $`{\displaystyle \frac{\omega _1^2+\omega _2^2+\sqrt{(\omega _1^2\omega _2^2)^2+4k^2}}{2}}`$
$`M^2`$ $`=`$ $`{\displaystyle \frac{\omega _1^2+\omega _2^2\sqrt{(\omega _1^2\omega _2^2)^2+4k^2}}{2}},`$ (3)
where $`\omega _{1,2}^2\omega _0^2[1+\varphi _{1,2}(t)]+k`$. Denoting the current state of the modulation by indices $`+`$ and $``$, notice that while $`P_{++}`$, $`M_{++}`$ and $`P_+=P_+`$ are always real, $`M_{}`$ becomes imaginary when $`A>1`$ while $`P_{}`$ and $`M_+=M_+`$ become imaginary when $`A>1+2k/\omega _0^2`$. Whether the system periodically alternates between harmonic behavior and one or more saddle nodes will thus depend on the parameter region and the phase difference $`\theta `$.

During a time interval $`\tau `$ with fixed $`\varphi _1,\varphi _2`$ we can relate $`𝑿(t+\tau )`$ to $`𝑿(t)`$ as $`𝑿(t+\tau )=\widehat{f}(\tau )𝑿(t)`$, where $`\widehat{f}`$ is the piecewise linear Floquet operator:

$$\widehat{f}(t)=\frac{1}{P^2M^2}\left(\begin{array}{cccc}m_1c_p+p_1c_m& m_1s_p+p_1s_m& k[c_p+c_m]& k[s_p+s_m]\\ m_1P^2s_pp_1M^2s_m& m_1c_p+p_1c_m& k[P^2s_pM^2s_m]& k[c_p+c_m]\\ k[c_p+c_m]& k[s_p+s_m]& m_2c_p+p_2c_m& m_2s_p+p_2s_m\\ k[P^2s_pM^2s_m]& k[c_p+c_m]& m_2P^2s_pp_2M^2s_m& m_2c_p+p_2c_m\end{array}\right)$$
(4)
where we make use of the shorthand notation $`c_p\mathrm{cos}(Pt)`$, $`s_pP^1\mathrm{sin}(Pt)`$, $`c_m\mathrm{cos}(Mt)`$, $`s_mM^1\mathrm{sin}(Mt)`$, $`m_{1,2}M^2\omega _{1,2}^2`$ and $`p_{1,2}P^2\omega _{1,2}^2`$. The Floquet operator is then finally obtained as the product of piecewise Floquet operators with arguments that depend on the phase difference between the modulations:

$$\widehat{F}(T)=\widehat{f}_+\left(\frac{T}{2}\frac{\theta }{\pi }\right)\widehat{f}_{}\left(\frac{T}{2}\left[1\frac{\theta }{\pi }\right]\right)\widehat{f}_+\left(\frac{T}{2}\frac{\theta }{\pi }\right)\widehat{f}_{++}\left(\frac{T}{2}\left[1\frac{\theta }{\pi }\right]\right).$$
(5)
The stability properties of the coupled system are determined by the magnitudes of the four eigenvalues $`\lambda _1,\mathrm{},\lambda _4`$ of the Floquet operator. Note that if $`\theta =0`$ or $`\theta =\pi `$, the above expression reduces to the product of only two matrices (as it should) since $`\widehat{f}(0)`$ is just the identity matrix. The addition of damping can be dealt with by noticing that if $`𝑿(t)`$ is a solution for $`\gamma =0`$, then $`𝒀(t)e^{\gamma t/2}\stackrel{~}{𝑿}(t)`$ is a solution for $`\gamma 0`$, where $`\stackrel{~}{𝑿}(t)=𝑿(t)|_{\omega _{1,2}^2\omega _{1,2}^2\gamma ^2/4}`$.

## 3 Results

In the absence of damping, Liouville’s theorem insures that $`det\widehat{f}=1`$, which is confirmed by cumbersome calculations and in turn implies $`\lambda _1\lambda _2\lambda _3\lambda _4=det\widehat{F}=1`$. In contrast with the single oscillator (see appendix), here this is not a sufficient condition to guarantee that the parametric resonance bifurcations occur at $`\lambda _j=\pm 1`$ : for coupled oscillators the largest eigenvalue can cross the unit circle $`|\lambda _j|=1`$ and hence enter a region of instability (parametric resonance) in other directions in the complex plane. The fourth order characteristic polynomial needs to be solved numerically to determine the instability boundaries; this is a computationally inexpensive procedure since we have an analytic expression for the Floquet operator.

### 3.1 In-phase modulations

The case $`\theta =0`$ corresponding to in-phase modulation $`\varphi _2(t)=\varphi _1(t)`$ simplifies sufficiently to yield analytic results for the stability boundaries. This is the simplest situation to be studied, since it is possible to reduce the original four-dimensional problem in this case to two known two-dimensional systems. This decoupling is accomplished by changing the coordinate system to the reference frame of the center of mass. Defining $`x(x_1+x_2)/2`$ and $`\rho (x_1x_2)/2`$, one obtains from Eq. (2)

$`\ddot{x}`$ $`=`$ $`\omega _0^2[1+\varphi _1(t)]x\gamma \dot{x}`$ (6)
$`\ddot{\rho }`$ $`=`$ $`\left\{\omega _0^2[1+\varphi _1(t)]+2k\right\}\rho \gamma \dot{\rho }.`$ (7)
Each of these equations is precisely that of a single parametric oscillator, one with frequency $`\omega _0`$ and modulation amplitude $`A`$, the other with rescaled parameters

$`\omega _0^{}`$ $``$ $`\sqrt{\omega _0^2+2k}`$
$`A^{}`$ $``$ $`A(\omega _0/\omega _0^{})^2.`$ (8)
We can thus use the known results for the single oscillator , for which a closed expression exists for the boundary of the instability regions in the $`(r,A)`$ plane (see the appendix).

For an undamped parametric oscillator with parameters $`(r,A)`$ parametric instabilities occur inside “tongues” that open up from integer and half-integer values of $`r`$ as discussed in the appendix and shown by the dashed curves in the first panel of Fig. 1. For the coupled problem in which we have transformed the problem to two independent parametric oscillators with different parameters, the instability regions are given by the overlap of the sets of tongues arising from each independent oscillator, one set emerging from integer and half-integer values of $`r`$ (“$`r`$-instability regions”), and the other from

$$r^{}(k)r\sqrt{1+2k/\omega _0^2}=\frac{\omega _0^{}}{\omega _p}$$
(9)
(“$`r^{}`$-instability regions”). When the coupling $`k`$ between the oscillators is small then the two sets of tongues almost overlap and one obtains in the undamped case the solid curves in the first panel, which show the instability boundaries when $`k/\omega _0^2=0.02`$. Still in the undamped problem, when the coupling $`k`$ is large the two sets of tongues occur at different parameter scales. The boundary diagram for $`k/\omega _0^2=3`$ is shown in lower left panel of the figure. In either case, since now we have two independent sets of instability regions the effect of the coupling $`k`$ for $`\theta =0`$ has been to enlarge the parametric resonance regime relative to that of two uncoupled oscillators.

This last result is completely general: the effect of coupling for in-phase parametric modulations is to enlarge the parametric resonance regime relative to the uncoupled case regardless of the form of the periodic modulation. The details of the phase boundaries will of course depend on the specific modulation function.

Damping, even in an isolated parametric oscillator, destroys much of the intricate boundary structure and increases the regions of stability. This is also the case for coupled oscillators, where increasing the damping tends to blur out the effects of coupling. These behaviors can be seen in the right column of Fig. 1. The regions of instability are now restricted to larger amplitudes $`A`$. Note the similarity between these two figures, which involve different couplings but now with substantial damping, $`\gamma /\omega _0=0.3`$. Note also that in the lower panel of this column (large coupling) the very first instability wedge is an $`r^{}`$-instability region while the other portions of the diagram (as well as the unstable regions in the upper weak-coupling case) include both $`r`$ and $`r^{}`$ instabilities.

For the particular case $`\theta =0`$ these results allow us to say something about the interesting problem of asymptotic synchronization. In regions where the center of mass coordinate $`x`$ is unstable but the relative coordinate $`\rho `$ is stable ($`r`$-instability regions that do not overlap with $`r^{}`$-instability regions) the coupled oscillators are synchronized if $`\gamma 0`$ ($`x_1=x_2`$), that is, the two oscillators move together about the origin with ever increasing amplitude. Conversely, if $`x`$ is stable but $`\rho `$ is unstable ($`r^{}`$-instability regions that do not overlap with $`r`$-instability regions), with $`\gamma 0`$ the oscillators become “antisynchronized” ($`x_1=x_2`$), that is, the two oscillators oscillate with ever increasing amplitude but in opposite directions, crossing one another each time they pass through the origin. Antisynchronization becomes more difficult to achieve with increasing coupling. If $`\gamma =0`$ the strict equalities $`x_1=x_2`$ or $`x_1=x_2`$ no longer hold in the non-overlapping regions, but the difference between $`x_1`$ and $`x_2`$ or $`x_2`$ is oscillatory and remains bounded. A simultaneous instability of both $`x`$ and $`\rho `$ involves an unstable center of mass coordinate and unbounded oscillations of each oscillator about this unbounded mean, which in turn involves a more complicated phase relation between the motions of the two oscillators. We return to this issue later.

### 3.2 Out-of-phase modulations

A wealth of very intricate results arises when $`\theta 0`$. In contrast with the $`\theta =0`$ case, it is now no longer clear how to break down the problem into simpler independent components (even for $`\theta =\pi `$) and, in particular, there is no longer a transparent way to relate the results to those of single parametric oscillators.

Perhaps the most striking feature of the $`\theta `$-dependence is its sensitivity: the unstable regions in the $`(r,A)`$ plane can change very abruptly (yet continuously) with the phase difference. As one might expect, this sensitivity is primarily modulated by the coupling $`k`$, as exemplified in Fig. 2 in the absence of damping. The left and right columns of Fig. 2 show the instability regions in the $`(r,A)`$ plane for $`k/\omega _0^2=0.05`$ and $`k/\omega _0^2=0.5`$, respectively. The evolution of the boundaries is shown as $`\theta `$ increases vertically along the columns (in order to save space, this and subsequent figures consistently omit labels, using the same scale as Fig. 1 in the $`(r,A)`$ plane, with the origin at the bottom left corner and tic marks separated by half a unit on both axes). While the most visible effects of the phase difference seem to concentrate on the low $`A`$ region for the smaller coupling (left column), the larger coupling (right column) induces a richer behavior, with (in)stabilities arising also for larger values of $`A`$ as $`\theta `$ changes. The phase difference can therefore create regions of stability and instability which are absent for in-phase modulations.

In order to address the effect of damping in this scenario, we have chosen to fix the coupling $`k`$ and tile several $`(r,A)`$ instability plots in such a way that $`\theta `$ varies in the vertical direction while $`\gamma `$ varies in the horizontal direction. Fig. 3 shows such a panel for a relatively small coupling $`k/\omega _0^2=0.12`$. One notices that the $`\theta `$ dependence is smoothed out as damping increases, so that almost no structure is visible along the rightmost column ($`\gamma /\omega _0=0.3`$). Fig. 4 shows the results for a larger coupling $`k/\omega _0^2=0.6`$. Note that the same tendency is observed: the instability boundaries become gradually less sensitive to $`\theta `$ as $`\gamma `$ increases. However, the rightmost column of Fig. 4 now shows much more structure than that of Fig. 3, indicating that for the same value of $`\gamma `$, the $`(r,A)`$ plane with larger $`k`$ shows a stronger $`\theta `$ dependence.

Figs. 3 and 4 suggest that $`\theta =\pi `$ gives rise to particularly stable behavior. This is in agreement with the vague intuitive notion that anti-phase modulations should be less prone to resonance than in-phase ones. In order to verify the extent to which this intuition is correct, a complementary view of these phenomena can be obtained by projecting the instability regions in the $`(k,\theta )`$ plane instead. This allows us to start with uncoupled oscillators ($`k=0`$) and observe how a given (in)stability evolves as $`k`$ and $`\theta `$ change. Indeed, it turns out to be possible to understand much of the behavior of the coupled system in terms of the behavior of the uncoupled system. To produce the representative results shown in Figs. 5 and 6, we fix several $`(r,A)`$ points in the single oscillator stability boundary diagram shown as black circles in the first panel of Fig. 1 and study the way in which variations in $`k`$ and $`\theta `$ affect these particular states. Two of the points in Fig. 1, $`(r,A)=(0.8,1.8)`$ and $`(1.25,0.6)`$ are stable states for the single oscillator (the first black dot touches the stability boundary in the figure only because it has been drawn large enough to render it visible; the point is well within the stable region). The other two points, $`(r,A)=(\mathrm{1.9.0.9})`$ and $`(2.25,1.3)`$ lead to unstable behavior of the single oscillator.

The first thing to be noted in Fig. 5 is that the horizontal axis has been rescaled in order to reveal the relevance of the variable $`r^{}`$ (see Eq. (9)). The top panels focus on $`r=1.9`$ and $`A=0.9`$, which is in the parametrically resonant regime for the uncoupled system. The small-$`r^{}`$ portion of the figure thus represents an unstable regime. One notes that for $`r^{}(k)`$ slightly above 2 the system becomes stable within a $`\theta `$ interval centered around $`\pi `$. Increasing $`k`$ a little further, this stability region then evolves in a very complex pattern, which includes reentrant “holes” of instability. After a distinguishable gap of instability, a somewhat simpler band of stability arises around $`r^{}3.5`$ starting at $`\theta =\pi `$. This band continues to larger values of $`k`$, with its outermost boundaries presenting a relatively simpler envelope than the low-$`k`$ pattern. The interesting point to be emphasized is that the band is perforated by gaps of instability most of which are centered precisely at integer and half-integer values of $`r^{}`$. This result is perhaps anticipated by the fact that the frequency $`\omega _0^{}`$ appears as the effective average diagonal frequency in the mean field equation of motion in ). The gaps are eventually closed by increasing the damping (top right panel), which also broadens the stability band and simplifies its dependence on the phase difference. The bottom panels in Fig. 5 are for $`r=2.25`$ and $`A=1.3`$, which again is in the resonant regime of the single oscillator. One notices the same pattern: for lower values of $`k`$, a complex shape emerges in the $`\theta `$ dependence of the stable regions. For sufficiently large values of $`k`$, a band of stability arises which has instability gaps basically centered at integer and half-integer values of $`r^{}`$. Damping (bottom right panel) causes gaps to disappear, creating a uniform region of stability centered around $`\theta =\pi `$.

Fig. 6 presents what may be regarded as the opposite situation, namely, when the original uncoupled system is stable. The top panels show the results for $`r=1.25`$ and $`A=0.6`$. Notice that the behavior is much simpler in this case, with the original stability being disturbed mostly around integer and half-integer values of $`r^{}`$ in the absence of damping (top left), with a relatively weak $`\theta `$ dependence. The effect of damping (top right) is to suppress most of these instability regions, yielding a predominantly stable $`(k,\theta )`$ plane. The bottom panels show results for an interesting intermediate situation: even though the uncoupled system is stable for $`r=0.8`$ and $`A=1.8`$, this point lies in a narrow corridor between two instability regions in the $`(r,A)`$ plane (see Fig. 1). One would therefore expect instabilities to arise more easily, and the immediate question is how the resulting diagram might reconcile the structures observed in Fig. 5 and the top panels of Fig. 6. The answer lies in a very rich structure in the $`(k,\theta )`$ plane (bottom left panel): initially, very small coupling induces an instability for all $`\theta `$. The now predominantly unstable system evolves in a manner similar to those of Fig. 5, a stability region being created around $`\theta =\pi `$ as $`k`$ increases, with string-like gaps of instability around the usual values of $`r^{}`$. The difference is that there is now also a second band of stability centered around $`\theta =0`$, presenting gaps at the same $`r^{}`$ positions. For $`k`$ sufficiently large, those two stability bands merge and one is left with a structure similar to that of the top panels: a predominantly stable region permeated by strings of instability. The bottom right panel of Fig. 6 shows the effect of damping, which significantly reduces the instability gaps thus greatly simplifying the picture.

Therefore the notion that modulations operating with a phase difference $`\theta =\pi `$ are less prone to parametric resonance is thus seen to be essentially correct. Framing this statement more carefully, our results show that for given $`r`$ and $`A`$, sufficiently large values of $`k/\omega _0^2`$ and $`\gamma /\omega _0`$ are able to induce, in the $`(k,\theta )`$ plane, a band of stability centered around $`\theta =\pi `$ even if the uncoupled oscillators are individually unstable. The width of this band can eventually comprise the whole $`2\pi `$ interval if the uncoupled system is originally stable.

## 4 Collective parametric instability

Bena and Van den Broeck  studied the stability boundaries of $`N`$ parametrically modulated oscillators $`\{x_i\}`$ each coupled to all the others by the same coupling constant $`2k/N`$ (in our notation). The phases $`\{\theta _i\}`$ are initially chosen at random from a uniform distribution in the interval $`[0,2\pi ]`$, remaining quenched thereafter. With square-block modulation the system is exactly solvable in the limit $`N\mathrm{}`$, where the mean-field solution becomes exact. The mean field equation is

$$\ddot{x}=\omega _0^2[1+\varphi _\theta (t)]x\gamma \dot{x}2k(xx)$$
(10)
where $`x`$ is the displacement of any oscillator in the chain, $`\varphi _\theta (t)`$ is the periodic modulation with phase $`\theta `$, and $`xN^1_{i=1}^Nx_i`$ is the mean displacement to be determined self-consistently. Bena and Van den Broeck note that the exact solution of this equation is

$$\left(\begin{array}{c}x(t)\\ \dot{x}(t)\end{array}\right)=𝐆_\theta (t)\left(\begin{array}{c}x(0)\\ \dot{x}(0)\end{array}\right)+2k𝐆_\theta (t)_0^t𝑑\tau 𝐆_\theta (\tau )^1\left(\begin{array}{c}0\\ x(\tau )\end{array}\right).$$
(11)
The propagator $`𝐆(t)`$ is known explicitly. Indeed, at $`t=T`$ it is $`𝐆(T)=e^{\gamma T}\widehat{F}(T)`$ where $`\widehat{F}`$ is precisely the single-oscillator Floquet operator given in Eqs. (A.9) and (A.10), but now with the frequencies shifted by the coupling constant

$$\omega _\pm =\sqrt{\omega _0^2(1\pm A)+2k\gamma ^2/4}.$$
(12)
Note that this is exactly the same as the propagator associated with the $`\rho `$ variable of Eq. (7) in the two-oscillator in-phase modulation problem, that is, the propagator associated with a single oscillator of frequency $`\omega _0^{}`$.

Bena and Van den Broeck identify two sorts of instabilities. One, which they call the “usual parametric resonance,” arises from the divergence associated with eigenvalues of $`𝐆`$ of magnitude greater than unity, that is, with the unbounded growth of the first term in Eq. (11) which in turn signals the unbounded growth of the amplitude of any typical oscillator in the chain. The stability boundaries associated with this type of instability are given precisely by Eq. (A.12) and are shown for the parameter choices indicated in the caption as the dotted curves in the top row panels of Fig. 7. In our in-phase two-oscillator parlance these are exactly the boundaries of the “$`r^{}`$-instability regions” defined in terms of the shifted frequency $`\omega _0^{}`$ (cf. Eq. (3.1)). The “usual” regions shrink in width and move toward lower $`r`$ and larger $`A`$ with increasing coupling $`k`$, a behavior already exhibited in the context of the in-phase two-oscillator results of Fig. 1. Indeed, this instability is beyond the scale of the figures in the large-coupling bottom row panels. The other type of instability, which they call a “collective instability,” is associated with the divergence of the mean $`x`$ and hence of the second term in Eq. (11). The collective instability boundaries are shown as dashed curves in all panels of Fig. 7. Note that the two types of instabilities may occur simultaneously, as seen in the instability region overlap evident in the top row panels of the figure. We return to this point below. Note also that with increasing coupling the system becomes increasingly stable, as one might expect, and that the unstable behavior becomes primarily collective.

We wish to explore whether our two-oscillator model is able to capture at least some of the behavior of the infinite system. In particular, we would like to investigate whether features of the collective instabilities of the mean field model are apparent in a system of only two oscillators with a fixed relative modulation phase. To make the comparison we also show in Fig. 7 the stability boundary for a single uncoupled oscillator (thin solid lines in left column) and for coupled oscillators with relative modulation phase $`\theta =\pi `$.

We make the following assertions: the two-oscillator system with any value of $`\theta `$ captures features of the overlap region of “usual” and “collective” instabilities. The purely “usual” regions are captured most accurately by the $`\theta =0`$ system, and the purely “collective” regime is best captured by the two-oscillator system with $`\theta =\pi `$. It is therefore this latter system that most fully captures (with unexpected detail) the principal features of collective behavior of the mean field model, and does so with increasing accuracy as the coupling between oscillators increases. We support these assertions, particularly the last one which is the one of most interest to us, with the results shown in Figs. 7 and 8.

Clearly, the $`\theta =0`$ system captures the full “usual” instability regime exactly since, as already stated, the “usual” instability is exactly the same as the “$`r^{}`$-instability.” This identity is not restricted to the square wave modulation but holds for any periodic modulation. However, the $`\theta =0`$ system does not capture the “collective” instability since the “$`r`$-instability” condition is that associated with a single parametric oscillator with the unshifted frequency $`\omega _0`$. Thus, for example, in the first panel of Fig. 7 the $`\theta =0`$ instability boundaries can be constructed from the combination of the dotted regimes and the thin solid lines (compare with the right lower panel of Fig. 1), whereas those of the mean field system include the same dotted regimes but now the very different dashed regions. The left lower panel shows an even greater difference between the $`\theta =0`$ two-oscillator model (whose “$`r`$-instability” boundaries are independent of $`k`$) and the mean field model (where the boundaries of the collective instability are sensitively dependent on $`k`$). Two coupled parametric oscillators with relative modulation phase $`\theta =0`$ therefore do not capture the collective features of the mean field model.

Our assertion that the two-oscillator system with any $`\theta `$ contains elements of the “usual” instabilities in the mean field model is simply a restatement of our earlier observation that “$`r^{}`$-instabilities” continue to appear even when one moves away from $`\theta =0`$ and all the way to $`\theta =\pi `$.

Consider now the two coupled oscillators with $`\theta =\pi `$. The stability boundaries are shown by the thick solid lines in the panels in the right column of Fig. 7. In the upper panel we observe that the first tongue approximates the region of “usual” mean field instability (in the non-overlapping regime) and that the remainder captures the collective instability boundary features surprisingly well (although it does miss the gap). We particularly point to the excellent fit of the leftmost boundary of this region. The agreement between the two models is even more dramatic in the lower panel, which corresponds to stronger coupling $`k`$. Again, the details are surprisingly well matched and the leftmost boundary of the region is captured essentially exactly.

To further support our analysis, and to gain a clearer understanding of the difference between “usual” and “collective” instabilities (which are both evidently already present in our two-oscillator system although the notion of a collective effect is not obvious in such a small system), we consider the motions that might characterize the instabilities. In the mean field system we conjecture that in the non-overlapping “usual” instability regime the mean is zero, $`x=0`$, but each oscillator oscillates about zero with ever increasing amplitude. The motion in the non-overlapping “collective” instability regime may involve an ever increasing mean with each oscillator oscillating about this moving mean with finite amplitude. This description is in accord with that of Bena and Van den Broeck . The overlap regions may involve both an increasing mean and oscillations of ever increasing amplitude about the moving mean. We have not ascertained these conjectures in the mean field system, but present results for the two-oscillator system that support this description.

Figure 8 shows trajectories for the two-oscillator $`\theta =\pi `$ system at the two points marked on the upper right panel of Fig. 7. The trajectories shown are those of each of the two oscillators as well as the mean trajectory. The upper panel is for parameter values in the unstable region that is not in the “collective” regime. It is tempting to associate this with the non-overlapping “usual” instability of the mean field model, an association that requires some caution. The figure indicates that not only does each oscillator and also the mean oscillate about zero, but all the trajectories, including the mean, appear to diverge. This behavior is that envisioned in our earlier discussion of the $`\theta =0`$ two-oscillator system in the regime where “$`r`$-instabilities” and “$`r^{}`$-instabilities” overlap, and is an indication that features of both kinds of instabilities persist even at $`\theta =\pi `$. We conjecture that “$`r`$-instabilities” are finite size divergences not present in the mean field model (and hence we would not expect to see this particular type of trajectory there). In the mean field limit only the “$`r^{}`$-instability” contributions in this regime persist, becoming the non-overlapping “usual” instability portion of the phase boundary portrait.

The lower panel of Fig. 8 is for parameter values in the non-overlapping “collective” regime. The mean indeed increases, and each oscillator oscillates about this increasing mean. It is interesting that the main feature of the “collective” instability, namely, an increasing mean displacement with individual oscillators oscillating about this mean, is already clearly captured by the two-oscillator model.

## 5 Conclusions

We have presented results for a model of two coupled parametric oscillators that can be solved for any value of the phase difference $`\theta `$ between the modulations. The other parameters of the problem are the oscillator frequency $`\omega _0`$, the modulation amplitude $`A`$ and period $`T=2\pi /\omega _p`$, the coupling $`k`$ and the damping $`\gamma `$.

We found that in-phase ($`\theta =0`$) coupled oscillator motion can be separated into two independent contributions. One involves the two oscillators moving together about zero (“$`r`$-instabilities”). The instability boundaries for this motion are identical to those of a single parametric oscillator of frequency $`\omega _0`$ and are independent of coupling $`k`$ since the spring connecting the oscillators is never disturbed. The other involves the two oscillators moving about zero but exactly in anti-phase with one another (“$`r^{}`$-instabilities”). The stability boundaries are sensitive to $`k`$ for these motions: the system becomes more stable with increasing coupling. We noted that these latter instabilities are exactly those identified as “usual” instabilities in the mean field model  and that they contribute to the instability boundaries in our two-oscillator model for any $`\theta `$, not just for $`\theta =0`$. We also showed that damping shrinks the instability regimes and smoothes the stability boundaries.

We showed that a change in $`\theta `$ can substantially modify the regions of parametric instability as projected in the usual $`(r=\omega _0/\omega _p,A)`$ plane, and that these changes are strongly affected by the coupling between the oscillators. In general, increasing $`\theta `$ up to $`\pi `$ provides greater stability but also leads to more intricate stability boundaries. An increase in $`k`$ and/or $`\gamma `$ also leads to increased stability. Alternatively, it is possible to understand the coupled system in terms of the uncoupled system by projecting the instability regions onto the $`(k,\theta )`$ plane. This rendition is particularly useful to show that $`\pi `$-centered bands of stability arise for sufficiently large $`k`$ and $`\gamma `$. We have thus identified all the trends of behavior in the two-oscillator model as each of the parameters is varied.

Our most interesting insights arise from a comparison of the two-oscillator results with the mean field model  consisting of $`N\mathrm{}`$ parametrically modulated oscillators all coupled to one another. In this latter system two types of instabilities have been identified: “usual” and “collective.” The “usual” instabilities are exactly our “$`r^{}`$-instabilities,” that is, the stability boundaries are identical. The interesting result is that the “collective” instabilities already appear in the two-oscillator model with $`\theta =\pi `$. This statement is based on the great similarity of the stability boundaries of the mean field and two-oscillator systems, especially with increasing coupling, and the specific features of the oscillator trajectories that typify the motions in each of these unstable regimes. It is perhaps surprising that a two-oscillator model can capture so much of the mean field collective behavior, and suggests that collective resonance in the latter may be dominated by phases quenched around $`\pi `$.

APPENDIX A - Parametric Resonance for a Single Oscillator
A brief review of parametric resonance for a single oscillator provides results called upon in the body of the paper for the coupled system . The equation of motion is

$$\ddot{x}+\omega ^2(t)x+\gamma \dot{x}=0.$$
(A.1)
The frequency $`\omega (t)`$ is a periodic function with period $`T2\pi /\omega _p`$. The simplest model of such a parametric modulation is perhaps the square wave

$$\omega ^2(t)=\omega _0^2\left[1+A\text{sgn}(\mathrm{sin}(\omega _pt))\right].$$
(A.2)
In the undamped case ($`\gamma =0`$) the oscillator frequency thus switches periodically between the two values $`\omega _\pm =\omega _0\sqrt{1\pm A}`$. If $`A1`$, the system jumps between two oscillatory regimes, each of which is stable. If $`A>1`$, then $`\omega _{}`$ becomes purely imaginary, corresponding to a saddle node behaviour in the phase space $`(x,\dot{x})`$. In this case, the system jumps between an unstable saddle and a stable oscillatory regime. As we shall see below, both instability and stability can occur for either $`A1`$ or $`A>1`$. With damping, the relevant frequencies are shifted (see below). The stability of the system depends on the interplay between the amplitude $`A`$ of the modulation and its frequency $`\omega _p`$ as well as the system parameters $`\omega _0`$ and $`\gamma `$.

The linearity of the equations leads naturally to the application of Floquet theory for the solution of the problem . A change of variables to

$$y(t)=e^{\gamma t/2}x(t)$$
(A.3)
yields the undamped equation

$$\ddot{y}+\stackrel{~}{\omega }^2(t)y=0$$
(A.4)
where

$$\stackrel{~}{\omega }^2(t)=\omega ^2(t)\frac{\gamma ^2}{4}.$$
(A.5)
For the square wave modulation we then have

$$\stackrel{~}{\omega }=\{\begin{array}{c}\stackrel{~}{\omega }_+\sqrt{\omega _0^2(1+A)\gamma ^2/4}\hfill \\ \stackrel{~}{\omega }_{}\sqrt{\omega _0^2(1A)\gamma ^2/4}\hfill \end{array}$$
(A.6)
Due to the periodicity of $`\stackrel{~}{\omega }(t)`$ the two independent solutions of Eq. (A.4) must also have period $`T`$. One can write for the original variable

$$𝐗(t+T)=e^{\gamma T}\widehat{F}(T)𝐗(t),$$
(A.7)
where

$$𝐗\left(\begin{array}{c}x\\ \dot{x}\end{array}\right)$$
(A.8)
and $`\widehat{F}`$ is the so-called Floquet operator that (aside from damping) propagates the system in phase space for one period:

$$\widehat{F}(T)=\widehat{f}_{}\left(T/2\right)\widehat{f}_+\left(T/2\right).$$
(A.9)
where

$$\widehat{f}_\pm (t)=\left(\begin{array}{cc}\mathrm{cos}\stackrel{~}{\omega }_\pm t+\frac{\gamma }{2\stackrel{~}{\omega }_\pm }\mathrm{sin}\stackrel{~}{\omega }_\pm t& \frac{1}{\stackrel{~}{\omega }_\pm }\mathrm{sin}\stackrel{~}{\omega }_\pm t\\ \frac{\omega _\pm ^2}{\stackrel{~}{\omega }_\pm }\mathrm{sin}\stackrel{~}{\omega }_\pm t& \mathrm{cos}\stackrel{~}{\omega }_\pm t\frac{\gamma }{2\stackrel{~}{\omega }_\pm }\mathrm{sin}\stackrel{~}{\omega }_\pm t\end{array}\right)$$
(A.10)
Equation (A.7) implies that

$$𝐗(t+nT)=e^{n\gamma T}\widehat{F}^n(T)𝐗(t),$$
(A.11)
so the long term behavior of the system is governed by the eigenvalues of $`e^{\gamma T}\widehat{F}`$. If any eigenvalue is larger than 1 in absolute value, there is parametric resonance and the system is unstable.

One can show that the two eigenvalues of $`e^{\gamma T}\widehat{F}(T)`$ obey the relation $`\lambda _1\lambda _2=1`$, reflecting the incompressible piecewise Hamiltonian flow. Therefore the eigenvalues are complex conjugates in the stable regions and real in the unstable regions. At the boundaries both eigenvalues equal either $`+1`$ or $`1`$, in which case the system oscillates with period $`T`$ or $`2T`$, respectively. Therefore the trace of the Floquet operator is $`\pm 2`$ at the stability boundaries, leading to the following boundary equation:

$$\mathrm{cos}\left(\frac{\stackrel{~}{\omega }_{}T}{2}\right)\mathrm{cos}\left(\frac{\stackrel{~}{\omega }_+T}{2}\right)\left(\frac{\stackrel{~}{\omega }_+^2+\stackrel{~}{\omega }_{}^2}{2\stackrel{~}{\omega }_+\stackrel{~}{\omega }_{}}\right)\mathrm{sin}\left(\frac{\stackrel{~}{\omega }_{}T}{2}\right)\mathrm{sin}\left(\frac{\stackrel{~}{\omega }_+T}{2}\right)=\pm \mathrm{cosh}\left(\frac{\gamma T}{2}\right).$$
(A.12)
In the absence of damping the solutions depend only on $`A`$ and the ratio $`r\omega _0/\omega _p`$, and these particular results are shown by the dashed curves in the first panel of Fig. 1. The regions of instability appear as “tongues” starting from integer and half-integer values of $`r`$ at small modulation amplitudes. Note that the stability boundaries surrounding the tongues that emerge from integer values of $`r`$ correspond to both Floquet eigenvalues being equal to $`1`$, while those surrounding tongues that emerge from half-integer values of $`r`$ are associated with the Floquet eigenvalues equal to $`1`$. Note that no abrupt transition (nor a visible transition of any kind) is seen at the line $`A=1`$ which marks the transition between an oscillating-oscillating behavior and a saddle-oscillating one. As one would expect, on the other hand, the system is nearly always unstable for sufficiently small $`\omega _p`$ and sufficiently large $`A`$ (upper right corner of the figure).

#### Acknowledgements:

We would like to thank C. Van den Broeck and R. Kawai for stimulating discussions. This work was supported in part by the National Science Foundation under grant No. PHY-9970699.