text:
# Mechanisms for electron transport in atomic-scale one-dimensional wires: soliton and polaron effects

## I Introduction

How will an electronic system conduct in the limit where it becomes entirely one-dimensional? This became an issue in the late 1980s, with the fabrication of semiconductor structures that possessed only a single conducting “channel”, and which exhibited the phenomenon of quantized conductance . With nanoscale fabrication techniques, the class of effectively one-dimensional systems has expanded to include: atomic-scale structures on surfaces produced by scanning probe lithography , self-assembly , or directed assembly ; break junctions ; carbon nanotubes ; and conjugated organic molecules . The measurement of the transport properties of individual examples of such systems is now possible; for example, conduction in individual molecules has recently been measured by scanning probes  or bulk electrodes ; transport in nanotubes has been measured using lithographed electrodes . Similar measurements may also soon become possible for atomic-scale wires on surfaces.

Many of these recent examples approach the ideal of a system whose atomic—as well as electronic—structure is truly one-dimensional. As such they may be expected to exhibit enhanced effects of correlations during electron transport. Electron-electron correlations have consequences including the possible formation of Luttinger liquids ; on the other hand, the effects produced by electron-atom correlations can include the Peierls transition , which opens a gap at the Fermi energy and renders one-dimensional metallic systems semiconducting. If these quasi-one-dimensional systems are to have applications in nanotechnology, it is vital to understand whether currents can pass through them, despite the possibility of a Peierls transition. In this paper we address this question, in the context of molecular wires; specifically, we consider finite-length poly-acetylene chains described by a fully quantum version of the Su–Schrieffer–Heeger (SSH) model .

We have chosen to study molecular conductors here, rather than any of the other types, because electron-lattice coupling in bulk samples is relatively well studied . The Peierls transition manifests itself in an alternation of single and double bonds (i.e., dimerisation) along the molecule’s length, and bulk electron transport is dominated by mobile intrinsic defects formed when electrons or holes are injected into this dimerised structure. For our purposes the most important defects are polarons (which can be thought of as a local reduction in the dimerisation around an injected charge) and charged solitons (topological defects in the bond length alternation). The soliton has an associated electronic state at mid-gap and the dimerisation changes sign, passing through zero in the neighbourhood of the defect .

In order to understand molecular charge transport on nanometre lengthscales, we must develop a theory of coherent transport that accounts for polaron and soliton formation. We took the first steps recently when we showed that the tunneling of carriers through molecular wires at low temperature is dramatically enhanced by the formation of virtual polarons . Therefore, previous treatments of tunneling in molecules based on elastic scattering significantly underestimate the conductance . However, the boundary conditions used in our previous work prevented solitons from forming in the molecule, as the dimerisation was constrained to be equal at both ends. In the present paper we study chains with an odd number of monomers; the ground state of such a system possesses a single soliton in the centre . We are thus able to study for the first time the role of the moving soliton, with its associated mid-gap state, and the competition between polarons and solitons in coherent transport.

## II Physical model

The model for the molecular wires includes delocalized $`\pi `$-electrons interacting with quantum phonons. The electronic Hamiltonian is expressed in the basis of the one-electron eigenstates (labelled by $`n,m`$) of the reference system. The e-ph coupling is linear in the phonon displacements and induces transitions between electronic states. The molecule Hamiltonian is

$$H=\underset{n}{}ϵ_nc_n^{}c_n+\underset{\lambda }{}\mathrm{}\omega _\lambda a_\lambda ^{}a_\lambda +\underset{\lambda ,n,m}{}\gamma _{\lambda nm}(a_\lambda ^{}+a_\lambda )c_n^{}c_m,$$
(1)
where $`c_n^{}`$ ($`c_n`$) creates (annihilates) an electron in the $`n`$-th electronic state with energy $`ϵ_n`$ and $`a_\lambda ^{}`$ ($`a_\lambda `$) creates (annihilates) a quantum of energy $`\mathrm{}\omega _\lambda `$ in the eigenmode of vibration $`\lambda `$ of the isolated molecule. The values of $`ϵ_n`$, $`\omega _\lambda `$ and $`\gamma _{\lambda nm}`$ are calculated from the ground state of the neutral molecule (described by the SSH model ) containing an odd number $`N`$ of monomers. The vibrational eigenmodes $`V_\lambda `$ and frequencies of the molecule are calculated within the harmonic approximation .

To calculate the transport properties through the wire, the left and right ends of the molecule (atomic sites $`i=1`$ and $`i=N`$) are connected to metallic leads via hopping integrals $`v_{L,R}`$ respectively. The leads are modeled as one-dimensional semi-infinite tight-binding chains with on-site energy $`ϵ_{L,R}`$ and inter-site hopping matrix elements $`\beta _{L,R}`$ (with no $`e`$-ph coupling). The scattering states $`|\mathrm{\Psi }`$ for a single incoming charge carrier are expanded inside the molecule onto the eigenstates $`|n,\{n_\lambda \}=c_n^{}_\lambda (a_\lambda ^{})^{n_\lambda }/\sqrt{n_\lambda !}|0`$, of the non-interacting $`e`$-ph system  A similar expansion is used in the leads, $`c_n^{}`$ is then replaced by the charge creation operator $`d_i^{}`$ on site $`i`$.. $`n_\lambda `$ is the occupation number of the phonon mode $`\lambda `$ and $`|0`$ is the vacuum state (representing the neutral ground state of the whole system, with a definite number of electrons in each part). The single added carrier can be anywhere in the system (on the left, the right, or in the molecule) and interacts with phonons only when inside the molecule. We only consider current-carrying states in which a single electron is added to the ground state; we expect this assumption to hold since the interval between electron transmission is much greater than the transit time .

We now briefly outline the calculation technique. The transport problem is solved by mapping the many body (one-electron/many-bosons) problem onto a single-electron problem with many scattering channels . Each channel represents the different scattering processes by which the electron might exchange energy with the phonons. For an initial phonon distribution $`b\{m_\lambda \}`$ and an incoming electron from the left lead, the outgoing channels in the left and right leads are associated with energy-dependent reflection $`r_{ab}(ϵ)`$ and transmission $`t_{ab}(ϵ)`$ coefficients respectively ($`a\{n_\lambda \}`$ is the phonon distribution after scattering in each outgoing channel). In the leads, the scattering states take the asympotic form of propagating Bloch waves with amplitudes $`r_{ab}`$ (reflection) and $`t_{ab}`$ (transmission). The wave vectors are given by the dispersion relations in each channel. For example, the energy $`ϵ_{in}`$ of the incoming electron from the left lead is related to the wave vector $`k_b^L`$ by the tight-binding-like dispersion relation $`ϵ_{in}=ϵ_L+2\beta _L\mathrm{cos}k_b^L`$. For the final energy $`ϵ_{fin}`$ of the electron transmitted to the right, one has: $`ϵ_{fin}=ϵ_R+2\beta _R\mathrm{cos}k_a^R`$. One can then project out the leads from the problem and work in the molecular wire subspace to solve for the value of the scattering state $`|\mathrm{\Psi }(E)`$. This state is obtained from propagating the source term $`|s(E)`$ (incoming electron from the left) via the effective Green’s function $`G(E)`$ defined in the molecular wire subspace: $`|\mathrm{\Psi }(E)=G(E)|s(E)`$, where $`G`$ is given by $`G(E)=[EH\mathrm{\Sigma }_L(E)\mathrm{\Sigma }_R(E)]^1`$, $`H`$ is the molecular wire Hamiltonian defined in Eq. (1) and $`\mathrm{\Sigma }_{L,R}(E)`$ are complex potentials arising from embedding the molecular wire spectrum into the continuum of states associated with the leads. Overall energy is conserved, so $`ϵ_{in}`$ and $`ϵ_{fin}`$ are related by: $`E=ϵ_{in}+_\lambda m_\lambda \mathrm{}\omega _\lambda =ϵ_{fin}+_\lambda n_\lambda \mathrm{}\omega _\lambda `$ §§§ In this paper, we consider the limit of low temperatures: the initial phonon distribution is $`b\{m_\lambda \}=\{0\}`$ where all optic phonon modes are in the ground state (in this limit, the injection energy equals the total energy $`ϵ_{in}=E=ϵ_{fin}+_\lambda n_\lambda \mathrm{}\omega _\lambda `$). It is a very good approximation even at room temperature for all optic modes except the soliton translation. For the soliton translation the condition $`kT\mathrm{}\omega `$ restricts us to $`T230K`$..

The linear system $`|\mathrm{\Psi }=G|s`$ is solved for a finite size basis set by truncating the phonon subspace (i.e considering the lowest occupation numbers up to $`n_{\mathrm{occ}}^{\mathrm{max}}`$ in each mode). Furthermore the electron is coupled inside the molecule to a finite number $`N_{\mathrm{ph}}`$ of the most relevant phonon modes (see next section). A detailed analysis of the validity of such approximations can be found in Ref. . From the solution $`|\mathrm{\Psi }=G|s`$, one can obtain the reflection $`r_{ab}`$ and transmission $`t_{ab}`$ coefficients for all the channels and hence the currents flowing through the wire. One can also calculate the expection values of any correlation functions between the electron and phonon degrees of freedom.

We define the transmission probability $`T_{ab}=|t_{ab}(ϵ_{in})|^2\beta _R\mathrm{sin}k_a^R(ϵ_{fin})/(\beta _L\mathrm{sin}k_b^L(ϵ_{in}))`$  $`T_{ab}`$ contains as usual the ratio of the velocity in the outgoing channel to that in the incoming channel.. In our model, $`|t_{ab}|^2`$ is proportional to $`|i=N,a|\mathrm{\Psi }|^2=|N,a|G|s|^2`$ and the source term is proportional to the velocity of the incoming electron . Then, working in the real-space representation, the inelastic transmission probability can be rewritten in the following usual form  :

$`T_{ab}(ϵ_{fin},ϵ_{in})=4{\displaystyle \frac{v_L^2}{\beta _L}}\mathrm{sin}k_b^L(ϵ_{in}){\displaystyle \frac{v_R^2}{\beta _R}}\mathrm{sin}k_a^R(ϵ_{fin})`$ (2)
$`\times |i=N|G_{ab}(E)|i=1|^2,`$ (3)
where $`N|G_{ab}(E)|1`$ is the matrix element of the Green’s function $`G`$ taken between the left side $`i=1`$ and the right side $`i=N`$ of the molecule and the phonon configurations before ($`b`$) and after ($`a`$) scattering  The factors $`v_{L,R}^2/\beta _{L,R}\mathrm{sin}k_{b,a}^{L,R}`$ in Eq.(3) are related to the imaginary parts of the potentials $`\mathrm{\Sigma }_{L,R}`$..

## III Results

First we describe the vibrational modes needed to obtain the relevant distortions of the molecule. We have shown that polaron formation inside the molecule after charge injection is due to the coupling to the long-wavelength optic modes . For odd-number chains where a soliton defect is present, we have to include in addition the modes responsible for the motion and the deformation of the soliton. We therefore include the soliton translation (the ‘Goldstone mode’), the ‘amplitude mode’ related to the deformation of the soliton width, and higher-order deformation modes (for example the so-called third mode) . Figure 1 shows the phonons considered for a molecule of length $`N=99`$. As an illustrative example, we also show in Fig. 2 the dimerisation $`d_j`$ for neutral isolated even- and odd-length molecules. The dimerisation is obtained, in terms of atomic displacements $`u_j`$, from the staggered difference between adjacent bond lengths: $`d_j=(1)^j(u_{j+1}2u_j+u_{j1})`$. For even $`N`$, the neutral molecule is perfectly dimerized: the dimerisation $`d_c`$ is constant in the middle of the chain despite the end effects. For odd $`N`$, a soliton appears in the middle of the chain, acting as a domain wall separating the two domains of opposite sign of dimerisation.

Now we turn on the lattice deformations induced by electron propagation in the (odd $`N`$) molecule connected to the electrodes \**\**\** Calculations were performed for an electron coupled to the $`N_{\mathrm{ph}}=5`$ lowest frequencies modes shown in Fig. 1 and $`n_{\mathrm{occ}}^{\mathrm{max}}=4`$. Calculations were also done for other sets of parameters. The same qualitative physics remain when the results are converged versus the basis set size.. The lattice deformations induced by the tunneling electron are given by the correlation function $`\delta _\lambda ^{[i]}`$ between the electron density $`P_i=c_i^{}c_i`$ on site $`i`$ and the displacement $`\mathrm{\Delta }_\lambda =(a_\lambda +a_\lambda ^{})\sqrt{\mathrm{}/2M\omega _\lambda }`$ of the mode $`\lambda `$ as $`\delta _\lambda ^{[i]}=P_i^{}\mathrm{\Delta }_\lambda P_i/P_i`$. We also define the distortion on site $`j`$ due to the displaced modes when the electron is on site $`i`$ as $`u_j^{[i]}=_\lambda \delta _\lambda ^{[i]}V_\lambda (j)`$. The dimerisation pattern $`d_j^{[i]}`$ is then calculated from $`u_j^{[i]}`$.

We plot on Fig. 3 the absolute value of the dimerisation $`d_j^{[i]}`$ from which the constant dimerisation $`d_c`$ and the end effects have been subtracted (the corresponding pattern for the isolated molecule is shown on Fig. 2). Such a choice permits us to represent the defects with more contrast. For an injection energy $`E=0`$ at mid-gap, the corresponding dimerisation is shown on Fig. 3(a). The bright feature around the middle of the chain represents the soliton. In the absence of $`e`$-ph coupling, this would correspond to charge injection in resonance with the one-electron soliton level (see the transmission curve below). The soliton would then be immobile, its position fixed in the middle of the chain and its width unchanged (i.e. straight vertical feature around $`j=50`$). However due to the $`e`$-ph coupling, the soliton defect becomes mobile and its width varies while the electron progates. Similar behaviour is observed for other injection energies around mid-gap. As $`E`$ approaches the valence band edge the mechanisms become different, as shown on Fig. 3(b). In parallel with the soliton delocalization and deformation, one observes the formation of a polaron. This is characterized by the extra bright feature around the first diagonal in Fig. 3(b). In principle, the (virtual) polaron corresponds to a local reduction of the dimerisation around the (tunneling) electron . However it appears as an increase (positive number) since $`|d_j^{[i]}|`$ is plotted on Fig. 3(b). Furthermore, the corresponding dimerisation is not simply the superposition of the polaron and the soliton—the two defects interact strongly together. For example, the dimerisation pattern for an electron at position $`i25`$ is characteristic of a merging of both the polaron and soliton. Such mechanisms, leading to strong lattice distortions, would be expected to affect the transmission through the molecular wires—effects that we now consider.

We define an effective total transmission probability $`T(E)`$ arising from the contribution of the different outgoing channels as $`T(E)_aT_{ab}(ϵ_{fin},ϵ_{in})\delta (ϵ_{fin}+_\lambda n_\lambda \mathrm{}\omega _\lambda ϵ_{in})`$ (see Fig. 4). Without $`e`$-ph coupling, the transmission curve presents the usual features: resonances at energies corresponding to the one-electron levels of the molecule with almost perfect transmission, and strong suppression of the transmission in the band-gap of the molecule for which propagation occurs by tunneling. The soliton resonance (at mid-gap) is much narrower than the other resonances. This is to be expected since the corresponding one-electron (soliton) state is much more localized than the other states. With $`e`$-ph interaction, we have already shown that a polaron can be formed in the molecule and that the soliton position and width are strongly modified. Such mechanisms affect the transmission in the following way: (i) the delocalization and deformation of the soliton broaden but lower the resonance peak at mid-gap, (ii) for larger injection energies, the formation of the polaron effectively reduces the apparent band-gap, and a polaron resonance peak appears (at $`E0.42`$ eV for the $`N=99`$ wire) inside the original gap. This effective band-gap reduction correlatively increases the transmission in the tunneling regime (i.e. for injection energies in the gap).

## IV Conclusion

The results presented above demonstrate the complexity of electronic transport through one-dimensional atomic-scale wires. We have illustrated the mechanisms of electron transport in molecular wires when $`e`$-ph interactions are included. Owing to $`e`$-ph coupling, polarons can be formed inside the molecule. Polaron propagation is the main mechanism of transport through perfectly dimerised (semiconducting) molecules. The presence of a mid-gap state associated to a soliton defect in the middle of the molecule involves different mechanisms for the transport. For injection energies around the mid-gap state, the delocalisation and deformation of the soliton is the main mechanism for electron transfer. With larger injection energies (still inside the energy band-gap of the molecule), a virtual polaron can be formed. The transport is then associated with more complicated mechanisms involving the interaction of both polaron and soliton. However in most cases, the effective reduction of the band-gap due to polaron formation, the delocalization and deformation of the soliton and the polaron-soliton interaction increase the transmission in the tunneling regime. Although we studied a model derived from organic molecules, there is good reason to believe that the $`e`$-ph coupling in other one-dimensional nanoscale wires is likely to give rise to similar phenomena. In the case of carriers injected into dangling-bond lines on the Si(001) surface , polaron states which are in many respects analogous to those in molecular wires are also formed . Specifically, we expect the present results to be valid for any atomic-scale wire in which there is a degeneracy between two different values of some order parameter (the dimerisation in the molecular case, the surface dimer buckling in the dangling bond line) which are related by a discrete symmetry and strongly coupled to the electrons.

###### Acknowledgements.

We acknowledge financial support from the U.K. EPSRC (Grant GR/M09193).


format: text

text:
# Symmetry-breaking instability in a prototypical driven granular gas

## I I. Introduction

Granular materials play an important role in industrial applications, geophysics and astrophysics. They are also of a great general interest to physicists, as each of the “phases” of granular materials: solid, liquid and gas have unusual properties that distinguish them from their classic (atomic) counterparts Jaeger ; Duran . We will consider granular gas (or rapid granular flow) and focus on a variant of clustering instability, a striking tendency of granular gases to form dense clusters. Clustering instability was first discussed in the context of a freely ”cooling” granular gas Hopkins ; Goldhirsch ; McNamara . Following these works, related clustering phenomena were investigated in driven granular gases as well, both in experiments Kudrolli ; Urbach  and in particle simulations Kadanoff2 ; Grossman ; Esipov .

Granular clustering results from energy losses by inelastic collisions, and it is a manifestation of thermal condensation instability, also encountered in other fields, for example in gases and plasmas that cool by their own radiation plasma . Since the discovery of the clustering instability, the validity of granular hydrodynamics (see Ref. Campbell  for a review) has been under scrutiny Kadanoff1 . In contrast to the clustering in a freely ”cooling” granular gas, where one deals with a complex time-dependent process, steady states are achievable in driven granular systems. One of the simplest settings of this type is driving the granulate by a side wall at a zero gravity. Therefore, an ensemble of inelastically colliding hard spheres, confined in a box and driven by one or two “thermal” walls has served as a prototypical driven granular system Kudrolli ; Kadanoff2 ; Grossman ; Brey ; Tobochnik ; Kudrolli2 ; LMS . Steady states of this system have served as test beds for granular hydrodynamics and its violations. The first analysis of this system in the physical literature was performed, in one dimension, by Kadanoff et al. Kadanoff2 . The nearly elastic particles were constrained to move on a straight line with energy input from the boundaries. Particle simulations Kadanoff2  showed that, for typical initial conditions, the system evolves to a state where the particles are separated into two groups. Almost all particles form a cluster in a small region of space, where they move with very small velocities, while a very few remaining particles move with high velocities. Clearly, this steady state cannot be described by granular hydrodynamics (actually, by any coarse-grained theory). Therefore, the results of Kadanoff et al. Kadanoff2  brought into question the validity of granular hydrodynamics in general.

This question was addressed in two subsequent theoretical works Grossman ; Esipov  that dealt with similar systems in two dimensions. Esipov and Pöschel Esipov  investigated an ensemble of nearly elastic hard disks in a circular box with the circumference serving as thermal wall. Grossman et al. Grossman  considered a rectangular box, one side of which served as thermal wall. Particle simulations Grossman ; Esipov  showed granular clusters: dense and “cold” regions of granulate developing away from the thermal wall. In terms of thecoarse-grained particle density, these steady-state clusters had simple shapes: azimuthally uniform (circular state) Esipov  and laterally uniform (strip state) Grossman . Grossman et al. also showed that, for nearly elastic collisions, the strip state is describable by a steady-state solution of granular hydrodynamic equations. The empiric constitutive relations suggested by Grossman et al. used simple interpolations between the low-density limit, where the constitutive relations are derivable systematically dilute , and high-density limit where, close to the dense packing, free volume arguments can be used.

The results of Refs. Grossman ; Esipov  showed that the anomaly observed in the one-dimensional setting Kadanoff2  does not persist in higher dimensions. Clustered states qualitatively similar to those of Grossman et al. were observed in experiment of Kudrolli et al. Kudrolli  who investigated a system of spherical particles in a box, rolling on a smooth surface and driven by a rapidly vibrating side wall. The number of particles served as the control parameter in Ref. Kudrolli . A dense cluster of the strip type was observed when the number of particles was big enough, in much the same way as in Ref. Grossman . The basic physics of the strip state is simple and can be explained by the following hydrodynamic argument. Because of the inelastic collisions the particle random motion slows down (that is, the granular temperature decreases) with the increase of the distance from the driving wall. To maintain the momentum balance, the granular density should increase with this distance. When the total number of particles is big enough (the rest of parameters being the same), the density contrast becomes large, and the enhanced density region away from the driving wall is observed as the strip state.

The prototypical system exhibits many interesting phenomena of non-hydrodynamic nature. These include inelastic collapse Esipov ; Tobochnik , possible lack of scale separation Grossman , non-Gaussianity in the particle velocity distribution Grossman ; Kudrolli2 , normal stress difference and pressure nonuniformities Brey , etc. For nearly elastic collisions, however, granular hydrodynamics was shown to yield an accurate quantitative description in the dilute limit Brey , and reasonably accurate description for moderate and high granular densities Grossman . Of course, the nearly elastic limit is quite restrictive for most of realistic granular flows. Still, this limit is conceptually important just because granular hydrodynamics can be used there. Granular hydrodynamics has a great predictive power and helps to identify important collective phenomena (shear flows and vortices, shocks, different modes of clustering flows etc.) that are difficult, if not impossible, to identify and predict in the language of individual particles. Once identified, these phenomena can then be investigated in experiment and simulations in more general (not necessarily hydrodynamic) formulations.

Therefore, granular hydrodynamics provide a leading-order approach to a big ensemble of nearly-elastically colliding grains. This approach has been taken recently by Livne et al. LMS  who employed granular hydrodynamics for a stability analysis of the strip state described above. The analysis revealed a spontaneous symmetry-breaking instability of the strip state with respect to perturbations along the strip. Well within the instability region, the clustered states with broken symmetry, found by a numerical solution of the steady state hydrodynamic equation, are strongly localized in the lateral direction: most of the particles are located in dense “islands” LMS . These results indicate that the prototypical system can show a non-trivial behavior even in the leading-order, hydrodynamic limit. Indeed, this systems can be put into the list of pattern-forming systems far from equilibrium Cross .

The present work focuses on a more detailed stability analysis of this system. Our first objective is to check to what extent the symmetry-breaking instability predicted in Ref. LMS  is sensitive to the constitutive relations. Livne et al. LMS  employed the empiric relations suggested by Grossman et al. Grossman . Here we shall use the better known Jenkins-Richman (JR) relations Jenkins . While the relations of Grossman et al. are more accurate for high densities (even including those close to the dense packing limit), the JR relations should work better at low and intermediate densities. We shall see, however, that the marginal stability curves, obtained with these two sets of relations, are not much different from each other. This implies that the symmetry-breaking instability is robust. Our second objective is to get more insight into the marginal stability problem and, where possible, to obtain analytic results. We shall show that the marginal problem is equivalent to an eigenvalue problem of quantum mechanics. An important finding here is a universal behavior of the marginal stability curves in the limit when the density perturbations are strongly localized at the elastic wall opposite to the thermal wall. In the dilute limit, we obtain analytically the short- and long-wavelength asymptotics of the marginal stability curves and density eigenfunctions. We also give the physical interpretation to the symmetry-breaking instability and to the density borders of the instability region.

The rest of the paper is organized as follows. In Sec. 2 we formulate the model and briefly describe the strip state: the simplest steady state of the prototypical system. Section 3 presents marginal stability analysis of the strip state and compares the results obtained for two different sets of constitutive relations. More results on marginal stability, including some analytic results in the dilute limit, are presented in Sec. 4. Section 5 includes a discussion and summary.

## II II. Prototypical system and strip state

.

The prototypical driven granular system in two dimensions include inelastically colliding hard disks of diameter $`d`$ and mass $`m=1`$, moving in a box with dimensions $`L_x\times L_y`$. The gravity force is zero. Collisions of disks with the walls $`x=0,y=0`$ and $`y=L_y`$ are assumed elastic. The wall $`x=L_x`$ is ”thermal” wall: upon collision a particle is assigned a random velocity taken from a Gaussian ensemble with temperature $`T_0`$. Energy input at the thermal wall balances the energy dissipation due to inter-particle collisions, so the system can reach a steady state. We shall parameterize the inelasticity of the particle collisions by a constant normal restitution coefficient $`r`$ and work in the nearly elastic limit: $`1r^21`$. In this limit, the Navier-Stokes granular hydrodynamics is expected to be sufficiently accurate in a system with a big number of particles and small Knudsen number. The possible steady states of the system are described by the steady state versions of the momentum and energy balance equations:

$$p=const\text{and}(\kappa T)=I,$$
(1)
where $`p`$ is the granular pressure, $`T`$ is the granular temperature, $`\kappa `$ is the thermal conductivity and $`I`$ is the rate of energy losses by collisions. We assume that the number density $`n`$ is not too big: $`n/n_c0.5`$, where $`n_c=2/(\sqrt{3}d^2)`$ is the (hexagonal) dense packing density. This assumption enables us to employ the constitutive relations derived by Jenkins and Richman Jenkins . For the steady state problem, the required constitutive relations include the equation of state $`p=p(n,T)`$ and relations for $`\kappa `$ and $`I`$ in terms of $`n`$ and $`T`$.

Let us introduce scaled coordinates: $`𝐫/L_x𝐫`$. In the new coordinates the box dimensions are $`1\times \mathrm{\Delta }`$, where $`\mathrm{\Delta }=L_y/L_x`$ is the aspect ratio. Introducing the normalized inverse density $`z(x,y)=n_c/n(x,y)`$, one can rewrite the energy balance equation in Eq. (1) in terms of $`z(x,y)`$:

$$(F(z)z)=\eta Q(z),$$
(2)
where $`F(z)=A(z)B(z)`$,

$`A(z)`$ $`=`$ $`{\displaystyle \frac{G\left[1+\frac{9\pi }{16}\left(1+\frac{2}{3G}\right)^2\right]}{z^{1/2}(1+2G)^{5/2}}},`$
$`B(z)`$ $`=`$ $`1+2G+{\displaystyle \frac{\pi }{\sqrt{3}}}{\displaystyle \frac{z(z+\frac{\pi }{16\sqrt{3}})}{(z\frac{\pi }{2\sqrt{3}})^3}},`$
$`Q(z)`$ $`=`$ $`{\displaystyle \frac{6}{\pi }}{\displaystyle \frac{z^{1/2}G}{(1+2G)^{3/2}}},`$
$`G(z)`$ $`=`$ $`{\displaystyle \frac{\pi }{2\sqrt{3}}}{\displaystyle \frac{z\frac{7\pi }{32\sqrt{3}}}{(z\frac{\pi }{2\sqrt{3}})^2}},`$ (3)
and $`\eta =(2\pi /3)(1r)(L_x/d)^2`$. Notice that, for an arbitrary small but finite inelasticity $`1r`$, the dimensionless parameter $`\eta `$ can be made arbitrary large, if the system size $`L_x`$ is large enough. The parameter $`\eta `$ differs from the parameter $``$ used by Livne et al. LMS  only by a numerical factor of order unity. Of most interest are regimes where $`\eta 1`$, see below.

The boundary conditions for Eq. (2) are determined by the properties of the particle-wall interactions. At the elastic walls $`x=0`$, $`y=0`$ and $`y=\mathrm{\Delta }`$ we should prescribe a zero normal component of the heat flux. In terms of the inverse density $`z`$ we have $`_nz=0`$ at these three walls. Here index $`n`$ denotes the gradient component normal to the wall. The constant temperature at the “thermal” wall $`x=1`$ yields the condition $`z(x=1,y)/y=0`$. To make the formulation of the problem complete, one more condition is needed. In experiment or particle simulations, the number of particles $`N`$ is fixed. This yields a normalization condition:

$$\frac{1}{\mathrm{\Delta }}_0^1_0^\mathrm{\Delta }\frac{dxdy}{z(x,y)}=f,$$
(4)
where $`f=n/n_c`$ is the area fraction of the grains and $`n=N/(L_xL_y)`$ is the average number density of the grains.

Equations (2)-(4) and the boundary conditions make a complete set. One can see that the governing parameters of this system are the scaled parameters $`\eta `$, $`f`$ and $`\mathrm{\Delta }`$. If the system is infinite in the $`y`$-direction, only two governing parameters: $`\eta `$ and $`f`$ remain. Notice that the steady-state density distributions are independent of the wall temperature $`T_0`$, in contrast to the similar problem with gravity, where the gravity acceleration, combined with $`T_0`$ and the (finite) system size in the direction of gravity, would make an additional governing parameter.

The laterally uniform steady state (strip state) corresponds to the one-dimensional ($`y`$-independent) solution $`z=Z(x)`$. It is described by the equations

$`(FZ^{})^{}=\eta Q,Z^{}|_{x=0}=0,`$
$`\text{and}{\displaystyle _0^1}Z^1(x)𝑑x=f,`$ (5)
where primes stand for the $`x`$ derivatives. For the strip state, the boundary condition at the wall $`x=1`$ drops out. This implies, in particular, that the density profile of the strip state is independent of the exact nature of the driving wall (thermal or vibrating wall) boundary . This degeneracy of the strip state is caused by the character of particle interaction: the hard-core potential does not introduce any characteristic energy Esipov . Notice that, instead of prescribing the grain area fraction $`f`$, one can prescribe the inverse density $`Z=Z_0`$ at $`x=0`$. This condition, combined with the no-flux condition at $`x=0`$ defines a Cauchy problem for $`Z(x)`$. Solving the Cauchy problem, one can then compute, from the last equation in Eq. (5), the respective value of $`f`$. At fixed $`\eta `$, there is a one-to-one correspondence between $`Z_0`$ and $`f`$. Therefore, an alternative parameterization of the strip state is given by the scaled numbers $`\eta `$ and $`Z_0`$. We shall see below that the same property keeps (and can be conveniently used) in the marginal stability problem. Figure 1 shows a typical example of the scaled density profile $`n(x)/n_c`$ of the strip state obtained by solving Eqs. (5) numerically.

## III III. Instability of the strip state: marginal stability curves

In general, the strip state is only one of the possible solutions of Eq. (2). Because of its nonlinearity, Eq. (2) may have additional solutions satisfying the same boundary conditions. When exist, these additional solutions are truly two-dimensional: the translational symmetry along $`y`$ is broken. An important class of these solutions bifurcate supercritically from the strip state LMS . Therefore, close to the bifurcation point, these solutions can be found by linearizing Eq. (2) around the strip state. A similar analysis was performed in Ref. LMS  for the constitutive relations of Grossman et al. Grossman . In the framework of a time-dependent hydrodynamic formulation, this analysis corresponds to marginal stability analysis of the strip state with respect to small perturbations along the strip LMS .

Substituting $`z(x,y)=Z(x)+\psi _k(x)\mathrm{cos}ky`$ and linearizing Eq. (2) with respect to the small correction $`\psi _k(x)\mathrm{cos}ky`$, we obtain:

$$\varphi ^{\prime \prime }\left(\frac{\eta Q_Z}{F}+k^2\right)\varphi =0.$$
(6)
Here $`\varphi (x)=F(x)\psi _k(x)`$, functions $`F`$ and $`Q`$ are evaluated at $`z=Z(x)`$, and subscript $`Z`$ means the $`z`$ derivative evaluated at $`z=Z(x)`$. The boundary conditions are

$$\varphi ^{}(x=0)=0\text{and}\varphi (x=1)=0.$$
(7)
Equation (6) coincides with the Schrödinger equation for an even wave function $`\varphi (x)`$ of a particle in the potential well

$$U(x)=\{\begin{array}{cc}\frac{\eta Q_Z}{F}\hfill & \text{if}|x|<1,\hfill \\ +\mathrm{}\hfill & \text{otherwise}.\hfill \end{array}$$
(8)
The quantity $`k^2`$ serves as the energy eigenvalue. The energy levels in the potential (8) are always discrete, and there is an infinite number of them. However, as the wave number $`k`$ should be real, only negative or zero energy levels are admissible. At fixed values of $`\eta `$ and $`f`$, the potential (8) admits at most one such energy level. The absence of negative energy levels implies that, in the vicinity of the strip state, there are no steady states different from it. The presence of a negative energy level corresponds to a “weakly two-dimensional” steady state, bifurcating from the strip state. We shall exploit the quantum-mechanical analogy more fully in Sec. 4. Here we report some numerical results. Figure 2 shows the marginal stability curves: the curves $`k=k(f)`$ at different values of $`\eta `$, computed numerically. In these computations, the parameter $`\eta `$ was taken large enough. The strip state is unstable below the respective curve and stable above the curve. Notice that, at fixed $`\eta `$, the instability is possible only within a finite interval of $`f`$: $`f_1(\eta )<f<f_2(\eta )`$. The same property was reported in Ref. LMS  for another set of constitutive relations. We shall give a physical explanation to this finding in Sec. 4. Notice (see also Ref. LMS ) that, at large $`\eta `$, the high-density stability border $`f_2`$ is quite small. The curves in Fig. 2 are actually plotted in scaled coordinates: $`k\eta ^{1/2}`$ versus $`f\eta ^{1/2}`$. It can be seen that, in the scaled coordinates, all the curves exit from the same point of the horizontal axis $`f\eta ^{1/2}`$. In addition, the maxima of all the curves are equal. These observations will be also explained in Sec. 4.

If the system is infinite in the lateral direction, $`\mathrm{\Delta }=\mathrm{}`$, while $`\eta `$ and $`f`$ are fixed, a whole continuum spectrum of wave numbers between $`k=0`$ and $`k=k(\eta ,f)`$ is unstable. Both in experiment, and in numerical simulations $`\mathrm{\Delta }`$ is finite. In this case $`k`$ becomes discrete because of the boundary conditions: $`k=m\pi /\mathrm{\Delta }`$, where $`m=1,2,\mathrm{}`$. For each $`m`$ we can find the critical value of the aspect ratio $`\mathrm{\Delta }`$ (let us call it $`\mathrm{\Delta }_m`$) such that for $`\mathrm{\Delta }>\mathrm{\Delta }_m`$ the strip state looses stability with respect to the $`m`$-mode. Obviously, $`\mathrm{\Delta }_m=m\mathrm{\Delta }_1`$, and $`\mathrm{\Delta }_1`$ is the lowest critical value for the symmetry-breaking instability. Figure 3 shows $`\mathrm{\Delta }_1`$ as a function of $`f`$ for different values of $`\eta `$. For fixed $`\eta `$ and $`f`$, the strip state is unstable above the curve. It is seen from Fig. 3 that, in order to observe the symmetry-breaking instability, one does not need to work with very large aspect ratios $`\mathrm{\Delta }`$: it is sufficient if the system is big enough, so that parameter $`\eta `$ is sufficiently large.

To what extent is the symmetry-breaking instability sensitive to the precise form of the constitutive relations? We compared the marginal stability curves $`\mathrm{\Delta }=\mathrm{\Delta }_1(f)`$ for different values of $`\eta `$ with the respective curves LMS  found for the constitutive relations of Grossman et al. Grossman . A typical example of this comparison is shown in Fig. 4. One can see that, qualitatively, the results are the same: the both curves describe a symmetry-breaking instability at a critical value of the aspect ratio that depends on the area fraction. In both cases, there are sharp low- and high-density borders of instability region. Therefore, we can conclude that the instability is robust and does not require a very special form of the constitutive relations. On the other hand, there is a noticeable (about 15%) difference in the exact positions of the marginal stability curves, so the instability provides a good quantitative test for constitutive relations of granular hydrodynamics.

## IV IV. More results on marginal stability

In this Section we investigate the marginal stability problem in more detail and obtain some analytic results in the dilute limit.

### IV.1 A. Localization and universality

Let us characterize the strip state by the scaled parameters $`\eta `$ and $`Z_0`$ and introduce a different rescaling of the coordinate: $`\overline{x}=x\eta ^{1/2}`$. In terms of the original, physical coordinate $`x_{ph}`$ the new rescaling is independent of the system size:

$$\overline{x}=\left(\frac{2\pi }{3}\right)^{1/2}\frac{(1r)^{1/2}x_{ph}}{d}.$$

Equations (5) for the strip state become

$`(FZ^{})^{}=Q,Z(\overline{x}=0)=Z_0`$
$`\text{and}Z^{}(\overline{x}=0)=0,`$ (9)
while the eigenvalue problem (6) and (7) reads

$`\varphi ^{\prime \prime }\left({\displaystyle \frac{Q_Z}{F}}+\overline{k}^2\right)\varphi =0,\varphi ^{}(\overline{x}=0)=0`$
$`\text{and}\varphi (\overline{x}=\eta ^{1/2})=0.`$ (10)
Now the primes denote the derivatives with respect to $`\overline{x}`$, while $`\overline{k}=k\eta ^{1/2}`$ is the new scaled wave number. Like the coordinate $`\overline{x}`$, the new scaled wave number $`\overline{k}`$ is independent of $`L_x`$:

$$\overline{k}=\left(\frac{3}{2\pi }\right)^{1/2}\frac{k_{ph}d}{(1r)^{1/2}},$$

where $`k_{ph}`$ is the physical wave number. The problem (9) and (10) is determined by two parameters: $`Z_0`$ and $`\eta `$. However, $`\eta `$ enters the rescaled equations only in one place: in the last boundary condition in Eq. (10) where it determines the scaled system size. If the wave function $`\varphi (\overline{x})`$ is strongly localized in the potential well $`U(\overline{x})`$ (correspondingly, the negative energy level is sufficiently deep), the results for $`\overline{k}`$ and $`\varphi (\overline{x})`$ become independent of $`\eta `$ at sufficiently large $`\eta `$. Indeed, in this case one can safely move the boundary $`\overline{x}=\eta ^{1/2}`$ to infinity. It is important that, in this case, the exact form of the boundary condition at the driving wall becomes insignificant, leading only to exponentially small corrections boundary . This universal “localization regime” was discovered in Ref. LMS  that employed a different set of constitutive relations. The criterion for localization can be obtained from the requirement that the localization length (which is of order $`\overline{k}^1`$) be much smaller than the (scaled) system size $`\eta ^{1/2}`$. In the physical units it corresponds to the short-wavelength limit of the bifurcating solution: $`k_{ph}L_x1`$. To fulfill this criterion, $`\overline{k}`$ should be far enough from the borders of the instability interval, where $`\overline{k}`$ vanishes. In the next subsection we will work in the dilute limit and rewrite this criterion in terms of the governing parameters of the problem.

Figure 5 shows the marginal stability curves $`\overline{k}=\overline{k}(Z_0)`$ for different values of $`\eta `$, obtained numerically. Instead of the borders $`f_1(\eta )`$ and $`f_2(\eta )`$ of the instability interval in terms of parameter $`f`$, the respective borders in terms of parameter $`Z_0`$ appear. One can see that, for large values of $`\eta `$, the marginal stability curves coincide in a wide region of $`Z_0`$ not too close to the borders of the instability interval. This region corresponds to strong localization. Figures 6-8 show the form of the potential (8) and the negative energy level $`k^2`$ in three characteristic cases (in these figures we returned to the rescaling of the coordinates and wavenumber by the system size $`L_x`$). Figure 6 corresponds to the region of parameters where the energy level is deep and eigenfunction is localized. Figures 7 and 8 correspond to the parameter regions close to the low- and high-density borders of the instability, respectively. There is no localization here. Notice the qualitative change in the form of the potential near the high-density stability border. Figure 9 shows the respective eigenfunctions in these three cases.

### IV.2 B. Marginal stability borders and physics of the instability

The low- and high-density stability borders $`f_1(\eta )`$ and $`f_2(\eta )`$ (or respective borders in terms of $`Z_0`$) are determined by the zero-eigenvalue ($`\overline{k}=0`$) solution of Eq. (10). This solution can be found if we know the strip solution $`Z(\overline{x};Z_0)`$ of Eq. (9). Indeed, it is easy to check that function $`\varphi _0(\overline{x})=FZ/Z_0`$ is a solution of Eq. (10) with $`\overline{k}=0`$, satisfying the boundary condition $`\varphi ^{}(\overline{x}=0)`$. Employing the second boundary condition $`\varphi (\overline{x}=\eta ^{1/2})=0`$ we obtain $`Z_1/Z_0=0`$, where $`Z_1=Z(\overline{x}=\eta ^{1/2};Z_0)`$. For a given $`\eta `$, this equation is an algebraic equation for $`Z_0`$. Our numerical results imply that this equation has only two solutions corresponding to the low- and high-density instability borders. The instability borders have a clear physical meaning which sheds light on the physics of the instability. Let us consider the granular pressure $`p=nT(1+2G)`$ Jenkins  of the strip state, and introduce a scaled pressure

$$P=\frac{p}{n_cT_0}=\frac{1+2G}{Z}\frac{T}{T_0}.$$

As $`P`$ is independent of the coordinates, we can compute it at the thermal wall $`\overline{x}=\eta ^{1/2}`$. Here $`T=T_0`$ and $`Z=Z_1`$, so we arrive at

$$P=\frac{1+2G(Z_1)}{Z_1}=P(\eta ,f).$$

Now let us compute the derivative $`P/f`$ at a constant $`\eta `$:

$$\frac{P}{f}=\frac{P}{Z_1}\frac{Z_1}{Z_0}\frac{Z_0}{f}.$$

One can easily check that the first and third multipliers in the right hand side of this relation are always negative. Therefore, the sign of $`P/f`$ is determined by the sign of $`Z_1/Z_0`$. As we have seen, the marginal stability borders are determined by equation $`Z_1/Z_0=0`$. Therefore, the steady-state pressure has its extremum points points exactly at the points $`f_1`$ and $`f_2`$. Figure 10 shows an example of the dependence $`P=P(f)`$ at a constant $`\eta `$ for the strip state, found numerically. One can see, that $`P`$ is a decreasing function of $`f`$ within the instability interval $`(f_1,f_2)`$, and an increasing function of $`f`$ outside the interval. The physical interpretation of these results is clear. The presence of the anomalous (falling) part of the $`P(f)`$ curve indicates instability, and it is caused by the destabilizing role of collisional heat losses Goldhirsch ; plasma . We can say that, on the interval $`(f_1,f_2)`$, the granulate has negative lateral compressibility. At $`f<f_1`$ the heat losses are too small to cause instability. The presence of the high-density border $`f_2`$ is caused by the finite-density corrections to the constitutive relations (that is, by the finite size of the particles). This is in contrast to radiative condensations in gases and plasmas, where such a stabilizing effect would be absent. Now consider a small density modulation of the strip state with a very long wavelength: $`k0`$. For this perturbation, the stabilizing effect of the lateral heat conduction vanishes, and the negative compressibility makes the strip state unstable on the interval $`(f_1,f_2)`$. For a non-zero $`k`$, the lateral heat conduction has a stabilizing effect. Therefore, a density modulation with too a short lateral wavelength should be stable, as Fig. 2 indeed shows.

### IV.3 C. Dilute limit

In this subsection we shall work in the dilute limit and derive several analytic results. We shall see that, at large $`\eta `$, the dilute limit faithfully reproduces the low-density parts of the marginal stability curves.

#### IV.3.1 1. Strip state and marginal stability problem

In the dilute limit, $`Z1`$, Eq. (9) for the strip state becomes

$$(Z^{3/2})^{\prime \prime }=3Z^{1/2},$$
(11)
where the primes denote the $`\overline{x}`$ derivatives. The boundary conditions are $`Z^{}(\overline{x}=0)=0`$ and $`Z(\overline{x}=0)=Z_0`$, where $`Z_0`$ is related to $`f`$ and $`\eta `$ by the normalization condition $`_0^\sqrt{\eta }Z^1𝑑\overline{x}=f\eta ^{1/2}`$. The solution of this problem is elementary:

$$\overline{x}=\frac{Z_0}{2}\left(\text{arccosh}\sqrt{\zeta }+\sqrt{\zeta ^2\zeta }\right),$$
(12)
where $`\zeta =Z/Z_0`$ and

$$Z_0=\frac{2\eta ^{1/2}}{f\eta ^{1/2}+\frac{1}{2}\mathrm{sinh}(2f\eta ^{1/2})}.$$
(13)
(Returning for a moment to the old rescaling of the coordinate, $`x_{ph}/L_xx`$, one can see that the density profile (12) is determined by a single parameter: $`\xi =f\eta ^{1/2}`$.) Eq. (10) for $`\varphi (\overline{x})=(\sqrt{3}/2)Z^{1/2}(\overline{x})\psi _k(\overline{x})`$ takes the form

$$\varphi ^{\prime \prime }\left(\kappa \frac{1}{\zeta ^2}\right)\varphi =0,$$
(14)
where $`\kappa =\overline{k}^2Z_0^2`$, and $`\zeta =\zeta (\overline{x})`$ is given, in an implicit form, by Eq. (12). The boundary conditions for Eq. (14) remain the same as in Eq. (10). As $`\zeta (\overline{x})`$ is a monotonic function of $`\overline{x}`$, we can change the independent variable in Eq. (14) from $`\overline{x}`$ to $`\zeta `$. The resulting equation for $`\varphi (\zeta )`$ is

$$4(\zeta 1)\zeta \varphi ^{\prime \prime }+2\varphi ^{}+(1\kappa \zeta ^2)\varphi =0,$$
(15)
where the primes now denote the $`\zeta `$ derivatives. The function $`\varphi (\zeta )`$ is defined on the interval $`1\zeta \zeta _1`$, where $`\zeta _1=Z_1/Z_0=\mathrm{cosh}^2\xi `$. One boundary condition is $`\varphi (\zeta =\zeta _1)=0`$ to which we may add an arbitrary normalization condition $`\varphi (\zeta =1)=1`$. An additional boundary condition, $`\varphi ^{}(\zeta =1)=(\kappa 1)/2`$, can be found from Eq. (15) itself, after substituting there $`\zeta =1`$. This eigenvalue problem includes a single parameter $`\xi `$, while $`\kappa `$ serves as the eigenvalue.

We have been unable to solve Eq. (15) analytically for a general $`\kappa `$. Still, several important asymptotics can be obtained.

#### IV.3.2 2. Zero-energy state and stability border $`f_1`$

For $`\overline{k}=0`$ Eq. (15) can be solved analytically:

$$\varphi (\zeta ,\overline{k}=0)\varphi _0(\zeta )=\sqrt{\zeta }\sqrt{\zeta 1}\text{arccosh}\sqrt{\zeta },$$
(16)
In other words, we impose a zero eigenvalue $`\overline{k}=0`$ and find the low-density stability border $`f_1=f_1(\eta )`$ from the boundary condition $`\varphi _0(\zeta _1)=0`$. We obtain an algebraic equation $`\mathrm{coth}(\xi _1)=\xi _1`$ for $`\xi _1=f_1\eta ^{1/2}`$. Its solution is $`\xi _1=1.19968\mathrm{}`$. This result explains why all marginal stability curves shown in Fig. 2 depart (almost) from the same point at the low-density side. Figure 11 compares the scaled quantity $`f_1\eta ^{1/2}`$ at different $`\eta `$, found numerically from Eq. (6), with this analytic prediction (a constant). The agreement is very good for large $`\eta `$. As $`\eta `$ goes down, $`f_1`$ increases and the dilute approximation starts to deteriorate.

Figure 12 compares the analytic result for the zero-energy eigenfunction $`\varphi _0(x)`$, given by Eqs. (16) and (12), with a numerical solution of Eq. (6) for $`\eta =10^4`$. The coordinate $`x`$ in Fig. 12 is rescaled by $`L_x`$. The analytic and numeric results are obtained for slightly different values of $`f_1`$ (see Fig. 11). One can see that the agreement is excellent.

#### IV.3.3 3. Short-wavelength limit: localization and universality

In the short-wavelength limit the system boundary $`\zeta =\zeta _1=\mathrm{cosh}^2\xi `$ can be moved to infinity. This requires a strong inequality $`\xi =f\eta ^{1/2}1`$. In this limit, the eigenvalue problem (15) does not include any parameter. The eigenvalue $`\kappa `$ should therefore be a number of order of unity, hence $`\overline{k}=A/Z_0`$, with constant $`A`$ of order of unity. The constant can be found numerically: $`A0.525`$. This simple result represents the low-density limit of the “universal” marginal stability curve, corresponding to strong localization. Figure 13 shows this asymptotics for $`\eta =10^5`$. One can see excellent agreement for large enough $`Z_0`$, but not too close to the higher-$`Z_0`$ (low-density) instability border. Near the instability border $`\xi `$ becomes of order of unity, and localization breaks down.

Returning to the parameters $`\eta `$ and $`f`$, and to the wave number $`k=k_{ph}L_x`$, we can rewrite the asymptotics $`\overline{k}=A/Z_0`$ as

$$k=\frac{A}{2}\left[f\eta ^{1/2}+\frac{1}{2}\mathrm{sinh}(2f\eta ^{1/2})\right].$$
(17)
This asymptotics is shown in Fig. 14.

The asymptotics (17) is valid when $`\xi =f\eta ^{1/2}1`$. It is easy to check that this criterion coincides, in the dilute limit, with the localization criterion $`\overline{k}^1\eta ^{1/2}`$ discussed in Sec. 4A. On the other hand, the parameter $`\xi `$ should not be too large, so that the dilute limit condition $`Z_01`$ is still satisfied, see Eq. (13). These two criteria can be rewritten as a strong double inequality for $`f`$:

$$\frac{1}{\eta ^{1/2}}f\frac{\mathrm{log}(8\eta ^{1/2})}{2\eta ^{1/2}},$$

that can be satisfied only for extremely large $`\eta `$.

#### IV.3.4 4. Long-wavelength limit: perturbation theory

Close to the low-density stability border, $`ff_1f_1=1.19968\mathrm{}\eta ^{1/2}`$ we can assume that $`\kappa 1`$ and solve Eq. (15) perturbatively. In the physical units, this strong inequality corresponds to the long-wavelength limit: $`k_{ph}L_x1`$. In its turn, the dilute limit requires $`\eta 1`$. We substitute in Eq. (15) $`\varphi (\zeta )=\varphi _0(\zeta )+\kappa \mathrm{\Omega }(\zeta )`$, where $`\varphi _0`$ is given by Eq. (16). Neglecting the $`\kappa ^2`$-term in Eq. (15), we obtain

$$4(\zeta 1)\zeta \mathrm{\Omega }^{\prime \prime }+2\mathrm{\Omega }^{}+\mathrm{\Omega }=\zeta ^2\varphi _0(\zeta ),$$
(18)
The normalization and boundary conditions at $`\zeta =1`$ are $`\mathrm{\Omega }(\zeta =1)=0`$ and $`\mathrm{\Omega }^{}(\zeta =1)=1/2`$, respectively. The latter condition follows from Eq. (18) itself. Equation (18) can be solved analytically. With the account of the two boundary conditions we obtain

$$\mathrm{\Omega }(\zeta )=\frac{1}{8}\mathrm{\Phi }_2(\zeta )+\mathrm{\Phi }_2(\zeta )I_1(\zeta )\mathrm{\Phi }_1(\zeta )I_2(\zeta ),$$
(19)
where $`\mathrm{\Phi }_1(\zeta )=(\zeta 1)^{1/2}`$ and $`\mathrm{\Phi }_2(\zeta )=2\varphi _0(\zeta )`$. $`I_1`$ and $`I_2`$ are indefinite integrals:

$$I_1=^\zeta \frac{\mathrm{\Phi }_1G_1}{G_2W}𝑑\zeta \text{and}I_2=^\zeta \frac{\mathrm{\Phi }_2G_1}{G_2W}𝑑\zeta ,$$
(20)
where $`G_1=\zeta ^2\varphi _0/4,G_2=\zeta (\zeta 1)`$, $`W=\mathrm{\Phi }_1\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}\mathrm{\Phi }_2`$ and the primes denote the $`\zeta `$ derivatives. Integrals $`I_1`$ and $`I_2`$ can be evaluated in elementary functions, but the results are too cumbersome to be presented here. The additional boundary condition \[at $`\zeta =Z_1/Z_0\zeta _1`$\] reads $`\varphi (\zeta _1)=\varphi _0(\zeta _1)+\kappa \mathrm{\Omega }(\zeta _1)=0`$ which yields

$$\kappa (\xi )=\varphi _0(\zeta _1)/\mathrm{\Omega }(\zeta _1),$$
(21)
where $`\zeta _1=\mathrm{cosh}^2\xi `$. At the low-density stability border $`f=f_1`$ we have $`\xi =\xi _1=1.19968\mathrm{}.`$ In the perturbative treatment, one should expand $`\kappa (\xi )`$ near $`\xi =\xi _1`$ up to the linear term $`\xi \xi _1`$. The zero-order term vanishes, and we obtain

$$\kappa (\xi )=\frac{\mathrm{\Phi }_2^{}(\xi _1)}{2\mathrm{\Phi }_1(\xi _1)I_2(\xi _1)}\left(\xi \xi _1\right).$$

In the physical variables we have

$$k_{ph}L_x=\frac{\eta ^{1/2}}{Z_0(\xi _1)}\left(\frac{\mathrm{\Phi }_2^{}(\xi _1)}{2\mathrm{\Phi }_1(\xi _1)I_2(\xi _1)}\right)^{1/2}\left(\xi \xi _1\right)^{1/2}.$$
(22)
As $`Z_0`$ is proportional to $`\eta ^{1/2}`$, the right hand side of Eq. (22) is actually independent of $`\eta `$. Evaluating the integral $`I_2(\xi _1)`$, we obtain $`I_2(\xi _1)=0.883381\mathrm{}`$. The final result is

$$k_{ph}L_x=2.5115\mathrm{}\left(\xi \xi _1\right)^{1/2}.$$
(23)
Alternatively, we obtain

$`\overline{k}(Z_0)`$ $`=`$ $`\left({\displaystyle \frac{\xi _1}{Z_0(\xi _1)}}\right)^{3/2}\left({\displaystyle \frac{Z_0Z_0(\xi _1)}{2I_2(\xi _1)}}\right)^{1/2}`$ (24)
$`=`$ $`0.000485\mathrm{}(Z_0(\xi _1)Z_0)^{1/2}.`$
The asymptotics (24) is depicted in Fig. 13. Close to the higher-$`Z_0`$ (low-density) stability border it shows good agreement with the marginal stability curve found numerically.

## V V. Summary and Discussion

We determined the criteria for the spontaneous symmetry-breaking instability of the laterally uniform granular cluster (strip state) in a prototypical driven granular gas. Working in the limit of nearly elastic particle collisions and low or moderate densities, we employed granular hydrodynamics with the Jenkins-Richman constitutive relations Jenkins . The instability of the strip state can be interpreted in terms of negative compressibility of the granulate in the lateral direction. An important limit is found, when the marginal stability curves are independent of the details of the boundary condition at the driving wall. In this regime the density perturbation is exponentially localized at the elastic wall opposite to the driving wall. Working in the dilute limit, we obtained some analytic asymptotics of the marginal stability curves.

The results of this work show that the symmetry-breaking instability predicted in Ref. LMS  is robust and does not require very special constitutive relations. The marginal stability curves obtained in this work are quite similar to those obtained earlier LMS  for a different set of constitutive relations (see Fig. 4). There are some quantitative differences, however. Therefore, the instability provides a sensitive test to the accuracy of constitutive relations.

This work was focused on the criteria of instability of the strip state. In systems sufficiently long in the lateral direction, instability occurs in a whole range of wave numbers $`k`$ (below the respective marginal stability curve). Correspondingly, multiple steady state solutions with different $`k`$ are possible. In a laterally infinite system, these solutions are periodic in the lateral coordinate. A finite system selects a finite number of wavelengths LMS . An important issue that was not addressed in this work is selection: what is the wavelength of the resulting symmetry-broken cluster in an infinite, or long enough, system? The selection has dynamical nature; this important issue is addressed elsewhere LMS2 .

Recently, the predicted symmetry-breaking instability has been observed in particle simulations Schwager . We hope it will be investigated in experiment, too. The experimental setting can be of the type used by Kudrolli et al. Kudrolli ; Kudrolli2 : a system of steel spheres, rolling on a smooth surface and driven by a rapidly vibrating side wall. The present work (see also Ref. LMS ) provides the region of parameters where the instability can be observed. An important issue is to eliminate the static friction between the particles and surface that occurs far enough from the driving wall. In experiment, this is achieved by slightly inclining the system, so that a very small gravity appears Kudrolli ; Kudrolli2 . As the result, the strip state moves down, toward the driving wall Kudrolli . The model problem investigated in the present work does not include gravity. We expect, however, that the symmetry-breaking instability should persist for a non-zero gravity. In fact, a similar instability has already been observed in particle simulations of a dilute two-dimensional granular bed fluidized by a rapidly vibrating bottom plate Sunthar . Under conditions of the simulations Sunthar  there was no direct, mechanical coupling between the bottom plate vibration and collective granular motions. Therefore, the vibrofluidized system, investigated in Ref. Sunthar , is similar (though not identical) to the model system driven by a thermal wall. As gravity introduces an additional scaled parameter, the phase diagram of this type of systems should be more complicated. For example, it is already known that, at some values of the scaled parameters, steady “thermal” convection (steady state of a different type) develops both in vibrofluidized systems Sunthar ; Wildman  and in systems driven by a “thermal” wall Ramirez ; He . Granular hydrodynamics will be instrumental in delineating the phase diagrams of these systems in the limit of nearly elastic collisions.

Finally, when inelasticity of the particle collisions is not small, the normal stress difference, possible lack of scale separation and non-Gaussianity in the velocity distribution may become important. The potential role of these effects in the symmetry-breaking instability should be the subject of further investigations.

## ACKNOWLEDGEMENTS

We are very grateful to P.V. Sasorov and J.P. Gollub for useful discussions. This research was supported by the Israel Science Foundation, founded by the Israel Academy of Sciences and Humanities.


text:
# V407 Vul: a direct impact accretor

## 1 Introduction

Accreting white dwarfs in binary stars fall into two groups according to the magnetic field of the white dwarf. The flow of material close to strongly magnetic white dwarfs is controlled by the magnetic field, and matter is channelled onto one or both magnetic poles. In this case, accretion energy is released as the material crashes into the white dwarf, emitting copious X-rays and optical cyclotron emission. The X-ray and cyclotron emission are modulated on the spin period of the white dwarf as the accreting poles rotate into and out of our view. Moreover, the cyclotron emission is circularly and linearly polarized. In the non-magnetic case, by contrast, the accreting material, having some initial angular momentum, forms a disc. The disc material accretes onto the white dwarf via an equatorial boundary layer where the kinetic energy of the gas is released. The resulting radiation is not significantly modulated, and there is no polarized cyclotron radiation. Pulsed X-ray emission from white dwarf binaries has therefore been regarded as a secure indication that the accretor is magnetic.

When X-ray pulsations on a period of $`9.5\mathrm{min}`$ were discovered from the star V407 Vul (=RX J1914.4+2457), they were immediately interpreted in terms of an accreting, magnetic white dwarf spinning on the same period (Motch et al. 1996). Similar spin periods are commonly seen in the “intermediate polar” class of cataclysmic variable star in which a relatively weakly-magnetized white dwarf spins faster than the binary orbit. However, in these stars one also sees other periods, related to the orbital period and the “beat” period between the spin and orbital periods. In V407 Vul, however, both X-ray and optical data show just the one period of $`9.5\mathrm{min}`$. Moreover, the X-ray flux from V407 Vul drops to zero between pulses, which is difficult to account for on an intermediate polar model. This led Cropper et al. (1998) to suggest, instead, that V407 Vul is a “polar”, in which the white dwarf’s magnetic field locks its spin to its companion star’s orbit. If so, then V407 Vul has an orbital period of only $`9.5\mathrm{min}`$, the shortest known for any binary system. Such a period implies that the donor is a helium-rich degenerate star, making V407 Vul the first magnetic member of the AM CVn stars. The AM CVn stars are a select group of eight mass-transferring binary stars (other than V407 Vul), which have periods ranging from $`17\mathrm{min}`$ (AM CVn, Nelemans, Steeghs & Groot 2001) to $`65\mathrm{min}`$ (CE 315, Ruiz et al. 2001). As for V407 Vul, their short periods imply that the donor stars are hydrogen-deficient, and indeed no hydrogen appears in their spectra. They are thought to form from initially detached double white dwarf systems, from systems with helium-star donors or from mass transfer initiated when a $`1\mathrm{M}_{}`$ donor starts to transfer mass to a white dwarf at the end of its core hydrogen burning (Nelemans et al. 2001; Podsiadlowski, Han & Rappaport 2001). Double white dwarfs that fail to become AM CVn stars are possible Type Ia supernova progenitors.

A problem with the polar model is that V407 Vul shows no optical polarization (Ramsay et al. 2000), possible on the polar model only if the white dwarf has either a very strong or relatively weak field (while remaining synchronized). In an effort to explain this, Wu et al. (2001) proposed a model in which the spin of the magnetic white is not synchronised with the orbit leading to dissipation of electric currents in the donor which produces unpolarized optical flux. However, the dissipation also leads to synchronization on a short timescale, which makes the chance of such a configuration low.

In this paper we present an alternative model for V407 Vul, in which the white dwarf need not be magnetic, but the X-rays will still be strongly modulated. V407 Vul may thus be the first example of a new class of X-ray emitting binary star. We start by describing the observational characteristics of V407 Vul that need explaining.

## 2 Observed features of V407 Vul

While there have been relatively few observations of V407 Vul, any model of the system must satisfy the following constraints:

1. The $`9.5\mathrm{min}`$ X-ray pulsations. The X-ray bright phase occupies half of the pulsation period; for the other half of the cycle, the X-ray flux is undetectable (Cropper et al. 1998). No other periodic signals are seen in X-rays.
2. The $`9.5\mathrm{min}`$ optical pulsations. Again, no other periodic signals are seen. The optical pulsations have a peak-to-peak amplitude of $`0.07`$ magnitudes. The optical flux peaks 0.4 cycles before the X-ray flux (Ramsay et al. 2000).
3. The X-ray spectrum. This is soft and can be fitted with an absorbed black-body spectrum of temperature $`40`$$`55\mathrm{eV}`$ (Motch et al. 1996; Wu et al. 2001).
4. The lack of optical polarization. Ramsay et al. (2000) measured $`0.3`$% circular polarization, a level consistent with zero given the systematic uncertainties.
5. The optical spectrum. Unfortunately this has not yet been published, but it is reported to be devoid of emission lines (Wu et al. 2001). Without having seen it, it is hard to judge the importance of this, but we will consider it as another possible constraint to satisfy.
6. The distance. V407 Vul is heavily absorbed from which Ramsay et al. (2000) obtain $`d>100\text{pc}`$. They further find that $`d<400\text{pc}`$ from a de-reddened $`I`$-band magnitude of $`15.5`$, although the “reddening” is deduced from the X-ray column, and may be too high since the colours end up being too blue even for a Rayleigh-Jeans spectrum. Ramsay et al. (2000) used $`A_V=5.6`$, but also state that colours deduced from $`A_V=4`$ do fit black-bodies, from which we estimate that the de-reddened $`I`$-magnitude lies in the range $`15.5<I<16.1`$.

It was the singly periodic signal and the on/off nature of the X-ray light-curve that led Cropper et al. (1998) to suggest their polar model. The absence of both polarization and optical line emission presented difficulties that led Wu et al. (2001) to develop their unipolar inductor model by analogy with the Jupiter-Io system. As outlined in the introduction, in their model, the magnetic white dwarf (the accretor in our model) is slightly asynchronous with the binary orbit and the resulting electric field drives currents that run between the two stars, leading to energy dissipation in both of them. Dissipation on the magnetic white dwarf powers the X-ray emission, while dissipation at the donor plus irradiation powers the optical flux. The irradiation is predicted to be comparable to the ohmic dissipation, which they suggest can explain the relatively weak optical modulation and the lack of line emission. Since the non-magnetic star dominates the optical emission, the absence of polarization is also explained.

The major problem with Wu et al.’s model is that it is short-lived: they estimate that it will last only $`1000\text{yr}`$. Even if all AM CVn systems pass through this stage, there would only be 1–10 such systems in our Galaxy at any one time, according to the formation rates of Nelemans et al. (2001) and Podsiadlowski et al. (2001). The chance of finding one within our neighbourhood is therefore small. Magnetic systems may in fact comprise only a small fraction of the total, making the probability of finding such a system very low. It is therefore worth searching for longer-lived models.

## 3 A new model for V407 Vul

Our idea is simple: in very close binary systems, the mass transfer stream can plough straight into the accretor even in the absence of the magnetic field. This happens if the minimum distance of the ballistic gas stream from the centre of mass of the accretor is smaller than the accretor’s radius, which is most famously the case in Algol binary stars, where the accretors are main-sequence stars. We propose that V407 Vul is the first instance of Algol-like direct impact in the case of a white dwarf accretor.

In Fig. 1 we show that the direct impact is possible using typical system parameters (section 4). The figure shows the ballistic path of the mass transfer stream calculated for an accretor mass of $`M_1=0.5\mathrm{M}_{}`$ and a donor mass of $`M_2=0.1\mathrm{M}_{}`$. The separation of the binary is fixed through Kepler’s laws by the masses and the orbital period, $`P=9.5\mathrm{min}`$, while we have used Nauenberg’s (1972) analytic formula for the radius of the white dwarf. If the $`9.5\mathrm{min}`$ pulsations truly reflect V407 Vul’s orbital period, direct stream impact is possible.

This model immediately explains the absence of polarization and the X-ray pulsations. There is no polarization because the white dwarf is non-magnetic, while the impact point is fixed in the rotating frame of the binary and will produce one pulse per orbit. Our model nicely explains the absence of X-ray flux for half the cycle since the impact will be on the equator of the accretor; this is a key problem with intermediate polar models of V407 Vul, as Cropper et al. (1998) point out.

There are two possible explanations for the optical pulsations. First, we do not expect the white dwarf to be synchronized with the orbit, but instead to spin rapidly as the result of the accretion of material of high specific angular momentum. The surface of the white dwarf will therefore move under the point of impact, and one expects a heated area to trail downstream from the spot. If this area is responsible for the optical pulsations, they would then peak at an earlier phase than the X-rays, as observed. The amplitude of the optical pulsations in this case depends upon the extent of the heated region and the brightness of the unheated photosphere of the white dwarf. In order to dominate the optical flux, given the lower limit on the distance of $`100\text{pc}`$ (Ramsay et al. 2000) and the unreddened $`I`$ magnitude between $`15.5`$ and $`16.1`$ (section 2), we require $`M_I<10.5`$$`11.1`$. For a $`0.5\mathrm{M}_{}`$ white dwarf this implies an effective temperature $`T_{\mathrm{eff}}>15`$,$`000`$$`21`$,$`000\text{K}`$ (Bergeron, Wesemael & Beauchamp 1995), which is nothing out of the ordinary. Indeed, given an accretion rate of $`10^8\mathrm{M}_{}\text{yr}^1`$ (section 5), we would expect a much higher temperature from compressional heating alone: Townsley & Bildsten (2002) give a surface temperature of $`27`$,$`000\text{K}`$ from compressional heating for a $`0.6\mathrm{M}_{}`$ white dwarf accreting at $`10^9\mathrm{M}_{}\text{yr}^1`$.

Alternatively, the optical pulsations may come from the heated face of the donor star. This could be heated by irradiation from the impact spot (but see the next section), or by photospheric emission from the white dwarf. Although the latter is normally ignored, once again the very close orbit of V407 Vul means that it is not just possible, but likely. The geometry of the impact site means that the heated face of the donor star will also produce a peak optical flux in advance of the X-ray flux, as observed, although the impact site would have to be on the side facing the donor to explain the $`0.4`$ cycle shift seen. The temperature of the heated face of the donor due only to photospheric emission from the accretor is $`T_\mathrm{h}(R_1/a)^{1/2}T_10.33T_1`$, where $`R_1`$ and $`T_1`$ are the radius and temperature of the accretor. Given that the donor is about 3 times the size of the accretor, and the temperature of the accretor estimated above, the donor could easily produce the optical pulsations.

Of the observational facts above, we are left to explain the reported lack of optical line emission and the softness of the X-ray spectrum.

### 3.1 Optical Line Emission

The direct-impact model provides a beautiful way to avoid X-ray irradiation of the donor star, and therefore any associated emission lines, altogether, because it is possible for the impact point to be hidden from the donor star. Fig. 1 illustrates just such a case; we will examine the parameter constraints needed for this to be the case in the next section. Nevertheless, it is not clear whether it is really necessary to exclude irradiation of the donor because any emission lines would execute high-amplitude sinusoidal motion on a period of only $`9.5\mathrm{min}`$ and would be smeared in wavelength, especially given that V407 Vul is optically faint (the orbital velocity of the donor in Fig. 1 is $`800\mathrm{km}\mathrm{s}^1`$). Given the hydrogen deficiency and the proximity of the two stars, one would only expect to see ionized helium emission. Since V407 Vul is only easily observable in the $`I`$-band, it is not clear whether any emission lines should have been detected, even if irradiation is taking place. There is, of course, some irradiation of the stream, but only where it nears the impact site. By this time it is moving very fast, and will be highly ionized, and since emission from this part is hard to detect even in much brighter and longer period polars, we do not regard this as a problem.

### 3.2 The Soft X-ray Spectrum

Direct impact accretion differs significantly from magnetically-confined accretion. First, the accretion stream will be much narrower in the direct impact case. The reason is that in the magnetic case, threading occurs over a range of radii leading to impact in a relatively large arc close to one or both poles of the white dwarf. On the other hand, as Lubow & Shu (1976) and Lubow (1989) showed, as the ballistic stream nears its closest approach to the centre of mass of the accretor, it becomes comparable in width to the hydrostatic scale height that an accretion disc would have at the same radius. At the white dwarf this is a width of order $`(kTR_1^3/GM_1m_p)^{1/2}5\times 10^3R_1`$ for $`T=20`$,$`000\text{K}`$ (stream from the heated face of donor) and $`M_1=0.5\mathrm{M}_{}`$. Therefore, accretion will take place over a fraction $`f0.01`$% of the white dwarf’s surface.

Together with an accretion rate of $`\dot{M}10^8\mathrm{M}_{}\text{yr}^1`$ (section 5), some 10–100 times higher than in polars, the narrow stream width makes the accretion in V407 Vul comparable to the blob accretion model (Kuijpers & Pringle 1982; Frank, King & Lasota 1988; King 2000), in which the soft X-ray excess of polars is explained by the stream breaking into dense blobs which are able to penetrate below the photosphere and therefore become thermalized, giving soft X-ray emission. This happens for blob densities $`\rho >10^7\mathrm{g}\mathrm{cm}^3`$ (Frank et al. 1988). For a stream width of $`10^4\mathrm{R}_{}`$, an accretion rate of $`\dot{M}=10^8\mathrm{M}_{}\text{yr}^1`$, and a free-fall velocity $`v_{ff}=4\times 10^8\mathrm{cm}\mathrm{s}^1`$, we find $`\rho 3\times 10^5\mathrm{g}\mathrm{cm}^3`$, well in excess of Frank et al.’s limit.

In magnetically-confined accretion, after the accretion shock, the material continues to flow along the field lines and can only come to a halt through cyclotron and bremsstrahlung cooling. This leads to the shock forming some height above the point at which the ram pressure of the stream matches the pressure in the white dwarf’s atmosphere. In the direct impact case, there is nothing to stop the stream expanding sideways following the shock which will cause adiabatic cooling and lead to the shock forming much closer to the point where the pressures match. This increases the chance that the shock will be buried. The final, and perhaps most important, difference is that the white dwarf moves rapidly beneath the impact site, which will sweep the heated material downstream. We suggest that it is this which spreads the emission over a large enough area to cause the spectrum to be very soft as observed. We can crudely estimate the cooling time as the thermal timescale of a column of gas heated to the observed temperature of $`T50\mathrm{eV}`$. The column density $`\mathrm{\Sigma }`$ is fixed by the penetration depth according to

$$\mathrm{\Sigma }g=P_{\mathrm{ram}}=\rho v_{\mathrm{ff}}^2,$$
(1)
where $`g=GM_1/R_1^2`$ is the surface gravity, $`\rho `$ is the stream density and $`v_{\mathrm{ff}}`$ is the free-fall velocity ($`v_{\mathrm{ff}}^2=2GM_1/R_1`$), which is comparable to that of the stream. Given an energy of $`3kT/2`$ per particle, and a loss rate of $`\sigma T^4`$, we obtain a cooling time $`t_c`$, in seconds, of

$$t_c1.0\left(\frac{\rho }{3\times 10^5\mathrm{g}\mathrm{cm}^3}\right)\left(\frac{R_1}{0.01\mathrm{R}_{}}\right)\left(\frac{T}{50\mathrm{eV}}\right)^3$$
(2)
Since the white dwarf spin period may be as short as 10 seconds, this cooling time is the right order of magnitude to lead to significant spreading of the emission. All that is needed is that the area is spread so that $`f0.1`$%, for the energy liberated by accretion at the rate $`10^8\mathrm{M}_{}\text{yr}^1`$ to be radiated away with a temperature of $`T50\mathrm{eV}`$, as observed.

## 4 Parameter Constraints

We now calculate the limits upon the parameters required for our model to be viable, fixing the orbital period at $`P=9.5\mathrm{min}`$. We have first the orbital separation in terms of the total mass of the binary

$$a=0.148\left(\frac{M_1+M_2}{\mathrm{M}_{}}\right)^{1/3}\mathrm{R}_{}.$$
(3)
For the radii of the white dwarfs we use Nauenberg’s 1972 formula

$$\frac{R}{\mathrm{R}_{}}=0.0112\left[\left(\frac{M}{1.433\mathrm{M}_{}}\right)^{2/3}\left(\frac{M}{1.433\mathrm{M}_{}}\right)^{2/3}\right]^{1/2}.$$
(4)
This is a lower limit, since Nauenberg’s relation is based upon zero temperature and no rotation, and departures from these assumptions act to increase the radius for a given mass. The radius of closest approach to the centre of the white dwarf relative to the separation of the binary is a function of mass ratio, for which Nelemans et al. (2001) give a fit based upon the calculations of Lubow & Shu (1975), and confirmed by our own integrations.

Next, we have the well-known relation between the mean density of the Roche lobe filling donor and the orbital period, which gives $`\overline{\rho }=4.3\times 10^3\mathrm{g}\mathrm{cm}^3`$ in the case of V407 Vul. Using Nauenberg’s formula, this corresponds to a mass of $`M_2=0.077\mathrm{M}_{}`$. More generally this is a lower limit because any hydrogen content or semi-degeneracy reduces the density for a given mass, and since the density increases monotonically with mass, the mass must increase to match a fixed density.

Two other constraints come from considerations of the stability of mass transfer which set upper limits to the mass ratio. To ensure that the Roche lobe does not shrink faster than the donor in the case of conservative mass transfer, the mass ratio must satisfy

$$q<\frac{5}{6}+\frac{1}{2}\frac{d\mathrm{ln}R_2}{d\mathrm{ln}M_2},$$

e.g. Nelemans et al. (2001). Nauenberg’s formula gives

$$\frac{d\mathrm{ln}R_2}{d\mathrm{ln}M_2}=\frac{1}{3}\left(\frac{1+x^{4/3}}{1x^{4/3}}\right),$$

where $`x=M_1/1.433`$. This strictly only applies to the limiting case of $`M_2=0.077\mathrm{M}_{}`$, but we will take it to apply to higher mass cases too. This is a conservative assumption, because radiative stars shrink on mass loss. A stricter constraint applies if accretion occurs directly onto the white dwarf and the angular momentum of the stream is *not* transferred back from the white dwarf to the orbit (i.e. tidal and magnetic coupling are ineffective). The limit then becomes

$$q<\frac{5}{6}+\frac{1}{2}\frac{d\mathrm{ln}R_2}{d\mathrm{ln}M_2}\sqrt{(1+q)r_h},$$
(5)
(Nelemans et al. 2001), where $`r_h`$ is the circularization radius as a fraction of $`a`$, which is a function of mass ratio for which we use the approximation of Verbunt & Rappaport (1988), accounting for their inverted definition of $`q`$.

Finally, if the reported absence of emission lines is significant (see section 3 for a discussion of this), there is another upper limit on the mass of the donor for a given accretor mass, required in order for the impact spot to be hidden from the donor star. We carried out our own integrations in the Roche potential to derive this upper limit as a function of accretor mass.

These constraints give the allowed parameter space indicated in Fig. 2 which shows that direct stream impact in a stable binary is possible if the masses of the two stars lie in either the large triangle ABC for the standard stability limit, or in the smaller triangle A’BC’ if the more rigorous stability criterion, Eq. 5, applies. This smaller region is also consistent with the requirement that the donor not be irradiated (curved dashed line in Fig. 2). Moreover, this happens for $`M_10.5\mathrm{M}_{}`$, a typical mass for a white dwarf. Thus this model is feasible without the need for fine tuning of the parameters. If the donor is degenerate, then, for no irradiation, the parameter space is restricted to $`M_2=0.077\mathrm{M}_{}`$, $`M_1=0.42`$$`0.49\mathrm{M}_{}`$.

## 5 Future Evolution

V407 Vul’s destiny is to become a standard member of the AM CVn systems. As the donor’s mass decreases, even if the accretor fails to increase in mass, the increase in the orbital period and the decrease in mass ratio, both act to make it less likely that the stream will hit the accretor directly. Eventually, the stream will orbit the accretor and a disc will form. Depending upon the precise parameters, it will do so after the loss of $`0.01`$$`0.1\mathrm{M}_{}`$ from the donor. The mass transfer rate for an orbital period of $`9.5\mathrm{min}`$ is $`10^8\mathrm{M}_{}\text{yr}^1`$ (Nelemans et al. 2001), so this phase will last $`10^6`$$`10^7\text{yr}`$. If all AM CVn systems pass through such a phase, then, according to the formation rates of Nelemans et al. (2001) and Podsiadlowski et al. (2001), there should be $`10^3`$$`10^5`$ systems like V407 Vul in the Galaxy. The nearest of these should be between $`100`$ and $`400\text{pc}`$ away, consistent with estimates of the distance to V407 Vul. This agreement favours the detatched double white dwarf route for the formation of AM CVn’s over either the helium star or the semi-degenerate routes, since the latter rarely reach orbital periods as short as 10 minutes (Nelemans et al. 2001; Podsiadlowski et al. 2001), whereas Nelemans et al. (2001) predict that most double white dwarfs will pass through a direct-impact phase if they avoid merging.

## 6 Gravitational Waves

V407 Vul is a very promising source of gravitational waves. Following Meliani, de Araujo & Aguiar (2000), who calculate strain amplitudes for many cataclysmic variable stars, and using $`M_1=0.5\mathrm{M}_{}`$ and $`M_2=0.1\mathrm{M}_{}`$, V407 Vul will produce a strain amplitude $`h=1.2\times 10^{21}(d/100\text{pc})^1`$ at Earth. This is comparable to the largest values listed by Meliani et al. (2000). In addition, the very short period of V407 Vul will lead to gravitational wave emission at a frequency of $`f_{\mathrm{gw}}=2/P=3.5\times 10^3\mathrm{Hz}`$, where the background noise of an instrument such as LISA is predicted to be much lower than at the lower frequencies characteristic of other cataclysmic variable stars. V407 Vul is thus one of the brightest prospects for space-based gravitational wave detection.

## 7 Conclusions

We have shown that plausible masses for the two components of V407 Vul suggest that the mass transfer stream in this system strikes the accreting white dwarf directly. The resulting spot explains the X-ray pulses observed from this system, without the need for a magnetic white dwarf. This is consistent with the lack of polarization from the system, and with the complete disappearance of X-ray flux in between pulses. Our model has a lifetime of $`10^6`$ to $`10^7\text{yr}`$ compared with the $`10^3\text{yr}`$ of Wu et al.’s unipolar inductor model, making it easier to reconcile with estimated formation rates for the AM CVn systems, as long as a substantial fraction of these systems pass through a direct impact-phase, as predicted by Nelemans et al. (2001) on the basis of a double white dwarf origin for these systems. V407 Vul is the first example of a new class of stream-fed, non-magnetic white dwarf accretors.


text:
# Degradation of LaMnO3-y surface layer in LaMnO3-y/ metal interface

## Acknowledgement

This work was supported by the Slovak Grant Agency for Science (Grants No.2/7199/20 and 1/7072/20) and in part by the European Commission (project GRT-CT-2000-05001 MULTIMETOX)