The discrete self-trapping equation

Physica 16D (1985) 318-338 North-Holland, Amsterdam

THE DISCRETE SELF-TRAPPING EQUATION J.C. EILBECK'~, P.S. LOMDAHL and A.C. SCOTT Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, N M 87545, USA

Received 15 October 1984

A simple system of ordinary differential equations is introduced which has applications to the dynamics of small molecules, molecular crystals, self-trapping in amorphous semiconductors, and globular proteins. Analytical, numerical and perturbation methods are used to study the properties of stationary solutions. General solution trajectories can be either sinusoidal, periodic, quasiperiodic or chaotic.

1. Introduction Self-trapping is now widely recognized as a generic mechanism in the physics of condensed matter. In broad terms self-trapping involves the complementary interaction of two cause and effect relationships: i) the localization of some conserved quantity influences its surrounding structure; and ii) this structural change reduces the energy of the conserved quantity preventing its dispersion. In this way one obtains a localized "lump" of the conserved quantity that is dynamically stable and often mobile. The concept of self-trapping was introduced a half century ago in a note by Landau on the motion of an electron in a crystal lattice .[1]. He suggested that an effect of a localized electron would be to polarize the crystal which, in turn, would lower its energy. Landau's suggestion was discussed in detail by Pekar [2] (who seems to have coined the term "polaron" for the localized electron plus lattice distortion), by Frthlich [3], and by Holstein [4]. Since 1970 the polaron has been studied by a number of authors [5]. Other examples of self-trapping involve electromagnetic en~'Permanent address: Department of Mathematics, HeriotWatt University, Riccarton, Edinburgh EH14 4AS, U.K.

ergy in a plasma [6], fight energy in an optical fiber [7], and hydrodynamic energy in a wave tank [8]. Our primary motivation for undertaking a systematic study of self-trapping has been to generalize the theory of solitons on the alpha-helix structure of protein which has been proposed by Davydov [9]. In 1973 he suggested that the amide-I (CO stretching) vibrations of pcptide groups located along the alpha-helix might become selftrapped through interaction with low frequency phonons. Numerical studies indicate that the physical parameters of natural protein are appropriate to support such a self-trapped state [10], and experimental confirmation has been obtained from infrared absorption and Raman scattering experiments on crystals of the synthetic polypeptide acetanilide [11]. A qualitatively similar concept which has evolved from biochemistry is referred to as a "conformon". As proposed by Volkenstein [12], the conformon is essentially a polaron; while the conformon proposed by Green and Ji [13] is quite similar to Davydov's soliton. When one deals with natural protein, it is desirable to avoid the simplifying assumption that self-trapped vibrational energy is propagating as a soliton along the regular structure of an alphahelix. The real context is the atomic structure of a particular globular protein as determined by X-ray

0167-2789/85/$03.30 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

J.c. Eilbeck et at/The discreteself-trapping equation

319

diffraction [14]. In this context one is interested in following the temporal evolution of probability amplitudes for vibrational quanta on each amino acid (or peptide group) for a particular protein. This evolution can be described by the following set of ordinary differential equations which we call the discrete self-trapping (DST) equation [15]:

convenient to set the largest component(s) of M equal to unity and carry the quotient as the parameter e. With 3, = 0, (1.1) reduces to the linear system

(id - too )A + 3,D(IAI2)A +

which has been studied in detail by Anderson [16]. His motivation was to describe the diffusion of conduction electrons (or holes) in amorphous semiconductors and the diffusion of spin states in spin glasses; thus the order (n) was considered to be a large number. The main result of these studies is that "localized states" occur if the probabilistic spread of the diagonal components of M is greater than the largest off-diagonal components. Our interest here is to study localization that arises from nonlinear interactions; so we remove "Anderson localization" by assuming all the diagonal elements of M to be zero. Then the unitary transformation A ~ A exp(-ito0t ) reduces (1.1) to the standard form of the DST equation

eMA

=

0.

(1.1)

Here

(1.2)

A = col(A t, A 2 . . . . . A . )

is a complex n-component vector each component of which represents the probability amplitude of finding some conserved quantity on the n th subunit of the structure. The diagonal matrix

B(IAI2)~_lA112IA212

(1.3) Ia,I 2

i-~-

)

~oo A + e M A = 0

iA" + 3,O(I,,112),,I + e M A = 0.

which appears in the second term of (1.1) represents the tendency of A to self-trap through a nonlinear interaction with the adjacent structure. The strength of this interaction is specified by the parameter 3, which is proportional to the square of the coupling between the conserved quantity and the structure divided by a structural spring constant [15]. For self-trapping of vibrational energy on a protein, the off-diagonal components of the n × n matrix

M = [ mjk ]

(d

(1.4)

are dipole-dipole interaction energies [14, 15]. For self-trapping of electric charge in a (possibly disordered) crystal, these components are energies of electronic overlap integrals [4]. In either case we expect mjk = mkj SO it is reasonable to assume the matrix M to be real and symmetric. This assumption will be made throughout the paper. It is

(1.5)

(1.6)

Recently it has been suggested that the density of states for a linear eigenvalue spectrum of a particular protein can be estimated from the fractal dimension of that protein [17], and the term fracton has been coined for those eigenvectors with length scales sufficiently short to be insensitive to the boundaries [18]. In this paper we avoid such considerations by noting that complete knowledge of the linear eigenvalue spectrum is readily obtained directly from the structure of M. In addition to the above mentioned motivations for studying the DST equation, we have recently become aware that this equation describes the nonlinear dynamics of small polyatomic molecules such as water, ammonia, methane, acetylene and benzene [19]. In these applications, nonlinearity enters to describe the "softness" of XH stretching vibrations, and (1.1) correctly describes the strong overtone series observed for such molecules [20]

J.C. Eilbeck et al./ The discrete self-trapping equation

320

and the transition to local mode behavior as the total energy is increased [21]. In a related development Collins has recently suggested that solitons may propagate around the ring of benzene [22]. In the next section are presented some general properties of (1.6) including conservation laws. Analytical and numerical calculations of stationary (or single frequency) solutions are presented in section 3 for n = 1, 2, 3 and 4 paying particular attention to bifurcations with respect to the parameters 3, and e. Our motivation for studying small values of n is to obtain precise results that can be generalized to systems of significantly larger order, but the above mentioned application of (1.1) as a model for small polyatomic molecules makes the results for low order systems intrinsically interesting. In section 4 we present a perturbation expansion for 13'1<< I~l and we count the number of real stationary solutions for lel << 13'1. The behavior of periodic, quasiperiodic and chaotic solutions is described in section 5. It is shown here that solution trajectories can be chaotic for n > 3 but not for n < 2. A final section lists our conclusions and a linear stability theory for stationary solutions is presented in an appendix.

2. General properties

The DST equation (1.6) can be derived from the Lagrangian

we note that both /:/= 0

(2.4)

and ~r = 0

(2.5)

when the { Aj } are solutions of DST. Stationary solutions of DST are defined as having the time dependence A ( t ) = , e x p (itot),

(2.6)

where ¢ is independent of time. Such solutions are of particular interest when One seeks the spectrum of absorption from an harmonic field. In this case we note that the harmonic frequency is not proportional to the energy as it would be if DST were linear. The relation is

o: = - d H / d N .

(2.7)

This follows if one notes that (1.6) can be written as Hamilton's equation

iA 1 = BH/OA?

(2.8)

and considers H to be evaluated as a function of N on a family of stationary solutions. Then (2.8) becomes - toA2 = ( d H / d N ) ( ON/OAy) which, in turn, implies (2.7). For a detailed analysis of (1.6) it is sometimes convenient to write A in terms of the magnitudes and phases of its components, thus

L = E [½ i( A T A j - AjA~ ) + ½v[AsI" ] J

A -- col (axe i01, az ei02..... aneW"),

+ e )-'. mjkA~.A k,

(2.9)

(2.1) where the aj are assumed to be real. Then (1.6) requires

j,k

thus the Hamiltonian (energy) is H = - ½ v ~ I A j l 4 - e ~ mjiATa ,. j

k--1

and

Defining the number as N = E[As.I 2 J

(2.10a)

(2.2)

j,k

(¢

(2.3) k-1

(2 .lOb)

J. C Eilbeck et al./ The discrete self-trapping equation

3. Stationary solutions In most cases no analytic forms for solutions of the DST equation (1.6) are known and even a qualitative classification of the numerical solutions is difficult. However there is an important subclass of solutions that can easily be classified and often expressed in simple form. These are the stationary solutions defined in (2.6) where #, is an n-vector that is time independent. Substituting (2.6) into (1.6), we get a nonlinear, algebraic, eigenvalue problem for to and ~, - t o ~ + -tD(l~lZ)~ + where

D(I,I 2) ffi

eM~:

0,

(3.1)

[ 01 1~212

.

(3.2)

I%12

0

When ~, = 0, (3.1) reduces to the linear eigenvalue problem eM~ = to# and the eigenvalues to and eigenvectors ~ can easily be calculated. In the case e = 0, we have ~ito ='Yl~bil2t~i" S i n c e to must be constant, each ~i must satisfy ~ ffi 0 or I,/,il2 ffi to/~,,. If ~ is real, this reduces to ~ = 0, + ~/w/T • Note the scale invariance: if ~ is a solution of (3.1) then so is ~ e x p ( i a ) for any a. In between these two limiting cases "t ffi 0 or e = 0 (or alternatively e-o oo or ~, -o oo) it is often possible to find analytic solutions of (3.1) connecting the solutions in each

321

limit. When analytic solutions are not known, numerical path continuation methods [15, 23] enable solution curves to be generated efficiently. The linear stability of each stationary solution can be calculated from the eigenvalues of a certain linear system (see [24] and appendix). This calculation gives either necessary conditions for stability or sufficient conditions for instability. For convenience we refer to solutions satisfying the necessary conditions for stability as stable: however such solutions may be unstable due to nonlinear instabilities or resonance conditions. In many cases the stability of such solutions has been checked by numerical integration. Since there are many solutions to describe, especially for the larger values of n, it is convenient to adopt a shorthand description for each solution branch. If we consider only the case cb real (subject to a gauge transformation), in the limit ~, ~ oo each ~,~ --, 0 or + vrN-/K, where K is the number of nonzero 'hi- To represent these three different limits we use the s y m b o l s . , 1', and ~ respectively. Thus a solution with n ffi 3, and ~1 = - g ' 2 N1/-ff~, ~3 --' 0, as 3' --> oo, would be written as (1" J, • ). Since all the cases we consider below are invariant under a cyclic permutation and the gauge invariance discussed above, we can without loss of generality arrange the ordering so that ~1 is in the 1' state. In addition, those ~ values which are identically zero for the whole branch will be denoted by 0, whereas complex ~ values will be denoted by *.

Table I Stationary solutions of (3.1) for the case n ffi 2 Mode description

(1'~) (1' -)

Mode amplitudes

~0

H

,/h ~. ,~ ffiv/N-/2

"~N+~

- ~3'N 2 -~N

~"/N - e

- ~gN 2 +eN

qh = -%ffil/~-/2

J.C. Eilbeck et al. / The discrete self-trapping equation

322

1)

The case n = 1

In this case there is only one solution (T) of (3.1) since M = 0. This is = v~-,

~ = N~,.

(3.3)

This stable solution harmonic oscillator.

2)

corresponds

to a single

The case n = 2

The solution curves for these modes are shown in fig. 1. At y = 0 there are two branches, (T J,) and (T T ) corresponding to the two eigenvalues of M. At 7 = 2to/e the (1' T) branch bifurcates to give a doubly degenerate (T • ) branch. The (T $ ) and (1' • ) branches are stable, whereas the ($1' ) branch is stable below "t= 2o~/e and unstable above.

3)

In view of our restrictions on M, the only case to consider is

M=[ O1 01]"

The case n = 3

For the remainder of this section we limit our consideration to the special case

(3.4) N = 1

A straightforward analysis of (3.1) leads to the family of stationary solutions displayed in table I. These are the only solutions up to the gauge invariance: ¢ ~ Cexp(ia). Note that (2.7) can be used to check the expressions for o~ and H.

(3.5)

which implies EI¢~I 2 = 1. With this choice of normalization, the scaling ~k= "tl/2#P is convenient for numerical calculations; also e can be scaled out by a change of the time variable, so that (3.1) can be

t.O Z

0

-2

0

2

4

6

8

Fig. 1. The relations between 7, ~, e and N for stationary solutions of (1.6) with n = 2 and M defined in (3.4). Stable solutions are indicated by full lines and unstable solutions by dotted lines. The large dot is a bifurcation point. The mode symbols [e.g ( T T ), ( T $ ), etc.] are explained in the text.

J.C. Eilbeck et aL / The discrete self-trapping equation

burgh). This is

solved in the form -w~k + O(l~12)~ + M~, = 0.

(3.6)

to+l=

Once @ has been calculated, 7 can be found from the relationship 7 = Eilffil 2. Consider the interaction matrix

M =

[ya ] 0

323

3f3- s i n ( 0 + vr/6) sin 30 ' 2 (-d + 1 sin(0 - rr/3),

2 Lk2 = $3 = ~ -vr~ + 1 sin0, •

(3.7)

1

2 2 7 = t~ 2 "t- ~ 2 "k t~3 ,

Physically this corresponds to all three sites interacting with equal strengths (e.g. an equilateral triangle configuration as for the N H stretch vibrations of ammonia). The system is invariant to any site p e r m u t a t i o n so the distinct asymptotic solutions are ( 1' • • ), (1' 1' • ), (1' $ 0), (1' 1' 1'), and (1' 1' $). In addition there is one, and only one c o m p l e x solution. T h e b r a n c h (1' 1' 1') has the analytical solution ~bl = ~ b 2=~b 3 = l / v ~ - , to=7/3+2e, H=-7/6 - 2 e . T h e b r a n c h (1' J,0) has the solution ~1 =

using the scaling of (3.6). For 0 < 0 < ~t/6 this gives the b r a n c h 1' • • (to --* oo as 0 -* 0, to -- 7 / 2 for 0 - - ~ t / 6 ) . T h e branch 1' 1`. is obtained for ~t/6 < 0 < ~t/3 (to --, oo as 0 --, or~3), and the b r a n c h $1' 1' is given for 2~r/3 < 0 < 5~r/6 (to -~ oo as 0 --, 2~r/3, to = - 1 for 0 = 5 u / 6 ) . T h e complex branch has the analytic solution ep1 = 1 / v ~ , ~2 = ~ ' = 1/v~- exp(2~ri/3), to -7/3-e, H=-7/6+e. To see that this is the only possible complex solution, consider the equations (3.1) for this choice of M and n = 3,

--t~2 = 1//v/~'

- w~x + 7[dP112d?l + e(qb2 + dp3) = 0,

(3.8a)

- to% + r1%12,/,2 + e(¢1 + % ) = o,

(3.8b)

-- t o ~ 3 + 7[~312~b3 -F e(dpl + ~b2) ---- 0 .

(3.8c)

¢~3 --~ 0, to = 7 / 2

-- e, H = - 7 / 4

+ e.

T h e other three real branches have a remarkable analytic solution in parametric form (private communication, D.C. Heggie, University of Edin-

/./

/s

Y

/"

///~°

~0.

,,¢.

o

-2

0

2

I

I

4

6

•

I

8

63 Fig. 2. N o t a t i o n as in fig. 1 for n = 3, M defined in' (3.7), N = 1, a n d e = 1.

10

324

J. C Eilbeck et aL / The discrete self-trapping equation

W i t h o u t loss of generality we can take ¢1 to be real. F r o m (3.8a) this means that I m ( ¢ 2 ) = - I m ( ¢ 3 ) . Multiplying (3.8b) by ¢~, (3.8c) by ¢~', and subtracting we get an equation whose imaginary part gives x 2 = x 2, where x 2 = Re(C2) and x 3 = Re(C3). It is easy to show t h a t the choice x 2 = - x 3 leads to x 1 -- 0, so this solution can be t r a n s f o r m e d to one with ¢2 and 03 real. The other choice, x 2 = x 3, gives ¢ 2 = ¢ ~ '. Examining the i m a g i n a r y p a r t of eq. (3.8b) multiplied by ¢~ gives x 2 = - ¢ 1 / 2 . Finally, subtracting (3.8c) from (3.8b) and then subtracting (3.8a) gives I m ( ¢ 2 ) = + 3 ~ 1 / 2 , which is the solution already given. Solution curves for this case are displayed in fig. 2, with unstable solutions shown by dotted lines. T h e r e is a bifurcation point at ( t o , " / ) = (7e/2, 9e/2). T h e ( 1' • • ) branch becomes stable at the point at which 0"//0to = 0. The ( 1' $ 0) branch destabilizes at to/e = 3.5385 .... (where the r.h.s, is the only root of x 3 - 8 x - 16 = 0) and stabilizes again at to/e = 8 (the two branches whose projections cross at this point do not touch). Stability calculations on the complex solution show that this b r a n c h is stable for all " / < 40, the limit of the numerical calculations.

4) The case n = 4 W e consider only two cases. The first corres p o n d s to four sites interacting with equal strengths, e.g. a 3-D regular tetrahedral structure such as methane. The second case corresponds to the four sites on a ring with only nearest neighbor interactions of equal strength. i. Equal interactions. Here the dispersion matrix is

IVl =

°111] 1 1 1

0 1 1

1 0 1

"

(3.9)

are two distinct complex solutions: ($ * $ *) for which Cj = ½ exp(0'*r/2); j = 1,2,3,4; to = V / 4 e; a n d ( 1' * * 0) for which Cj = ( 1 / v ~ - ) exp (2/j'1r/3); j = 1, 2, 3; Ca = 0; to = "//3 - e. N o other solution b r a n c h e s are known. Several of the real branches have analytic solutions;

(1" 1' I' T), ¢ 1 = ¢ 2 = ¢ 3 = ¢ 4 = 1 / 2

(1' ~ 1' ~), ¢1 = - ¢ 2 = ¢ 3

(1"1'--),

= -¢4,

to = "//4 - e;

( 1" ,~00),

¢1 = - - ¢ 2 = 1 / ~ - ,

('~ T " " ),

to = "//2 - e; Cx = 02 = ½[1 + (1 - 64~2/"/2)1/2] 1/2, ¢2 = 04 = ½[1 - (1 - 64e2/"/2)1/2] 1/2,

¢3=04=0,

to = "//2 + e; (1' 1' • ,~),

¢1 = 0 2 =

0,=

7

(1 >/", g-

3"/]

-

¢4r--__.-- 21£[' (13/ 2- -~ -

+

+(1

,

"}-

to = ( " / - 2e)/3. There are bifurcation points on the ( 1' 1' 1' 1' ) curve at " / = 8e, on the (1' $00) curve at " / = 2e, and on the ( $1' • $ ) curve at " / = 7e/2. A full picture o f the solution curves and bifurcation points with stability information for the real solutions is given in fig. 3. The solution p a t h s in the region - 1 < to < + 1 are shown in detail in fig. 4. N u m e r i c a l linear stability calculations show that the (1' $00) b r a n c h is stable for to >_. 10.8 and 0 < to < 4.34, the ( 1' • • • ) branch is stable for to > 4.98. Tests on the stability of the complex solutions suggest that the (1' * $ *) branch is always stable, whereas the ($ * * 0) branch is unstable for - e < to < 4.16e a n d stable above this region. ii. Nearest neighbor interactions. In this case the dispersion matrix is

This system is invariant to any permutation of sites, so the distinct asymptotic solutions are

(1"'"),

,

to = "//4 + 3e;

(1"~00), (1'1'1"),

(l'l'Tt),

( 1' 1' • $ ),( 1' $ ~ J, ), and ( 1' ~ 1' J, ). In addition there

[ 0 1 0 1 1 M=

1

0

1

0

1

0

1

0

1

0

0

1 "

(3.10)

J.C. Eilbeck et al.// The discrete self-trapping equation

325

o

S

......................

....

..........

~x~> .

o

o

-2

I

I

I

I

2

6

10

14

60

Fig. 3. Notation as in fig. 1 for n = 4 ,

M defined in (3.9), N = l ,

ix .

¢q-

o

I

0 C,O

Fig. 4. Detail of fig. 3 for - 1 N ~o _~ 1.

and ~ = 1 .

326

J.C Eilbeck et a l l The discrete self-trapping equation

This system is invariant to any cyclic permutation of sites, so the distinct real asymptotic solutions a r e ( / ' 1' 1' 1 ' ) , ( t I' I' "),(T 1' • "), (I' . . "),(T • I' "), (1' ~, 1' ~), (1' "~ J, ~), ( t O , O ) , (T 1" T ~), (T ~ • "), ( 1' 1' • 4 ), ( 1' • 1' 4 ). In addition we know one complex solution: ( 1' * $ * ), ~ = ½(1, i, - 1, - i), co = 3'/4. T h e following real branches have analytic solutions:

(1' T T t), 0 1 = 0 2 = 0 3 = 0 4 = ½ , co = 3'/4 + 2e; (t

~ T ~),

I#1 = - - 0 2 = 0 3

= --04=½,

co = 3'/4 - 2E; (I'O,LO),

0t = -03 = 1/~-, 02 = 04 = 0, co = 3'/2;

(T 1' * J,), 0 1 = % = - % = co = 3'/4;

-04=½,

(1' " 1' " ), 01 = 03 = ½[1 + (1 - 64e2/3"2)1/2]1/2, 02 = 04 = ½[1 - (1 - 64e2/3'2)1/2] 1/2, co = 3'/2; 01 = 0 2 = ½[1 + (1 - 16e2/y2)l/2] 1/2,

(I'T..),

03 = 04 = ½[1 - (1 - 16e2//3'2)1/2] 1/2, co = 3'/2 + e; 01 = - 0 2 -- ½[1 + (1 - 16e2/3'2)1/2] 1/2,

(1'~..),

03 = - 0 4 = - ½[1 - (1 - 1682//y2)1/2] 1/2, co = 3'/2 - e.

A n analytic but implicit solution for the ( 1' 1' • $ ) branch has also been found (private c o m m u n i c a t i o n , D.C. Heggie). This is 01 = ] s i n ( x 1 -

]~r),

03=~sin(xl+½,r),

04=~sin(x2-}~r),

cos2xl + cos2x2 = ½ and =

- 3e/(4

5)

Very large n

Some general results for the case n >> t have been o b t a i n e d f r o m the perturbation theory described in the following section. Also in ref. 15 we have presented stationary solutions of D S T for n = 100 in the context of self-trapping on crystalline acetanilide [11]. It is sometimes possible to use the solutions discussed above to construct solutions for (3.1) with larger n. As an example consider the dispersion matrix

02 = ] s i n ( x 2 + ]~r),

where

co = y / 3

0 . 6 5 9 0 5 8 . . . < Xx < 2 , r / 3 (on the b r a n c h - ½ c o s - l ( ¼ ) < X2 < 0), then co varies in the range 2 < co < oo. There are bifurcation points on the ( $ 1 ' 1' 1') b r a n c h at 3' = 4 e , 8~, on the (1'" T ' ) b r a n c h at 3' = 602e, on the (1' 1' $ $) b r a n c h at 3' = 4e, and o n the ( 1' ~ • • ) b r a n c h at 3' = 6e. T h e solution curves and stability information for the real solutions are plotted in fig. 5. N o t e that the ( 1' 0 ~ 0) curve (starting at (co, 7) = (0, 0)) overlaps the ( 1' • T • ) curve for co > 4, but the two curves do not meet. The former solution is stable for co > 4 but the latter solution is unstable for co > 4. Similarly the unstable solutions (1' 1' ~, ~) a n d (1' 1' I' £) overlap approximately (without touching) for co > 8. N u m e r i c a l calculations suggest that the complex solution (1' * J, * ) is stable for all co.

sin xtsin

X

2)"

This solution is explicit if X1 is used as a parameter; if X1 varies in the range ½ c o s - t ( ~ ) =

M=

[!100 11 0

1

0

0

1

0

1

0

0

1

0

0

•

.

0

0

(3.11)

,

1

0

.x.

for which (1.1) is a discrete form of the nonlinear SchriSdinger equation [25] with periodic b o u n d a r y conditions• I n this case we can generate solutions of period 1 and 2 (in the index j of the dependent variable 0j) for arbitrary n. T h e period 1 solution is 0y = 1 / f n ; j = 1 . . . . . n; a n d co = y / n + 2e. The equation for co holds if n >_ 3. If n = 1, then co = 3'. If n = 2, then co = y / 2 +e.

J.C Eilbeck et al./ The discreteself-trapping equation

/

O

327

'0,-"

ix. o

o.

o

-2

2

6

10

14

6) Fig. 5. Notation as in fig. 1 for n = 4 , M defined in (3.10), N = l ,

The simplest period 2 solution (n even) is ¢i = -~i+l=l/¢n, i=1,3 ..... n-l, ~0=7/n-2e, n > 4. For n = 2 the equation for ~ becomes ~0 = y / 2 - e. The second period 2 solution (n even) corresponding to the unsymmetric solution in table I is ~i = [(1 + s ) / n ] 1/2, ~i+1 = [(1 - s ) / n ] 1/2, i=1,3 ..... n-l, s = ( 1 - 4 e 2 n 2 / y 2 ) 1/2, o~=

and e = l .

present a perturbation expansion about this limiting case for stationary solutions. This expansion holds only for I'/I << lel and reduces to the linear eigenfunctions in zero order. We also consider the parameter range Id "~ I'll and determine the number of real stationary solutions.

27/n. This solution bifurcates from the period 1 solution given above at y = 2en. These results hold for n > 4; if n = 2 then replace e by e / 2 to regain the solution given in table I. Other solutions of period 3 and higher can be constructed in a similar way. It is important to note that the stability of these solutions may change as n is increased: for example the n = 2 ( 1 ' . ) solution is stable for all of the values of 7 for which this solution exists, but the corresponding period 2 solution for n = 4 ( 1' • 1' • ) is unstable for all values of 7.

4. Perturbation theory It was noted in the introduction that, with 7 = 0, D S T reduces to a linear system. In this section we

1) Perturbation expansion for I~'1 << Id F o r stationary solutions, (1.6) reduces to (3.1) which we rewrite in the form [-YD(I¢,l 2) + eM - o~l] 4) = 0,

(4.1)

where I is the identity matrix. We assume that eM is nonsingular and has n distinct eigenvalues: ~0~°) < o~°) < oo(3 °) < --- < ~0{,°), with corresponding eigenvectors #~o),~o), #~0),..., 4,,(0). Choose one eigenvalue and eigenvector and drop the subscripts for typographical convenience. The stationary solution that reduces to this choice for 7 = 0 can be written as a power series in

J.C. Eilbeck et aL / The discrete self-trapping equation

328

7.

Thus

respect to the inner product

(u, ~) = ~--~.uj*60j,

1# = 1#(0) Jr 71#(1) + 7 21#(2) Jr "" • , o~=o~

(0)+Tw

(1)+72w

(2)+

(4.2)

'''.

Substituting (4.2) into (4.1) and equating terms with equal powers of 7 leads to the following set of equations: [ eM - ~(°)1] 1#(o) = 0,

(4.3.0)

[eM - 0~(°)111# 0) = ~(1)1#(o)_

D(11#(o)12)1#(o),

we observe that L is self-adjoint and its null vector is 1#(o). Thus the solvability condition (1#(o), F1 ) = 0

= (1#(o), i)1#(o)) (1#(o), 1#(o))

[ ~M - ~(°)1] 1#(2) = ~0(2)1#(0)+ ~0)1#0) _ D(11#(°)l 2) 1#(1) 0

_

(

°)* +

0

) (I) ( o ) .

(4.' 4.'

]

+~o).~(o)~. ~.

• - - etc • - •

[eM - a~(°)l] 1#(") =Fm ( (o(,-- 1) ..... (o0), 1#(m-- 1) . . . . . 1#(0)). (4.3.m)

Consider solving these equations in sequence. By assumptions on ~0(°) and 1#(0),(4.3.0) is satisfied. Proceeding to (4.3.1), it is convenient to write it in the form (4.3.1')

where L = [~M - ~(°)1]

(4.4)

and

F1 - ~(1)1#(o) _ D(11#(°)12)1#(o).

~jl@~°)l a - EiI@)O)i------~ .

(4.8)

4(°),

(4.3.2)

L1# (1) - F1,

(4.7)

is satisfied for

(4.3.1)

' ( d~(I).~.(0)* .l,(1)*.t.(O) '~ ~k "1"1 'f'l + ~¢'1 "~'1 /

(4.6)

J

(4.5)

Solution of (4.3.1) requires care because the matrix L is singular. To insure the existence of a solution, F x must be orthogonal to the null space of the adjoint of L. Defining orthogonality with

With this value of w(1), (4.3.1) can be solved for 1#(1). The solvability condition for (4.3.2) is then (1#(o), F2)= 0 which determines ~o(2) and permits (4.3.2) to be solved for 1#(2). This process can be continued indefinitely to determine any coefficient in the power series of (4.2). For sufficiently small 7 we expect these series to converge. In many of the simpler cases considered in the previous section, this series converges after only one step; i.e. ~0(J) and 1#(J) are zero for j > 1. These results can be summarized as follows: i) Assume eM is nonsingular with distinct eigenvalues. Assume also the 7 lies in the range - 70 < 7 < +70 where 7o > 0. For 70/lel sufficiently small, (1.6) has n stationary solutions. ii) For each of these stationary solutions, the corresponding function ~0(7) has the slope d~o dy

~o(1)

at

0, where

-

-

=

7 =

(4.9) ~o (I)

is given in (4.8).

2) The number of real stationary solutions for 171

lel <<

Assume e = 0. Stationary solutions of (1.6) are then determined by

[YD(11#(°)12) - (o(°)1] 1#(o) = 0,

(4.10)

J. C. Eilbeck et al. / The discrete self-trapping equation

where ~(o) and ~(o) nowcorrespond to the fimit e = O rather than y = 0 as in (4.2).Thus

(~lq,~°)l z - ~(°>) q,~°)= O,

(4.11)

for j = 1,2 . . . . . n. Either q~)0)= 0 or Iq,)°)l2 = ~0(°)/'r. Assuming the latter for the first K values of j and the normalization condition (3.5), we conclude that (4.12a)

a~ °) = y / K ,

Iq,)°)12 = a / K , Iq~°)l2 = 0,

for j = 1,2 ..... K, for j = K + 1 . . . . . n.

5. General solutions In this section we discuss general solution trajectories of DST (1.6) following the format of section 3.

1) The case n --~ 1 The general solution is identical to the stationary solution discussed in section 3.

(4.128) (4.12c)

Disregarding phase, there are n such solutions for K = 1, ½ n ( n - 1) solutions for K = 2 and in general n!

329

(4.13)

2) The case n = 2 In this case, DST is exactly integrable. To see this consider (2.10). Defining 0 ~ 0 i -- 0 2

(5.1)

and noting that a 2 + a 2 ffi N, (2.10) can be written al = e( N - a 2 ) X/2sin 0,

solutions. When e 4= 0, the phases of the (q,j) cannot be arbitrarily chosen but must satisfy the requirement that the right-hand sides in (2.10a) be zero. Furthermore the right-hand sides in (2.10b) must be identical for all j. For real solutions the first requirement is satisfied by making ( O j - O k ) an integer multiple of ~r for all j and k. This choice leaves the aj free for adjustment to make the right-hand sides equal in (2.10b). If K = 1, the phase is not significant. If K ffi 2, the phase of the second variable with respect to the first can be either 0 or ~r. In general there are 2 K- 1 significant phases. Thus the number of solutions for each value of K is (~.)2 N-1 and the total number of real stationary solutions is

(g)2 K-l= ½(3"-1).

(4.14)

Ki1

From the preceding section, not all of these are stable.

b=(2a2-N)

y-

_cos___o0 ) ax(N_a2)l/2

(5.2a)

'

(5.2b)

where the most general interaction .matrix is defined in (3.4). Solution trajectories of (5.2) are obtained in analytical form by writing the Hamiltonian (2.2) as a function of the "reduced variables" a 1 and 0. Thus H(al,0)=

- ½ Y ( N 2 + 2a 4 - 2 N a 2) - 2eal(N-a2)t/ZcosO,

(5.3)

where H is the constant energy of a particular trajectory. Several families of solution trajectories in the a t, 0 phase plane, which have been plotted from (5.3), are displayed in figs. 6-8. It is interesting to note that the fixed points in figs. 6-8 correspond to the stationary solutions discussed in section 3 displayed in table I. In

J.C. Eilbeck et aL / The discrete self-trapping equation

330 1.00

particular it is clear from fig. 6 that the stationary solution ( 1' 1' ) is a stable fixed point for "yN < 2e. For "yN > 2e (see figs. 7 and 8) it turns into an (unstable) saddle point and two new stable centers, corresponding to (T • ) and (. 1' ), appear.

0.75 -

O-

0.50

0.25

3)

0.00

I

0.00

1.57

I

3.!4

4.71

6.28

0 Fig. 6. G e n e r a l s o l u t i o n trajectories for n = 2 projected o n t o the a t , # p l a n e u s i n g (5.3) with N = 1, • = 1 and 7 = 0.5.

1.00

0.75 -

0~

0.50 -

0.25

The

case n = 3

We consider the interaction matrix defined in (3.7) which has eigenvalues (2,1,1). In the special case 7 = 0, DST is linear and therefore integrable. The power spectrum of a solution trajectory will display peaks at the eigenvalues of M. In the special case e = 0, DST reduces to three uncoupled oscillators each corresponding to the case n = 1 (see above). Thus DST is again exactly integrable. In the intermediate case, 7 ~: 0 and e ~ 0, we consider numerical solutions using a computer code (GLOP) which is described in [14]. In the following discussion we restrict our attention to solution trajectories that evolve from the initial conditions A d O ) = 1,

-

A 2 (0) = 0.001i, 0.00

a

I

0.00

1.57

31.14

4.71

Fig. 7. Same as fig. 6 except 7 = 2.2.

1.00

0.75-

0.50-

0.25

0.00

0.00

3'.

4'.71

0 Fig. 8. S a m e as fig. 6 except 7 = 5.

(0) = 0.

6.28

e

d

(5.4)

6.28

Here the small initialization on A 2 is to break symmetry. If A2(0 ) were zero, A2(/) would equal A3(t ) and the solution trajectories would be identical to the n = 2 case. Fig. 9a shows the projection of a solution trajectory on the Re(A1), I m ( A t ) plane for 7 = 0 and e -- 1. The corresponding power spectrum (fig. 9b) is calculated as the square of the Fourier transform of Re(A1). As expected this power spectrum has peaks at frequencies 2 and 1. Figs. 10a and 10b show corresponding calculations for 7 = 8 and e --- 0. Since N 1 is equal to 1, (3.5) oscillates harmonically with a frequency equal to 8. This is confirmed by the power spectrum in fig. 10b. A selection of intermediate cases are shown in figs. l l a - p . As 7 / e increases from zero (fig. l l a - d ) ,

J.C. Eilbel:k et al./ The discrete self-trapping equation

331

1.000"

b

a

0 500"

3"

0. 000"

.6"

/

-0 500-

-21°88o'

- 0 . 500'

.g"

oooo'

RE(At)

o soo'

~.ooo

-12

4'

8

12' FREQUENCY

16 • . t.m

Fig. 9. (a) Projection of solution trajectory on the R e ( A l ) , I m ( A 1 ) plane for 3' = 0 and e = 1. (b) Corresponding power-spectrum obtained as the square of the Fourier transform of Re(A~).

1. 0 0 0

b -3"

O. 5 0 0 "

,< -6"

O. 0 0 0 "

9j

-0. 500

"-'~.°B80'

-o. 50o'

o. ooo'

o. 5oo'

-12

,. ooo

,'

RE ( A 1 )

d FREQUENCY

,2'

,e • o •.m

Fig. 10. Same as fig. 9 except 7 = 8 and e = 0.

the power spectra become quasiharmonic and eventually chaotic (fig. 11e-1) in the manner originally suggested by Landau [26]. For y > 6.5e (fig. l l m - p ) the power spectrum again becomes quasiperiodic. 4) The

case n = 4

From numerical studies of this case we obtain solution trajectories and power spectra that are qualitatively similar to those in figs. 11.

5)

Very large n

As was noted in the introduction, our primary motivation for this study has been to understand the nonequilibrium dynamics of globular proteins for which n = 200. Using GLOP [14] we have obtained some results for crystalline acetanilide (a model protein) [11] for n = 100 which are presented in [15]. Detailed dynamical calculations for specific globular proteins are currently under way and will be published in due course.

J.C. Eilbeck et al./ The discrete self-trapping equation

332

0

b

C

d

Fig. 11. Selection of phaseplane trajectories and corresponding power-spectra obtained in the same way as for figs. 9 and 10. The scales on the axes are the same as in figs. 9 and 10. All cams have e = 1 while ~ is increasing: "t = 0.5(a), 1.0(b), 1.5(c), 2.0(d), 2.5(e), 3.0(0, 3.5(g), 4.0(h), 4.50), 5.0(j), 5.5(k), 6.0(1), 6.5(m), 7.0(n), 7.5(o), 8.0(p).

J.C. Eilbeck et aL // The discrete self-trapping equation

e

g

h

Fig. 11. (cont.)

333

334

J.C. Eilbeck et a l . / The discrete self-trapping equation

k

Fig. 11. (cont.)

J.C. Eilbeck et al./ The discrete self-trapping equation

Ill

n

J

0

P

Fig. 11. (concl.).

335

336

J.C. Eilbeck et al./ The discrete self-trapping equation

6. Conclusions The first message that we hope the reader takes from this paper is this: The discrete self-trapping equation (DST) is generic in the spirit of modem nonlinear science. To see this, consider (1.1). With the parameters ~, and e set to zero, it describes a system of n linear oscillators, each having the same frequency (too). Making e > 0 introduces dispersive interactions between these oscillators and making ~, > 0 introduces the nonlinear effects of self-trapping. The dynamic interaction between these two effects is displayed by solutions of DST. Like other generic equations (for example K o r t e w e g - d e Vries, sine-Gordon, nonlinear Schr6dinger, etc.), it is an important model in applied science with the following applications: i) Dynamics of small molecules [19]. The XH stretch vibrations in water (n = 2), ammonia (n = 3), methane (n = 4), and benzene (n = 6) are of particular interest, ii) Dynamics of molecular crystals [15]. Here the interest is in highly anharmonic hydrogen bonded crystals such as acetanilide, iii) Self-trapping in amorphous solids [16]. An important future research problem is to consider the augmentation of "Anderson localization" by anharmonic effects, iv) Globularprotein [9, 10, 13 14]. Here the aim is to solve one of the most fundamental problems in modem biology: How do functional proteins store and transport biological energy? Specific results presented here for stationary solutions include complete analytic descriptions for n < 2, bifurcation diagrams (including calculations of hnear stability) for n = 3 and 4, and some conclusions from perturbation theory for arbitrarily large n. DST is shown to be exactly integrable for n < 2 and for all n if V = 0 or e = 0. From our numerical studies of the general case (~, 4:0 and e ~ 0) for n = 3 and 4, the following picture emerges. As ~,/e is increased from zero, KAM islands become smaller and regions of chaotic instability grow. For "//e ~ 1, chaos dominates. Then, as ~,/e is further increased, the sea of KAM instability

recedes and new islands appear. A similar situation was reported in ref. 27.

Appendix A

Linear stability theory In this appendix we sketch the linear stability theory for the DST eqs. (1.1), following the results of [24]. We consider the stability of stationary solutions only, i.e.

A(t) = exp (itot)~,

(A.1)

where ~b is time independent, and to, ~ satisfy (2.3). Consider now an arbitrary perturbation u(t) to the solution (A.1), written in the form a ( t ) = exp (itot)(~ + u(t)).

(A.2)

Notice that we have added the perturbation in a frame rotating with the stationary solution: since u(t) is an arbitrary function of t, this does not involve any approximations, but it is a crucial step in simplifying the resulting mathematics. We first consider the case when # is real. Insert (A.2) into (1.1), use the fact that to, ~ satisfy (3.2), and drop terms quadratic and higher in the perturbation u(t) to get the linear equation satisfied by u(t) for small times

iil-o~u+

,/diag ( (i)2}(2u+

u*)+eMu=O.

(A.3)

Here diag ( f / ) is the diagonal matrix

diag { f,} =

Ii

A

° L

and u* is the complex conjugate of u r Due to the choice of perturbation (A.2), (A.3) is an autonomous set of equations. Now write u - - u 1 + iu2

J.C Eilbeck et a l l The discrete self-trapping equation

where u~ and //2 are real n-vectors, then (A.3) becomes d['q] d-7 u2 ~ M - ~1 +'yD3(~)

0

u2 '

337

where q~ = ~1 + i~2, etc. If all the eigenvalues of the coefficient matrix in (A.7) have real parts < 0 ( > 0) this is a necessary (sufficient) condition for stability (instability).

N o t e added in proof

(A.4) where I is the n × n unit matrix and Dr(q0 = diag { l~ 2 }. Denote the n x n blocks in the 2n x 2n matrix in (A.4) by C, d respectively so that (A.4) becomes (A.5)

d t [u2J = [d

N o t e that - C is the matrix multiplying q~ on the left-hand side of (3.2), and d is the Jacobian matrix of this equation. Eigenvalue of the coefficient matrix in (A.5) are given by det [ d e - h2I ] = O.

(A.6)

A necessary condition for the stability of the solutions (A.1) is that all the h i have real parts < 0. F r o m (A.6) this can only be achieved if all the h ] are purely real and negative, giving purely imaginary Xi. A sufficient condition for instability is that at least one X] is not purely real and negative, giving one root h i which has positive real parts. Hence a simple algebraic test (performed numerically if necessary) gives necessary conditions for stability or sufficient conditions for instability for the stationary solutions (A.1) of (1.1). Further details are given in [25]. The case when q~ is complex can be treated in a similar manner: in this case (A.4) becomes

~-

u2

to+l=

30- sin (8 + 00) sin 30

2 ¢5

+2 =

sin(0-

lk3 = to4 = 7 3 v~ + 1 sin0, 4 i~l

where 0o = t a n - x(3¢3) - ~r/3 = 0.33347 . . . . For 0 < 0 _< ~r/6 this gives the branch 1' -" •, for ~r/6 < 0 < ~ ' / 3 this gives the branch T $ ? - , and for 2~r/3 < 0 < ~r - 00 this gives the branch 1' T 1' $. This solution can be generalised to the case of arbitrary n >_ 3, in the equal interactions case (rnij = 1, i ~ j ; rn, = 0). In this case we have

to+l=

3(n 2 - 3n + 3)x/2 sin (0 + 0o) sin 3 0

qq = - - ~ v/-~ + 1 s i n ( 0 - -~'),

¢2=%

.....

q% = ~3 f ~ + 1 sin0,

i=l

eM - to| + ,fdiag ( 3,~ 1 + *i~ }

2~/diag { ~n~i2 }

After this paper was written, we found the analytic form for the remainder of the stationary solutions for the case 4.i, n = 4, equal interactions, M given by (3.9). These are, in parametric form,

u2 ' (A.7)

where 0o -- tan -1 [(n - 1 ) ¢ 3 / ( n - 3)] - ~t/3. In particular for 0 < 0 _< 7r/6 we get the stationary "soliton" mode ? ( . ) ' - 1 in which all of the energy becomes localized on one site as "t ~ oo.

338

J.C. Eilbeck et al./ The discrete self-trapping equation

References [1] L.D. Landau, Phys. Zeit. Sowjetunion 3 (1933) 664. [2] S. Pekar, J. Phys. U.S.S,R. 10 (1946) 341,347. [3] H. Fr6hlich, H. Pelzer and S. Zienau, Phil. Mag. 41 (1950) 221; H. Fr6hlich, Adv. in Phys. 3 (1954) 325. [4] T. Holstein, Ann. Phys. 8 (1959) 325, 343. [5] S. Fischer and S.A. Rice, J. Chem. Phys. 52 (1970) 2089. H. Haken and P. Reineker, Zeit. f. Phys. 249 (1972) 253. H. Sumi, J. Phys. Soc. Japan 32 (1972) 616. H. Haken and G. Strobe, Zeit. f. Phys. 262 (1973) 135. K. Tomioka et al., J. Chem. Phys. 59 (1973) 4157. W. Weidlich and W. Hendorfer, Zeit. f. Phys. 268 (1974) 133. G.C. Morris and M.G. Sceats, Chem. Phys. 3 (1974) 342. D.P. Craig and L.A. Dissado, Chem. Phys. 14 (1976) 89. D. Yarkony and R. Silbey, J. Chem. Phys. 65 (1976) 1042. D.P. Craig and L.A. Dissado, Chem. Phys. Lett. 44 (1976) 419. DP. Craig, L.A. Dissado, and S.H. Walmsley, Chem. Phys. Lett. 46 (1977) 191. D.K. Campbell, A.R. Bishop and K. Fesser, Phys. Rev. B26 (1982) 6862. H. De Raedt and A. Lagendijk, Phys. Rev. B80 (1984) 1671. [6] V.E. Zakharov, Zh. Eksp. Teor. Fiz. 62 (1972) 1745 [Sov. Phys. JETP 35 (1972) 908]. J. Gibbons, S.G. Thomhill, M.J. Wardrop, and D. ter Haar, J. Plasma Phys. 17 (1977) 153. V.E. Zakharov and A.B. Shabat, Zh. Eksp. Teor. Fiz. 61 (1971) 118 [Sov. Phys. JETP 34 (1972) 62]. [7] A. Hasegawa and F. Tappert, Appl. Phys. Lett. 23 (1973) 142. C.F. Mollenauer, R.H. Stolen and J.P. Gordon, Phys. Rev. Lett. 23 (1973) 142. A. Hasegawa and Y. Kodama, Proc. IEEE 69 (1981) 1145. [8] J. Wu, R. Koolian and I. Rudnick, Phys. Rev. Lett. 52 (1984) 1421. A. Larraza and S. Putterman, Theory of nonpropagating hydrodynamic solitons, Phys. Lett. A 103 (1984) 15. [9] A.S. Davydov and N.I. Kislukha, Phys. Stat. Sol. (b) 59 (1973) 465. A.S. Davydov, J. Theor. Biol. 38 (1973) 559; Physica Scripta 20 (1979) 387; Biology and Quantum Mechanics (Pergamon, Oxford, 1982); Sov. Phys. Usp. 25 (1982) 899 and references therein. [10] J.M. Hyman, D.W. McLaughlin and A.C. Scott, Physica 3D (1981) 23. A.C. Scott, Phys. Rev. A26 (1982) 578; 27 (1983) 2767; Physica Scripta 25 (1982) 651. L. MacNeil

and A.C. Scott, Physica Scripta 29 (1984) 284. [11] G. Careri, U. Buontempo, F. Carta, E. Gratton and A.C. Scott, Phys. Rev. Lett. 51 (1983) 304. G. Careri, U. Buontempo, F. Galluzzi, A.C. Scott, E. Gratton and E. Shyamsunder, Spectroscopic evidence for Davydov-like solitons in acetanilide, Phys. Rev. B30 (1984) 4689. [12] M.V. Volkenstein, J. Theor. Biol. 34 (1972) 193. [13] D.E. Green and S. Ji, Proc. Nat. Acad. Sci. (USA) 69 (1972) 726. [14] P.S. Lomdahl, Nonlinear dynamics of globular protein, Los Alamos Lab. Report LA-UR-83-2252; in: Nonlinear Electrodynamics in Biological System, eds. W.R. Adey and A.F. Lawrence (Plenum, New York, 1984) 143. [15] J.C. Eilbeck, P.S. Lomdahl and A.C. Scott, Soliton structure in crystalline acetanilide, Phys. Rev. B30 (1984) 4703. [16] P.W. Anderson, Phys. Rev. 109 (1958) 1492. D.J. Thouless, Phys. Rep. 67 (1980) 5, 13 (1974) 95. [17] H.J. Stapleton, J.P. Allen, C.P. Flynn, D.G. Stinson and S.R. Kurtz, Phys. Rev. Lett. 45 (1980) 1456. J.P. Allen, J.T. Colvin, D.G. Stinson, C.P. Flynn and H.J. Stapleton, Biophys. J. 38 (1982) 299. [18] S. Alexander, C. Laermans, R. Orbach and H.M. Rosenberg, Phys. Rev. B28 (1983) 4615. [19] A.C. Scott, P.S. Lomdahl and J.C. Eilbeck, Between the local mode and normal mode limits, Chem. Phys. Lett. 113 (1985) 29. [20] J.W. Ellis, Trans. Faraday Soc. 25 (1929) 888. [21] B.R. Henry, in: Vibrational Spectra and Structure, J.R. Durig, ed. (Elsevier, Amsterdam, 1981) 269. M.L. Sage and J. Jortner, Adv. Chem. Phys. 47 (1981) 293. [22] M.A. Collins, Adv. Chem. Phys. 53 (1983) 225. [23] D.W. Decker and H.B. Keller, Comm. Pure and Appl. Math. 34 (1981) 149 and references therein. M. Kubicek, A.C.M. Trans. Math. Soft. 2 (1976) 98. [24] J. Carr and J.C. Eilbeck, Stability of stationary solutions of discrete self-trapping equation (to appear in Phys. Lett.). [25] R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, and H.C. Morris, Solitons and Nonlinear Wave Equations (Academic Press, London 1982). [26] L.D. Landau, Dok. Akad. Nauk. SSSR 44 (1944) 339; C.R. Acad. Sci. URSS 44 (1944) 311. [27] Genesis 7:11 through 8: 5.