text:
# Resonantly Forced Inhomogeneous Reaction-Diffusion Systems

## I Introduction

Patterns in reaction-diffusion systems where the kinetics are spatiotemporally modulated can display a variety of phenomena that are not found in homogeneous systems. Many reaction-diffusion processes of practical interest take place in inhomogeneous media or may be coupled to external processes that affect the kinetics in a non-uniform manner. A convenient experimental system for studying the effect of spatiotemporal modulations on pattern dynamics in spatially distributed systems is the ruthenium-catalyzed Belousov-Zhabotinsky (BZ) reaction.  This reaction is light-sensitive and thus the kinetics may be modulated by projecting a pattern of illumination of varying intensity onto the reaction medium.

Recent studies have made use of the light-sensitive BZ system to investigate wavefront propagation in systems with spatially disordered excitability. Kádár et al. studied stochastic resonance in a system with a periodically regenerated noise pattern. Sendiña-Nadal et al. studied percolation and roughening of wavefronts in systems with quenched spatially disordered excitability;  excitable spiral waves in the presence time-varying disorder were also studied. The dynamics of reaction-diffusion waves in inhomogeneous excitable media is thought to be relevant to cardiac fibrillation since electrical waves may be disrupted by irregularities in the heart muscle medium. Noise is thought to play a role in initiation and propagation of waves in neural tissue. Earlier investigations into spatially inhomogeneous excitable systems include a study of a BZ medium containing catalyst-coated resin beads which served as nucleation sites for wavefronts , numerical simulations of a spatially-distributed network of coupled excitable elements which exhibited spontaneous wave initiation, stochastic resonance and fragmentation of wavefronts,  and an excitable cellular automaton in which the refractory times of the elements were assigned randomly.  The effect of stochastic spatial inhomogeneities on other types of reaction-diffusion systems has been much less studied.

Periodically forced reaction-diffusion systems have also been investigated. Petrov et al. and Lin et al. subjected an oscillatory version of the light-sensitive BZ reaction to periodic spatially uniform illumination.  As the ratio of the forcing frequency to the natural frequency neared various resonances patterns were observed. In subsequent numerical studies the observed transition between labyrinthine and non-labyrinthine two-phased patterns at the 2:1 resonance was reproduced in the periodically forced Brusselator.  Belmonte et al. observed a transition from a stable spiral to turbulence in a BZ reaction when resonant forcing was applied.  Resonantly forced oscillatory systems have been investigated by numerical simulation as well as theoretically by Coullet and Emilsson,  Coullet et al.,  Elphick et al.,  and Chaté et al.  These systems exhibit a number of interesting pattern-forming phenomena.

Given the wide range of phenomena arising from spatial disorder in excitable systems and resonant forcing of oscillatory systems, one might expect that oscillatory systems with spatially inhomogeneous forcing have the potential to exhibit interesting new features. The research presented here explores qualitatively the phenomenology of such systems and characterizes some of the phenomena quantitatively. We restrict the focus of our study to systems in two spatial dimensions.

## II Periodically Forced Oscillatory Systems

Consider an externally forced oscillatory reacting system described by the ordinary differential equation,

$$\frac{d𝐜(t)}{dt}=𝐑(𝐜(t);𝐚,𝐛(t)),$$
(1)
where $`𝐜(t)`$ is a vector containing the concentrations of reagents. The reaction rates are described by the nonlinear vector function $`𝐑`$ which depends on a collection of parameters $`𝐚`$, such as rate constants and constant concentrations of pool chemicals, as well as parameters $`𝐛(t)`$ which comprise the periodic forcing and are of the form $`𝐛(t)=𝜼_0\mathrm{\Phi }(\omega _\mathrm{f}t)`$ with $`𝜼_0`$ the constant forcing amplitude and $`\mathrm{\Phi }`$ a $`2\pi `$-periodic function giving the form of the forcing. As mentioned above, such forcing may be implemented by periodic illumination of the system if the reaction is light sensitive.

If $`𝐛(t)=\mathrm{𝟎}`$ we suppose the unforced reacting system has a stable limit cycle $`𝐜_0(t)`$ with period $`T_0=2\pi /\omega _0`$. In such a system there exists an infinite number of limit cycle solutions, $`𝐜_0^{}(t)=𝐜_0(t+\mathrm{\Delta }t)`$ which differ from $`𝐜_0(t)`$ only by an arbitrary phase shift $`2\pi \mathrm{\Delta }t/T_0`$. Limit cycle attractors are neutrally stable to phase perturbations corresponding to translations along the orbit. A system following a limit cycle will, in general, have undergone a phase shift when it returns to the limit cycle after experiencing a small perturbation.

These characteristics of the unforced oscillator may be contrasted with those of the forced oscillator.  If $`\omega _\mathrm{f}/\omega _0`$ is sufficiently close to an irreducible ratio of integers $`n/m`$, and if the forcing amplitude is sufficiently large, then the oscillations may become entrained to the external forcing and the system possesses $`n`$ stable limit cycle solutions of period $`T=nT_\mathrm{f}=2n\pi /\omega _\mathrm{f}mT_0`$ which are mapped into each other under phase shifts $`tt+kT/n`$ for $`k=0,1,2,\mathrm{}`$. A system following one of these limit cycles will return to it with no phase shift after a small perturbation. This discrete, finite collection of limit cycles may be contrasted with the infinite and continuous collection of limit cycles in the unforced case. The entrained resonantly forced oscillator is a system with $`n`$ stable states, defined by the phase of the oscillations rather than by the system’s location in phase space.

The general form of an oscillatory reaction-diffusion system with spatially inhomogeneous periodic forcing is

$$\frac{c(𝐫,t)}{t}=𝐑(𝐜(𝐫,t);𝐚,𝐛(𝐫,t))+𝐃^2𝐜(𝐫,t),$$
(2)
where $`𝐃`$ is a diagonal matrix of diffusion coefficients. The parameters responsible for the periodic forcing $`𝐛(𝐫,t)`$ now depend on space as well as time and are of the form $`𝐛(𝐫,t)=𝜼(𝐫,t)\mathrm{\Phi }(\omega _\mathrm{f}t)`$. The random variable $`𝜼(𝐫,t)`$ accounts for the fact that the forcing amplitude may vary stochastically in space and time.

In a spatially distributed system with spatially uniform forcing, $`𝐛(𝐫,t)=𝐛(t)=𝜼_0\mathrm{\Phi }(\omega _\mathrm{f}t)`$, the diffusive coupling and the stability of the $`n`$ limit cycles to phase perturbations leads to the formation of domains of the different phase locked states. At the domain walls separating the different phase locked states the phase of the oscillations shifts by an amount determined by the character of the phase locking. In two dimensions, $`n`$-armed spiral waves may form, given suitable initial conditions, if the domain walls have non-zero velocity. The core of the spiral wave is a phase singularity at which all $`n`$ states meet and around which the phase advances by $`2m\pi `$. The rotation of these phase locked spirals is a result of the propagation of the phase dislocations comprising the domain walls; thus, they rotate much more slowly than spiral waves in the unforced system, where the spiral rotation frequency is equal to the frequency of the local oscillations.

In this article we investigate systems subjected to periodic forcing with a stochastic component, either quenched disorder, $`𝜼(𝐫,t)=𝜼(𝐫)`$, where the periodically applied perturbation has a random distribution in space which does not change in the course of the evolution, and systems with dynamic disorder, $`𝜼(𝐫,t)`$, where the spatial distribution of the perturbation may change with time. The reduction of a resonantly forced oscillatory system to a complex Ginzburg-Landau (CGL) equation by Gamabaudo  and by Elphick et al  may be simply extended to systems with quenched disorder and consequently we have employed the CGL equations in our studies of this case. For time-varying disorder, one must reconsider the derivation of the appropriate CGL equation because of the presence of another time scale associated with the noise process; hence, we have chosen instead to study the FitzHugh-Nagumo system as a typical example of an oscillatory reaction-diffusion system.

In such cases a new length, $`\mathrm{}_s`$, related to the spatial correlation range of the noise distribution enters the problem. The behavior of the system will depend on the magnitude of this length relative to that of other important lengths in the system, such as the diffusion length $`\mathrm{}_D`$ and the typical width of a domain wall or propagating front, $`w_d`$. If $`\mathrm{}_s`$ is sufficiently small compared to both of these lengths the system will appear effectively homogeneous and behave as if it were subject to periodic forcing with amplitude determined by the spatial average $`\overline{𝜼}`$ of $`𝜼(𝐫)`$. If $`\mathrm{}_s`$ is very large compared to these lengths then the system dynamics may be simply represented in terms of the dynamics of a collection of uniformly forced patches. The interesting regime is when these length scales are comparable and our studies focus on these cases.

## III Forced CGL with quenched disorder

In the vicinity of the Hopf bifurcation point a periodically forced oscillatory reaction-diffusion system may be reduced to its normal form, the forced complex Ginzburg-Landau equation (FCGL) . This reduction is valid provided the system is near an $`n:m`$ resonance, where $`n=1,2,3,4`$. Such a reduction may also be carried out for the case of quenched disorder, $`𝜼(𝐫,t)=𝜼(𝐫)`$, and the FCGL takes the form,

$`{\displaystyle \frac{A(𝐫,t)}{t}}`$ $`=`$ $`(\mu +i\nu )A(1+i\beta )|A|^2A`$ (4)
$`+\gamma (𝐫)\overline{A}^{n1}+(1+i\alpha )^2A,`$
where $`\gamma (𝐫)`$ accounts for the spatial variation of the forcing amplitude. The complex amplitude $`A(𝐫,t)`$ describes slow modulations in the frequency and amplitude of oscillations of the original system in Poincaré planes taken at the forcing period; $`\overline{A}`$ denotes the complex conjugate of $`A`$. Equation (4) with $`\gamma (𝐫)\gamma _0`$, a constant, has been used as a model for an oscillatory reaction-diffusion system with spatially uniform resonant forcing.

Consider the normal form of a single resonantly forced oscillator obtained by omitting the diffusion terms and spatial dependence in $`\gamma (𝐫)`$,

$$\frac{dA(t)}{dt}=(\mu +i\nu )A(1+i\beta )|A|^2A+\gamma \overline{A}^{n1}.$$
(5)
For Eq. (5) a critical $`\gamma _c`$ exists such that for $`\gamma <\gamma _c`$ the equation exhibits a stable limit cycle solution, while for $`\gamma \gamma _c`$ there are $`n`$ stable fixed points. These correspond to the $`n`$ stable limit cycle solutions of the original system and can be mapped into each other by phase shifts $`AAe^{i2\pi /n}`$. More specifically, the fixed points can be found from the stationary solutions of Eq. (5) by expressing the complex amplitude in the form $`A=Re^{i\varphi }`$. The modulus, $`R_0`$, of the non-zero fixed points of Eq. (5) depends on $`\gamma `$ according to

$$\gamma ^2=\frac{(R_0^2\mu )^2+(\nu \beta R_0^2)^2}{R_0^{2(n2)}}.$$
(6)
For some values of $`n`$ and $`\gamma `$ Eq. (6) permits multiple $`R_0`$ values, some of which may correspond to unstable fixed points. Figure 1 shows a $`\gamma `$ vs. $`R_0`$ curve for $`n=3`$; the upper branch corresponds to the stable fixed points while the lower branch describes the nonzero unstable fixed points. The parameter values used to construct this figure are $`\mu =1`$, $`\nu =0`$, $`|\beta |=0.6`$. Note the presence of a critical forcing amplitude $`\gamma _c0.58`$ below which phase locking does not occur. For $`n=3`$, Eq. (6) gives $`\gamma _c=\left[2\left((1+\beta ^2)(\mu ^2+\nu ^2)\right)^{\text{“sevenrm1}/\text{“sevenrm2}}2(\mu +\beta \nu )\right]^{\text{“sevenrm1}/\text{“sevenrm2}}`$.

As an example of quenched disorder in resonantly forced oscillatory systems the $`\gamma (𝐫)`$ fields were taken to be dichotomous random variables. Two values for the forcing amplitude, $`\gamma _1`$ and $`\gamma _2`$, were chosen. The two-dimensional system was partitioned into square cells and the value of $`\gamma (𝐫)`$ in each cell was chosen to be either $`\gamma _1`$, with probability $`p`$, or $`\gamma _2`$ with probability $`q=1p`$. More precisely, if the noise cells have dimension $`s\times s`$ and the system’s dimensions are $`W\times L=sN_W\times sN_L`$ then

$$\gamma (𝐫)=\underset{i=1}{\overset{N_W}{}}\underset{j=1}{\overset{N_L}{}}\xi _{ij}\mathrm{\Theta }_{ij}(𝐫),$$
(7)
where

$$\xi _{ij}=\{\begin{array}{cc}\hfill \gamma _1\text{with probability}& p\hfill \\ \hfill \gamma _2\text{with probability}& q=1p,\hfill \end{array}$$
(8)
and

$`\mathrm{\Theta }_{ij}(𝐫)`$ $`=`$ $`\mathrm{\Theta }_{ij}(x,y)`$ (9)
$`=`$ $`\theta (x(i1)s)\theta (isx)`$ (11)
$`\times \theta (y(j1)s)\theta (jsy),`$
where $`\theta `$ is the Heaviside function and $`(ij)`$ are the discrete coordinates of a noise cell. For quenched disorder this distribution is fixed for all time. The cell size, the values of $`\gamma _1`$ and $`\gamma _2`$ and the seeding probabilities $`p`$ and $`q`$ are the relevant parameters to consider. The probabilities $`p`$ and $`q`$ were typically, but not necessarily, independent of position; exceptions will be noted as they occur. This $`\gamma (𝐫)`$ field has the mean value $`\overline{\gamma }=p\gamma _1+q\gamma _2`$, and spatial autocorrelation

$`C(𝐫^{})`$ $`=`$ $`{\displaystyle \frac{\delta \gamma (𝐫^{}+𝐫^{\prime \prime })\delta \gamma (𝐫^{\prime \prime })}{\delta \gamma (𝐫^{\prime \prime })\delta \gamma (𝐫^{\prime \prime })}},`$ (12)
$`=`$ $`\{\begin{array}{cc}\left(1\frac{|x^{}|}{s}\right)\left(1\frac{|y^{}|}{s}\right)& \mathrm{i}f|x^{}|s\mathrm{a}nd|y^{}|s\\ & \\ 0& \mathrm{o}therwise\end{array}`$ (16)
where $`\delta \gamma (𝐫)=\gamma (𝐫)\overline{\gamma }`$, $`𝐫^{}=(x^{},y^{})`$ and the average $``$ is taken over all $`𝐫^{\prime \prime }`$.

The studies described in this paper investigate patterns at the 3:1 resonance. For such systems there are two interesting cases for resonant forcing with a dichotomous $`\gamma (𝐫)`$ field. In the first case both $`\gamma _1`$ and $`\gamma _2`$ lie above the phase locking threshold, $`\gamma _c`$, so that all regions of the medium are entrained to the forcing. In the second case only one of $`\gamma _1`$ or $`\gamma _2`$ lies above the threshold and the medium consists of a mixture of entrained and non-entrained regions. In the simulations described below, including those in Sec. IV, numerical integration was performed using explicit forward differencing and a second-order discrete Laplacian.

### A All sites phase locked; front roughening

When both $`\gamma _1`$ and $`\gamma _2`$ lie above the phase locking threshold $`\gamma _c`$ then all regions of the medium are tristable. The system possesses three quasi-homogeneous stationary states in which the complex amplitude fluctuates about an average value. Domain walls separating these phase locked states are in general non-stationary in the $`3:1`$ resonant regime since all three phases are inequivalent, although their velocity may pass through zero as parameters are tuned. Initially planar fronts in these inhomogeneously forced systems roughen as they propagate. Figure 2 shows an example of a rough interface separating two of the three phases.

Front roughening is also observed if $`\gamma _1<\gamma _c<\gamma _2`$, (i.e. when the medium consists of a mixture of entrained and non-entrained regions) but is difficult to study because spontaneously nucleated patterns interfere with propagating fronts. This case will be examined below.

The propagating fronts in this system experience local velocity fluctuations arising from spatial variations in the $`\gamma (𝐫)`$ field. Diffusion will tend to eliminate front roughness generated by random fluctuations in the $`\gamma (𝐫)`$ field; consequently, the front dynamics should obey the Kardar-Parisi-Zhang (KPZ) equation,

$$\frac{h(x,t)}{t}=\overline{v}+D\frac{^2h}{x^2}+\frac{\lambda }{2}\left(\frac{h}{x}\right)^2+\zeta (x,t),$$
(17)
where $`h(x,t)`$ is the front position, $`\overline{v}`$ is the average front velocity, $`D`$ and $`\lambda `$ are phenomenological coefficients and $`\zeta (x,t)`$ is Gaussian white noise with zero mean and correlation $`\zeta (x^{},t^{})\zeta (x^{\prime \prime },t^{\prime \prime })=2\mathrm{\Gamma }\delta (x^{}x^{\prime \prime })\delta (t^{}t^{\prime \prime })`$. In such a circumstance the interface width $`w(t)=(L^1_i(h(x_i,t)\overline{h}(t))^2\mathrm{\Delta }x)^{\text{“sevenrm1}/\text{“sevenrm2}}`$ increases with time as $`w(t)t^{\widehat{\beta }}`$ for short times; the average width of a saturated front, $`w_\mathrm{s}`$, scales with system size as $`w_\mathrm{s}L^{\widehat{\alpha }}`$, where $`\widehat{\alpha }=1/2`$ and $`\widehat{\beta }=1/3`$.

We have verified these KPZ scaling properties for a FCGL system with parameter values ($`\alpha =1`$, $`\beta =0.6`$, $`\mu +i\nu =1`$) and the forcing field parameters ($`\gamma _1=0.60`$, $`\gamma _2=1`$, $`p=0.50`$) for which the critical forcing amplitude is $`\gamma _c0.58`$ (cf. Fig. 1). The front properties were measured in a frame moving with the front. The system dimensions were $`W=100`$ and $`L=100`$, $`200`$, $`400`$ and $`800`$. The noise grain size was $`s\times s=W/25\times W/25`$. The boundary conditions were periodic along the boundaries perpendicular to the front, and no-flux along the boundaries perpendicular to the direction of front motion.

The saturated front width $`w_\mathrm{s}`$ was found to scale with system size as $`w_\mathrm{s}L^{1/2}`$. Plotting $`w(t)/L^{\widehat{\alpha }}`$ against $`t/L^{\widehat{\alpha }/\widehat{\beta }}`$ collapses the $`w(t)`$ versus $`t`$ data for four different $`L`$ onto a single curve when the KPZ values $`\widehat{\alpha }=1/2`$, $`\widehat{\beta }=1/3`$ are used, as shown in Fig. 3. Here $``$ denotes an average over realizations.

The width temporal autocorrelation function

$$C_w(t)=\frac{\delta w(t+t^{})\delta w(t^{})}{\delta w(t^{})\delta w(t^{})},$$
(18)
was calculated for fronts in the saturated regime at four system sizes $`L=100`$, $`200`$, $`400`$ and $`800`$. Here $``$ represents a time average within the saturated regime. Temporal correlations were found to decay to zero, with the rate of decay decreasing as system size increases (Fig. 4). This is expected, since the larger system size allows larger fluctuations in the front profile which require longer times to form and decay.

To investigate the presence of spatial correlations in saturated front profiles, the height spatial autocorrelation function

$$C_h(x)=\frac{\delta h(x+x^{})\delta h(x^{})}{\delta h(x^{})\delta h(x^{})},$$
(19)
was calculated for the various system sizes. In Eq. (19) the inner $``$ refer to averaging over $`x^{}`$ in a single front profile and the outer $``$ refers to averaging over different saturated front profiles. Results for $`L=100`$ and $`L=800`$ are shown (Fig. 5).

In the case of a front free of spatial structure, i.e. a periodic random walk that returns to its initial position in $`N`$ steps, the height spatial autocorrelation function is of the form $`C_h^R(x)=1(2𝒟/L)x(Lx)`$, where $`𝒟`$ is a phenomenological diffusion coefficient.  The deviations of the fit of the numerically determined $`C_h(x)`$ from a quadratic function indicate the existence of spatial correlations.

### B A medium with phase locked and oscillatory sites; spontaneous nucleation of target patterns

When $`\gamma _1<\gamma _c`$ and $`\gamma _2>\gamma _c`$, the system consists of a mixture of tristable regions and oscillatory regions. If the density of oscillatory sites, $`p`$, is low the diffusive coupling maintains the value of $`A`$ within the oscillatory regions near that in the adjacent tristable regions. The medium behaves essentially like a tristable medium and supports three quasi-homogeneous stable stationary states, travelling kink-like domain walls and three-armed spiral waves. The oscillatory sites provide a spatial inhomogeneity that leads to roughening of domain walls and to fluctuations of the concentrations within domains.

With increasing $`p`$, target patterns are observed (Fig. 6). They consist of concentric, approximately circular domain walls moving outwards from a central region (a “pacemaker”) where they are initiated periodically. Within each ring of the target the concentration is quasi-homogeneous. Note that all images in Fig. 6 are from realizations with identical parameters, however a range of wavelengths is observed.

The probability of a realization possessing a target pattern was measured as a function of $`p`$ and system size (Fig. 7). It was found to approach zero for low $`p`$, one for high $`p`$, and to increase rapidly around some critical $`p_c`$. Figure 7 shows results for the parameter values ($`\mu +i\nu =1`$, $`\alpha =1`$, $`\beta =0.6`$) and forcing field parameters ($`\gamma _1=0`$, $`\gamma _2=1`$); qualitatively similar results were also obtained for ($`\mu +i\nu =1`$, $`\alpha =1`$, $`\beta =0.6`$, $`\gamma _1=0`$, $`\gamma _2=1`$) and ($`\mu +i\nu =0`$, $`\alpha =2`$, $`\beta =0.82`$, $`\gamma _1=0`$, $`\gamma _2=1`$). For larger system sizes, the probability of occurrence of a target pattern is higher, consistent with a fixed probability per unit area of the medium nucleating a target pattern.

The occurrence of target patterns may be explained by supposing that diffusion provides a mechanism for averaging $`\gamma (𝐫)`$ over some length scale, and that regions of the medium behave qualitatively like a uniform medium with the same local average $`\overline{\gamma }(𝐫)`$. Thus the medium is locally tristable for a high $`\overline{\gamma }(𝐫)`$, but if $`\overline{\gamma }(𝐫)`$ is sufficiently low it will be oscillatory. As the density of oscillatory sites increases, it becomes increasingly likely that there will exist regions with a local $`\gamma `$ sufficiently low to exhibit oscillatory dynamics.

We can determine a profile of the average $`\gamma `$ field around the average pacemaker. We define $`\overline{\gamma }_\mathrm{p}`$ as the average value of $`\gamma `$ a distance $`R`$ from the center of a pacemaker. This is given by

$$\overline{\gamma }_\mathrm{p}(R)=\frac{1}{2\pi R}\frac{d}{dR}\pi R^2\overline{\gamma }_\mathrm{d}(R),$$
(20)
where

$$\overline{\gamma }_\mathrm{d}(R)=\frac{1}{\pi R^2}_{|𝐫𝐫_0|R}\gamma (𝐫)𝑑𝐫,$$
(21)
is the average value of $`\gamma `$ over a disk of radius $`R`$, averaged over all pacemakers in all realizations with a given set of parameters. In practice, $`\overline{\gamma }_\mathrm{p}`$ is calculated as the discrete derivative

$$\overline{\gamma }_\mathrm{p}(R)=\frac{\pi (R+\mathrm{\Delta }R)^2\overline{\gamma }_\mathrm{d}(R+\mathrm{\Delta }R)\pi R^2\overline{\gamma }_\mathrm{d}(R)}{\pi (R+\mathrm{\Delta }R)^2\pi R^2}.$$
(22)
Figure 8 shows $`\overline{\gamma }_\mathrm{p}`$ and $`\overline{\gamma }_\mathrm{d}`$ for $`p=0.30`$, system size $`L\times L=100\times 100`$, noise grain size $`s\times s=L/200\times L/200`$, and other parameter values equal to those used to measure the nucleation probability curves shown in Fig. 7. Qualitatively similar plots, in which $`\overline{\gamma }_\mathrm{p}`$ and $`\overline{\gamma }_\mathrm{d}`$ increase from a low value at $`R=0`$ to the mean-field value, $`\overline{\gamma }`$, at high $`R`$ were observed for other values of $`p`$. Also shown in Fig. 8 are $`\overline{\gamma }_\mathrm{p}`$ calculated from the neighborhoods of points chosen randomly from the realizations rather than from points centered on pacemakers, and the mean-field value, $`\overline{\gamma }`$.

Consider the following simple model of pacemaker formation. We divide the system into $`𝒩`$ “sites” of size $`\mathrm{\Delta }x\times \mathrm{\Delta }x`$ such that the value of $`\gamma `$ is constant over a site; i.e., the sites may be noise domains or subdivisions of noise domains. We assume that, on average, pacemakers have a radius of $``$ and contain $`𝒩_p=\pi (/\mathrm{\Delta }x)^2`$ sites. Thus we divide the system into $`𝒩/𝒩_p`$ independent domains which are potential pacemakers. We assume that such a domain is a pacemaker if the average value of $`\gamma `$ within it is less than some threshold for oscillatory behavior, $`\gamma ^{}`$. If $`N_1`$ is the number of sites on which $`\gamma =\gamma _1`$ and $`N_2=𝒩_pN_1`$ is the number of $`\gamma _2`$ sites, then the potential pacemaker is a pacemaker when

$$\frac{\gamma _1N_1+\gamma _2N_2}{𝒩_p}=\frac{\gamma _1N_1+\gamma _2(𝒩_pN_2)}{𝒩_p}\gamma ^{},$$
(23)
that is, when

$$N_1\frac{\gamma _2\gamma ^{}}{\gamma _2\gamma _1}𝒩_p.$$
(24)
Defining

$$N^{}=\frac{\gamma _2\gamma ^{}}{\gamma _2\gamma _1}𝒩_p,$$
(25)
the probability that a domain is a pacemaker is

$$𝒫=\underset{k=N^{}}{\overset{𝒩_p}{}}\left(\genfrac{}{}{0pt}{}{𝒩_p}{k}\right)p^kq^{𝒩_pk},$$
(26)
and it is not a pacemaker with probability $`𝒬=1𝒫`$. The probability that the entire system possesses no pacemakers is $`𝒬^{𝒩/𝒩_p}`$ and it possesses one or more pacemakers with probability

$$P(𝒩,p)=1𝒬^{𝒩/𝒩_p}.$$
(27)
The parameters $`𝒩_p`$ and $`\gamma ^{}`$ uniquely determine the $`P(𝒩,p)`$ surface. The values giving the best fit to the experimental data were $`\gamma ^{}=0.65`$ and $`𝒩_p=180`$ which implies $`3.8`$. These values are compared with the numerically determined $`\overline{\gamma }_\mathrm{d}(R)`$ curve in Fig. 8. The simple model predicts that the average pacemaker should have $`\overline{\gamma }_\mathrm{d}()=\gamma ^{}`$; the point $`(,\gamma ^{})`$ lies close to, but not on, the $`\overline{\gamma }_\mathrm{d}`$ curve.

In order to provide further insight into the criteria necessary for a region to act as a pacemaker, a series of studies was carried out in which the $`\gamma `$ field consisted of a disk of radius $`R`$ sites with a density $`p_{\mathrm{i}n}`$ of oscillatory sites. This disk was embedded in a field with density of oscillatory sites $`p_{\mathrm{o}ut}`$. In all cases $`p_{\mathrm{i}n}`$ was greater than $`p_{\mathrm{o}ut}`$, so that the disk could act as a pacemaker. Multiple realizations of the evolution were simulated at various values of $`R`$, $`p_{\mathrm{i}n}`$ and $`p_{\mathrm{o}ut}`$, and the fraction of realizations Pr(br) in which the central disk emitted target waves into the surrounding medium (i.e. in which “breakout” occurred) was measured.

Figure 9 shows Pr(br) as a function of $`R`$ and $`p_{\mathrm{i}n}`$ for $`p_{\mathrm{o}ut}=0.10`$. As one would expect, Pr(br) increases as $`R`$ increases and decreases as $`p_{\mathrm{i}n}`$ decreases, with the exception that at $`p_{\mathrm{i}n}=1`$ the probability of breakout increases with $`R`$ until $`R=3`$ and then decreases for $`3<R<5`$, after which it increases. (cf. Figs. 9 and 10). For $`p_{\mathrm{i}n}=0.95`$, when there is a small fraction of tristable sites in the inner disk, the magnitude of the decrease is substantially reduced (cf. Fig. 10). Apart from the anomalous region at low $`p_{\mathrm{i}n}`$ the trends that Pr(br) increases with increasing $`R`$ and decreases with decreasing $`p_{\mathrm{i}n}`$ are consistent with the notion that pacemakers form when the local density of oscillatory sites is high. The behavior for $`p_{\mathrm{o}ut}=0.15`$ was qualitatively similar.

The anomalous decrease between $`R=3`$ and $`R=5`$ is possibly related to the fact that the anomalous behavior begins when $`R\mathrm{}_D`$, where $`\mathrm{}_D`$ is the diffusion length of the unforced system ($`\gamma (𝐫)\gamma _1=0`$). We find the diffusion length $`\mathrm{}_D=\sqrt{D\tau }3.24`$ by taking the diffusion coefficient $`D`$ to be unity, and the characteristic time $`\tau `$ to be equal to $`2\pi /\beta \mu `$, the period of homogeneous oscillations in the unforced system for the parameters $`\beta =0.6`$, $`\mu =1`$ used in these studies. Parts of the system separated by a length greater than $`\mathrm{}_D`$ evolve independently over time intervals less than $`\tau `$. For large $`R`$ one observes that pacemaker nucleation occurs locally on the boundary of the disk. As the disk perimeter increases with $`R`$ so does the probability of forming a local pacemaker on the disk boundary.

In addition to the target patterns discussed above, one may also observe spiral waves if the initial conditions contain a phase defect. An example of such a spiral is shown in Fig. 11. It was formed in a realization with $`\overline{\gamma }=0.50<\gamma _c`$, hence the medium may be thought of as oscillatory, exhibiting three-fold symmetric relaxational local dynamics, rather than as tristable.

## IV Systems with time-varying spatial disorder

We now consider situations where the spatial distribution of forcing amplitudes varies in time. For this purpose we examine the behavior of the spatially distributed FitzHugh-Nagumo (FHN) system

$`{\displaystyle \frac{u(𝐫,t)}{t}}`$ $`=`$ $`uu^3v+D_u^2u`$ (28)
$`{\displaystyle \frac{v(𝐫,t)}{t}}`$ $`=`$ $`ϵ\left(uav+b(𝐫,t)\right)+D_v^2v,`$ (29)
subject to such time-varying noise distributions. Here $`u(𝐫,t)`$ and $`v(𝐫,t)`$ are the “concentrations”, $`a`$ and $`ϵ`$ are constant parameters, $`D_u`$ and $`D_v`$ are the diffusion coefficients and $`b(𝐫,t)`$ is a forcing function of the form $`\eta (𝐫,t)\mathrm{cos}\omega _\mathrm{f}t`$ with the forcing amplitude $`\eta (𝐫,t)`$ a random variable. We have chosen this form to make explicit the periodic component of the external forcing whose amplitude is a periodically or randomly updated spatial random variable. In the applications described below $`\eta (𝐫,t)`$ is updated on a time scale that is some multiple of the forcing period.

Before considering Eq. (29) we discuss the associated system of ordinary differential equations

$`{\displaystyle \frac{d}{dt}}u(t)`$ $`=`$ $`uu^3v`$ (30)
$`{\displaystyle \frac{d}{dt}}v(t)`$ $`=`$ $`ϵ(uav+b),`$ (31)
in which $`b`$ is a constant. If $`0<a<1`$ the system possesses a single fixed point. If, in addition, $`b=0`$, then the fixed point is unstable and the system exhibits a stable limit cycle whenever $`aϵ<1`$. In this article we consider only systems with $`0<a<1`$ and $`aϵ<1`$. Figure 12 shows the limit cycle for a system with $`a=0.3`$, $`ϵ=0.1`$ and $`b=0`$. Also shown are the cubic $`u`$-nullcline and the linear $`v`$-nullcline for different $`b`$ values. The effect of varying $`b`$ is to shift the $`v`$-nullcline relative to the $`u`$-nullcline. As $`|b|`$ increases from zero the limit cycle contracts in the phase plane, eventually collapsing to a stable fixed point at a Hopf bifurcation that occurs when $`|b|=b_H=\left(1(2a+a^2ϵ)/3\right)\left((1aϵ)/3\right)^{\text{“sevenrm1}/\text{“sevenrm2}}`$. For $`|b|>b_H`$ the system exhibits excitable kinetics. Figure 12 shows $`v`$-nullclines corresponding to $`b=\pm 0.46`$, which lies just inside the excitable regime.

Given this background, we consider the dynamics of Eq. (29) as an example of an oscillatory reaction-diffusion system with periodic forcing. The forcing amplitude field $`\eta (𝐫,t)`$ was similar to the $`\gamma (𝐫)`$ fields used in the quenched disorder studies described earlier in that the system was divided into squares which were randomly assigned one of the two forcing intensities $`\eta _1`$ and $`\eta _2`$. This disordered forcing amplitude field was periodically updated, with the new values of the amplitude in the spatial distribution drawn from the same dichotomous distribution of amplitudes. The updating period was taken to be $`nT_\mathrm{f}`$ time units, which is the period of the corresponding $`n:m`$ entrained ordinary differential equation. (The investigations described here consider only the case where the updating is on resonance, i.e. where the interval between updates is $`nT_\mathrm{f}`$. We shall not consider the case of periodic updating where the updating is off-resonance.)

The $`nT_\mathrm{f}`$-periodic updating of the forcing amplitude field may have a phase offset relative to the $`T_\mathrm{f}`$-periodic forcing $`\mathrm{cos}\omega _\mathrm{f}t`$. We describe this offset with the parameter $`\sigma `$ which ranges from 0 to 1 and specifies the phase offset in units of $`nT_\mathrm{f}`$, i.e. as a fraction of a period. Thus, updates occur at times

$$\tau _k=\left((k1)+\sigma \right)nT_\mathrm{f}\mathrm{f}ork=1,2,3\mathrm{},$$
(32)
in addition to the initial specification of the random forcing field at $`\tau _0=0`$. To formalize the foregoing: the studies were carried out with $`b(𝐫,t)=\eta (𝐫,t)\mathrm{cos}\omega _\mathrm{f}t`$, where

$$\eta (𝐫,t)=\underset{k=0}{\overset{\mathrm{}}{}}\underset{i=1}{\overset{N_W}{}}\underset{j=1}{\overset{N_L}{}}\xi _{ij}^k\theta (t\tau _k)\theta (\tau _{k+1}t)\mathrm{\Theta }_{ij}(𝐫),$$
(33)
where $`\theta `$ is the Heaviside function, $`\mathrm{\Theta }_{ij}(𝐫)`$ is the characteristic function selecting the square with discrete coordinates $`(ij)`$ as in Eq. (11), and

$$\xi _{ij}^k=\{\begin{array}{cc}\hfill \eta _1\text{with probability}& p\hfill \\ \hfill \eta _2\text{with probability}& q=1p.\hfill \end{array}$$
(34)
The spatial average and time average are equal and given by $`\overline{\eta }(t)=p\eta _1+q\eta _2=\eta (𝐫)`$. The space-time autocorrelation function is

$`C(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{\delta \eta (𝐫^{}+𝐫,t^{}+t)\delta \eta (𝐫^{},t^{})}{\delta \eta (𝐫^{},t^{})\delta \eta (𝐫^{},t^{})}}`$ (35)
$`=`$ $`\{\begin{array}{cc}(1\frac{|x|}{s})(1\frac{|y|}{s})(1\frac{|t|}{nT_\mathrm{f}})& \mathrm{i}f|x|s,|y|s\mathrm{a}nd|t|nT_\mathrm{f}\\ 0& \mathrm{o}therwise\end{array}.`$ (38)
The investigations described in this section considered systems at the $`n:m=3:1`$ resonance, using the parameters $`a=0.3`$, $`ϵ=0.1`$, $`\omega _\mathrm{f}/\omega _0=3.05`$, $`D_u=D_v=0.25`$ with the forcing field parameters $`\eta _1=0`$, $`\eta _2=0.92`$, $`p=q=0.50`$ and noise grain size $`s\times s=4\times 4`$. The corresponding mean field system, Eq. (29) with $`b(𝐫,t)=\eta _0\mathrm{cos}\omega _\mathrm{f}t`$, $`\eta _0=\overline{\eta }=0.46`$, lies in the entrained regime and admits three-armed phase locked spiral waves (Fig. 13).

Figure 14 shows an example of a spiral wave in the FHN system with quenched disorder, analogous to that considered in Sec. III B for the CGL equation. For this case the $`\eta (𝐫,t)`$ field may be described within Eqs. (33)-(34) if we take $`\tau _1=+\mathrm{}`$. Substantial front roughening is apparent. This system exhibits a phenomenon not seen in the studies in the inhomogenously forced CGL described in Sec. III. In the uniformly forced FHN, with $`\eta (𝐫,t)\eta _0`$ and other parameters the same, the front velocity passes through zero as $`\eta _0`$ is varied. Variations in the effective local $`\eta `$ values, combined with front curvature effects, result in frequent local pinning of the fronts. The fronts may be depinned through coupling to mobile portions of the front, or by the perturbation provided by a following front as it approaches near the pinned front. Thus, the fronts move with an irregular stop-start motion that is controlled by the pinning and depinning events. It seems unlikely that the resulting fronts could obey KPZ scaling; indeed, inspection of Fig. 14 suggests the front profile is unlikely to even remain a single-valued function of position. Realizations of spirals in a system with smaller size eventually reached a stationary configuration where fronts were pinned everywhere along their lengths.

For the FHN system with time-varying disorder and the aforementioned parameters, three quasi-homogeneous states were observed, similar to the behavior seen in the CGL with quenched disorder. The existence of noise-update events breaks the symmetry between the different entrained states of the system. Similarly, domain walls are now no longer equivalent and may travel at different velocities. Consider the 3:1 forced system and arbitrarily label the phases 1, 2, and 3. In the following discussion, a  front means a domain wall between phases 3 and 1, with phase 3 on the left; its opposite front is . There are three front types: , ,  (and their opposites). The velocities of these fronts were measured as function of $`\sigma `$ (Fig. 15). It suffices to measure the velocities for $`5/6<\sigma <1`$, the values for other $`\sigma `$ follow from the system’s symmetry under $`tt+T_\mathrm{f}`$ and $`(u,v,t)(u,v,t+T_\mathrm{f}/2)`$ . All three fronts move to the left (positive velocity).

Depending on $`\sigma `$, their velocities rank as $`v_{12}>v_{23}>v_{31}`$ or as $`v_{23}>v_{12}>v_{31}`$. We note that for all $`\sigma `$, $`v_{12}>v_{31}`$. Thus, if a system has initial conditions consisting of two plane fronts, \[31 . . . 12\] (where . . . represents a region of phase 1), the  front will move faster than the  front and the distance between the two will decrease. Eventually, the  front will closely approach the  front and a new stable propagating front consisting of a thin layer of phase 1 connecting phases 3 and 2 will result. We term such a front a compound front and denote it . Similarly, for $`\sigma `$ values where $`v_{23}>v_{12}>v_{31}`$, the compound front  exists and can be obtained from the starting configuration \[12 . . . 23\].

As one might expect, compound fronts cannot be made from a slower moving front following a faster moving front. For example,  is not stable; it splits into \[23 . . . 31\] because $`v_{23}>v_{31}`$.

The velocities of these compound fronts were measured as a funtion of $`\sigma `$ and, within our numerical accuracy, were found to lie between the velocities of the two simple fronts from which they were derived (i.e. $`v_{12}>v_{312}>v_{31}`$ and $`v_{23}>v_{123}>v_{12}`$ ) (Fig. 15). Depending on $`\sigma `$ either $`v_{312}>v_{23}`$ or $`v_{23}>v_{312}`$. For $`v_{312}>v_{23}`$ one expects the travelling pulse solution  to be stable, and this is indeed the case, while for $`v_{23}>v_{312}`$ the pulse  is unstable (it splits into  and ) while the pulse  is stable. The pulse velocities were measured, they are essentially the same as the velocity of the faster moving component of the pulse.

The existence of stable pulse solutions joing two domains of the same phase raises the question of whether a “one-armed” spiral whose arms consist of the pulse can exist. Figure 16 shows a stable spiral for $`\sigma =11/12`$, a regime where the velocity ordering is $`v_{12}>v_{312}>v_{23}>v_{31}`$. It is not a “one-armed” spiral with a  pulse front but could be viewed as a two-armed spiral with arms consisting of phases 3 and 2, with fronts of type  and . Since the  front velocity is greater than that of the  front, one expects that as the waves travel outward phase 3 will shrink and phase 2 will grow, and far from the core the waves will become a train of  pulses.

Figure 17 shows a spiral for $`\sigma =1`$. In this regime the velocity ordering is $`v_{23}>v_{123}>v_{3123}>v_{12}>v_{312}>v_{31}`$. The stable pulse is  rather than . Far from the core we expect the waves to become a train of  pulses and in Fig. 17 this can indeed be seen to happen.

The motion of the spiral core was recorded for a realization of the dynamics with $`\sigma =11/12`$. The trajectory of the core, $`𝐫_\mathrm{c}(t)=(x_\mathrm{c}(t),y_\mathrm{c}(t))`$, describes a “noisy flower pattern” (Fig. 18). Both the periodic looping motion and distortions of the simple flower pattern due to the noise are evident. A plot of $`|𝐫_\mathrm{c}|^2`$ vs. $`t`$ shows periodic behavior with period $`\mathrm{17\hspace{0.33em}000}`$, which is also the mean period of rotation of the spiral. In the mean field system with $`\eta (𝐫,t)\overline{\eta }`$ the core is stationary. It is also stationary for the uniformly forced systems with $`\eta (𝐫,t)\eta _1`$ and $`\eta (𝐫,t)\eta _2`$, the two extreme values of $`\eta `$ in the dichotomous noise process. Consequently, the core motion is a result of the time-varying spatial disorder of the forcing amplitude field.

We have also investigated time-varying noise where updates occurred at Poisson-distributed intervals instead of periodically. The Poisson distribution used was $`\mathrm{P}r(t\mathrm{\Delta }t_{\mathrm{}}t+dt)=(1/\overline{t})e^{t/\overline{t}}dt`$ where $`\overline{t}=\mathrm{\Delta }t_{\mathrm{}}=nT_\mathrm{f}`$. With this choice of $`\overline{t}`$ the mean time between updates is the same as in the on-resonance periodic updating case discussed previously and corresponds to one period of the entrained system. Thus, if $`\mathrm{\Delta }t_{\mathrm{}}`$ are chosen from this distribution then updates of the forcing field occur at times

$$\tau _k=\underset{\mathrm{}=1}{\overset{k}{}}\mathrm{\Delta }t_{\mathrm{}}$$
(39)
in addition to the initial specification of the forcing field at $`\tau _0`$. With these definitions of $`\tau _k`$ in place of Eq. (32), $`\eta (𝐫,t)`$ is as given in Eqs. (33)-(34).

We expect that in this system the three phases will be equivalent on average on time scales longer than the average interval between $`\eta (𝐫)`$ field updates. The observed spiral shown in Fig. 19 confirms this equivalence. The three arms seen in the figure are approximately equivalent and, when the animation of the dynamics is viewed, this equivalence is preserved in time.

## V Discussion

We have explored the phenomenology of resonantly forced oscillatory reaction-diffusion systems subject to both quenched and time-varying disorder in the forcing amplitude field. Noteworthy phenomena found when there is quenched disorder are front roughening and spontaneous nucleation of target patterns. Spontaneous nucleation of target patterns arises because diffusion effectively causes averaging of $`\gamma (𝐫)`$ ocally over some length scale; hence, the medium locally behaves like a uniform system with the same $`\overline{\gamma }`$. Alternatively but equivalently, we may describe the dynamics as arising from competition between two regimes selected by the dichotomous forcing values which lie on either side of a bifurcation point. In some spatial regions one regime dominates the dynamics while in other regions the second regime dominates.

For time-varying disorder, the symmetry of the system’s three states is broken, and thus the velocities of different domain walls differ. As a result, travelling front structures other than simple kink-like fronts may exist. These are the compound fronts and pulses. The velocities of the domain walls depend on the phase of the forcing field updating, therefore the updating parameter $`\sigma `$ selects between regimes in which different sets of compound fronts and pulses exist. This asymmetry among states also leads to spiral waves with inequivalent arms.

Meandering core motion with a noisy flower-like trajectory is seen with time-varying disorder. This core motion arises purely from inhomogeneities in $`\eta (𝐫,t)`$ since none of the uniformly forced systems with $`\eta (𝐫,t)\overline{\eta },\eta _1,\eta _2`$ exhibit core meandering.

Some of the effects described herein may arise from the spatial inhomogeneity of the forcing, and not necessarily from its disordered nature. One may also consider regular patterns of inhomogeneous forcing and investigations of this type are in progress. Many of the effects found in this study, such as target pattern nucleation, KPZ front roughening, and front dynamics for time-varying disorder arise from the stochastic nature of the forcing field and would not be found in a system with a regular pattern of forcing.

There are opportunities for further exploration of spatially disordered resonantly forced systems; for example, the effects of noise on various phenomena or bifurcations known in spatially uniform resonantly forced systems have not been studied. Additionally, one may investigate oscillatory systems possessing richer limit cycle dynamics than the relatively uncomplicated examples considered in this paper.

Systems of the types described here could be readily realized in an experimental setup. Standard methodology for investigation of spatial disorder in the light-sensitive excitable BZ reaction involves use of a computer-controlled video projector to project a precisely controllable spatiotemporal pattern of illumination intensity onto the reaction medium. Thus, the only necessary changes are the use of the light-sensitive BZ in the oscillatory regime and reprogramming of the projector to provide a periodic illumination signal incorporating appropriate stochastic spatiotemporal modulation of the light intensity.

## Acknowledgements

This work was supported in part by the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# References

Summation Formulas for the product of the q-Kummer Functions from $`E_q(2)`$

H. Ahmedov1 and I. H. Duru2,1

1.Feza Gürsey Institute, P.O. Box 6, 81220, Çengelköy, Istanbul, Turkey 111E–mail : hagi@gursey.gov.tr and duru@gursey.gov.tr.

2.Trakya University, Mathematics Department, P.O. Box 126, Edirne, Turkey.

Abstract

Using the representation of $`E_q(2)`$ on the non-commutative space $`zz^{}qz^{}z=\sigma `$; $`q<1`$, $`\sigma >0`$ summation formulas for the product of two, three and four q-Kummer functions are derived.

1. Introduction

Properties of manifolds can be investigated by means of their automorphism groups. Non-commutative spaces too are studied similarly. For example the quantum groups $`E_q(2)`$ and $`SU_q(2)`$ are the symmetry groups of the quantum plane and the quantum sphere respectively -. The group representation theory gives the possibility to construct the complete set of orthogonal functions on these deformed spaces. For example the Hahn-Exton q-Bessel and q-Legendre functions appears as the matrix elements of the unitary representations of $`E_q(2)`$ - and $`SU_q(2)`$   which are the complete set of orthogonal functions on the quantum plane and the quantum sphere respectively. Using group theoretical methods the invariant distance and the Green functions have also been written on the quantum sphere  and the quantum plane .

In recent works we have studied the non-commutative space $`[z,z^{}]=\sigma `$ ( i. e. the space generated by the Heisenberg algebra) by means of its automorphism groups $`E(2)`$ and $`SU(1,1)`$ , . The basis in this noncommutative space where irreducible representations of $`E(2)`$ are realized were found to be the Kummer functions which involves the coordinates $`z,z^{}`$ not as their arguments but as indices. That study enable us to obtain generic summation formulas involving Kummer and Bessel functions. For $`SU(1,1)`$ case the basis is given in terms of the hypergeometric functions having the non-commutative coordinates $`z`$ and $`z^{}`$ as the parameters. Again we derived generic summation formulas involving hypergeometric and Jacobi functions. This analysis enable us to construct different complete sets of orthogonal functions on the non-commutative space. Both studies also provide new group theoretical interpretations for the already known relations involving special functions.

Motivated by the outcomes of the above mentioned studies, in the present work we consider the two parametric deformation of the plane which is the $``$-algebra $`P_q^\sigma `$ generated by $`z`$ and $`z^{}`$ with

$$zz^{}qz^{}z=\sigma ,q<1,\sigma >0$$
(1)
which possesses the symmetry of the group $`E_q(2)`$. In $`\sigma 0`$ limit it becomes the usual quantum plane. In $`q1`$ limit it turns to the algebra of functions on the Heisenberg algebra. This study allows us to obtain many identities involving several Hahn-Exton q-Bessel and Moak q-Laguerre functions which are the special forms of the q-Kummer functions. Note that previously some formulas involving q-Laguerre functions were derived by making use of the representation theory of the q-oscillator algebra -. Some relations involving the basic Bessel and Laguerre functions were also considered in .

In Section 2 we realize $`E_q(2)`$ as the automorphism group of the non-commutative space $`P_q^\sigma `$.

In Section 3 we construct the basis in $`P_q^\sigma `$ where the irreducible representations of $`E_q(2)`$ are realized.

Section 4 is devoted to the generic summation formulas for the product of two, three and four q-Kummer functions.

In Section 5 some simple examples are presented.

2. $`E_q(2)`$ as the symmetry group of $`P_q^\sigma `$

Euclidean group transformations of $`z,z^{}`$ plane are given by

$`\delta (z)`$ $`=`$ $`B+Az,`$ (2)
$`\delta (z^{})`$ $`=`$ $`B^{}+A^{}z^{}.`$ (3)
The conditions

$$\delta (z)\delta (z^{})q\delta (z^{})\delta (z)=1$$
(4)
and

$$(\delta (z))^{}=\delta (z^{})$$
(5)
imply the $`E_q(2)`$ relations

$$BB^{}=qB^{}B,AB=qBA,AB^{}=qB^{}A,A^{}=A^1.$$
(6)
Note that $`z,z^{}`$ commute with $`B,B^{}`$ and $`A`$. (Starting with(4) we employ $`\sigma =1`$. When we need to calculate $`\sigma 0`$ limit we replace $`z,z^{}`$ by $`\frac{z}{\sigma },\frac{z^{}}{\sigma }`$). Formulas

$`z|n,j`$ $`=`$ $`\sqrt{(n)_q}|n1,j,`$
$`z|n,j`$ $`=`$ $`\sqrt{(n+1)_q}|n+1,j`$
$`B|n,j`$ $`=`$ $`q^{\frac{j}{2}}|n,j1,`$ (7)
$`B^{}|n,_j`$ $`=`$ $`q^{\frac{j+1}{2}}|n,j+1,`$
$`A|n,j`$ $`=`$ $`|n,j2,`$
where

$$(n)_q=\frac{1q^n}{1q}$$

give the solution of (1) and (6) in some suitable domain of the Hilbert space $`X`$ with the basis $`\{|n,j\}`$, $`n=0,1,2,\mathrm{}`$ and0, $`j𝐙`$. Let us define in $`X`$ a new basis such that

$`\delta (z)|n,j^{}`$ $`=`$ $`\sqrt{(n)_q}|n1,j^{}`$ (8)
$`\delta (z)|n,j^{}`$ $`=`$ $`\sqrt{(n)_q}|n+1,j^{}.`$ (9)
Due to

$$ze_q^{xz^{}}=xe_q^{xz^{}}+e_q^{qxz^{}}z,$$
(10)
with

$$e_q^x=\underset{k=0}{\overset{\mathrm{}}{}}\frac{x^k}{(k)_q!}$$
(11)
being the q-deformed exponential function, we conclude that

$$|0,j^{}=e_q^{A^{}Ba^{}}\sqrt{e_q^{B^{}B}}|0,j$$
(12)
is the ground state of the new basis:

$$\delta (z)|0,j^{}=0.$$
(13)
Applying the creation operator $`(\delta (z))^{}`$ on this state we can generate the desired basis in $`X`$:

$$|n,j^{}=\frac{(\delta (z^{}))^n}{\sqrt{(n)_q!}}|0,j^{}.$$
(14)
We also have

$$\delta (z)=UzU^{},$$
(15)
where $`U`$ is the unitary operator in

$$|n,j^{}=U|n,j.$$
(16)
Thus $`\delta `$ defines the homomorphic map of $`P_q^\sigma `$ into the $``$-algebra generated by $`A,B,z`$ and their adjoints. Applying twice this map we get

$$\delta ^{}(\delta (z))=B^{}+A^{}\delta (z)=B^{\prime \prime }+A^{\prime \prime }z,$$
(17)
where

$$B^{\prime \prime }=B^{}+A^{}B,A^{\prime \prime }=A^{}A,$$
(18)
is the group multiplication in $`E_q(2)`$ ( the operator $`B^{\prime \prime }`$ has the same analytic properties as $`B`$ ). $`\delta `$ defines the representation of $`E_q(2)`$ in $`P_q^\sigma `$.

Before closing this section we give the explicit formula for the matrix representation of $`U`$:

$$U_{(mi)(nj)}=mi|n,j^{}.$$
(19)
For $`|n,j=|n|j`$ we first define

$$U_{mn}=n|\delta (z^{}))^ne_q^{A^{}Bz^{}}|0\sqrt{\frac{e_q^{B^{}B}}{(n)_q!}}$$
(20)
which is the function of $`B`$, $`B^{}`$, $`A`$ and $`A^{}`$. Then

$$U_{(mi)(nj)}=i|U_{mn}|j.$$
(21)
After some algebra we get

$$U_{mn}=A^mB^{nm}\mathrm{\Phi }_{mn}(\eta )\mathrm{for}nm$$
(22)
and

$$U_{mn}=q^{\frac{(mn)(mn1)}{2}}A^m(B)^{mn}\mathrm{\Phi }_{nm}(\eta )\mathrm{for}mn,$$
(23)
where $`\eta ^2B^{}B`$ and

$$\mathrm{\Phi }_{mn}(\eta )=\sqrt{\frac{(n)_q!}{(m)_q!}}\frac{\sqrt{e_q^{\eta ^2}}}{(nm)_q!}\mathrm{\Phi }^q(q^m,q^{1+nm};q^{n+1}\eta ^2).$$
(24)
Here

$$\mathrm{\Phi }^q(a,b;x)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{q^{\frac{k(k1)}{2}}(a;q)_k}{(q;q)_k(b;q)_k}((1q)x)^k$$
(25)
which in $`q1`$ limit reduces to the Kummer function:

$$\underset{q1}{lim}\mathrm{\Phi }^q(q^c,q^d;x)=\mathrm{\Phi }(c,d;x).$$
(26)
We call it the q-Kummer function. The functions $`\mathrm{\Phi }_{nm}`$ can also be expressed in terms of Moak’s q-Laguerre polynomials

$$L_n^{q(\alpha )}(x)=\frac{(q^{1+\alpha };q)_n}{(q;q)_n}\mathrm{\Phi }^q(q^n,q^{1+\alpha };q^{1+\alpha +n}x)$$
(27)
as

$$\mathrm{\Phi }_{mn}(\eta )=\sqrt{e_q^{\eta ^2}\frac{(m)_q!}{(n)_q!}}L_m^{q(nm)}(\eta ^2),\mathrm{for}nm.$$
(28)
3. Irreducible representations of $`E_q(2)`$ in $`P_q^\sigma `$

The irreducible representation of the deformed enveloping algebra $`U_q(e(2))`$

$$P^{}P=qPP^{},KP=qPK,P^{}K=qKP^{}$$
(29)
defined by the weight $`\lambda 𝐑`$ in the non-commutative space $`P_q^\sigma `$ is given by

$$R(K)D_j^\lambda (z,z^{})=q^jD_j^\lambda (z,z^{})$$
(30)
$$R(P)D_j^\lambda (z,z^{})=\lambda q^{\frac{j}{2}}D_{j1}^\lambda (z,z^{}),$$
(31)
$$R(P^{})D_j^\lambda (z,z^{})=\lambda q^{\frac{j+1}{2}}D_{j+1}^\lambda (z,z^{}),$$
(32)
where

$`R(K)z^nz^m`$ $`=`$ $`q^{mn}z^nz^m`$
$`R(P)z^nz^m`$ $`=`$ $`iq^n(m)_qz^nz^{m1}`$ (33)
$`R(P^{})z^nz^m`$ $`=`$ $`iq^{n+1}(n)_qz^{n1}z^m`$
defines the right realization of $`U_q(e(2))`$. (30) implies

$$D_j^\lambda (z,z^{})=\{\begin{array}{c}f_j^\lambda (\zeta )z^j\mathrm{for}j0\\ z^jf_j^\lambda (\zeta ),\mathrm{for}j0\end{array}$$
(34)
where

$$\zeta 1(1q)z^{}z.$$
(35)
Inserting the ansatz (34) in (31) and (32) we get

$$f_j^\lambda (\zeta )=\frac{q^{\frac{j^2}{4}}(i\lambda )^j}{(j)_q!}\mathrm{\Phi }^q(\zeta ^1,q^{j+1};q^{j+1}\lambda ^2\zeta ).$$
(36)
In the derivation of (36) we used

$$z^nz^n=()^nq^{\frac{n(1n)}{2}}(\zeta ^1;q)_n\zeta ^n.$$
(37)
By means of the universal T-matrix we can exponentiate (30) and get

$$\delta (D_j^\lambda )=\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}t_{ji}^\lambda D_i^\lambda ,$$
(38)
where

$`t_{ij}^\lambda `$ $`=`$ $`q^{\frac{i^2j^2}{4}}{\displaystyle \frac{(i\lambda B)^{ij}A^j}{(ij)_q!}}\mathrm{\Phi }^q(0,q^{1+ij};(q1)q^{1j}(\lambda \eta )^2)\mathrm{for}ij`$ (39)
$`t_{ij}^\lambda `$ $`=`$ $`q^{\frac{i^2j^2}{4}}{\displaystyle \frac{A^j(i\lambda B^{})^{ji}}{(ji)_q!}}\mathrm{\Phi }^q(0,q^{1+ji};(q1)q^{1i}(\lambda \eta )^2)\mathrm{for}ji`$ (40)
are the matrix elements of the irreducible representations of $`E_q(2)`$ , . We can express them in terms of the Hahn-Exton q-Bessel functions

$$J_k^q(x)=\frac{x^k}{(k)_q!}\mathrm{\Phi }^q(0;q^{1+k}|q;(q1)qx^2)$$
(41)
as

$`t_{ij}^\lambda `$ $`=`$ $`(\sqrt{1}q^{\frac{1}{4}})^{ij}V^{ij}A^jJ_{ij}^q(q^{\frac{j}{2}}\lambda \eta )\mathrm{for}ij`$ (42)
$`t_{ij}^\lambda `$ $`=`$ $`(\sqrt{1}q^{\frac{1}{4}})^{ji}V^{ij}A^jJ_{ji}^q(q^{\frac{i}{2}}\lambda \eta )\mathrm{for}ji,`$ (43)
where $`V`$ is the unitary operator defined by $`B=V\eta `$.

In $`\sigma 0`$ limit the non-commutative space $`P_q^\sigma `$ becomes the quantum plane $`E_q(2)/U(1)`$ generated by $`B`$, $`B^{}`$:

$$\underset{\sigma 0}{lim}D_j^{\sqrt{\sigma }\lambda }(\frac{B}{\sqrt{\sigma }},\frac{B^{}}{\sqrt{\sigma }})=t_{j0}^\lambda .$$
(44)
In $`q1`$ limit $`P_q^\sigma `$ becomes the noncommutative space generated by the Heisenberg algebra :

$$\underset{q1}{lim}D_j^\lambda (z,z^{})=\mathrm{\Phi }(zz^{};1+j;\lambda ^2)\frac{(i\lambda z)^j}{j!}.$$
(45)
In $`\sigma 0`$ and $`q1`$ limit we arrive at the complex plane $`E(2)/U(1)`$:

$$\underset{\sigma 0}{lim}\underset{q1}{lim}D_j^{\sqrt{\sigma }\lambda }(\frac{re^{i\psi }}{\sqrt{\sigma }},\frac{re^{i\psi }}{\sqrt{\sigma }})=i^je^{ij\psi }J_j(\lambda r).$$
(46)
4. Summation formulas for the q-Kummer functions

(15) and (38) imply

$$UD_j^\lambda U^{}=\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}t_{ji}^\lambda D_i^\lambda $$
(47)
or

$$UD_j^\lambda =\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}t_{ji}^\lambda D_i^\lambda U,$$
(48)
$$D_j^\lambda =\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}t_{ji}^\lambda U^{}D_i^\lambda U.$$
(49)
The above formulas define the summation of products of two, three and four q-Kummer functions. Sandwiching (47), (48) and (49) between the states $`m|`$ and $`|n`$ and using

$$(D_j^\lambda )_{mn}=\sqrt{\frac{(n)_q!}{(m)_q!}}f_j^\lambda (q^m)\delta _{j,nm},\mathrm{for}j0,$$
(50)
$$(D_j^\lambda )_{mn}=(D_j^\lambda )_{nm}$$
(51)
we get

$`{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}(D_j^\lambda )_{ss+j}U_{ms}U_{s+jn}^{}`$ $`=`$ $`(D_{nm}^\lambda )_{mn}t_{jnm}^\lambda \mathrm{for}j0`$ (52)
$`{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}(D_j^\lambda )_{sjs}U_{msj}U_{sn}^{}`$ $`=`$ $`(D_{nm}^\lambda )_{mn}t_{jnm}^\lambda \mathrm{for}j<0,`$ (53)
$`{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}t_{jsm}^\lambda (D_{sm}^\lambda )_{ms}U_{sn}`$ $`=`$ $`(D_j^\lambda )_{njn}U_{mnj}\mathrm{for}nj`$ (54)
$`{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}t_{jsm}^\lambda (D_{sm}^\lambda )_{ms}U_{sn}`$ $`=`$ $`0\mathrm{for}n<j`$ (55)
and

$$\underset{s,l=0}{\overset{\mathrm{}}{}}t_{jls}^\lambda (D_{ls}^\lambda )_{sl}U_{ms}^{}U_{ln}=(D_j^\lambda )_{mn}\delta _{jnm}.$$
(56)
In the above formulas $`U^{}`$s and $`t^{}`$s are given in terms of the q-Kummer functions of the operator $`\eta =B^{}B`$ (see (22), (23) and (39).

In the coming section we give some simple examples.

5. Examples

A. For $`n=m=0,j0`$ (52) implies

$$\underset{s=0}{\overset{\mathrm{}}{}}\frac{q^{\frac{s(1s)}{2}}\eta ^{2s}}{(s)_q!}\mathrm{\Phi }^q(q^s,q^{1+j};q^{1+j+s}\lambda ^2)=\frac{(i\lambda \eta )^j(j)_q!}{\sqrt{e_q^{\eta ^2}e_q^{q^j\eta ^2}}}J_j^q(\lambda \eta )$$
(57)
which in $`q1`$ limit gives  (page 1038, Eq. 3 of 8.975)

$$\underset{s=0}{\overset{\mathrm{}}{}}\frac{\eta ^{2s}}{s!}\mathrm{\Phi }(s,1+j;\lambda ^2)=j!(\eta \lambda )^je^{\eta ^2}J_j(\lambda \eta )$$
(58)
B. For $`n=0`$ and $`kjm`$ (54) implies

$$\underset{s=0}{\overset{\mathrm{}}{}}q^{\frac{s(s+m2k)}{2}}C_{sm}\eta ^sJ_{s+km}^q(q^{\frac{s+k}{2}}\lambda \eta )=\frac{q^{\frac{k(km)}{2}}\lambda ^k\eta ^{km}}{(m)_q!(km)_q!}\mathrm{\Phi }^q(q^m,q^{1+km};q^{1+k}\eta ^2),$$
(59)
where

$$C_{sm}=\frac{()^s\lambda ^{ms}}{(s)_q!(ms)_q!}\mathrm{\Phi }_q(q^s,q^{1+ms};q^{1+m}\lambda ^2),\mathrm{for}ms.$$
(60)
For $`sm`$ one has to replace $`m,s`$ with $`s,m`$ in the right hand side of the above expression. When $`m=0`$ we have

$$\underset{s=0}{\overset{\mathrm{}}{}}q^{\frac{s^2}{2}+sk}\frac{(\lambda \eta )^s}{(s)_q!}J_{s+k}^q(q^{\frac{s+k}{2}}\lambda \eta )=q^{\frac{k^2}{2}}\frac{(\lambda \eta )^k}{(k)_q!}$$
(61)
which is the quantum analogue of a known formula  (page 974, Eq 1 of 8.515).

C. For $`j=\lambda =0`$ (56) implies the unitarity condition for the operator $`U`$:

$$\underset{s=0}{\overset{\mathrm{}}{}}(U_{sm})^{}U_{sn}=\delta _{nm},$$
(62)
where we used

$$U_{ms}^{}=(U_{sm})^{}.$$
(63)
For $`n=m`$ with $`x=\eta ^2`$ we have

$$\underset{s=0}{\overset{n1}{}}\frac{q^{\frac{(ns)(ns+1)}{2}}(n)_q!x^{ns}}{(s)_q!}(L_s^{q(sn)}(x))^2+\underset{s=n}{\overset{\mathrm{}}{}}\frac{q^{\frac{(ns)(ns+1)}{2}}(s)_q!x^{sn}}{(n)_q!}(L_n^{q(ns)}(x))^2=e_{q^1}^x$$
(64)
In deriving the above examples one frequently uses the identities

$$B^kB^k=q^{\frac{k(k+1)}{2}}\eta ^{2k},B^kB^k=q^{\frac{k(1k)}{2}}\eta ^{2k}.$$
(65)


text:
# Theory of nuclear spin conversion in ethylene

## I Introduction

Nuclear spin isomers of molecules were discovered nearly 70 years ago when ortho and para hydrogen isomers were separated for the first time. Although it was realized already at the time of this discovery that many other symmetrical molecules have nuclear spin isomers too, until recently almost nothing was known about isomers of molecules heavier than hydrogen. This happened because of the lack of practical separation methods. The separation method based on deep cooling was applicable to hydrogen (and deuterium) but failed in the case of heavier molecules. Presently a few new methods for spin isomer separation have been proposed and successfully tested which has advanced this field significantly. Yet we are at the very early stage of this research. The list of molecules for which the nuclear spin isomers have been separated is rather short: H2 , CH3F , H2O , CH$`{}_{2}{}^{}0`$ , Li2 , H$`{}_{3}{}^{}{}_{}{}^{+}`$ , and C2H4 .

The CH3F molecules occupies a special place in this list because it is the only polyatomic molecule for which the isomer conversion mechanism has been identified. The CH3F spin isomers appeared to be extremely stable surviving $`10^9`$ gas collisions. Nevertheless, the isomer conversion was found to be based on ortho-para state mixing induced by tiny intramolecular hyperfine interactions and interruption of this mixing by collisions. This specific type of relaxation was proposed to refer to as quantum relaxation. It is essential for the physics of molecular spin isomers to understand which processes are responsible for the spin conversion in other molecules.

Recently first separation of nuclear spin isomers of ethylene molecules (13CCH4) has been performed . The conversion rate has been determined and it was shown that the rate increases proportional to the gas pressure, $`\gamma /P510^4`$ s-1/Torr. An understanding of mechanism of the conversion in ethylene needs theoretical investigation of the process.

In the present paper we perform first theoretical analysis of the spin conversion in ethylene molecules. The purpose is to verify the consistency of the experimental data on isomer conversion in ethylene with the conversion by quantum relaxation. We consider the same isotope species, 13CCH4, for which the experiment was done . An essential difference of the present study from the well-understood case of CH3F molecules is that now we have to work with an asymmetric-top molecule. This appeared to be much more complicated than the case of symmetric-top molecules. We would like to point out that spin conversion in asymmetric tops were considered theoretically previously for the CH2O and H20 molecules .

## II Isomer conversion by quantum relaxation

Although the spin conversion in molecules by quantum relaxation is explained in other publications we will describe briefly its essence here for the convenience of the reader. At room temperature almost all molecules are situated in their ground electronic and vibrational states. Suppose that these states are separated into two subspaces which are the nuclear spin ortho and para states, as is shown in Fig. 1 for the particular case of the 13CCH4 molecules. Note that some molecules may have more than two nuclear spin isomer forms, e.g., methane, or normal ethylene. On the other hand, number of molecules, e.g., H2, CH3F, and 13CCH4 too, have just two type of spin species.

The relaxation process which we are going to consider has two main ingredients. First, the ortho and para quantum states of the test molecules are not completely independent. There is small intramolecular perturbation, $`\widehat{V}`$, which is able to mix the ortho and para states. This perturbation can mix, in general, many ortho-para level pairs. But in Fig. 1 the mixing is present just for one pair of states as will be the case for 13CCH4, see below. Second, the test molecules are embedded into an environment which is able to induce fast relaxation inside the ortho, or para subspace, but is not able to produce direct transitions between these subspaces. This implies that the relevant cross-section, $`\sigma (ortho|para)=0`$. The relaxation by direct transitions is the process opposite to the quantum relaxation.

The isomer conversion by quantum relaxation consists in the following. Suppose that at the instant $`t=0`$ a test molecule is placed into the ortho subspace. Due to collisions with surrounding particles, the test molecule starts to perform fast migration along rotational states inside the ortho subspace. This is the familiar rotational relaxation. This running up and down along the ladder of the ortho states continues until the molecule jumps in the state $`\alpha `$ which is mixed by an intramolecular perturbation with the energetically close para state $`\alpha ^{}`$. During the free flight after this collision, para state $`\alpha ^{}`$ will be admixed to the ortho state $`\alpha `$. Consequently, the next collision can transfer the molecule in other para states and thus localizes it inside the para subspace. Such a mechanism of spin isomers conversion was proposed in the theoretical work by Curl et al  (see also ).

The quantum relaxation of spin isomers can be quantitatively described in the framework of the kinetic equation for density matrix . One needs to split the molecular Hamiltonian into two parts

$$\widehat{H}=\widehat{H}_0+\widehat{V},$$
(1)
where the main part of the Hamiltonian, $`\widehat{H}_0`$, has pure ortho and para states as the eigenstates; the perturbation $`\widehat{V}`$ mixes the ortho and para states. If at initial instant the nonequilibrium concentration of, say, ortho molecules $`\delta \rho _o(t=0)`$ was created, the system will relax then exponentially, $`\delta \rho _o(t)=\delta \rho _o(0)e^{\gamma t}`$, with the rate

$$\gamma =\underset{a^{}p,ao}{}\frac{2\mathrm{\Gamma }F(a^{}|a)}{\mathrm{\Gamma }^2+\omega _{a^{}a}^2}\left(W_p(\alpha ^{})+W_o(\alpha )\right);F(a^{}|a)\underset{\nu ^{}p,\nu o}{}|V_{\alpha ^{}\alpha }|^2.$$
(2)
The sets of quantum numbers $`\alpha ^{}\{a^{},\nu ^{}\}`$ and $`\alpha \{a,\nu \}`$ consist of the degenerate quantum numbers $`\nu ^{}`$, $`\nu `$ and the quantum numbers $`a^{}`$, $`a`$ which determine the energy of the states. In (2) $`\mathrm{\Gamma }`$ is the decay rate of the off-diagonal density matrix element $`\rho _{\alpha ^{}\alpha }(\alpha ^{}para;\alpha ortho)`$ assumed here to be equal for all ortho-para level pairs; $`\omega _{a^{}a}`$ is the gap between the states $`a^{}`$ and $`a`$; $`W_p(\alpha ^{})`$ and $`W_o(\alpha )`$ are the Boltzmann factors of the corresponding states. For the following it is convenient to introduce the strength of mixing, $`F(a^{}|a)`$, which sums the intramolecular couplings over all degenerate states.

## III Rotational states of ethylene

Ethylene is a popular object in high resolution infrared spectroscopy. This fortunate circumstance made available rather accurate data on molecular parameters and position of molecular rotational levels. The molecular structure, numbering of hydrogen atoms and orientation of molecular system of coordinates are presented in Fig. 2. 13CCH4 is a plane, prolate, nearly symmetric top having the symmetry point group C2v. The characters of the group operations and the irreducible representations are given in the Table 1. We give in the Table 1 also two isomorphic groups, molecular symmetry group C$`{}_{2v}{}^{}(M)\text{[10]}`$ and the point group D2. The bond lengths and angles which were used in our calculation are as follows: $`r_{CH}=1.087\AA `$, $`r_{CC}=1.339\AA `$, the angle $`\alpha _{HCH}`$=117.40 .

The fact that the ethylene molecule is an asymmetric top complicates in two aspect the theoretical analysis of the isomer conversion. First, the energy levels and wave functions for asymmetric tops can be found only numerically but not analytically as is possible for symmetric-top molecules. Second, the rotational quantum number $`k`$ (projection of molecular angular momentum on the molecular symmetry axis) becomes an approximate quantum number in asymmetric tops. Consequently, one should calculate numerically more ortho-para matrix elements than in the conversion of symmetric-top molecules, see, e.g., .

Rotational states of 13CCH4 can be determined using the Hamiltonian of Watson  and the set of molecular parameters from Refs. . The rotational Hamiltonian up to sextic order terms has the form

$`\widehat{H}_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}(B+C)𝐉^2+(A{\displaystyle \frac{1}{2}}(B+C))J_z^2\mathrm{\Delta }_J(𝐉^2)^2\mathrm{\Delta }_{JK}𝐉^2J_z^2\mathrm{\Delta }_KJ_z^4`$ (5)
$`+H_J(𝐉^2)^3+H_{JK}(𝐉^2)^2J_z^2+H_{KJ}𝐉^2J_z^4+H_KJ_z^6`$
$`+{\displaystyle \frac{1}{4}}(BC)F_0\delta _J𝐉^2F_0\delta _KF_2+h_J(𝐉^2)^2F_0+h_{JK}𝐉^2F_2+h_KF_4`$
where $`𝐉`$, $`J_x`$, $`J_y`$, and $`J_z`$ are the molecular angular momentum operator and its projections on the molecular axes. The $`B`$, $`C`$, and $`A`$ are the parameters of a rigid top which characterize the rotation around $`x`$, $`y`$, and $`z`$ molecular axes, respectively (see, Fig. 2). The rest of parameters account for the centrifugal distortion effects . In Eq. (5) the notation was used

$$F_nJ_z^n(J_x^2J_y^2)+(J_x^2J_y^2)J_z^n.$$
(6)
The Hamiltonian (5) can be diagonalized for each value of $`J`$ and $`M`$ in the basis of the symmetric-top quantum states $`|J,k,M>`$ ($`JkJ`$), where $`J`$, $`k`$ and $`M`$ are the quantum numbers of the angular momentum, its projection on the molecular symmetry axis and on the laboratory quantization axis, respectively. The Hamiltonian (5) has diagonal in $`k`$ matrix elements due to the first two lines and matrix elements having $`|\mathrm{\Delta }k|=2`$. The calculations of the rotational eigen states can be simplified if one uses the Wang basis  with $`K=|k|`$:

$`|\alpha ,p>`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|\alpha >+(1)^{J+K+p}|\overline{\alpha }>\right];0<KJ,`$ (7)
$`|\alpha _0,p>`$ $`=`$ $`{\displaystyle \frac{1+(1)^{J+p}}{2}}|\alpha _0>;K=0.`$ (8)
Here $`p=0,1`$ and the sets of quantum numbers are $`\alpha \{J,K,M\}`$; $`|\overline{\alpha }>\{J,K,M\}`$; $`\alpha _0\{J,K=0,M\}`$. Depending on the parity of $`J`$, $`K`$ and $`p`$, the states (8) generate 4 different irreducible representations of the molecular symmetry group D2 (Table 1). The molecular Hamiltonian is full symmetric (symmetry A1). Consequently, the matrix elements between the states of different symmetry disappear. Thus diagonalization of the total Hamiltonian in the basis of (8) is reduced to the diagonalization of four independent submatrices, each for the states of particular symmetry. The rotational states of asymmetric top can be expanded over the basis states (8)

$$|\beta ,p>=\underset{K}{}A_K|\alpha ,p>,$$
(9)
where $`A_K`$ stands for the expansion coefficients. The summation index, $`K`$, is given explicitly in (9), although $`A_K`$ depends on the other quantum numbers as well. All coefficients in the expansion (9) are real numbers because the Hamiltonian (5) is symmetric in the basis $`|\alpha ,p>`$.

There are a few schemes for practical classification of the rotational states of asymmetric tops . In this paper we will use the scheme which is somewhat better adapted to the description of nuclear spin isomers. We will designate the rotational states of asymmetric top by indicating $`p`$, $`J`$ and prescribing the allowed $`K`$ values to the eigenstates keeping both in ascending order. For example, the eigenstate having $`p=0,J=20`$, being the third in ascending order will be designated by $`𝒦=4`$, because the allowed $`K`$ in the expansion (9) are $`K=0,2,4\mathrm{}`$20. (Note the difference between the two characters $`K`$ and $`𝒦`$). It gives unambiguous notation of rotational states for each of the four species A1, A2 ($`𝒦`$-even) and B1, B2 ($`𝒦`$-odd). This classification establishes the connection with the prolate symmetric top for which $`𝒦=K`$. One should remember that physical meaning of the quantum number $`𝒦`$ is clearly limited. To illustrate this we consider as an example the rotational states of A1 and A2 symmetry ($`𝒦`$-even) for $`J=20`$. The upper panel in Fig. 3 shows that the energy of the rotational states is not determined solely by $`𝒦`$ as it would be for a rigid symmetric top. The graph shows the difference in energy between rotational states A1 and A2 having the same $`𝒦`$-number but different quantum numbers $`p`$. The states A1, $`𝒦=0`$ and A2, $`𝒦=0`$ are omitted from the graph because the latter state does not exist. As it is seen from the data the splitting is significant for low values of $`𝒦`$ and rapidly disappears as $`𝒦`$ increases. Note, that the A1 and A2 states for a rigid symmetric-top molecule would be degenerate for all $`𝒦`$.

The low panel in Fig. 3 illustrates the property of the eigen state expansions over the basis states (8). In this panel the squared magnitude of the two $`A_K`$ coefficients in each eigenstate of the symmetry A1 ($`J=20`$) is given. The first one (”$`𝒦`$-term” in the Fig. 3) is $`A_K^2`$ having $`K=𝒦`$. These coefficients appeared to be the biggest coefficients in each expansion. This adds some physical insight to the proposed classification scheme. The ”second term” is the second biggest coefficient in the expansion. Again, at low $`𝒦`$ values there are more than one significant terms in the expansion (9). As $`𝒦`$ grows, the contribution from just one term becomes predominant. Note that the same graph for a symmetric-top molecule would show only one significant term, $`A_{K=𝒦}`$, in the expansions of all eigenstates.

## IV Ortho and para isomers of ethylene

The total molecular wave function is a product of spatial wave function and spin wave function. The nuclear spin states in 13CCH4 are of A1 and B1 species, having the statistical weights 20 and 12, respectively. The ortho spin states (symmetry A1) have the two hydrogen pairs (H1-H2 and H3-H4) either both in triplet state, or both in singlet state. The para species (symmetry B1) have one pair of protons in singlet state but the other pair in triplet state. Each symmetry operation of the C2v group interchange in 13CCH4 even number of protons. Consequently, the total wave function must be unchanged under all operations, thus belonging to the representation A1. In order to have the total wave function being of species A1 one has to have the spatial wave function being of symmetry A1 and B1 as the spin wave functions are. This implies that the rotational states A1 and B1 should be only positive (even in parity) but the rotational states A2 and B2 should be only negative (odd in parity).

Summarizing this discussion we can write the total states in the 13CCH4 molecules. The ortho states can be presented as

$$|\mu >=|\beta ,p>|i_{12},i_{34},i_C>;i_{12},i_{34}=t,s;𝒦\text{even}.$$
(10)
Here $`i_{12}`$, $`i_{34}`$, and $`i_C`$ designate the spin states of the two pairs of protons and the nucleus 13C (Fig. 2). For ortho molecules the spin states of the proton pairs should be either both triplet, $`|t,t>`$, or both singlet states, $`|s,s>`$.

The para states can be presented as

$$|\mu ^{}>=|\beta ^{},p^{}>|i_{12}^{},i_{34}^{},i_C^{}>;𝒦^{}\text{odd}.$$
(11)
For para molecules one pair of proton is in singlet, but the other in triplet state, thus the proton spin states are $`|t^{},s^{}>`$, or $`|s^{},t^{}>`$. In (10), (11) and further we use unprimed parameters for ortho species and primed parameters for para species; $`p=0,1`$ indicates positive, or negative sign of the state, respectively. Eqs. (10), (11) imply that the ortho-para state mixing in ethylene needs coupling of states having different parity of $`𝒦`$.

The relative position of ortho and para states is important for the calculation of the conversion rate (see Eq. (2)). In the case of a symmetric-top molecule the ortho-para energy gaps can be expressed as a polynomial of $`J`$ and $`K`$. This gives “accidental” (isolated) resonances between the ortho and para states in symmetric top. In the case of an asymmetric top one has a phenomenon which can be called collapse of ortho and para states. It appears as a progressive decrease of the ortho-para energy gaps between the states of identical rotational momenta, $`J`$, as $`J`$ increases. This is illustrated in Fig. 4, which shows the gaps between the 13CCH4 ortho states, $`𝒦=0`$, and the para states, $`𝒦^{}=1`$, both having positive sign ($`p=0`$) and the same $`J`$. If such sequence of ortho-para level pairs would exist for symmetric tops it would dominate the conversion. In the case of asymmetric tops the situation is different. We will see below that this sequence of closed ortho-para level pairs do not produce a significant contribution to the isomer conversion in 13CCH4.

The Boltzmann factors $`W_o(\alpha )`$ and $`W_p(\alpha ^{})`$ in Eq. (2) determines the relative population of rotational states in the ortho and para families,

$$\rho _\alpha =\rho _oW_o(\alpha );\rho _\alpha ^{}=\rho _pW_p(\alpha ^{}),$$
(12)
where $`\rho _o`$ and $`\rho _p`$ are the total densities of ortho and para molecules, respectively. The Eqs. (12) imply the equilibrium distributions inside the ortho and para subspaces. This is fulfilled with high accuracy even if the ratio $`\rho _o/\rho _p`$ is out of equilibrium because rotational relaxation is on many orders of magnitude faster than the ortho-para conversion. The partition functions for ortho and para molecules at room temperature (T=295 K) are found to be equal to

$$Z_{ortho}=2.6610^4;Z_{para}=1.6010^4.$$
(13)
In the calculation of these partition functions we took into account the degeneracy over $`M`$, nuclear spins, including also spin of nucleus $`{}_{}{}^{13}C`$, parity of states, as well as the restrictions imposed by quantum statistics.

## V Ortho-para state mixing

In the present paper we will consider the ortho-para conversion induced by the two hyperfine interactions, viz., magnetic dipole-dipole interaction between the molecular nuclei (spin-spin interaction, $`\widehat{V}_{SS}`$) and the nuclear spin-rotation interaction, $`\widehat{V}_{SR}`$. Thus the total intramolecular perturbation able to mix the ortho and para states is

$$\widehat{V}=\widehat{V}_{SS}+\widehat{V}_{SR}.$$
(14)
All matrix elements of $`\widehat{V}`$ are diagonal in parity $`p`$.

### A Nuclear spin-spin coupling

The ortho-para conversion in molecules induced by nuclear spin-spin coupling were investigated in  (other references see in the review ). The spin-spin Hamiltonian for the two magnetic dipoles $`\mu _1`$ and $`\mu _2`$ separated by the distance $`r`$ has the form

$`\widehat{V}_{12}`$ $`=`$ $`P_{12}𝐓^{(12)}{}_{}{}^{}\widehat{𝐈}_{}^{(1)}\widehat{𝐈}^{(2)};`$ (15)
$`T_{ij}^{(12)}`$ $`=`$ $`\delta _{ij}3n_in_j;P_{12}=\mu _1\mu _2/r^3I^{(1)}I^{(2)}h,`$ (16)
where $`\widehat{𝐈}^{(1)}`$ and $`\widehat{𝐈}^{(2)}`$ are the spin operators of the particles 1 and 2, respectively; n is the unit vector directed along r; $`i`$ and $`j`$ are the Cartesian indices. The second rank tensor T(12) in the Eq. (16) represents a spatial part of the spin-spin interaction. The second rank tensor $`\widehat{𝐈}^{(1)}\widehat{𝐈}^{(2)}`$ acts on spin variables.

The total spin-spin interaction in 13CCH4 is composed from the interactions between all pairs of molecular nuclei. One can show that the spin-spin interactions between the protons 1-3, 2-4, 1-2, and 3-4 have the spatial part which can mix the quantum states only if they have $`K`$ numbers of the same parity. Consequently, these terms do not contribute to the ortho-para state mixing in ethylene. The complete nuclear spin-spin interaction in 13CCH4 able to mix the ortho and para states reads

$$V_{SS}=V_{SS}^{(14)}+V_{SS}^{(23)}+V_{SS}^{(C1)}+V_{SS}^{(C2)}+V_{SS}^{(C3)}+V_{SS}^{(C4)}.$$
(17)
Here the upper indices refer to the hydrogen nuclei in the molecule 1…4, and to the nucleus 13C, which has spin equal 1/2 (see Fig. 2). Calculation of the spin-spin matrix elements can be simplified by a few observations which we discuss in more detail for the perturbations $`V_{SS}^{(14)}`$ and $`V_{SS}^{(23)}`$. First, one can proof by applying proper symmetry operation that:

$$<\mu ^{}|V_{SS}^{(14)}|\mu >=<\mu ^{}|V_{SS}^{(23)}|\mu >.$$
(18)
Next, the two, out of four, matrix elements between the ortho and para spin states vanish,

$$<t^{},s^{}|V_{SS}^{(14)}|s,s>=<s^{},t^{}|V_{SS}^{(14)}|s,s>=0,$$
(19)
because in these two cases one has the matrix elements of a vector (spin operator) between the states both having zero spin. On the other hand, the remaining two matrix elements are equal to each other,

$$<t^{},s^{}|V_{SS}^{(14)}|t,t>=<s^{},t^{}|V_{SS}^{(14)}|t,t>.$$
(20)
To summarize, one can conclude that the variety of matrix elements for the operators $`V_{SS}^{(14)}`$ and $`V_{SS}^{(23)}`$ is reduced to one matrix element, e.g., $`<\beta ^{},p,t^{},s^{}|V_{SS}^{(14)}|\beta ,p,t,t>`$. Obviously, the carbon spin state, which is omitted in this expression for simplicity, should be unchanged in this matrix element.

Further, we write this matrix element using an expansion over symmetric-top states (9)

$$<\beta ^{},p,t^{},s^{}|V_{SS}^{(14)}|\beta ,p,t,t>=<t^{},s^{}|\left[\underset{K^{},K}{}A_K^{}^{}A_K<\alpha ^{},p|V_{SS}^{(14)}|\alpha ,p>\right]|t,t>.$$
(21)
This expression reduces the calculation of the spin-spin matrix elements in asymmetric tops to the calculation of symmetric-top matrix elements. Solution for the latter can be found in , which allows to express the strength of mixing in ethylene by $`\widehat{V}_{SS}`$ as

$`F_{SS}(a^{}|a)`$ $`=`$ $`(2J^{}+1)(2J+1)𝒯^2|{\displaystyle \underset{K>0,q}{}}qA_{K+q}^{}A_K\left(\begin{array}{ccc}J^{}& 2& \hfill J\\ Kq& q& \hfill K\end{array}\right)`$ (27)
$`+{\displaystyle \frac{1+(1)^{J+p}}{\sqrt{2}}}A_1^{}A_0\left(\begin{array}{ccc}\hfill J^{}& 2& \hfill J\\ \hfill 1& 1& \hfill 0\end{array}\right)|^2.`$
Here $`q=\pm 1`$; (: : :) stands for the 3j-symbol and the notation was used

$$𝒯^2=2|P_{14}𝒯_{2,1}^{(14)}|^2+2|P_{C1}𝒯_{2,1}^{(C1)}|^2+2|P_{C3}𝒯_{2,1}^{(C3)}|^2;𝒯=46.5\text{kHz},$$
(28)
where $`P`$-factors are equivalent to the similar factor in (16); $`𝒯_{2,1}^{(mn)}`$ are the spherical components of the corresponding T-tensor calculated in the molecular frame. The numerical value of $`𝒯`$ was found using the molecular structure from . The selection rules for the spin-spin mixing in ethylene read

$$\mathrm{\Delta }p=0;|\mathrm{\Delta }J|2,$$
(29)
and in addition, parity of $`𝒦^{}`$ and $`𝒦`$ is opposite.

### B Spin-rotation coupling

Ortho-para conversion in molecules induced by spin-rotation coupling were studied in  (more references can be found in ). Nuclear spin-rotation coupling in molecules is due to magnetic fields produced by the molecular electrical currents. The spin-rotation perturbation can be presented as

$$\widehat{V}_{SR}=\frac{1}{2}(\underset{i}{}\widehat{𝐈}^{(i)}𝐂^{(i)}\widehat{𝐉}+h.c.);i=1,2,3,4.$$
(30)
For the spin-rotation perturbation relevant to the ortho-para mixing, the index $`i`$ should refer only to the hydrogen nuclei.

The calculation of the second rank spin-rotation tensor C is rather complicated problem which is not solved yet completely. We will use the following estimation of C. First we note that contribution to the spin-rotation coupling arising from the electric fields in the molecule has been shown to be very small compared to the part having “magnetic” origin and can be neglected . Then, we split the C-tensor arising from magnetic fields of the moving charges into two parts,

$$𝐂={}_{}{}^{e}𝐂+{}_{}{}^{n}𝐂,$$
(31)
produced by the molecular electrons and nuclei, respectively. The nuclear contribution, nC, is obtained as a first order average in the vibrational ground state, and depends only on the nuclear coordinates. Contrary, the electron part, eC, is a second order series expansion which involves the full electronic spectrum. Fortunately, $`xz`$ and $`zx`$ components which are responsible for the mixing of states having $`\mathrm{\Delta }K=1`$ in C2H4 are vanishing by symmetry requirements .

Therefore, the C-tensor, effective in ortho-para conversion, can be written (in Hz) for the $`i`$-th as

$`\text{C}^{(i)}`$ $`=`$ $`{\displaystyle \underset{ki}{}}b_k\left[(𝐫_k𝐑_k)\mathrm{𝟏}𝐫_k𝐑_k\right]𝐁;`$ (32)
$`b_k`$ $`=`$ $`2\mu _pq_k/c\mathrm{}R_k^3,`$ (33)
where Rk is the radial vector from the proton H(i) to the charge $`k`$; rk is the radial from the center of mass to the particle $`k`$; $`q_k`$ are the nuclei’ charges; B is the inverse matrix of inertia moment. B is a diagonal matrix having the elements $`B_{xx}=58.6`$ GHz, $`B_{yy}=48.7`$ GHz, and $`B_{zz}=291.9`$ GHz. Index $`k`$ runs here over all nuclei in the molecule except the proton $`i`$.

The spherical components of the spin-rotation tensor of the rank $`l`$ ($`l=1,2`$) for the $`i`$-th proton calculated in the molecular frame, $`𝒞_{l,q}^{(i)}`$, can be determined using Eq. (33). For the ethylene molecular structure from the Ref.  and bare nuclei’ charges these components are

$`𝒞_{2,1}^{(1)}=3.8\text{kHz};𝒞_{1,1}^{(1)}=2.9\text{kHz};`$ (34)
$`𝒞_{2,1}^{(3)}=4.1\text{kHz};𝒞_{1,1}^{(3)}=3.2\text{kHz}.`$ (35)
The difference between absolute values of $`𝒞_{l,q}^{(1)}`$ and $`𝒞_{l,q}^{(3)}`$ appears because of the shift of the molecular center of mass caused by the bigger mass of 13C in comparison with 12C.

Similar to the previous section, one can reduce the calculation of the matrix elements of the spin-rotation coupling in asymmetric tops to the calculation of the symmetric-top matrix elements. The latter can be found in . Thus one has the the strength of mixing due to the spin-rotation coupling in ethylene,

$`F_{SR}(a^{}|a)`$ $`=`$ $`2(2J^{}+1)(2J+1){\displaystyle \underset{i}{}}|{\displaystyle \underset{K>0,q}{}}A_{K+q}^{}A_K\mathrm{\Phi }(i;J^{},K+q|J,K)`$ (37)
$`+{\displaystyle \frac{1+(1)^{J+p}}{\sqrt{2}}}A_1^{}A_0\mathrm{\Phi }(i;J^{},1|J,0)|^2.`$
Here $`q=\pm 1`$; $`i=1,3`$ denotes the hydrogen nuclei 1 and 3. Note, that the protons 2 and 4 were also taken into account because they produce mixing equal to that of protons 1 and 3. In (37) the notation was used

$`\mathrm{\Phi }(i;J^{},K^{}|J,K)`$ $`=`$ $`{\displaystyle \underset{l}{}}\sqrt{2l+1}𝒞_{l,q}^{(i)}\left(\begin{array}{ccc}J^{}& l& J\\ K^{}& q& K\end{array}\right)\times `$ (45)
$`\left[y(J)(1)^l\left\{\begin{array}{ccc}\hfill J^{}& J& \hfill l\\ \hfill 1& 1& \hfill J\end{array}\right\}+y(J^{})\left\{\begin{array}{ccc}\hfill J& J^{}& \hfill l\\ \hfill 1& 1& \hfill J^{}\end{array}\right\}\right],`$
where {: : :} stands for the 6j-symbol; $`y(J)=\sqrt{J(J+1)(2J+1)}`$. The selection rules for the ortho-para mixing by spin-rotation perturbation in ethylene read

$$\mathrm{\Delta }p=0;|\mathrm{\Delta }J|1.$$
(46)
And again, parity of $`𝒦^{}`$ and $`𝒦`$ is opposite.

## VI Conversion rates and discussion

The decoherence rate, $`\mathrm{\Gamma }`$, is another important parameter for the isomer conversion by quantum relaxation (see Eq. (2)). This parameter is rather difficult to calculate or to measure. Due to its physical meaning the value of $`\mathrm{\Gamma }`$ is close to the rotational state population decay, which determines the line broadening in microwave rotational spectra. Unfortunately, the ethylene does not have pure rotational spectra because it has no permanent electric dipole moment. Or, to be precise, very small dipole moment in the case of 13CCH4.

In the present calculations we take the estimation of $`\mathrm{\Gamma }`$ from the line broadening in infrared spectra, for which the information is available. Recently the line broadening was accurately measured for the C2H4 molecules embedded in nitrogen . One should not expect big difference for the line broadening in the isotope species C2H4 and 13CCH4. From the data  one can estimate the decoherence rate as

$$\mathrm{\Gamma }/P=210^7\text{s}^1\text{/Torr}.$$
(47)
This rate is 10 times smaller than the $`\mathrm{\Gamma }`$ determined for the isomer conversion in CH3F . The decrease of $`\mathrm{\Gamma }`$ is because ethylene is not a polar molecule and thus has its rotational relaxation slower. Probably, the estimation (47) is too low for the case of pure ethylene gas because ethylene has polarizability larger than nitrogen. For the buffer gases having large polarizability, like Kr, or SF6, nearly three times bigger line broadening was found . Nevertheless, we will use the estimation (47) which can be considered as a low limit for $`\mathrm{\Gamma }`$.

Now we are ready to calculate the conversion rate in ethylene. The final expression for the conversion rate, $`\gamma `$, is again the Eq. (2) having the strength of mixing, $`F(a^{}|a)=F_{SS}(a^{}|a)+F_{SR}(a^{}|a)`$. The results of the calculations are given in the Table 2. From these data we note that there is no appreciable contribution to the conversion from the collapsing ortho-para level pairs presented in Fig. 4. Another observation is that the contribution from the mixing of states having $`|\mathrm{\Delta }𝒦|>1`$ is also small. This is the consequence of the fact that the ethylene molecule is rather close to a symmetric top. The main contributions to the conversion come from just two level pairs, both having the same quantum numbers ($`J^{}`$,$`𝒦^{}`$)–($`J`$,$`𝒦`$), but different $`p`$. The most important ortho-para level pair (1,27,7)–(1,28,6) has the frequency gap $``$1.0 GHz. This level pair is indicated in Fig. 1 from which one can see that this pair of states is situated at rather high energies. Nearly 10% of the conversion rate is due to the pair (0,27,7)–(0,28,6). This level pair has the same strength of mixing as the first pair but three times bigger level splitting, $`\omega `$. The properties of the states ($`1,J=28,𝒦=6`$) and ($`1,J=27,𝒦=7`$) are illustrated in Fig. 5 which presents the expansion coefficients, $`A_K`$, for these states. One can note that the biggest coefficient in both expansions has the value of $`K`$ equal to the value of $`𝒦`$ chosen for the designation of the eigenstate in our classification scheme.

The total conversion rate in 13CCH4 which is the sum over all ortho-para level pairs having $`J^{}`$ and $`J`$ up to 40, was found equal

$$\gamma /P=2.710^4\text{s}^1/\text{Torr}.$$
(48)
This rate is close to the experimental value $`(5.2\pm 0.8)10^4`$ s-1/Torr .

The calculations of the conversion rates were repeated using less accurate molecular parameters from . We obtained essentially the same results. The conversion is again determined by the same two level pairs and the value of the total conversion rate is $`3.510^4`$ s-1/Torr, which is close to the more precise value (48).

Our model allows to determine the pressure dependence of the conversion rate in 13CCH4 at the conditions of the experiment . By comparing the level splitting of the most important ortho-para level pair (1,27,7)–(1,28,6), which is $``$1.0 GHz, with the decoherence rate (47) one can conclude that the case corresponds to the limit $`\mathrm{\Gamma }\omega `$. From Eq. (2) one can see that in this limit linear grows of the conversion rate, $`\gamma `$, versus gas pressure should take place. (Note that $`\mathrm{\Gamma }`$ is proportional to the gas pressure.) This dependence was indeed observed in the experiment .

Magnitude of the calculated conversion rate relies heavily on the determination of the ortho-para level splitting. The gap of the most important level pair appeared to be rather large, $``$1.0 GHz, and is unlikely that it has been determined with significant error. This can be proven by comparing this splitting with the one calculated using less accurate molecular parameters . These parameters give the gap 0.914 GHz, which is close to the value given in the Table 2.

The two sources of the ortho-para mixing have been analysed. The spin-spin interaction between the molecular nuclei was possible to calculate rather accurately. The uncertainty of this perturbation is mainly due to small errors in the knowledge of the molecular structure. On the other hand, our analysis has shown that the spin-spin coupling contributes less than 1% to the ethylene conversion and thus negligible. There are a few reasons for this reduction. First, the spin-spin tensor has relatively small value due to the big distances between interacting nuclei and small magnitude of the magnetic moment of the carbon nucleus. Second, due to high symmetry of the molecular structure interactions between many pairs of protons do not contribute to the ortho-para mixing. Moreover, there are no close ortho-para level pairs at small $`J`$, where spin-spin coupling can compete with the spin-rotation coupling.

We have taken into account also the spin-rotation perturbation. The magnitude of the spin-rotation tensor was estimated and it was demonstrated that the main contribution to the ethylene conversion comes from the spin-rotation mixing. Strong contribution from the spin-rotation coupling is due to the fact that the strength of the mixing $`F_{SR}`$ grows as $`J^3`$ in comparison with the slower ($`J^2`$) grows of the spin-spin strength of mixing $`F_{SS}`$. This makes the “accidental” resonance at big $`J`$ (1,27,7)–(1,28,6) being so important. An extra cause of the efficiency of the spin-rotation coupling is that the spin-rotation tensors of first and second rank contribute both to the conversion. On the other hand, not all close ortho-para level pairs contribute to the conversion. In the particular case of the collapsing ortho-para levels (Fig. 4) there is no strong mixing mainly because of interferences between many symmetric-top components with different quantum numbers $`K`$, and also because for these $`J^{}=J`$ pairs the first rank spin-rotation tensor does not contribute.

Error in the estimation of $`\mathrm{\Gamma }`$ introduces systematic uncertainty into the final result. But it does not undermine the model as a whole. First, one can turn around the approach and determine which decoherence rate $`\mathrm{\Gamma }`$ should be chosen in order to fit the experimental data. In order to reproduce the experimental rate of conversion one should take the decoherence rate, $`\mathrm{\Gamma }/P=3.910^7`$ s-1/Torr, which is less than two times different from our rough estimation (47). Of course, such approach will be justified only when the spin-rotation perturbation in ethylene will be known from an independent source. There is another point which is worth to mention. Recently it was proposed a new approach to the isomer conversion in which the knowledge of $`\mathrm{\Gamma }`$ is not necessary at all . It is based on a fast linear sweep of the molecular levels through the ortho-para resonance. Nuclear spin conversion in such experimental arrangement does not depend on $`\mathrm{\Gamma }`$ and is dependent on the strength of ortho-para mixing and the population of the mixed states.

Finally, it seems interesting to emphasize the difference in the conversion rates produced by the two pairs of identical rotational quantum numbers (27,7)–(28,6) but different in parity. This peculiarity arises from the asymmetric top properties, where energy of states depend on parity. This in turn has dynamical consequence that the spin isomers equilibrate first in one parity manifold (here the odd one).

## VII Conclusions

The theoretical model for the nuclear spin conversion in ethylene (13CCH4) has been developed. For the first time a theory of the spin isomer conversion in asymmetric tops was possible to compare with experimental data on the isomer conversion. We have found that the two experimental results : the magnitude of the ethylene (13CCH4) isomer conversion rate and its pressure dependence, are consistent with the spin conversion governed by quantum relaxation. The ortho-para state mixing is performed in this molecule mainly by the coupling between the protons’ spins and the molecular rotation. We have identified also the two pairs of ortho-para states which are almost completely determine the spin conversion in 13CCH4

## Acknowledgments

The authors are indebted to A. Fayt for the possibility to use the latest set of the ethylene molecular parameters prior to publication and to M. Irac-Astaud for stimulating discussions.

Table 1. The character table for the symmetry groups C2v, C2v(M), D2 and classification of the basis states (8).

| C$`{}_{2v}{}^{}(M)`$ | E | (12)(34) | E | (12)(34) | ortho | ortho | para |
| --- | --- | --- | --- | --- | --- | --- | --- |
| C2v | E | C2 | $`\sigma _v`$ | $`\sigma _v^{}`$ |  |  |  |
| D2 | E | R$`{}_{z}{}^{}{}_{}{}^{\pi }`$ | R$`{}_{y}{}^{}{}_{}{}^{\pi }`$ | R$`{}_{x}{}^{}{}_{}{}^{\pi }`$ | $`K`$-even | $`K`$=0 | $`K`$-odd |
| A1 | 1 | 1 | 1 | 1 | $`p`$=0 | $`J,p`$-even, |  |
| B2 | 1 | -1 | -1 | 1 |  |  | $`p`$=1 |
| A2 | 1 | 1 | -1 | -1 | $`p`$=1 | $`J,p`$-odd, |  |
| B1 | 1 | -1 | 1 | -1 |  |  | $`p`$=0 |

Table 2. The most important ortho-para levels and their contributions to the spin conversion in ethylene. In calculations the molecular parameters from  were used. The total rate combines the contribution from all ortho-para level pairs having $`J40`$.

| Level pair | Energy | $`\omega /2\pi `$ | $`F_{SS}`$ | $`F_{SR}`$ | $`\gamma /P`$ |
| --- | --- | --- | --- | --- | --- |
| $`p^{},J^{},𝒦^{}`$-$`p,J,𝒦`$ | (cm-1) | (MHz) | (MHz2) | (MHz2) | (10-4 s-1/Torr) |
| 1,27,7–1,28,6 | 871.53 | 1006.7 | 4.2$`10^2`$ | 4.4 | 2.49 |
| 0,27,7–0,28,6 | 871.69 | -3435.5 | 4.2$`10^2`$ | 4.4 | 0.21 |
| 0,4,3–0,6,2 | 53.72 | -3180.6 | 1.1$`10^3`$ | 0 | 3.4$`10^3`$ |
| 1,26,3–1,24,6 | 680.68 | 944.3 | 3.4$`10^4`$ | 0 | 5.6$`10^4`$ |
| 0,20,1–0,20,0 | 356.02 | 1894.4 | 7.3$`10^8`$ | 2.4$`10^6`$ | 4.8$`10^6`$ |
| Total rate |  |  |  |  | 2.70 |