text:
# Bifurcations in annular electroconvection with an imposed shear

## I Introduction

Phenomological amplitude equation models are important tools in the study of symmetry-breaking, pattern-forming bifurcations. They have proven particularly useful at the weakly nonlinear level. While their broad application stems from the common symmetries that are shared by many different systems, it is the coefficients or specific parameters in amplitude equations that distinguish each system. It is an interesting question to examine how these coefficients, especially those of the nonlinear terms, change as the properties of the system are varied. There are only a few experimental systems, such as Rayleigh-Bénard convection, Taylor vortex flow, and electrohydrodynamic convection in nematic liquid crystals, for which these coefficients have been quantitatively determined.

Three-dimensional, nonlinear systems are prone to develop complicated spatial and temporal patterns even when only weakly nonequilibrium. The spatio-temporal patterns are often the result of successive symmetry breaking bifurcations. It is useful to study pattern formation in low-dimensional systems which are close to equilibrium but have little symmetry so that there are only a very limited set of symmetry-breaking bifurcations available. In general, one seeks the most complex dynamics that can be realized in as simple and restricted a system as possible. In our system, annular electroconvection, we exploit the strict two-dimensionality of a submicron smectic A liquid crystal film. The lower dimensionality greatly reduces the variety of possible pattern states and so makes it easier to experimentally study the rich nonlinear properties of the basic pattern.

In this paper, we report an experimental study of the properties of the bifurcations in annular electroconvection with a variable circular Couette shear. Our primary focus is to investigate the variation of the cubic nonlinear coefficient of an amplitude equation as the experimental parameters are systematically changed. This coefficient determines the saturation behaviour of the convective flow velocity. We begin by establishing that the experimental system is well-described by a theoretical model that reduces to an amplitude equation near the onset of the pattern-forming instability. We test the adequacy of the theoretical model by direct comparisons between the experimental results and theoretical predictions at the linear level. We then proceed to describe how the coefficient of the lowest order nonlinear term depends on the experimental parameters. Interestingly, we find that the coefficient can be made to change sign and thus that the primary bifurcation can be tuned between subcritical and supercritical. In this paper, we show from symmetry considerations that the amplitude equation that describes the sheared case is the complex Ginzberg-Landau equation. Elsewhere, we have derived this result directly from the full electrohydrodynamic equations using a multiple-scales analysis. Calculation of the parameter dependence of the coefficients of this equation is beyond the scope of this paper, but in future, we expect to be able to compare the experimental results presented here with theory at the weakly nonlinear level.

Our system consists of a freely suspended annular smectic A film which can be driven out of equilibrium by electrical forces and can be independently subjected to a shear flow as shown in Fig. 1. It was previously described in Refs.  and . The film was suspended between concentric circular electrodes and was driven by a voltage difference between the electrodes. The inner electrode could be rotated about its axis, thereby imposing a circular Couette shear. The flow pattern that emerges when the film is driven sufficiently hard is referred to as electroconvection and consists of an array of counter-rotating vortices. Our primary probe of the electroconvection amplitude was the excess electrical current carried convectively by the flow. The annular geometry and electrical nature of the experiment make it possible to study radially forced convection under steady shear, a combination more often associated with geophysical flows than with small laboratory systems. On the other hand, the accurately two-dimensional nature of the flow makes it far simpler than most planetary analogs, and thus more amenable to theoretical investigation.

The driving force originates from the interaction of the radial electrical field and the surface charge density that develops on the film’s free surfaces. Below the onset of convection, the film’s electrical state is independent of the imposed shear in the film. It is thus straightforward to superpose shear flows with the radial driving, as these are characterized by independent control parameters. The simplest such shear flow is steady azimuthal or circular Couette flow. The annular geometry is naturally periodic so that the shear flow is closed and leads to a net mean flow around the annulus.

The addition of a shear flow to the annular electroconvection alters the symmetries of the base state of electroconvection. When shear is absent, the base state is invariant under azimuthal rotation and reflection in any vertical plane containing the rotation axis through the center of the annulus. The electroconvective state is then stationary and appears with the spontaneous breaking of the azimuthal invariance. When sheared, the reflection symmetry of the base state is not present. When driven, the pattern once again breaks the azimuthal symmetry, but since the reflection symmetry is absent due to shear, the pattern is free to travel azimuthally in the direction of the mean flow. In addition, as we show in detail below, the shear flow alters the primary bifurcation, making it hysteretic.

The equilibrium properties of freely suspended liquid crystal films have been extensively studied; see Ref.  for a recent review. Smectic liquid crystals consist of layers of orientationally ordered long molecules which readily form suspended films. In smectic A, the average orientation of the long axis of the molecules, and hence the optic axis, is normal to the layer plane. The layers are of uniform thickness and within each layer the distribution of molecules is isotropic. Smectic A exhibits two-dimensional isotropic fluid properties in the layer plane while flows perpendicular to the layers are strongly inhibited. Other material properties such as the electrical conductivity and dielectric permitivity are also isotropic in the plane of the layers. Uniform suspended smectic films are always an integer number of smectic layers thick and while they flow readily, they seldom change thickness;

they are robustly two-dimensional. Unlike soap solutions, smectic liquid crystals have very small vapor pressures, hence the film can be enclosed in an evacuated environment. The reduction in ambient pressure leads to a proportional reduction in the air drag to which the film is subjected. This system thus has several attractive features for the study of bifurcations in simple nonlinear systems that may be described by amplitude equation models. Our study greatly extends and complements previous work on electroconvection in isotropic liquids, nematics, smectic A and higher smectic phases.

The paper is organized as follows. In Section II A we describe the results of a linear stability analysis of our system. In the process we define the various dimensionless parameters that characterize the experiment. In Section II B, we use symmetry arguments to justify the form of the relevant amplitude equation. Our description of the experimental apparatus and protocol follows in Section III. In Section IV A we outline our data analysis procedure. We present experimental results that pertain to linear analysis in Section IV B. An experimental determination of the cubic coefficient of the amplitude equation is found in Sections IV C and IV D. Section V discusses the relationships between our results and those for other similar systems and section VI presents a brief summary and conclusion.

## II Theoretical Background

In this section, we briefly review the linear and weakly nonlinear theory that is relevant to our data analysis and results.

A linear stability analysis of annular electroconvection with circular Couette flow was reported in Ref. . The theory is constructed with a number of simplifying assumptions. The fluid film is treated as an annular sheet of inner radius $`r_i`$, outer radius $`r_o`$ and thickness $`s`$. A schematic is shown in Fig. 1. An important geometric parameter is the radius ratio $`\alpha =r_i/r_o`$. The film width $`d=r_or_i`$ is assumed to be much greater than $`s`$. The viscous fluid is assumed to flow only in two-dimensions, and be incompressible, Ohmic and of uniform electrical conductivity. We denote the fluid density by $`\rho `$, its molecular viscosity by $`\eta `$ and its electrical conductivity by $`\sigma `$. Only the charges due to free surfaces are included and bulk dielectric effects inside the film are neglected. The electrodes are assumed to be of negligible thickness, and fill the rest of the plane not occupied by the annular film. Outside the plane of the film, we assume an empty space of dielectric permittivity $`ϵ_0`$. Couette shear in the theoretical model is imposed by specifying the appropriate velocity boundary conditions on the edges of the film. Linear theory predicts the position of the onset of convection and the degree to which onset is suppressed by the shear. These are discussed in detail in Ref.  and are further analyzed in Section IV B below.

From very general symmetry considerations, we can deduce the form of an amplitude equation which is valid in the weakly nonlinear regime just above onset. Here, a complex-valued amplitude describes the slowly-varying magnitude and phase of the spatially periodic physical fields, for example the fluid velocity. We find an expression of the well-known Landau form. The coefficients of the amplitude equation can in principle be calculated from the underlying electrohydrodynamic equations; a calculation of this kind is in progress and will be reported separately. The magnitude of the complex amplitude is directly related to the quantity we measured experimentally, the total current carried by the film. We can thus fit current-voltage data to the real part of the amplitude equation to experimentally determine its coefficients. This is done in Sections IV C and IV D below.

### A Outline of linear stability theory

The theoretical treatment of Ref.  examined the linear stability of the two-dimensional annular fluid film when it is subjected to a voltage $`V`$ at the inner electrode of the annulus while the outer electrode is held at ground potential. In addition, the theory allowed for the rotation of the inner edge of the annulus at angular frequency $`\omega `$, while the outer edge is held fixed. The rotation drives a Couette shear flow in the film below the onset of convection. It can be shown that any arbitrary rotation of the inner and outer edges can be transformed to a rotation of only the inner edge by moving to the appropriate rotating coordinates. Such a coordinate transformation has no effect on the dynamics as long as the flow is strictly two-dimensional. The theoretical model is completely specified by the radius ratio $`\alpha `$ and three additional dimensionless parameters:

$$\frac{ϵ_0^2V^2}{\sigma \eta s^2},𝒫\frac{ϵ_0\eta }{\rho \sigma sd},\mathrm{and}e\frac{\rho \omega r_id}{\eta }.$$
(1)
The control parameter $``$ is a measure of the external driving force which is proportional to the square of the applied voltage $`V`$ while $`𝒫`$ is a ratio of the time scales of electrical and viscous dissipation processes in the film. The Reynolds number $`e`$ is a measure of the strength of the applied shear, which is regarded as a second control parameter. It has been established that the instability leads to a one-dimensional pattern of $`m`$ vortex pairs where $`m`$ is the azimuthal mode number. Linear stability predicts the value of the critical control parameter $`_c`$ and the critical mode number $`m_c`$ at which the film becomes marginally unstable. $`_c`$ and $`m_c`$ are in general functions of $`\alpha `$, $`𝒫`$, and $`e`$. A detailed discussion of the calculated values of $`_c`$ and $`m_c`$ under various conditions is given in Ref. . When $`e=0`$ (i.e. with no applied shear), it was found that $`_c`$ and $`m_c`$ are independent of $`𝒫`$. In the following, these critical values for zero shear will be denoted $`_c^0`$ and $`m_c^0`$.

When $`e>0`$, it was found that $`_c>_{c}^{}{}_{}{}^{0}`$, for any $`\alpha `$ and $`𝒫`$, i.e. that the shear always suppresses the onset of convection. It is convenient to measure the relative degree of suppression for a given $`\alpha `$ and $`𝒫`$ in terms of the reduced quantity

$$\stackrel{~}{ϵ}(\alpha ,e,𝒫)=\left[\frac{_c(\alpha ,e,𝒫)}{_c^0(\alpha )}\right]1,$$
(2)
which is a monotonically increasing function of $`e`$ for all $`𝒫`$. In Section IV B below, we extend the discussion in Ref  by comparing these results to data for several different $`\alpha `$, and for a wide range of $`𝒫`$.

### B Nonlinear theory: The amplitude equation

The amplitude equation appropriate to our system follows from symmetry considerations. We defer a detailed calculation of its numerical coefficients to future work.

In the absence of shear, the base state solution of the electrohydrodynamic equations is invariant under azimuthal rotations and reflections through any plane perpendicular to the annular plane and containing its center. Just above onset, we assume that all the physical fields are nonaxisymmetric and proportional to a rapidly-varying spatial oscillation $`e^{im\theta }`$, where $`m`$ is an azimuthal mode number. This requires that the differential equation for the slowly-varying complex amplitude $`A_m`$ of mode $`m`$, be unchanged under the transformations

$$\theta \theta +\theta ^{},A_mA_me^{im\theta ^{}},$$
(3)
and

$$\theta \theta ,A_m\overline{A}_m.$$
(4)
In the above, $`\theta `$ is an azimuthal angle, and the overbar denotes complex conjugation. The most general amplitude

equation that is invariant under these symmetry operations has the Landau form,

$`\tau _tA_m`$ $`=`$ $`ϵA_mg|A_m|^2A_m+h|A_m|^4A_m\mathrm{},`$ (5)
where $`\tau `$, $`g`$, and $`h`$ are real-valued coefficients. The small parameter $`ϵ`$ is a reduced control parameter given by

$$ϵ=\left[\frac{}{_c}\right]1,$$
(6)
and is a measure of the distance from threshold. The amplitude equation is accurate for small $`ϵ`$ and describes a bifurcation from the $`A_m0`$ state ($`ϵ<0`$) to the $`A_m0`$ state with $`2m`$ vortices ($`ϵ>0`$). A sharp bifurcation occurs at $`ϵ=0`$.

If the fluid is subjected to a circular Couette shear, the base state is only invariant under azimuthal rotations, Eqn. 3. In this case, the general amplitude equation takes the complex Landau form

$`\tau (_tia_{Im})A_m`$ $`=`$ $`ϵ(1+ic_0)A_mg(1+ic_2)|A_m|^2A_m`$ (7)
$`+`$ $`h(1+ic_3)|A_m|^4A_m\mathrm{},`$ (8)
where $`a_{Im}`$ is the imaginary part of the eigenvalue of the unstable mode $`m`$ at onset and the coefficients $`c_0`$, $`c_2`$, $`c_3`$ are real-valued. Eqn. 8, which here is simply deduced from symmetry considerations, has been rigorously derived for annular electroconvection with shear from the basic electrohydrodynamic equations. In general, it describes a Hopf bifurcation to a pattern with a complex amplitude $`A_m`$ which travels in one azimuthal direction.

In both Eqns. 5 and 8, we have omitted terms involving the azimuthal gradients of $`A_m`$, which if included would give the so-called Ginzburg-Landau equation. While such terms are potentially quite interesting, they do not appear to be directly relevant to interpreting our current-voltage data. Our annular system is azimuthally periodic and thus lacks lateral boundaries at which $`A_m0`$ where gradients would be important. Thus, any gradients in $`A_m`$ would have to arise from purely dynamical effects, such as Eckhaus or other phase instabilities.

From current-voltage data, we can extract the reduced Nusselt number $`n`$, which is a global dimensionless quantity defined as the ratio of current carried by convection to that by conduction. The reduced Nusselt number $`n`$ is related to the electric Nusselt number $`𝒩`$ by $`n=𝒩1`$. One can show that the amplitude can be scaled so that $`n=|A_m|^2`$.

To solve for the real and imaginary parts of the amplitude equation, let $`A_m(t)=A(t)e^{i\mathrm{\Phi }(t)}`$, where $`A(t)`$ and $`\mathrm{\Phi }(t)`$ are real. Substitution into Eqn. 8, gives

$`\tau _tA`$ $`=`$ $`ϵAgA^3hA^5\mathrm{},`$ (9)
$`\tau (_t\mathrm{\Phi }a_{Im})`$ $`=`$ $`ϵc_0gc_2A^2+\mathrm{}.`$ (10)
We will show in Section IV A how the raw current-voltage data can be transformed into measurements of $`ϵ`$ and $`n`$, and hence be used to measure the real amplitude $`A`$ via a nonlinear fit procedure. We focus on extracting the coefficient $`g`$ of the cubic nonlinearity. The sign of $`g`$ determines whether the pitchfork bifurcation to electroconvection is supercritical (“forward”, $`g>0`$ ) or subcritical (“backward”, $`g<0`$ ). A tricritical bifurcation occurs when $`g=0`$. The results for $`g`$ are discussed in Sections IV C and IV D.

## III Experimental Apparatus

In this Section we describe the liquid crystal sample and the experimental apparatus. We performed a large number of precise current-voltage measurements on electroconvecting annular films with a wide range of radius ratios and applied shears. In overall structure, the apparatus consisted of the annular electrodes which were housed in a vacuum chamber, an electrometer circuit and a stepper motor, both under computer control. The chamber also allowed films to be drawn under vacuum and inspected under reflected white light to check thickness uniformity. A schematic is shown in Fig. 2.

### A The liquid crystal

As in previous experiments, our study used smectic A octylcyanobiphenyl (8CB). 8CB has a smectic A phase between $`21^{}\mathrm{C}`$ and $`33.5^{}\mathrm{C}`$. The electrical conductivity of “pure” 8CB is due to several ionic impurities of varied and unknown concentrations. In order to control the nature of the ionic species contributing to the electrical conductivity, we doped the 8CB with tetracyanoquinodimethane (TCNQ), a material believed to form charge transfer complexes with the host, so that the dominant impurity species was the dopant. To prepare the doped material, TCNQ was dissolved in acetonitrile and added to the 8CB sample. The acetonitrile was then evaporated in a vacuum oven while warming the mixture so that the 8CB was in its isotropic phase. The samples used had concentrations of $`2.96\times 10^4`$, $`1.11\times 10^4`$ and $`7.62\times 10^5`$, of TCNQ by weight. Experiments with significantly higher or lower dopant concentrations were found to be less useful due to irreproducible non-ohmic behaviour below the onset of convection.

### B the annulus

The annular electrodes were constructed out of stainless steel. The inner electrode was a circular disk of diameter $`2r_i`$. The outer electrode was a circular plate of diameter $`9.00`$ cm with a central hole of diameter $`2r_o`$. The outer electrode was $`0.73\pm 0.01`$ mm thick. By using several pairs of inner and outer electrodes, we conducted experiments at six different radius ratios between $`\alpha =0.33`$ and $`\alpha =0.80`$. We used inner electrodes with radii ranging between $`r_i=3.60\pm 0.01`$ mm and $`r_i=5.26\pm 0.01`$ mm. The radii of the outer electrodes were between $`r_o=5.57\pm 0.01`$ mm and $`r_o=11.25\pm 0.01`$ mm.

The annulus was housed in a vacuum chamber which was evacuated by a rotary pump. The air was evacuated slowly so as to prevent vigorous air flows that may cause the film to rupture. The surrounding air was pumped down to an ambient pressure between $`0.15.0`$ torr. At these pressures, the mean free paths of $`N_2`$ and $`O_2`$ are in the range $`0.50.05`$ mm. This was comparable to the film width $`d`$, so that the air drag on the film was negligible. All experiments were performed at the ambient room temperature of $`23\pm 1^{}`$C, well below the smectic A-nematic transition at $`33.5^{}`$C for undoped 8CB.

The inner electrode was made to rotate about its axis by means of a high precision stepper motor operating at $`25600`$ steps per revolution. The motor was located outside the vacuum chamber and was connected through a rotating seal. The inner electrode was adjusted to rotate true to its axis to within $`50\mu `$m at angular frequencies up to $`\omega =6\pi `$ rad/s.

### C Drawing the film and determining its thickness

A motorized film-drawing assembly was housed in the vacuum chamber. The film was spread across the annular gap between the electrodes using a stainless steel razor blade inclined at $`25^{}`$. The blade was moved away from the annulus while voltage was applied, to avoid perturbing the electric fields near the film.

The thickness of the film was determined from the interference color of the film under reflected white light using a low power videomicroscope. Using standard colorimetric functions a color-thickness chart was calculated for 8CB. Since smectic films are constrained to be an integer number of smectic layers thick, the film color can be used to identify the film thickness measured in smectic layers. Each layer of smectic A 8CB is $`3.16`$ nm thick. Most of the experiments were performed with films between $`25`$ and $`85`$ layers thick. Over most of the middle of this range, the film thickness can be determined to within $`\pm 2`$ layers, while close to the ends of the range a more conservative estimate of $`\pm 5`$ layers was used. The colorimetric method worked well for film thicknesses where the strong interference color can be unambiguously matched to a color chart. Very thin films that appeared black and thicker, nearly white films were not used.

During the course of an experiment, the film color was monitored to check that it remained uniform in thickness to within $`\pm 1`$ layer. Films drawn with non-uniform thicknesses, left to themselves, tend to anneal to become uniform in thickness and hence uniform in color. The annealing process can be accelerated by electroconvecting and shearing the films.

Current-voltage runs during which the film thickness spontaneously became non-uniform were abandoned. A small number of visualization experiments were performed on films which had a nonuniformity in thickness of about $`\pm 2`$ smectic layers, or about $`5\%`$ of $`s`$. Usually these nonuniform films had two thicknesses and in reflection displayed two nearby interference colors. The advection of the thickness nonuniformities was used to visualize the flow pattern and thus to determine the azimuthal mode number $`m`$ by counting vortices. For simplicity, some of these observations were done at atmospheric pressure. This visualization method, while somewhat crude, was preferrable to using suspended smoke particles, which tend to perturb the conductivity of the film.

### D Current-voltage measurements

Since the current transported through the film was picoamperes in magnitude, particular care had to be exercised to avoid stray currents. The inner rotating electrode was electrically isolated from all other components by a highly insulating sleeve. The current was carried onto the rotating electrode by a set of silver-graphite slip rings inside the vacuum chamber. A Keithley model 6517 electrometer was used both as a voltage source and a picoammeter. The ‘high’ of the programmable dc voltage source of the electrometer was connected to the rotating inner electrode, which is the only part of the entire apparatus that is not at ground potential. The outer annulus was connected to the input of the electrometer which is held at very nearly ground potential. The outer annulus was electrically isolated from true ground by teflon washers which served to eliminate leakage currents which would otherwise be added to the signal. The rest of the apparatus was grounded. Electrical noise was reduced by shielding the electrodes and most of the

experimental components in a large Faraday cage which doubled as the vacuum chamber. Low noise triaxial cables and feedthroughs were used to connect the outer electrode to the electrometer.

An example of a current-voltage characteristic in a film without shear is shown in Fig. 3a. It consists of data obtained for incremental and decremental voltages. The current-voltage characteristic clearly shows two regions: one for voltages smaller than a critical voltage $`V_c`$ and one for voltages greater than $`V_c`$. The critical voltage that separates these two regions is in the vicinity of the kink in the current-voltage characteristic. In the regime $`VV_c`$, the current is linearly dependent on the voltage and the film is ohmic. Experiments in 8CB with significantly higher and lower concentrations of TCNQ show, at least initially, non-ohmic current-voltage characteristics even for $`VV_c`$, due to electrochemical effects.

We estimated $`V_c`$ before each run then defined an approximate reduced control parameter $`ϵ=(V/V_{c}^{}{}_{}{}^{est})^21`$. Current-voltage data was then obtained by making variable steps in voltage in such a way as to make equal increments in $`ϵ`$. We used several $`ϵ`$ step sizes, which became finer near threshold. Similarly, we used a variable waiting time after each step, which allowed for longer relaxation times closer to $`V_c`$. Very long relaxation times were not feasible due to the drift of the electrical conductivity. After the wait time, between $`100200`$ measurements of the current, each separated by $`25`$ milliseconds, were averaged. The error in the average was taken to be the standard deviation of the mean or $`1\%`$, the reading error of the picoammeter, whichever was greater. Each run consisited of incrementing $`ϵ`$ up to a predetermined maximum and then decrementing it to zero voltage again. Once the data had been acquired, more precise $`ϵ`$ values at each voltage were found using drift-corrected values of $`V_c`$ from a fit procedure described in the next section.

When the film was sheared by rotation of the inner electrode, it was allowed at least $`30`$ seconds after a change of shear rate to attain a steady state before the current-voltage characteristics were obtained. In all runs, the shear rate was established and then held fixed during a current-voltage sweep. Fig. 3b displays a representative current-voltage characteristic in a film under shear. The principle effects of the shear are to suppress the onset of convection and, eventually, to make the primary bifurcation hysteretic.

## IV Experimental Results

With the exception of some simple flow visualization which used slightly nonuniform films, all our results were obtained by fitting current-voltage data taken with uniform films. It is possible to fit the data in such way as to extract all the unknown material parameters in a nearly model-independent way, as well as to correct for small drifts in conductivity.

The first task was to determine the dimensionless parameters relevant to each current-voltage sweep: the radius ratio $`\alpha `$, the dimensionless ratio $`𝒫`$ and the Reynolds number $`e`$. We can then properly non-dimensionalize the current-voltage data and fit it to an amplitude equation model. The details of this fit are described in the section IV A below.

In Section IV B, we consider the features of the data that are predicted by the linear stability theory developed in Ref. . For zero shear, these are the critical mode number $`m_c^0`$, and the critical voltage $`V_c^0`$. We use the latter to fix the last unknown material parameter, the viscosity $`\eta `$. We also compare our data to the linear stability prediction for the shear-suppression of the onset of convection. This was done for six values of $`\alpha `$, and for wide ranges of $`𝒫`$ and $`e`$.

In the weakly nonlinear regime, the value of the coefficient of the cubic nonlinearity, $`g`$ in the amplitude equation, was of primary interest. In Section IV C we report results for $`g`$ at various $`\alpha `$ and $`𝒫`$ for $`e=0`$, that is, in the absence of shear. We find that $`g`$ can become slightly negative for small $`\alpha `$, so that the bifurcation is weakly backward. The result for $`g`$ in the limit $`e=0`$, $`\alpha 1`$ is compared to the theoretical result for rectangular films. Finally, in Section IV D we describe the results for $`g`$ in sheared films, where $`e>0`$. We find that increasing $`e`$ has a strong effect on the nature of the bifurcation, driving $`g`$ negative and making the bifurcation backward.

### A Fitting the current-voltage data

Except for a small drift, the film is ohmic below the onset of convection. We used fits in this range to determine some important material parameters in order to scale the data.

Each voltage-current measurement in the ohmic regime, $`(V,I)`$, constitutes an experimental determination of the film’s conductance, $`c=I/V`$. For a film of radius ratio $`\alpha `$, thickness $`s`$ and conductivity $`\sigma `$, the conductance is given by

$$c=\frac{2\pi \sigma s}{\mathrm{ln}(1/\alpha )}.$$
(11)
Interestingly, the conductance is independent of the size of the film, i.e. independent of $`r_i`$ or $`r_o`$. Since $`\alpha `$ is merely a geometrical parameter, measurements of $`c`$ are effectively measurements of $`\sigma s`$, eliminating the need to determine these separately.

Using Eqn. 11, $`𝒫`$, defined in Eqn. 1, can be expressed in terms of $`c`$ as

$$𝒫=\frac{ϵ_0\eta }{\rho \sigma sd}=\frac{2\pi ϵ_0\eta }{\rho (r_or_i)\mathrm{ln}(1/\alpha )}\frac{1}{c},$$
(12)
where $`d=r_or_i`$. Similarly, the Reynolds number $`e=\rho \omega r_id/\eta `$ can be calculated from the measured angular frequency $`\omega `$ of the inner electrode in rad/s, given $`\rho `$ and $`\eta `$. The density $`\rho `$ of 8CB at room temperature is $`1.0\times 10^3`$ kg/m<sup>3</sup>. In both $`𝒫`$ and $`e`$ , the only remaining undetermined material parameter is the viscosity $`\eta `$. We fixed this parameter using a combination of experimental and theoretical results, as described in Section IV B. The result was $`\eta =0.18\pm 0.03`$ kg/ms.

Thus, for each current-voltage sweep, we were able to determine $`𝒫`$ and $`e`$. The drift in conductivity caused a drift in $`c`$ and hence in $`𝒫`$ over the course of a sweep. Furthermore, $`c`$ is only known directly during the ohmic parts at the beginning and end of each sweep. We associated a mean $`𝒫`$ with each sweep by averaging over $`c`$ data before and after. This results in $`10\%`$ errors in $`𝒫`$ for any one sweep, and a slow, uncontrolled evolution of $`𝒫`$ from sweep to sweep.

The whole range of current-voltage data is fit to the amplitude equation model. Eqn. 8 gives the full, time-dependent complex coefficient amplitude equation required by symmetry. Only the real and steady state part of Eqn. 8, given by

$$ϵAgA^3hA^5+f=0,$$
(13)
is required to model the current-voltage data. Here we have truncated at the quintic order and augmented Eqn. 8 with a phenomenological field term $`f`$. This term models the rounding of the bifurcation due to nonideal systematic effects which are slightly symmetry-breaking, resulting in an “imperfect” bifurcation. For all of the fits, we found $`f1`$.

As discussed in section II B, the amplitude $`A`$ which appears in Eqn. 13 is related to the reduced Nusselt number $`n`$ by

$$n=\frac{I}{I_{cond}}1=\frac{I}{cV}1=A^2,$$
(14)
while the reduced control parameter $`ϵ`$ is given by

$$ϵ=\frac{}{_c}1=\left(\frac{V}{V_c}\right)^21.$$
(15)
Eqns 13-15 contain five parameters, $`c`$, $`V_c`$, $`g`$, $`h`$, and $`f`$. While $`c`$ has a nearly constant initial slope and $`V_c`$ is marked by a relatively obvious kink in the raw $`(I,V)`$ data, the amplitude $`A`$ is indirectly deduced via the pair of transformations Eqns. 15 and 14 which are nonlinear. This amplitude is in turn fit using Eqn. 13, which is again nonlinear. Thus, the parameters $`(g,h,f)`$ are rather distantly related to the raw $`(I,V)`$ data. If $`g<0`$, the bifurcation is hysteretic, and $`A`$ is multivalued over some ranges of $`ϵ`$. Also, the nature of the model necessarily involves several fit parameters which are not independent. Consequently, the determination of these parameters is much more difficult than $`V_c`$. As discussed in Sections IV C and IV D, $`(g,h,f)`$ can be influenced by small systematic effects in the data leading to scatter which is larger than the statistical uncertainties in the fits. Nevertheless, the general trends are robust.

The drift of the film’s conductivity, and hence its conductance $`c`$, introduces some additional complications. Recall that the onset of convection occurs when the control parameter $``$ equals a critical value $`_c`$ given by

$$_c=\frac{ϵ_0^2V_{c}^{}{}_{}{}^{2}}{\sigma \eta s^2}\mathrm{or}V_c=\frac{s}{ϵ_0}\sqrt{_c\sigma \eta }.$$
(16)
Since the electrical conductivity is not constant, the value of $`V_c`$ required to calculate $`ϵ`$ slowly changes during the course of an experiment. Combining Eqns. 11 and 16 yields

$$V_c(c)=\sqrt{\frac{s\eta _c\mathrm{ln}(1/\alpha )c}{2\pi ϵ_{0}^{}{}_{}{}^{2}}}.$$
(17)
We corrected for the drift in $`c`$ by tracking its value during the ohmic parts of the sweep and using a linear interpolation during the portion of the sweep where the film is convecting. In this way, each voltage and current measurement was non-dimensionalized with its own value of $`c`$ and $`V_c(c)`$ to produce fittable values of $`ϵ`$ and $`A`$. We also transformed the errors in the current measurements $`\mathrm{\Delta }I`$ into amplitude errors $`\mathrm{\Delta }A`$.

We approached the fitting problem in the following bootstrap fashion. It was easy to ascertain bounds on $`V_c`$ by inspecting the current-voltage characteristics, see for example Figs. 3a and b. Each current-voltage characteristic was scrutinized and two voltage intervals were chosen. The first interval contained the critical voltage at which the conducting state became unstable to the convecting state. The second interval contained the voltage at which the convecting state became marginally unstable to the conduction state. Guesses of $`V_c`$ in both these intervals were chosen at random using a uniform deviate random number generator. The conductances $`c`$ were then determined for the ohmic regimes and from a linear interpolation for the convecting portion, as described earlier. The raw current-voltage data was then transformed into $`ϵ`$ and $`A\pm \mathrm{\Delta }A`$ and fit to Eqn. 13 by varying only the three parameters $`g`$, $`h`$ and $`f`$ in a weighted Levenberg-Marquardt nonlinear fitting procedure which minimized $`\chi ^2`$. We then did a Monte-Carlo optimization of the randomly chosen $`V_c`$.. The best $`V_c`$ was taken to be the one that minimized $`\chi ^2`$ over the Monte Carlo sample. The uncertainty in $`V_c`$ was the standard deviation of the uniform deviate on the constraint interval. Corresponding to this best $`V_c`$ were the three best fit parameters

$`g`$, $`h`$ and $`f`$. The uncertainties in $`g`$, $`h`$ and $`f`$ were estimated by a Monte Carlo data decimation step with $`V_c`$ constrained to its best fit value.

For several reasons, it was necessary to restrict the fit to the neighborhood of $`ϵ0`$. In the first place, amplitude equation models are only rigorously valid in the limit $`ϵ1`$, although in practice they have been found to apply over a more extended range. Secondly, restricting the range of $`ϵ`$ reduced the impact of the residual, uncorrected component of the conductivity drift. In many cases, the drift effects were sufficiently large that only the data acquired with increasing $`ϵ`$ were fit.

Figure 4 shows a typical result of transforming current-voltage data to reduced Nusselt numbers $`n`$ and amplitudes $`A`$ and fitting to Eqn. 13. The transition from conduction to convection is continuous and the best fit parameter $`g>0`$, indicating that the bifurcation is supercritical. Figure 5 shows the analogous result for a film under a strong shear. For this case, the transition from conduction to convection is discontinuous and $`g<0`$; the bifurcation is subcritical. In all fits, we found $`h>0`$, and $`0<f1`$. We performed a large number of such fits and surveyed the dependence of $`g`$ on $`\alpha `$, $`𝒫`$, and $`e`$. In addition, the fits provided determinations of the critical voltage $`V_c`$ as a function of $`\alpha `$, $`𝒫`$, and $`e`$ via the bootstrap process described above. These results are discussed in Sections IV C and IV D below.

### B Tests of linear stability theory

This section compares the experimental measurements with the predictions of the linear stability theory given in Ref. .

The simplest feature of linear theory that can be compared with experiment concerns the critical mode number for zero shear, $`m_c^0`$. As mentioned in Section III, some experiments were performed in films with slight thickness nonuniformity. This permitted the qualitative visualization of the flow field and a quantitative determination of the mode number $`m`$ of the stationary pattern of vortices. Unless the bifurcation is strongly subcritical, we expect $`m=m_c^0`$ close to onset. The observed zero shear mode number was in excellent agreement with predictions of linear stability analysis. Table I summarizes the results. Rapid rotation of the vortex pattern and larger hysteresis prevented a systematic study of $`m`$ under shear. Qualitative observations confirm the general prediction that $`m_c(e>0)<m_c^0`$, i.e. that shear reduces the number of vortices.

The primary theoretical result of Ref.  is the prediction of the critical voltage $`V_c`$ required for the onset of electroconvection. $`V_c`$ is given for the general case by Eqn. 16. Even though the viscosity $`\eta `$ is not well-known, we can test the zero shear theory for various thicknesses $`s`$ and conductivities $`\sigma `$ using the following scheme. Denoting the zero shear value of $`V_c`$ by $`V_c^0`$, we write

$$(V_c^0)^2(\alpha )=\left[\frac{\sigma \eta s^2}{ϵ_{0}^{}{}_{}{}^{2}}\right]_c^0(\alpha ).$$
(18)
Using Eqn. 11, Eqn. 18 can be expressed more conveniently as

$$v^2=\frac{4\pi ^2ϵ_{0}^{}{}_{}{}^{2}\sigma (V_c^0)^2(\alpha )}{(\mathrm{ln}(1/\alpha ))^2_c^0(\alpha )}=\eta c^2.$$
(19)
Written in this way, Eqn. 19 expresses a proportionality between $`v^2`$ and $`c^2`$ in which the viscosity $`\eta `$ is the only unknown parameter. Consistency with this proportionality over a wide range of parameters serves as a test of the linear theory for $`V_c^0`$ as well as a determination of $`\eta `$. There is one caveat: this analysis is not entirely experimental but requires the theoretical value of $`_c^0(\alpha )`$. The quantity $`v`$, which we refer to as the scaled critical voltage, was found as follows. For each set of current voltage data, the fit procedure outlined in Section IV A was used to deduce a critical voltage $`V_c^0`$ and conductance $`c`$ at onset. The film thickness $`s`$ was deduced from the color of the film. Using the measured radius ratio $`\alpha `$ and Eqn. 11, $`\sigma `$ was calculated. Finally, the numerical result for $`_c^0(\alpha )`$ was used to find $`v^2`$ which was plotted against $`c^2`$. Figure 6 shows the results obtained from $`228`$ current-voltage runs at six different $`\alpha `$, numerous different $`s`$ and conductivities in the range $`5.9\times 10^8<\sigma <8.4\times 10^7\mathrm{\Omega }^1\mathrm{m}^1`$. Consequently, the range of $`𝒫`$ is very broad.

Within some scatter, Fig. 6 exhibits the predicted proportionality over several decades. Hence, the theory properly accounts for the scaling of the critical voltage with respect to the film thickness and radius ratio.

A single-parameter linear fit to $`v^2=\eta c^2`$ gives $`\eta =0.18\pm 0.03`$ kg/ms. This is a reasonable value for the viscosity; while it has not been independently measured, it is believed to be of order of $`0.1`$ kg/ms.

We now turn to the case of nonzero shear. The main prediction of the linear stability analysis is that the onset of electroconvection is suppressed by Couette shear. This suppression is both $`\alpha `$ and $`𝒫`$ dependent. Returning to Eqn. 2, the degree of suppression $`\stackrel{~}{ϵ}`$ is given by

$$\stackrel{~}{ϵ}(\alpha ,e,𝒫)=\left[\frac{_c(\alpha ,e,𝒫)}{_c^0(\alpha )}\right]1=\left(\frac{V_c(\alpha ,e)}{V_c^0(\alpha )}\right)^21.$$
(20)
We used the two equivalent expressions for $`\stackrel{~}{ϵ}`$ in Eqn. 20 to calculate the suppression theoretically and experimentally. $`V_c(\alpha ,e)`$ was found from the nonlinear fit, along with the conductance $`c`$ for that particular sweep. It was necessary to correct $`V_c^0(\alpha )`$ for the drift of $`c`$ in order to calculate $`\stackrel{~}{ϵ}`$. This was done using values of $`V_c^0(\alpha )`$ taken from zero shear sweeps performed before and after each sheared sweep. The variation of $`V_c^0(\alpha )`$ with $`c`$ was modelled with a linear fit which was used to find the drift-corrected value of $`V_c^0(\alpha )`$. In this way each experimental value of $`\stackrel{~}{ϵ}`$ could be be computed for a constant $`c`$. The uncertainty in $`\stackrel{~}{ϵ}`$ was dominated by the uncertainty

in the drift correction to $`V_c^0`$. The conductance $`c`$ for each $`\stackrel{~}{ϵ}`$ can be used to calculate its associated value of $`𝒫`$. Thus, we arrive at sets of experimental values of $`\stackrel{~}{ϵ}(\alpha ,e,𝒫)`$ for fixed $`\alpha `$ and $`e`$ and ranges of $`𝒫`$. These can be compared in detail to the theoretical results for suppression as a function of these parameters given in Ref. . Figures 7a and b show this comparison at $`\alpha =0.47`$ and $`\alpha =0.64`$ respectively. In each case the theoretical curves are calculated for the mean value of the range of $`𝒫`$ spanned by the data. For a discussion of the distinction between the local and nonlocal approximations, see Ref. .

Note that the ranges of $`e`$ for these two $`\alpha `$ are different by a factor of ten and that the suppressions are also very different. It is quite interesting that a rather mild shear ($`e3)`$ suppresses the onset of electroconvection so strongly that $`\stackrel{~}{ϵ}14`$ (i.e. $`_c=15_c^0`$). We have also studied the suppression at $`\alpha =0.33,0.56,0.60`$, and $`0.80`$. The results are similar to Fig. 7. The suppression data covers the range $`10<𝒫<130`$ and $`0e<3`$. In each case the agreement between theory and experiment is fair, with the theory slightly underestimating the suppression in most cases. The upper bound on $`e`$ could easily be extended using increased rotation rates, while smaller $`𝒫`$ would require larger, thicker films.

The degree of agreement shown in Figs. 7 is essentially independent of the value of $`\eta `$. Recall that $`\eta `$ was determined by a single parameter fit to Eqn. 19. Since the $`\eta `$ dependence in the $`e`$ scaling of both the theory and the experiment are proportional to $`1/\eta `$, any change in $`\eta `$ multiplies both by the same factor. This simply results in a rescaling of the $`e`$ axis in Figs. 7a and b, with no change in the quality of the agreement.

We close this section with some brief speculation on what might account for the scatter in Fig. 6 and for the remaining discrepancies between theory and experiment in the suppression curves.

The most likely sources of these discrepancies are the imperfect geometry of the film and its finite thickness. The 2D theoretical model assumed that the flow velocity is independent of position over the thickness of the film. This may be inexact for thicker films. Since the electrical forcing is localized near the free surfaces, it seems likely that the surfaces are preferentially driven in thick films so that the motion is not accurately 2D. The film’s edges are also imperfect because small wetting layers unavoidably exist on the circumferences of the inner and outer electrodes. These may produce electrical or velocity boundary conditions that are not exactly those assumed by the theory. In general, however, the linear stability theory works remarkably well considering its rather simple assumptions.

### C Coefficients of the Cubic Nonlinearity without Shear

In this section, we consider the experimental results in the nonlinear regime, beginning with the unsheared $`e=0`$ case. These are expressed in terms of the cubic Landau coefficient $`g`$ of the amplitude equation, Eqn. 13. Even with $`e`$ fixed at zero, we will show that this coefficient is an interesting and nontrivial function of the remaining dimensionless quantities, $`𝒫`$ and $`\alpha `$.

The radius ratio $`\alpha `$ is a fixed geometric parameter for each run, determined only by the choice

of the electrode dimensions. The $`𝒫`$ appropriate to each film was deduced after each run as described in Section IV A. Since $`𝒫`$ is proportional to the conductivity, which exhibits a slow drift, a wide range of $`𝒫`$ could be investigated. While they are independent in principle, $`\alpha `$ and $`𝒫`$ are not easy to separate experimentally. This practical constraint derives from Eqn. 12 in which $`𝒫d^1`$, where the width of the film $`d=r_or_i`$. In practice, $`d`$ is largest at small $`\alpha `$, so the runs with the smallest $`𝒫`$ occur for small $`\alpha `$.

The nonlinear coefficient $`g(\alpha ,𝒫)`$ was determined as described in Section IV A. Figures 8a and b show $`g`$ as a function of $`𝒫`$ for two different $`\alpha `$. The scatter that is manifest in these plots exceeds the statistical uncertainty of the fit. As discussed in Section IV A, the scatter originates from systematic effects due to the non-ideal features of the experiment. Nevertheless, the gross features of the $`𝒫`$ dependence of $`g`$ can still be extracted.

At $`\alpha =0.33`$, the measurements explored the range $`2<𝒫<8`$. These were the smallest $`𝒫`$ reached in the experiment. It is clear from Fig. 8a that $`g`$ increases with $`𝒫`$, and passes through zero for small enough $`𝒫`$. Thus, for $`\alpha =0.33`$, we find that the bifurcation from the conduction state to the $`m=4`$ vortex state is subcritical ($`g<0`$) for $`𝒫\stackrel{<}{}5`$ and supercritical ($`g>0`$) for $`𝒫\stackrel{>}{}5`$. Near $`𝒫\stackrel{}{=}5`$ we pass a tricritical point. It is interesting to observe that in each case the bifurcation involves the same two symmetry states, yet its subcriticality depends on $`𝒫`$.

For all the other, larger $`\alpha `$ investigated, $`g`$ was found to be independent of $`𝒫`$. This is illustrated in Fig. 8b for $`\alpha =0.64`$ and was also true for $`\alpha =0.47,0.56,0.60,`$ and $`0.80`$. However, for all of these cases, $`𝒫`$ was never less than $`10`$ and can be as large as several hundred. Thus, it is unclear whether the subcriticality at $`\alpha =0.33`$ is due to the smallness of $`𝒫`$ or the smallness of $`\alpha `$, or some combination. It would be interesting to determine whether $`g`$ also becomes negative for large $`\alpha `$ as $`𝒫`$ decreases. The regime of large $`\alpha `$ and small $`𝒫`$, while accessible in principle, would require significantly larger electrodes or much thicker films. All electroconvection experiments on freely suspended films in the rectangular geometry have reported supercritical bifurcations. However, these experiments were also at large $`𝒫`$ and so the possibility of a subcritical bifurcation in rectangular films cannot be excluded.

The case of large $`𝒫`$ is easier to understand. In the governing electrohydrodynamic equations (see Ref. ), the parameter analogous to the Prandtl number only appears as the inverse, $`𝒫^1`$, multiplying certain nonlinear terms. Hence, it is perhaps not surprizing that $`g`$ becomes independent of $`𝒫`$ for $`𝒫\stackrel{>}{}10`$, where these terms become negligible.

In order to examine how $`g`$ depends on $`\alpha `$, we removed the $`𝒫`$ dependence by averaging over data at various $`𝒫`$. For the five largest $`\alpha `$, $`g`$ is independent of $`𝒫`$ over broad ranges of $`𝒫`$. A weighted average of $`g`$ was obtained for these cases. For $`\alpha =0.33`$, where some $`𝒫`$ dependence was evident, only $`g`$ values in the narrow range $`2.1<𝒫<4.4`$ were averaged.

In this averaging the systematic scatter in $`g`$ was treated as random. It is thus likely that the true uncertainty in the average values of $`g`$ will be much larger than the standard deviation of the mean. The ultimate test of this procedure lies in the comparison of the averaged values of $`g`$ with theoretical predictions. As we describe below, the limited comparison we are able to make at present is not unfavourable. The line in Fig. 8b shows the average value of $`g`$ for the plotted data. All of the results for $`g(\alpha )`$ are tabulated in Table II and plotted in Fig. 9.

It is clear from Fig. 9 that, overall, $`g(\alpha )`$ increases with $`\alpha `$. As discussed above, the larger $`\alpha `$ data are also associated with larger $`𝒫`$. At present, no theoretical predicion is available for $`g`$ at arbitrary $`\alpha `$ and $`𝒫`$. However, we expect $`g`$ to approach a limiting value quite rapidly as $`\alpha 1`$. This limit corresponds to an unbounded lateral geometry in which the film is a strip of fluid suspended at its long parallel edges by two semi-infinite plate electrodes. In this limit the discrete azimuthal mode numbers $`m`$ are replaced by a continuous wavenumber $`k=m/\overline{r}`$, where $`\overline{r}=(r_i+r_o)/2`$ is the mean radius. From linear theory, we found that the critical parameters $`_c^0`$ and $`k_c^0=m_c^0/\overline{r}`$ at $`\alpha =0.80`$ are already very close to the limiting values for $`\alpha 1`$.

This limiting behaviour is reasonable given the large dimensionless circumference for large $`\alpha `$. The relevant aspect ratio $`\mathrm{\Lambda }`$ is the circumference at the mean film radius, $`\pi (r_i+r_o)`$, divided by the width of the film $`d=r_or_i`$, so that $`\mathrm{\Lambda }=\pi (1+\alpha )/(1\alpha )`$. At $`\alpha =0.80`$, $`\mathrm{\Lambda }28`$. Almost all the experiments performed in the rectangular geometry had $`\mathrm{\Lambda }\stackrel{<}{}10`$  and were well-modelled by the theory for an unbounded strip for which $`\mathrm{\Lambda }=\mathrm{}`$. It is thus resonable to expect $`g`$ at $`\alpha =0.80`$ to be close to its limiting value for $`\alpha =1`$. Similarly, as the large $`\alpha `$ data also involve very large values of $`𝒫`$, we may also employ the theory in the $`𝒫\mathrm{}`$ limit.

Weakly nonlinear analysis of unsheared electroconvection in the ‘plate’ electrode geometry was presented in Ref.  for the case $`𝒫=\mathrm{}`$. The result of that analysis, $`g=2.842`$ is shown in Fig. 8. It is in good agreement with a reasonable extrapolation of the data to the $`\alpha 1`$ limit.

When $`g<0`$, the onset of convection becomes hysteretic and the quintic term with coefficient $`h`$ in Eqn. 13 becomes significant. Since $`0<f1`$, the width of the hysteresis loop $`\delta ϵ`$ in $`ϵ`$ is given by $`\delta ϵ=g^2/4h`$. When $`g0`$, $`\delta ϵ=0`$. The hysteresis width with zero applied shear is plotted in Fig. 10 for $`\alpha =0.33`$ and various $`𝒫`$. Note that the hysteresis vanishes for $`𝒫\stackrel{>}{}5`$ and even when non-zero, is always rather small, $`\delta ϵ0.02`$. This is in contrast with the sheared case, discussed in the next Section, for which $`\delta ϵ0.1`$.

### D Coefficients of the cubic nonlinearity with shear

In the presence of shear, the coefficient $`g`$ of the cubic nonlinearity is strongly dependent on the Reynolds number $`e`$ of the Couette flow imposed in the base state. We find that the shear drives $`g`$ negative, producing a hysteretic Hopf bifurcation to convection in the form of travelling vortices. In this Section, we describe the $`e`$ dependence of $`g`$, which, as in the previous Section, is also a function of the radius ratio $`\alpha `$ and the dimensionless ratio $`𝒫`$.

Figure 11 shows a representative example of the $`e`$ dependence of $`g`$, in this case for $`\alpha =0.47`$. At $`e=0`$, the value of $`g`$ is found, as described in Section IV C, by averaging over a range of $`𝒫`$, which was always $`>10`$. In the case shown in Fig. 11, $`g(e=0)=1.64\pm 0.06`$ with a mean $`𝒫`$ of $`16.3`$. The plotted $`g`$ data for $`e>0`$ had $`13.3<𝒫<21.7`$.

Our main result is that $`e`$ plays the part of a second control parameter which allows the primary bifurcation

to electroconvection to be tuned between super- and subcritical. For $`\alpha =0.47`$, the bifurcation is supercritical at $`e=0`$ and weakens as the Reynolds number increases. The bifurcation becomes tricritical at $`e0.2`$ and is subcritical thereafter. The subcriticality at first deepens with increasing Reynolds number, until at $`e0.85`$ a minimum value of $`g\stackrel{}{=}3.7`$ is reached. For $`e>0.85`$, the bifurcation remains subcritical but $`g`$ becomes an increasing function of the Reynolds number. For the range of Reynolds numbers investigated the bifurcation does not become supercritical again.

There is some scatter in the data, nonetheless, the overall trends are clear. The systematic deviations are comparable to those in Figs. 8a and b. The results were qualitatively similar for $`\alpha =0.56,0.60,0.64,`$ and $`0.80`$. The tricritical Reynolds number at which $`g=0`$ will be denoted $`e_T`$. The value of $`e_T`$ varys for different $`\alpha `$ and $`𝒫`$, but the general features of the $`g`$ vs. $`e`$ curve are preserved over a wide range of parameters. Table III lists the values of $`e_T`$ under various conditions.

The variation of $`e_T`$ is probably due to the combination of changing both $`\alpha `$ and $`𝒫`$, rather than to the variation either parameter separately. These are difficult to experimentally disentangle and the parameter space is large. For zero shear, $`g`$ was found to be independent of $`𝒫`$ for $`𝒫\stackrel{>}{}10`$, but under shear there is insufficient data to draw many conclusions about the $`𝒫`$ dependence of $`g`$. If linear theory is any guide, we expect there will in general be a greater $`𝒫`$ dependence in the sheared case. It was established in Ref.  that the linear theory result for $`_c`$ was independent of $`𝒫`$ for $`e=0`$,

but weakly dependent on $`𝒫`$ for $`e0`$. It appears from the tabulated results that as $`𝒫`$ and $`\alpha `$ decrease, $`e_T`$ increases.

The minimum value of $`g`$ as a function of $`e`$ also vary for different $`\alpha `$ and $`𝒫`$. We will refer to the coordinates of this special point as $`g_{min}`$ and $`e_{min}`$. Table IV lists the minimum values observed. A look at Fig. 11 shows that the minimum value is only known up to the density of data and to within the scatter in the whole plot. The minimum is a local one and is observed within ranges of $`e`$ which had different maxima for different $`\alpha `$.

Since the parameter space for the data is defined by several parameters it becomes difficult to make meaningful comparisons of the experimental results when more than one parameter changes. However, one fortunate comparison can be gleaned from Table IV. At $`\alpha =0.47`$, the film had $`𝒫=15.3`$ while at $`\alpha =0.80`$, $`𝒫=12.0`$. Since the radius ratios are very different and the $`𝒫`$ are not, it is not unreasonable to directly compare the values of $`g_{min}`$ and $`e_{min}`$ for these cases. It is evident that the bifurcation is much more strongly subcritical at $`\alpha =0.80`$ than at $`\alpha =0.47`$, and the values of $`e_{min}`$ are also very different.

The coefficient of the quintic nonlinearity $`h`$ as function of $`e`$ is studied in Ref. . Figure 12 shows the the hysteresis width $`\delta ϵ=g^2/4h`$ as a function of $`e`$ for $`\alpha =0.47`$. The hysteresis when $`e>0`$ is much larger than that described in the previous Section which was found for $`e=0`$ at small $`\alpha `$ and $`𝒫`$.

The Nusselt number of the convection is independent of the Hopf frequency (i.e. the travelling rate of the vortex pattern). It would be interesting to measure this frequency, which is governed by the imaginary part of the complex Ginzberg-Landau equation, Eqn. 10, and depends on $`ϵ`$, $`\alpha `$, $`e`$ and $`𝒫`$. This would however require visualization of the flow pattern or some spatially localized probe of the amplitude.

## V Discussion

Because of its rather simple symmetry and forcing scheme, annular electroconvection under applied shear has many similarities to other systems and models. In this section, we discuss these similar systems with a view to putting annular electroconvection into a general context.

Electroconvection is obviously analogous to buoyancy-driven thermal convection. In fact, it can be shown that the linear stability problem for surface charge driven electroconvection in thin films reduces to the 2D Rayleigh-Bénard problem if the nonlocal coupling between fields and charges is neglected. Detailed calculations show that this “local” approximation is an accurate one, even when applied in the weakly nonlinear regime.

There have been some theoretical studies of 3D Rayleigh-Bénard convection (RBC) in the presence of a plane Couette shear flow. The theory assumed the usual theoretical geometry of a fluid layer confined between infinite, perfectly conducting horizontal planes. Our annular system approaches the plane one in the $`\alpha 1`$ limit. (The $`\alpha <1`$ annular geometry has finite extent and curvature; the implications of these are discussed below.) Linear stability analysis of a plane Couette base state to RBC reveals stability differences between transverse roll disturbances (with axes perpendicular to the shear flow), and longitudinal roll disturbances (with axes parallel to the shear flow). Longitudinal-roll disturbances have identical stability properties to unsheared RBC, and are always more unstable than the transverse-roll disturbances. In fact, longitudinal-roll disturbances have stability properties that are independent of any uni-directional shear flow along the axis of the rolls. In our 2D system, these more unstable longitudinal roll modes do not exist and the vortices we see correspond to truly 2D transverse rolls. Thus, we observe the analog of a state that would normally be preempted by longitudinal rolls if the geometry were not constrained.

According to linear theory, transverse-roll disturbances in unbounded, sheared RBC, like our vortices, exhibit suppression, or added stability due to the shear, under plane Poiseuille or plane Couette or any mixture of these two flows. The onset Rayleigh number for transverse rolls is a monotonically increasing function of the shear Reynolds number, similar to what we found for electroconvection. Furthermore, the critical wavenumber of the most unstable transverse disturbance was found to be a monotonically decreasing function of the shear Reynolds number, which is analogous to the reduction of $`m_c`$ by shear that we observed. Transverse rolls also travel under unidirectional plane Couette shear, again analogous to our travelling vortex state.

Whereas unbounded RBC with plane Couette shear cannot be realized experimentally, RBC has been studied experimentally and theoretically in narrow slots with open through-flows. Quasi-2D transverse rolls can be stabilized by wall effects in slots. The through-flow consists of a weak Poiseuille flow with a very small Reynolds number. Its effects on RBC are well understood. In brief, the onset of convection is again suppressed, but the first instability is convective (i.e. it grows only downstream of a localized perturbation), rather than absolute. The resulting convection pattern drifts in the direction of the through flow. It is interesting that the clear distinction between convective and absolute instability is blurred in annular electroconvection with shear, in which the ‘through’ flow loops back on itself. The annular geometry is naturally closed. The linear stability analysis of Ref.  only treated absolute instablity. In principle, the system might still be considered to be convectively unstable if localized perturbations grew only as they travelled azimuthally.

Taylor Vortex flow (TVF) is an extensively studied pattern-forming instability with geometric similarities to annular electroconvection. However, the instability leading to TVF depends crucially on axial disturbances to 3D Couette flow. Purely 2D circular Couette flow is, in fact, linearly stable. TVF is the result of an instability of a 3D shear flow, while what we have studied here is the effect of a shear flow on the electroconvective instability.

Agrait and Castellanos theoretically studied the effect of a 3D Couette shear on an electrohydrodynamic instability in TVF geometry. Their electrohydrodynamic system consisted of a nearly insulating fluid confined between metallic electrodes. Charge injection, a process by which charge carriers are created at the electrodes, occurs when strong electric fields are applied. The interaction of this volume charge density with the applied electric field leads to electroconvection instabilities. Agrait and Castellanos considered electroconvection due to a radial field with charge injection on either cylinder. Both cylinders were permitted to rotate to produce a general Couette shear. They found that shearing enhanced the instability, leading to a 3D flow that resembles TVF. This is in direct contrast to the stabilizing effects we observed in 2D.

Although shear and rigid rotation are conceptually quite different, it is interesting to compare their effects on instabilities. Again, we find crucial distinctions between 2D and 3D systems.

The added stability in sheared annular electroconvection is a consequence of the shear and not of rotation. Under rigid rotation, where the inner and outer electrodes are co-rotating, one can transform to rotating co-ordinates in which the electrodes are stationary. This transformation introduces a Coriolis term $`2\mathrm{\Omega }\widehat{𝐳}\times \stackrel{}{𝐯}=2\mathrm{\Omega }\psi `$ in the Navier-Stokes equation which may be absorbed into the pressure gradient term $`P`$ and eliminated. Thus, in a purely 2D system, rigid rotation and the non-rotating, unsheared case have identical stability. It also follows, since the transformation is general and the unsheared bifurcation is stationary, that the resulting nonlinear vortex pattern above onset must be stationary in the co-rotating frame.

This lack of dependence on rigid rotation may be contrasted with a large class of 3D and quasi-2D rotating Rayleigh-Bénard systems, where rotation produces added stability but the absence of strictly 2D flow results in a time-dependent (precessing) convection pattern in the co-rotating frame. The precession seen in these systems is analogous to the travelling of the sheared patterns that we observed only to the extent that rotation and shear break the same symmetry, which allows travelling patterns. Symmetry aside, the physical origins of added stability and precession are fundamentally different between the 2D sheared and 3D rotating cases.

In general, 3D rotating convecting systems may also support strictly 2D ‘Taylor column’ solutions which do not precess in the rotating frame and whose onset occurs at the same critical Rayleigh number as in the absence of rotation. The system most analogous to ours is the interesting but experimentally unrealizable situation of 2D RBC in a rotating annular geometry with purely radial gravity and heating. Theoretical studies found columnar solutions which are very closely analous to our vortices. One might hope to approach this limit in experiments using radially temperature gradients imposed between rapidly rotating concentric cylinders. Purely columnar solutions have not observed in any rotating RBC experiment because the boundary conditions at the top and bottom of the cylinder must be stress free, a requirement that cannot be attained in terrestrial RBC experiments. In contrast, two-dimensionality, stress free surface conditions and radial driving forces all arise naturally in the electroconvection of an annular suspended film. Thus, the vortices that occur in annular electroconvection without shear are accurate 2D analogs of ‘Taylor columns’.

The similarity between our system and many other better-studied systems and models suggests that many linear and nonlinear techniques developed for other problems can be fruitfully brought to bear on annular electroconvection. In addition, the system may allow the study of bifurcation scenarios that are not experimentally realizable in other similar systems.

.

## VI Conclusion

In this paper, we have reported a wide ranging experimental study of the primary bifurcation to electroconvection in sheared two-dimensional annular films. Our principal experimental probe was the excess current carried by the film due to convection, which can be directly related to the amplitude of the convective flow. In all, we examined annuli with six different radius ratios $`\alpha `$, with $`0.33\alpha 0.80`$. For these, the Reynolds number of the applied Couette shear varied between $`0e<3`$. The explored range for the dimensionless parameter $`𝒫`$, which varies with film thickness and conductivity, was $`1<𝒫<150`$.

We compared this data, in the first instance, with the predictions of linear theory. The data for the critical voltage $`V_c^0`$ for the onset of convection without shear could be combined for all six radius ratios and was shown to obey the scaling predicted by linear theory. This data could then be used in a slightly model-dependent way to fix the only remaining unknown material parameter, the fluid viscosity. We could then compare the experimentally measured suppression of the onset by shear to linear theory in an essentially pararmeter-free way. The agreement was satisfactory for a wide range of $`\alpha `$, $`𝒫`$, and $`e`$. Using a simple visualization scheme, we also confirmed that the azimuthal mode number $`m`$ near onset was close to the most unstable critical value $`m_c`$ predicted by linear theory.

Nonlinear fits to current-voltage data above the onset of convection were then used to infer the real part of the amplitude of convection. Under shear, the amplitude can be shown to obey a complex Ginzburg-Landau equation; its real part was fit to the time-independent steady state version of this equation in order to extract the coefficient $`g`$ of the cubic nonlinearity. When shear is absent, the primary bifurcation is a pitchfork while with shear, the bifurcation is a pitchfork Hopf. We determined the functional dependence of $`g`$ as $`\alpha `$, $`e`$ and $`𝒫`$ were varied. When shear is absent, it was found that for $`\alpha =0.47,0.56,0.60,0.64,0.80`$ and $`𝒫\stackrel{>}{}13`$ that $`g`$ was independent of $`𝒫`$. In addition, in this regime $`g>0`$, hence the bifurcation was supercritical. For $`\alpha =0.33`$, $`g`$ was found to be an increasing function of $`𝒫`$ for $`2<𝒫<8`$. More importantly, it was found that the bifurcation was subcritical($`g<0`$) for $`𝒫\stackrel{<}{}5`$ and supercritical($`g>0`$) for $`𝒫\stackrel{>}{}5`$. In overall trend, we found that $`g`$ is an increasing function of $`\alpha `$ for zero shear. The largest $`\alpha =0.80`$ was sufficiently close to the limiting case of $`\alpha 1`$ that we could extrapolate the measured values of $`g`$ for comparison to the value of $`g`$ calculated from weakly nonlinear theory for the corresponding ‘plate’ electrode geometry. This quantitative comparison gave reasonable agreement.

When shear was applied, it had a strong effect on the subcriticality of the primary bifurcation. Measurements of $`g`$ as a function of the shear Reynolds number $`e`$ revealed that for $`\alpha =0.47,0.56,0.60,0.64,0.80`$ there existed a tricritical Reynolds number $`e_T`$ below which $`g>0`$ and the onset is a supercritical Hopf bifurcation. For $`e>e_T`$, $`g<0`$ and the onset became a subcritical Hopf bifurcation with the degree of subcriticality a nontrivial function of $`e`$. Hence, the Reynolds number forms an interesting second control parameter in this system which can be used to vary the nature of the primary bifurcation.

The data presented in this paper represents only a sparse sampling of the full 3D parameter space of $`g(\alpha ,𝒫,e)`$. Within this space, there presumably exist continuous 2D loci on which $`g=0`$ and the primary bifurcation is tricritical. On other loci, the hysteresis $`\delta ϵ`$ is a local maximum. As a function of $`\alpha `$, linear theory shows that there exist special values of $`\alpha `$ at which two azimuthal mode numbers $`m`$ and $`m+1`$ are simultaneously linearly unstable. All of the $`g`$ results we have described pertain to the primary bifurcation; above onset we also observed numerous secondary bifurcations which take the form of hysteretic jumps in $`m`$. The location of these secondary bifurcations is a function of $`e`$. These considerations suggest that more-or-less discontinuous jumps in the behavior of $`g`$ may occur at parameter values across which the azimuthal mode number $`m`$ of the nonlinear pattern changes discontinuously. The rich phenomenology of even the primary bifurcation presents a significant challenge to the weakly nonlinear theory of this instability.

## acknowledgements

This work was supported by the Natural Science and Engineering Research Council of Canada and the U.S. Department of Energy under contract W-7405-ENG-36.


format: text

text:
# Do gravitational waves carry energy-momentum and angular momentum?

## I Introduction

In General Relativity (GR) the gravitational field $`\mathrm{\Gamma }_{kl}^i`$ does not possess any energy-momentum tensor but it only possesses the so-called “energy-momentum pseudotensors”. It is a consequence of the Einstein Equivalence Principle (EEP). So, many authors  doubt in reality of the energy-momentum and angular momentum transfer by gravitational waves. As an argumentation some of these authors use the fact that for the majority exact solutions of the vacuum Einstein equations which, as we think, can represent gravitational waves, pseudotensors globally vanish in some coordinates. In consequence, they give in these coordinates “no gravitational energy and no gravitational energy flux”. Other of these authors  try to use vanishing of the components $`{}_{g}{}^{}t_{}^{ok}`$ (or $`{}_{g}{}^{}t_{0}^{k}`$) of a gravitational energy-momentum pseudotensor $`{}_{g}{}^{}t_{}^{ik}`$ (or $`{}_{g}{}^{}t_{i}^{k}`$) as a coordinate condition coupled to the Einstein equations and then also get (in special coordinates), as they assert, “global vanishing of the pure gravitational energy and gravitational energy flux”.

However, such conclusions are physically uncorrect. Firstly, these authors neglect an important role of the four-velocity $`\stackrel{}{v},\stackrel{}{v}\stackrel{}{v}=1`$ of an observer O, which is measuring gravitational (or other) field, in definition of the energy density $`ϵ`$ and energy flux $`P^i`$ of the field. Namely, the correct definitions of the energy density and its flux for such observer are the following

$$ϵ=T^{ik}v_iv_k=T_{ik}v^iv^k,$$
(1)
$$P^i=\left(\delta _k^iv^iv_k\right)T^{kl}v_l;$$
(2)
not

$$ϵ=T^{00},P^\alpha =T^{0\alpha }.$$
(3)
$`i,k,l,\mathrm{},=0,1,2,3;\alpha ,\beta ,\gamma ,\delta ,\mathrm{},=1,2,3`$. $`T^{ik}`$ mean here the components of an energy-momentum tensor.

$`T^{00}`$ and $`T^{0\alpha }`$ give energy density $`ϵ`$ and energy flux $`P^i`$ for the observer O iff in the used global coordinates we have for the observer $`v^i=\delta _0^i,v_k=\delta _k^0`$.***This means physically that the observer O is at rest with respect to the used coordinates.

The above facts are also true for a gravitational energy-momentum pseudotensor $`{}_{g}{}^{}t_{}^{ik}`$ (or $`{}_{g}{}^{}t_{i}^{k}`$).

In consequence, even if $`{}_{g}{}^{}t_{}^{0k}=0`$ (or $`{}_{g}{}^{}t_{0}^{k}=0`$) globally in the used coordinates then, as one can easily see, not for all observers $`ϵ=P^i=0`$. All depends on the four-velocity $`\stackrel{}{v}`$ of the observer O.

Moreover, one coordinate condition of the kind $`{}_{g}{}^{}t_{}^{0k}=0`$ (or $`{}_{g}{}^{}t_{0}^{k}=0`$) does not suffice for the all variety of the gravitational energy-momentum pseudotensors. In fact, we need here an own coordinate condition for every gravitational energy-momentum pseudotensor.

Secondly, and it is the most important fault, these authors forget that gravitational energy-momentum (and gravitationam angular momentum) pseudotensors, as being functions of th Levi-Civita connection coefficients Physically, this connection plays the role of the total gravitational strenghts. describe energy-momentum of the total gravitational field which is a combination of the real gravitational field (for which $`R_{iklm}0`$) and inertial forces field (for which $`R_{iklm}=0`$). The inertial forces field is generated by the used coordinates. This is also a consequence of the EEP.

So, if in some coordinates $`{}_{g}{}^{}t_{}^{ik}`$ or ($`{}_{g}{}^{}t_{i}^{k}`$) (or only $`{}_{g}{}^{}t_{}^{ok}`$ or $`{}_{g}{}^{}t_{0}^{k}`$) globally vanish, this does not mean that in these coordinates we have no pure gravitational energy and no pure gravitational energy flux.Even if we confine to the observers which are at rest in the used global coordinates and for which $`{}_{g}{}^{}t_{}^{00}`$ (or$`{}_{g}{}^{}t_{0}^{0}`$) is “energy density” and $`{}_{g}{}^{}t_{}^{0\alpha }`$ (or $`{}_{g}{}^{}t_{0}^{\alpha }`$) give the “energy flux”. Simply, it means only that in such coordinates energy and energy flux of the real gravitational field cancel with the energy and energy flux of the inertial forces field.§§§Energy-momentum of the inertial forces field gives contribution to energy-momentum pseudotensors. Energy and energy-momentum flux of the real gravitational field which has $`R_{iklm}0`$ always exist and are non-null. In order to show this very important fact one can use the canonical superenergy and canonical angular supermomentum tensorsOr other superenergy and angular supermomentum tensors. for gravitational field \[10-14\].

The canonical superenergy tensor $`{}_{g}{}^{}S_{i}^{k}`$ and canonical angular supermomentum tensor $`{}_{g}{}^{}S_{}^{ikl}=()_gS^{kil}`$ are constructed in such a way that they extract covariant information about real gravitational field which is hidden in canonical energy-momentum and canonical angular momentum pseudotensors. Namely, these superenergy and angular supermomentum tensors are obtained from suitable gravitational pseudotensors by some kind of averaging and they are functions of the curvature tensor and its covariant derivatives only. So, they really describe only real gravitational field, which has $`R_{klm}^i0`$.

The paper is organized as follows. In Section II we remember canonical superenergy and supermomentum tensors for gravitational field in GR. In Section III we give a summary about applications of the canonical superenergy and supermomentum tensors for gravitational field to analysis of the gravitational waves and in Section IV we give some concluding remarks.

## II Constructive definition of the canonical superenergy tensor and canonical angular supermomentum tensor

We remember here the constructive, general definition of the superenergy tensor $`S_a^b(P)`$ applicable to gravitational field and to matter fields.

In the normal coordinates NC(P) \[see, e.g., 20,21,22\], we define

$$S_{(a)}^{(b)}(P):=()\underset{\mathrm{\Omega }P}{lim}\frac{\underset{\mathrm{\Omega }}{}\left[T_{(a)}^{(b)}(y)T_{(a)}^{(b)}(P)\right]𝑑\mathrm{\Omega }}{1/2\underset{\mathrm{\Omega }}{}\sigma (P;y)𝑑\mathrm{\Omega }},$$
(4)
where

$$T_{(a)}^{(b)}(y):=T_i^k(y)e_{(a)}^i(y)e_k^{(b)}(y),$$
(5)
$$T_{(a)}^{(b)}(P):=T_i^k(P)e_{(a)}^ie_k^{(b)}(P)=T_a^b(P)$$
(6)
are the so-called physical or tetrad components of the pseudotensor (or tensor) field $`T_i^k(y)`$ which describes energy-momentum. $`\{y^i\}`$ are the normal coordinates. $`e_{(a)}^i(y),e_k^{(b)}(y)`$ mean here an orthonormal tetrad $`e_{(a)}^i(P)=\delta _a^i`$ and its dual $`e_k^{(a)}(P)=\delta _k^a`$ parallelly propagated along geodesics throught $`P`$ ($`P`$ = origin of the NC(P)).

We have

$$e_{(a)}^i(y)e_i^{(b)}(y)=\delta _a^b.$$
(7)
$`\mathrm{\Omega }`$ is a sufficiently small ball

$$y^{0^2}+y^{1^2}+y^{2^2}+y^{3^2}R^2$$
(8)
for an auxilliary positive-definite metric $`h^{ik}:=2v^iv^kg^{ik}`$, i.e., $`\mathrm{\Omega }`$ can be given as

$$h_{ik}y^iy^kR^2.$$
(9)
$`v^i`$ denotes here the components of the four-velocity $`\stackrel{}{v}:v^lv_l=1`$ of an observer O which is at rest at the beginning P of the used normal coordinates NC(P).

The ball $`\mathrm{\Omega }`$ surrounds P: P is centre of the ball.

$$\sigma (P;y)\dot{=}\frac{1}{2}\left(y^{0^2}y^{1^2}y^{2^2}y^{3^2}\right)$$
(10)
is the two-points world function which has been introduced by J.L. Synge years ago . The symbol $`\dot{=}`$ means that an equation is only valid for special coordinates. The latter can be covariantly defined by the eikonal-like equation

$$g^{ik}_i\sigma _k\sigma =2\sigma $$
(11)
together with $`\sigma (P,P)=0,_i\sigma (P,P)=0`$.

The $`\mathrm{\Omega }`$ can be also given by the inequality

$$h^{ik}_i\sigma _k\sigma R^2.$$
(12)
Since the tetrad components and normal components are equal at the point P, we will write the components of any quantity attached to the point P without tetrad brackets , e.g., we will write $`S_a^b(P)`$ instead of $`S_{(a)}^{(b)}(P)`$ and so on.

If $`T_i^k`$ are the components of a symmetric energy-momentum tensor of matter then we will get from (4)

$${}_{m}{}^{}S_{a}^{b}(P)=\delta ^{mn}_{(m}{}_{n}{}^{})\widehat{T}_a^b.$$
(13)
Hat over a quantity denotes its value at P.

By using of the four-velocity $`\stackrel{}{v}`$ of a fictitious observer O being at rest at the beginning P of the used NC(P) and local metric $`\widehat{g}^{ab}\widehat{=}\eta ^{ab}`$ one can write this covariantly as

$${}_{m}{}^{}S_{a}^{b}(P;v^l)=\left(2\widehat{v}^l\widehat{v}^m\widehat{g}^{lm}\right)_{(l}_{m)}\widehat{T}_a^b.$$
(14)
The last formula gives us the canonical superenergy tensor for matter.

For the gravitational field $`\mathrm{\Gamma }_{kl}^i=\{_{kl}^i\}`$, after substituting into (4) $`T_i^k=_Et_i^k`$, where

$`{}_{E}{}^{}t_{i}^{k}`$ $`=`$ $`{\displaystyle \frac{c^4}{16\pi G}}\{\delta _i^kg^{ms}(\mathrm{\Gamma }_{mr}^l\mathrm{\Gamma }_{sl}^r\mathrm{\Gamma }_{ms}^r\mathrm{\Gamma }_{rl}^l)`$ (15)
$`+`$ $`g_{,i}^{ms}[\mathrm{\Gamma }_{ms}^k{\displaystyle \frac{1}{2}}(\mathrm{\Gamma }_{tp}^kg^{tp}\mathrm{\Gamma }_{tl}^lg^{kt})g_{ms}`$ (16)
$``$ $`{\displaystyle \frac{1}{2}}(\delta _s^k\mathrm{\Gamma }_{ml}^l+\delta _m^k\mathrm{\Gamma }_{sl}^l)]\}`$ (17)
is the canonical Einstein energy-momentum pseudotensor for the gravitational field, we will obtain

$${}_{g}{}^{}S_{a}^{b}(P;v^l)=(2\widehat{v}^l\widehat{v}^m\widehat{g}^{lm}){}_{E}{}^{}\widehat{T}_{a}^{b}{}_{lm}{}^{}.$$
(18)
$`{}_{E}{}^{}T_{a}^{b}_{lm}`$ $`=`$ $`{\displaystyle \frac{2\alpha }{9}}[B_{alm}^b+P_{alm}^b`$ (19)
$``$ $`{\displaystyle \frac{1}{2}}\delta _a^bR_m^{ijk}\left(R_{ijkl}+R_{ikjl}\right)+2\delta _a^b\beta ^2E_{(l|g}E_{|m)}^g`$ (20)
$``$ $`3\beta ^2E_{a(l|}E_{|m)}^b+2\beta R_{(ag)(l|}^bE_{|m)}^g].`$ (21)
In the last formula

$$\alpha =\frac{c^4}{16\pi G}=\frac{1}{2\beta }$$
(22)
and

$$E_i^k:=T_i^k\frac{1}{2}\delta _i^kT$$
(23)
is the modified energy-momentum of matter. On the other hand

$$B_{alm}^b:=2R_{(l|}^{bik}R_{aik|m)}\frac{1}{2}\delta _a^bR_l^{ijk}R_{ijkm}$$
(24)
are the components of the Bel-Robinson tensor, while the tensor

$$P_{alm}^b:=2R_{(l|}^{bik}R_{aki|m)}\frac{1}{2}\delta _a^bR_l^{ijk}R_{ikjm}$$
(25)
is closely related to the latter.

The tensor $`{}_{g}{}^{}S_{a}^{b}(P;v^l)`$ is the canonical superenergy tensor for gravitational field.We must emphasize that the canonical superenergy tensor $`{}_{g}{}^{}S_{a}^{b}(P;v^L)`$ is originally a function of the curvature components $`R_{iklm}`$ and tensor Ricci components $`R_{ik}`$. We have eliminated the all Ricci components by use of the Einstein equations $`R_{ik}=\beta E_{ik}`$.

In vacuum the tensor $`{}_{g}{}^{}S_{a}^{b}(P;v^l)`$ takes the simpler form

$`{}_{g}{}^{}S_{a}^{b}(P;v^l)`$ $`=`$ $`{\displaystyle \frac{8\alpha }{9}}(2\widehat{v}^l\widehat{v}^m\widehat{g}^{lm})[\widehat{R}_{(l|}^{b(ik)}\widehat{R}_{aik|m)}`$ (26)
$``$ $`{\displaystyle \frac{1}{2}}\delta _a^b\widehat{R}_{(l|}^{i(kp)}\widehat{R}_{ikp|m)}]`$ (27)
and the quadratic form $`{}_{g}{}^{}S_{ab}^{}(P;v^l)v^av^b`$, where $`v^av_a=1`$, is positive-definite.

The canonical angular supermomentum tensor we define in analogy to the definition of the canonical superenergy tensor. Namely, we define in the normal coordinates NC(P)

$$S^{(a)(b)(c)}(P)=S^{abc}(P):=()\underset{\mathrm{\Omega }P}{lim}\frac{\underset{\mathrm{\Omega }}{}\left[M^{(a)(b)(c)}(y)M^{(a)(b)(c)}(P)\right]𝑑\mathrm{\Omega }}{1/2\underset{\mathrm{\Omega }}{}\sigma (P;y)𝑑\mathrm{\Omega }},$$
(28)
where, as formerly,

$$M^{(a)(b)(c)}(y):=M^{ikl}(y)e_i^{(a)}(y)e_k^{(b)}(y)e_l^{(c)}(y),$$
(29)
$`M^{(a)(b)(c)}(P)`$ $`:=`$ $`M^{ikl}(P)e_i^{(a)}(P)e_k^{(b)}(P)e_l^{(c)}(P)`$ (30)
$`=`$ $`M^{ikl}(P)\delta _i^a\delta _k^b\delta _l^c=M^{abc}(P)`$ (31)
are physical (or tetrad) components of the field $`M^{ikl}(y)=()M^{kil}(y)`$ describing angular momentum.

For matter, as $`M^{ikl}(y)=()M^{kil}(y)`$, we take

$$M^{ikl}(y)=\sqrt{|g|}\left(y^iT^{kl}y^kT^{il}\right),$$
(32)
where $`T^{ik}=T^{ki}`$ are the components of the symmetric energy-momentum tensor of matter\**\**\**This tensor is the source in the Einstein equations. and $`\left\{y^i\right\}`$ denote normal coordinates.

The expression (26) gives us the total angular momentum densities, orbital and spin, because the dynamical tensor $`T^{ik}=T^{ki}`$ is obtained from canonical one by means of the Belinfante symmetrization procedure (and, therefore, includes material spin).

For gravitational field, as $`{}_{g}{}^{}M_{}^{ikl}(y)`$, we favorize and take the expression most closely related to the Einstein canonical energy-momentum complex (See Appendix). We will call this expression canonical also. Namely, as $`{}_{g}{}^{}M_{}^{ikl}(y)`$, we take the expression given by Bergmann and Thomson

$${}_{g}{}^{}M_{}^{ikl}(y)=_FU^{i[kl]}(y)_FU^{k[il]}(y)+\sqrt{|g|}\left(y_{BT}^it^{kl}y_{BT}^kt^{il}\right).$$
(33)
In the last formula

$${}_{F}{}^{}U_{}^{i[kl]}:=g^{im}{}_{F}{}^{}U_{m}^{[kl]},$$
(34)
where $`{}_{F}{}^{}U_{m}^{[kl]}`$ mean von Freud superpotentials and

$${}_{BT}{}^{}t_{}^{kl}:=g_E^{ki}t_i^l+\frac{g_{,p}^{mk}}{\sqrt{|g|}}{}_{F}{}^{}U_{m}^{[lp]}$$
(35)
is the Bergmann-Thomson gravitational energy-momentum pseudotensor .

One can interpret the Bergmann-Thomson expression (27) as a sum of the spinorial part

$$S^{ikl}:=_FU^{i[kl]}_FU^{k[il]}$$
(36)
and orbital part

$$O^{ikl}:=\sqrt{|g|}\left(y_{BT}^it^{kl}y_{BT}^kt^{il}\right)$$
(37)
of the gravitational angular momentum densities.

Substituting (26) into (23) we will get the canonical angular supermomentum tensor for matter

$`{}_{m}{}^{}S_{}^{abc}(P;v^l)`$ $`=`$ $`2[(2\widehat{v}^a\widehat{v}^p\widehat{g}^{ap})_p\widehat{T}^{bc}`$ (38)
$``$ $`(2\widehat{v}^b\widehat{v}^p\widehat{g}^{bp})_p\widehat{T}^{ac}].`$ (39)
If we substitute (27) into (23), then we will obtain the gravitational canonical angular supermomentum tensor††††††We have also eliminated in (33) the components of Ricci tensor with the help of the Einstein equations.

$`{}_{g}{}^{}S_{}^{abc}(P;v^l)`$ $`=`$ $`\alpha (2\widehat{v}^p\widehat{v}^t\widehat{g}^{pt})[\beta (\widehat{g}^{ac}\widehat{g}^{br}\widehat{g}^{bc}\widehat{g}^{ar})_{(t}\widehat{E}_{pr)}`$ (40)
$`+`$ $`2\widehat{g}^{ar}_{(t}\widehat{R}_p^{(b}{}_{r)}{}^{c)}2\widehat{g}^{br}_{(t}\widehat{R}_p^{(a}_{r)}^{c)}`$ (41)
$`+`$ $`2/3\widehat{g}^{bc}\left(_r\widehat{R}_{(t}^r{}_{p)}{}^{a}\beta _{(p}\widehat{E}_{t)}^a\right)`$ (42)
$``$ $`2/3\widehat{g}^{ac}\left(_r\widehat{R}_{(t}^r{}_{p)}{}^{b}\beta _{(p}\widehat{E}_{t)}^b\right)].`$ (43)
In vacuum, i.e., when $`T_{ik}=0E_{ik}:=T_{ik}1/2g_{ik}T=0`$, the gravitational canonical angular supermomentum tensor $`{}_{g}{}^{}S_{}^{abc}(P;v^l)=()_gS^{bac}(P;v^l)`$ simplifies to the form

$`{}_{g}{}^{}S_{}^{abc}(P;v^l)`$ $`=`$ $`2\alpha (2\widehat{v}^p\widehat{v}^t\widehat{g}^{pt})[\widehat{g}^{ar}_{(p}\widehat{R}_t^{(b}_{r)}^{c)}`$ (44)
$``$ $`\widehat{g}^{br}_{(p}\widehat{R}_t^{(a}{}_{r)}{}^{c)}].`$ (45)
It is interesting that the orbital part $`0^{ikl}=\sqrt{|g|}\left(y_{BT}^it^{kl}y_{BT}^kt^{il}\right)`$ of the $`{}_{g}{}^{}M_{}^{ikl}`$ gives no contribution to the tensor $`{}_{g}{}^{}S_{}^{abc}(P;v^l)`$. Only spin part $`S^{ikl}=_FU^{i[kl]}_FU^{k[il]}`$ gives non-zero contribution to this tensor. Notice also that the canonical angular supermomentum tensors $`{}_{g}{}^{}S_{}^{abc}(P;v^l)`$ and $`{}_{m}{}^{}S_{}^{abc}(P;v^l)`$, gravitation and matter, needn’t any radius vector for their own existing.

## III Summary about application of the superenergy and angular supermomentum tensors to analysis of the gravitational waves

One can easily check by a direct calculation that the canonical superenergy tensor $`{}_{g}{}^{}S_{i}^{k}(P;v^l)`$ and the canonical angular supermomentum tensor $`{}_{g}{}^{}S_{}^{ikl}(P;v^m)`$ give positive-definite superenergy density, non-vanishing superenergy flux and non-vanishing angular supermomentum flux for every known solution to the vacuum Einstein equations which represents a gravitational wave with $`R_{iklm}0`$. This is true in any admissible coordinates and it all was showed in our papers \[10,12–14\]. It results from the above fact that every real gravitational wave, which has $`R_{iklm}0`$, must carry gravitational energy-momentum and gravitational angular momentum. If not, then there would be a contradiction between “energy level” and “superenergy level” because the canonical superenergy and angular supermomentum tensors originated by some kind of averaging of the suitable pseudotensors. Notice also in this context that, as it follows from the definition (4), the $`ϵ(\mathrm{\Omega }):=()_gS_{ik}(P;v^l)v^iv^k\underset{\mathrm{\Omega }}{}\sigma (P;y)𝑑\mathrm{\Omega }`$ gives the approximate (relative with respect P) energy contained in the sufficiently small domain $`\mathrm{\Omega }`$ and $`P^i(\mathrm{\Omega }):=()1/2\left(\delta _k^iu^iu_k\right)_gS_l^k(P;v^m)v^l\underset{\mathrm{\Omega }}{}\sigma (P;y)𝑑\mathrm{\Omega }`$ gives Poynting’s vector of the (relative with respect P) energy flux contained in $`\mathrm{\Omega }`$. $`ϵ(\mathrm{\Omega })`$ and $`P^i(\mathrm{\Omega })`$ do not vanish in any admissible coordinates for a real gravitational wave which has $`R_{iklm}0`$.

So, the non-zero superenergy density and its non-zero flux really demand non-zero gravitational energy and its non-zero flux.

On the “energy-momentum level”, where we use pseudotensors, the very important fact that any gravitational wave which has $`R_{iklm}0`$ always transfer energy-momentum is camouflaged in some coordinates by energy-momentum of the inertial forces field.

The analogical considerations and conclusion which are based on the definition (23) are valid also for the angular momentum.

## IV Concluding remarks

If you want to get the correct information about energy-momentum and angular momentum of the real gravitational field by use pseudotensors, then you must use the pseudotensors in very special coordinates only (see eg. \[25-27\]). Namely, you must use pseudotensors in a coordinates in which $`\mathrm{\Gamma }_{kl}^i`$ describe only the real gravitational field. The examples of such coordinates are given, e.g., by global Bondi-Sachs coordinates for a closed system  or, in general, by normal coordinates NC(P) \[20-22, 28,29\].

In order to get information on gravitational energy-momentum and angular momentum in arbitrary admissible coordinates one must use covariant expressions which depend on curvature tensor. Our canonical superenergy and angular supermomentum tensors are exactly the quantities of such a kind. In application to gravitational radiation these quasilocal quantities unambiguously show that any gravitational waves which has $`R_{iklm}0`$ always transfer energy-momentum and angular momentum. So, the conclusions given in the papers like , where authors use pseudotensors only are uncorrect.‡‡‡‡‡‡In the paper  the Author uses also stationary Tolman integral to analysis of a dynamical system. Of course, this is unjustified physically.

## A The canonical energy-momentum complex in General Relativity

One can easily transform \**\**\** It is the most easily to do by use the formalism of the exterior differential forms. Einstein equations to the superpotential form

$$\sqrt{|g|}\left(T_i^k+_Et_i^k\right)={}_{F}{}^{}U_{i}^{[kl]}{}_{,l}{}^{}.$$
(A1)
In the last formula $`T^{ik}=T^{ki}`$ is the symmetric energy-momentum tensor for matter, $`{}_{E}{}^{}t_{i}^{k}`$ mean the components of the Einstein canonical energy-momentum pseudotensor for gravitational field and $`{}_{F}{}^{}U_{i}^{[kl]}`$ denote von Freud superpotentials. The sum

$$\sqrt{|g|}\left(T_i^k+_Et_i^k\right)=:_EK_i^k$$
(A2)
is called the canonical Einstein energy-momentum complex for matter and gravitation in General Relativity and it is usually denoted by $`{}_{E}{}^{}K_{i}^{k}`$.

In extended form we have

$`{}_{E}{}^{}t_{i}^{k}`$ $`=`$ $`{\displaystyle \frac{c^4}{16\pi G}}\{\delta _i^kg^{ms}(\mathrm{\Gamma }_{mr}^l\mathrm{\Gamma }_{sl}^r\mathrm{\Gamma }_{ms}^r\mathrm{\Gamma }_{rl}^l)`$ (A3)
$`+`$ $`g_{,i}^{ms}[\mathrm{\Gamma }_{ms}^k1/2(\mathrm{\Gamma }_{tp}^kg^{tp}\mathrm{\Gamma }_{tl}^lg^{kt})g_{ms}`$ (A4)
$``$ $`1/2(\delta _s^k\mathrm{\Gamma }_{ml}^l+\delta _m^k\mathrm{\Gamma }_{sl}^l)]\},`$ (A5)
$${}_{F}{}^{}U_{i}^{[kl]}=\frac{c^4}{16\pi G}\frac{1}{\sqrt{|g|}}g_{ia}\left[(g)\left(g^{ka}g^{lb}g^{la}g^{kb}\right)\right]_{,b}.$$
(A6)
From (A.1) there follow the local, differential energy-momentum conservation laws

$$\left[\sqrt{|g|}\left(T_i^k+_Et_i^k\right)\right]_{,k}=0,$$
(A7)
for matter and gravitation.