{'fraction_non_alphanumeric': 0.0374331550802139, 'fraction_numerical': 0.009893739031905926, 'mean_word_length': 4.154925564498896, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unstable decay and state selection II

The decay of unstable states when several metastable states are available for occupation is investigated using path-integral techniques. Specifically, a method is described which allows the probabilities with which the metastable states are occupied to be calculated by finding optimal paths, and fluctuations about them, in the weak noise limit. The method is illustrated on a system described by two coupled Langevin equations, which are found in the study of instabilities in fluid dynamics and superconductivity. The problem involves a subtle interplay between non-linearities and noise, and a naive approximation scheme which does not take this into account is shown to be unsatisfactory. The use of optimal paths is briefly reviewed and then applied to finding the conditional probability of ending up in one of the metastable states, having begun in the unstable state. There are several aspects of the calculation which distinguish it from most others involving optimal paths: (i) the paths do not begin and end on an attractor, and moreover, the final point is to a large extent arbitrary, (ii) the interplay between the fluctuations and the leading order contribution are at the heart of the method, and (iii) the final result involves quantities which are not exponentially small in the noise strength. This final result, which gives the probability of a particular state being selected in terms of the parameters of the dynamics, is remarkably simple and agrees well with the results of numerical simulations. The method should be applicable to similar problems in a number of other areas such as state selection in lasers, activationless chemical reactions and population dynamics in fluctuating environments.

the non-linear coupling between parallel rolls of different wave numbers may be derived [2], and when the effects of the other modes are incorporated [3], a set of coupled stochastic differential equations of the kind mentioned above are generated. In superconductivity, exactly the same equations as in the above example govern state selection in, for example, narrow superconducting rings. This is simply because the amplitude equation which governs the instabilities in the fluid dynamics example is nothing but the Ginzburg-Landau equation for a superconductor [4]. A mode truncation then gives exactly the same set of coupled differential equations acted upon by noise [5]. In chemical kinetics, it has become clear in the last decade or two, that there are many important chemical reactions in which a barrier to the formation of an excited state is not present. Examples of these activationless reactions [6] are the electronic relaxation of triphenylmethane dyes and barrierless electron transfer in solution. The potential discussed above is a reaction potential energy surface in this case [7]. In lasers, such as the ring dye laser [8], decay of an unstable mode to metastable ones can occur under the right operating conditions. The coordinates in this case are the mode amplitudes and Langevin-type equations are derived within the semiclassical theory of the laser [9]. In population dynamics, the Gause model of two competing species [10] again falls into this class. When the competition is played out with a fluctuating environment, the resulting stochastic differential equations once again fall into the generic class that we have been discussing [11]. It should, however, be noted that in this case there is no potential and moreover the noise is multiplicative.
Although many of these phenomena have been known for some time, the investigation of the state selection aspect has been hampered by the lack of a suitable calculational tool. In the case where a potential V exists, the noise is additive and only two modes, x and y, are considered, the equations take the forṁ where η x and η y are white noises of strength D. But, even in this, the simplest non-trivial case, these equations are difficult to study mainly because, as we shall discuss later, in the region where state selection occurs, the coupled nature of the equations, their non-linearity and the noise, are all important. However, as we have shown [1], there is a method that can take all these aspects into account, and that is the path-integral formulation of stochastic dynamics. This method succeeds where others fail because the equations can be represented as a path-integral without approximation; systematic approximation techniques for path-integrals developed over the years can then be used as the basis of a calculational scheme. Most of the previous theoretical work on this problem has been limited to systems with one degree of freedom, governed by the single equationẋ = −V ′ (x) + η(t), because many of the complexities mentioned above are not present in this one-dimensional case. Suzuki and co-workers developed a theory for the decay of an unstable state in one dimension in a series of papers [12], as did a number of other authors [13][14][15][16][17][18][19]. However, the one-dimensional theory, although much easier to deal with, has none of the subtleties inherent in state-selection in higher dimensions: the decay is simply either to the right or to the left. Some studies [20][21][22] purported to go beyond one dimension, but in fact considered spherically symmetric potentials, so that the problem could be reduced to a quasi-one-dimensional problem in the radial coordinate. Once again, the resultant structure is not rich enough to address the question of state selection. Other approaches, such as replacing the stochastic differential equations by the corresponding deterministic equations, but with random initial conditions [23], are unsatisfactory in other ways. Probably the investigation nearest to our own has been by Mangel [24], however his main interest was not state selection.
The plan of the paper is as follows. In Sec. II we discuss the assumptions underlying the picture of state selection which we have outlined and introduce a generic model which we use to describe our calculational scheme in detail. An important aspect of our method is the realization that, typically, a state is selected well before the system reaches this chosen state. This fact is used to simplify the problem in Sec. III, where it is also shown that naive calculational prescriptions, based on linearization of the initial dynamics, fail. A more systematic approach based on path-integral techniques is introduced in Sec. IV and the relevant optimal path is determined. This gives the leading order contribution in the limit D → 0. The next to leading order contributions are determined in Sec. V. The results of this analytic approach are given and compared to Monte Carlo simulations in Sec. VI and our conclusions are presented in Sec. VII. There are two appendixes. Appendix A describes some of the elementary approaches described in Sec. III in more detail. Appendix B contains technical aspects relating to the determination of the y-coordinate of the optimal path and the calculation of the action of that path.

II. GENERAL CONCEPTS
We have already given several examples of situations where state selection, following the decay of an unstable state, occurs. In this section we will explore the different ways in which the initial state of the system could arise, that is, the origin of the unstable state, and give an intuitive description of the decay and subsequent state selection, which will form the basis of our analytical approach.
An unstable state may arise by mechanisms which are either non-adiabatic or adiabatic. In the non-adiabatic case, the system, in an initially (meta)stable state, is transported very rapidly to an unstable state. The time scale for the transition from the stable to unstable state is much more rapid than the characteristic time scale defining the natural dynamics of the system. In this context the transition from the stable to the unstable state can be ignored and we simply characterize the system as having been prepared in an unstable state. Examples of this type were mentioned in Sec. I and include chemical reactions with no activation barrier and population dynamics in a fluctuating environment. Quasi-one dimensional superconductors [5] are another example. In the case of the chemical reactions, an ultra-fast light pulse (a femtosecond laser) pumps molecules into an excited state. In terms of potential surfaces, the molecules are initialized in an unstable state on the reactive surface. The subsequent reaction can be viewed as a nuclear rearrangement on this reactive potential energy surface. The evolution of the reaction, that is the relative preponderance of reactants and products, can be studied using other ultrafast lasers. In the population dynamics example, the population is initialized so as to consist of only a few individuals of each species. During the initial period the population of both species grow exponentially, independently of the other, until competition drives the population of one of the species to zero.
In the adiabatic case, an unstable state may arise when a (meta)stable state is transported through the instability over a time scale that is slow compared to the time scale characterizing the natural dynamics of the system. An example of this process is provided by a quasi one-dimensional superconductor as studied in Ref. [25]. In this case, the system is driven by a voltage source that accelerates the supercurrent. However, this acceleration cannot continue indefinitely and the system becomes unstable, a situation associated with the critical current of the superconductor. Once unstable, the relaxation process occurs via "phase slips" in which spatially localized regions of the wire temporarily lose their superconducting properties and carry "normal", i.e. non-superconducting, current. This process dissipates the excess energy and the system relaxes back to a locally stable state. Mesoscopic wires, i.e. wires that are not in the thermodynamic limit, have a finite number of discrete metastable stationary current carrying states. State selection in this case is characterized by the competition between these metastable current carrying states. The superconductivity instability described above is an example of the Eckhaus instability.
Eckhaus instabilities arise in many physical systems in addition to quasi one-dimensional superconductors, including fluids, nematic liquid crystals and lasers. In these systems, stationary one-dimensional periodic patterns are stable for a range of wave vectors, Q. For instance, in the case of the common Eckhaus instability, stationary solutions exist for Q 2 < 1, but are only stable if Q 2 < 1/3. By changing the control parameter slightly, Q-states with Q 2 < 1/3 may be shifted into the unstable Q 2 > 1/3 regime. The essential features of the subsequent changes which follow from this may be understood by performing a mode-decomposition of the relevant amplitude equation and keeping only the previously stable mode and the destabilizing modes with the largest growth rates [5].An adiabatic elimination of the unstable mode then leaves us with coupled ordinary differential equations for the amplitudes of the destabilizing modes. This adiabatic assumption is equivalent to assuming that the form of these differential equations does not change on the time-scale of state selection, i.e. in the time taken for one of the destabilizing modes to outcompete the others [25] An obvious question which now arises is the following: how do we know that the system is initialized in a state where the amplitudes of the competing modes, final products, etc are zero? Or expressed in the topographical terminology of Sec. I: How do we know that we begin exactly at the top of the hill? Here we are assuming that the small non-zero amplitude, required to initiate the growth, is induced by fluctuations i.e. by the noise. Once again, there are several possibilities. It may be that in every realization of the stochastic process the system starts near, but not necessarily at the origin. As long as there is no bias favoring any one of the competing modes, an average over many realizations of the process will give the same results as if the system had started at the origin in each case. Alternatively, the initial probability distribution may be so sharply peaked about the origin, that any small initial bias is irrelevant to the ultimate choice of state made by the system. Yet another possibility is that it is just not possible to start exactly at the origin (the example taken from population dynamics is a illustration of this). In these cases the sensitivity of the final result to changes in the initial starting point will have to be carefully investigated. In addition to all of these specific reasons, on general grounds beginning at the origin seems very natural, since it simply means that none of the final products or final modes are initially present.
Since there appear to be good reasons to initialize the system at the origin in many cases, it will be assumed in the subsequent theoretical development. In all of the examples given in Sec. I, the competing modes grow exponentially for an initial period, without any significant interaction with each other. This corresponds to a linear growth law of the formẋ l = α l x l , where x l is the amplitude and α l the growth rate of the lth mode. Clearly, the main interest from a state-selection point of view, is in the nature of the non-linear interaction terms between the modes. Since the purpose of this paper is to give a clear and intuitive introduction to our calculational scheme, we shall choose to illustrate our method on an example where only two competing modes are present, and where the interaction terms may be derived from a potential. In other words, the governing equations will have the form (1). Obviously, two is the minimum number of modes we need in order to discuss state selection, but having only two modes moving on a potential surface has the advantage that we can describe our method using the pictorial language of "hills and valleys" which we have already adopted on several occasions. Furthermore, since as we have already mentioned, the equations derived in both the fluid dynamics and superconductivity examples are identical [2,5] we will work with these. They areẋ where x and y are the two competing modes with growth rates α and β respectively. All five parameters characterizing the model (α, β, γ, δ and ǫ) are assumed to be positive. The noise terms η x (t) and η y (t) are taken to be Gaussian random variables with zero mean and where i and j are either x or y and D is the noise strength. Equations (2) and (3) fall into the class of those given by (1), since they are derivable from the potential We wish, however, to stress that the method is not restricted to systems with only two modes nor is the existence of a potential a prerequisite. This will become clear as the method is explored in subsequent sections. It is interesting to note that the simple generalization of the model in which the strengths of the noises η 1 (t) and η 2 (t) are both of the same order, but not equal, already leads to a more complicated situation. If we suppose that the strengths of these noises are D 1 and D 2 respectively, then we may introduce ζ i (t), i = 1, 2 such that the correlation function of these new noise terms is exactly (4), by writing η i (t) = ρ 1/2 i ζ i (t) where the ρ i are constants given by D i = ρ i D. This rescaling of the noise may, in turn, be eliminated from Eqns. (2) and (3) by rescaling x(t) and y(t) by defining new variables x 1 and x 2 via x(t) = ρ 1/2 1 x 1 (t) and y(t) = ρ 1/2 2 x 2 (t). New constants replacing δ and ǫ which absorb this change of scale can easily be defined, but the interaction terms γxy 2 and γyx 2 become γ 12 x 1 x 2 2 and γ 21 x 2 x 2 1 respectively, with γ 12 = γ 21 . Equations such as this are not derivable from a potential, nevertheless minor modications of our method are applicable to this case.
The potential V (x, y) for a particular choice of parameters is shown in Fig 1(a). Before discussing the dynamics of the system, we have to identify the extrema of the potential. These are given by the intersection of the two sets of curves ∂V /∂x = 0 and ∂V /∂y = 0 or, explicitly, by the intersection of {x = 0, δx 2 + γy 2 = α} and {y = 0, γx 2 + ǫy 2 = β}. These two ellipses, along with the x and y-axes, are shown in Fig 1(b). By comparing Figs 1(a) and Fig 1(b) we see that there is (i) a maximum at the origin (denoted by a filled circle), (ii) four minima at the intersection of the ellipses and the axes (denoted by open circles), and (iii) four saddle points at the intersection of the two ellipses (denoted by crosses). The precise positions of these extrema are: Saddles : (x, y) = 1 In fact, there are some conditions which need to be imposed on the parameters of the model if saddle points are to exist. These are that αγ − βδ, βγ − αǫ and γ 2 − δǫ should all have the same sign. Since we will normally assume that the stability parameters, δ and ǫ, which ensure that the potential is bounded below, are small, we will take βγ > αǫ and αγ > βδ, which implies that γ 2 > δǫ. The dynamics given by (1), (4) and (5) is equivalent to an overdamped particle moving on the potential surface shown in Fig 1(a) and which is also acted upon by white noise. In a typical realization, the particle will begin at the maximum and perform a random walk which is in the vicinity of the origin at early times, but which eventually explores an ever larger region. Eventually, the particle gets far enough from the origin that the deterministic dynamics, specified by V , begins to have a significant impact and the particle accelerates down towards the saddles and the minima. In this paper our main concern is state selection: we are not concerned with the approach the particle makes to a particular minimum, since by this stage this minimum will almost certainly be the selected state. For this not to be so, the particle would have to hop over a barrier -an extraordinarily rare event. In fact, as is clear from Fig  1(a), as soon as the particle has passed the x or y coordinate of one of the saddle points, it has effectively chosen the final state and its subsequent motion is of little interest to us. We therefore arrive at the conclusion that the saddle points are the major factor influencing state selection and by comparison the minima are of little consequence. This can be made more transparent if we imagine that δ and ǫ are very small, so that the minima (6) are now at very large values of x and y. By contrast, the positions of the saddles have changed much less: in the limit where δ and ǫ go to zero their x and y coordinates tend to the finite values β/γ and α/γ, whereas the minima tend to infinity. The topography now consists of long valleys, leading from the head of the valley between two saddles, all the way down to a minima at the end of the valley. As soon as the particle enters the valley it is extremely unlikely to escape over the sides, and so extremely unlikely it will do anything else but move to the minima at the end of that particular valley. This makes it very clear that it is the saddles, and not the minima, that we should focus on if we wish to understand state selection.
These ideas are well illustrated by plotting out a few Monte Carlo trajectories. Such simulations are extremely simple to carry out, and later we will compare our analytic expression for the probability of a particular state being selected with the proportion of runs which ended up in that state. For the moment, however, we are interested in the nature of individual trajectories. The five which are shown in Fig 2 are typical: the particle carries out Brownian motion about the origin, to a greater or lesser extent, and then selects a particular valley. Even after state selection, there may still be large deviations, but eventually the particle settles down near to the axis along which the valley runs.

III. THE REDUCED POTENTIAL
We have seen in the last section that the essential features of state selection become much clearer in the limit δ, ǫ → 0, when the saddle points are completely separated from the minima, which have moved off to infinity. For the rest of the paper we work in this limit in order to illustrate our method in the simplest possible way. We will call the problem defined in Sec. II when the limit is taken, the reduced problem. The reduced potential is then defined to be A plot of V R looks very similar to the plot of the full potential V in Fig 1(a), since this figure emphasizes the littlechanged central features of the potential. The main difference is that whereas in V the valley floors start to ascend again at large |x| and |y|, in V R they keep descending. However, as was made clear in Sec II, this behavior is of no interest to us, and the fact that V and V R differ in this way is immaterial. From (8), we see that V R has only two types of extrema: (i) a maximum at the origin, and (ii) four saddle points at (x, y) = (±X min , ±Y min ) where The reason for the subscript "min" is that the values given by (9) are the smallest values of |x| and |y| at which we can assume that state selection has taken place. A simple stability analysis near these saddles shows that the stable and unstable directions are at an angle of π/4 to the axes. For example, in the positive quadrant, if the particle approaches the saddle along the line of steepest descent y = x, it will have a tendency to move down the slopes which are to the right and left, rather than carry on through the saddle and straight up the incline ahead. Just as the plot of V R is similar to that of V if we cut it off near to the head of the valleys, so the Monte Carlo simulations are similar in both cases if we do not follow the trajectories too deeply into the valleys. In Fig. 3, the range of x and y used to display the trajectories are smaller than those in Fig. 2, so that the central region is emphasized. This makes the initial Brownian dynamics seem more obvious, but in reality the nature of the trajectories before and during state selection are essentially identical to those of the full problem shown in Fig. 2. In the next section, we will develop a calculational scheme, based on the path-integral formulation of the stochastic process, to determine the probabilities of entering the x-or y-valleys as a function of the model parameters α, β, γ and D. However, one might feel that such a sophisticated theory is not necessary: it should be possible to obtain a satisfactory theory by constructing an approximation based on the picture we have built up in this and the previous section. We will therefore end this section by constructing an example of such a theory and show that it is unsatisfactory for a variety of reasons.
A simple theory of state selection might have the following ingredients: 1. In the initial period the noise and linear growth of the modes dominate. Therefore, neglect the non-linear interaction between the modes (i.e. set γ = 0).
2. State selection will be specified in the following way. Define four sectors in the xy-plane by drawing lines through the origin and the saddle points. The particle will select the state corresponding to a particular sector if it is in that sector for x ≈ X min and y ≈ Y min .
A calculation based on these two assumptions is given in Appendix A. Even within the framework that these provide, there turn out to be many possible variants, each giving slightly different results. Leaving this aside for the moment, it is shown in Appendix A that in one of the simplest schemes along these lines, the probability of ending up in an x-valley is (see Eqn. (A14)) whereD This has the correct qualitative behavior: when α = β, that is, ρ = 0, N x = 1/2. As ρ increases, N x increases towards unity. A more careful calculation gives (Eqn. (A15)) which again has the correct qualitative features. We now compare these two formulas with the results of Monte Carlo simulations. We took β, γ and D to be fixed in the simulations, at values 1.0, 1.0 and 0.01 respectively, and varied α from 1.0 to 3.0 in steps of 0.2. For each value of α we performed a large number of runs where the particle started at the origin and ended up in one or other of the valleys (for more details see Sec. VI). We then simply counted the number of times that the particle ended up in a particular valley and expressed this as a fraction of the total number of runs made. The results are shown in Fig  It is clear that the results of the naive approach based on assumptions 1 and 2 above, worked out in Appendix A and given by (10) and (12) are not quantitatively correct. They mainly overestimate the probability of the particle ending up in the x-valley. In addition, Eqn. (12), which represents a more refined calculation, actually compares less favorably with the simulation results, than does (10). This clearly points to the unsatisfactory nature of this simpleminded scheme, as one would expect the more refined calculation to compare more favorably with the simulation results.
We do not wish to pursue techniques based on this approach further in this paper. It was included specifically to show that there is a delicate interplay between the non-linearity and the noise in (2) and (3) and any attempt to incorporate one or the other of these using an unsystematic approximation scheme, as above, is likely to give disappointing results. We do not rule out being able to construct a particular scheme of this type which will give reasonable results: there are enough possible variants that it may be possible to interpolate between these by introducing free parameters which can then be fitted to the Monte Carlo data. However, such an ad hoc scheme is unnecessary, since we will show in the remainder of the paper that a systematic approach exists which gives simple formulas which are in good agreement with the Monte Carlo data. The method is based on the use of optimal paths, and so we begin with a brief review of the path-integral formulation of Langevin dynamics.

IV. OPTIMAL PATHS
In this section we review the formulation of stochastic differential equations, such as (1), as functional integrals and obtain the dominant contribution to the conditional probability we wish to determine in the limit where the noise strength tends to zero.
The conditional probability that the system is in the state (x f , y f ) at time T , given it was initially in the state (0, 0) at t = 0 is where IC denotes the initial condition x(0) = 0, y(0) = 0 on the stochastic process and x(T ) and y(T ) are solutions of (2) and (3). The average in (13) is over Gaussian white noises η x (t) and η y (t) with zero mean and correlation function given by (4). In terms of functional integrals (13) equals where C is a normalization constant and the subscript η on x(T ) and y(T ) is to emphasize that they depend on η(T ) through (2) and (3). Using these equations to perform a functional change of variable from (η x , η y ) to (x, y) yields where J is the Jacobian of the transformation. Expressed as a path-integral where r = (x, y). The action S[ r] and the Jacobian J[ r] are functionals which are given by and [26] J For D → 0, the path-integral (16) is dominated by solutions of the Euler-Lagrange equations δS/δ r(t) = 0, which satisfy the boundary conditions r(0) = 0 and r(T ) = r f . Let the solution of least action be denoted by r c (t). Then writing r(t) = r c (t) + δ r(t), we have Scaling δ r by D 1/2 and performing the Gaussian functional integral yields The new overall constant, C ′ , is immaterial since, as we have seen in Sec. III, we are interested only in the probability of ending up in the x-valley or y-valley, and we normalize these probabilities according to N x + N y = 1. We will therefore omit the overall constant from now on.
While all the derivations in this section have been carried out for potential systems with two degrees of freedom, it should be clear that they generalize in an obvious way to systems of more than two degrees of freedom and those where no potential exists [26]. We may summarize the result of performing the functional steepest descent on (16) to next-to-leading order by where the leading order contribution, P (0) , is just the action of the optimal path r c (t) and the next-to-leading order contribution is Here L is the matrix formed from the second order functional derivative of the action functional evaluated at the optimal path: The result (21) is the starting point for our method: if we can determine the functions P (0) and P (1) , then we will have a form for the conditional probability valid when the noise is weak. This in turn will enable us to obtain a formula for the probability that either the x-valley or y-valley is selected, as a function of α, β, γ and D. We shall devote the rest of this section to the determination of P (0) and leave the calculation of P (1) to the next section. The case of interest to us in this paper is when V is the reduced potential (8). The Euler-Lagrange equations for this problem areẍ It is important to realize that there are two distinct dynamics associated with the problem under consideration. The first is the stochastic dynamics given by (1) with the potential (8). This was our starting point, the basis of the intuitive discussion of the dynamics given in Secs. II and III, and the dynamics of the Monte Carlo simulation. The second dynamics is the deterministic dynamics, given by (24) and (25), which describes the D → 0 limit of the stochastic dynamics. They are quite different and it is important not to carry over intuition from one to the other without careful consideration. From (17), where Since the last term in (26) is a constant, and consequently gives zero variation, the Euler-Lagrange equations (24) and (25) correspond to classical mechanics in the potential When considering the optimal paths such as r c (t), it is the potential U R , and not V R , that is relevant. A plot of U R (x, y) is shown in Fig. 5. We will, for concreteness, focus on paths which end in the positive x-valley; paths which end in the other valleys are treated in exactly the same way. The boundary conditions on these paths are x(0) = y(0) = 0, x(T ) = x f and y(T ) = y f , with y f ≪ x f . We therefore simplify (24) and (25) by keeping only terms which are linear in y, to obtain to lowest order,ẍ = α 2 x and (29) In our earlier treatment of this problem [1], we restricted attention to end-points lying on the x-axis: x f = X and y f = 0 (and correspondingly y f = Y and x f = 0 for paths ending in the positive y-valley) arguing that the relative contributions coming from paths ending on the axes should be approximately equal to the relative contributions coming from the paths ending anywhere in the valleys. In fact, as shown in Ref. [1], this turns out to be an excellent approximation, but it has the deficiency that X is not determined. In this paper, we go beyond this approximation, which allows us to explore its validity, but also, as we will see, the precise value of X within this improved treatment does not enter; the result is independent of X, provided that X is large enough. The quantity y f has a different character to X: it will be integrated over at a later stage, when the probability flux through the valley is determined. The solution of (29) which satisfies the boundary conditions x(0) = 0 and x(T ) = X is exactly the same as in Ref. [1], namely, The solution of (30) which satisfies the boundary conditions y(0) = 0 and y(T ) = y f (we reserve the notation Y for the final y value for the end-point of paths going along the y-valleys) is discussed in Appendix B. There it is shown that an excellent approximation to the solution of this equation (indistinguishable from the numerical solution) can be found by first solving the equation for small t and then for large t, and matching the two at some intermediate matching time t m . Comparison with the numerical solution shows that the approximation remains good for a wide range of values of t m from about 0.1T to 0.6T or 0.7T . In order to be able to obtain a simple form for the solution, we have assumed that T is large, in the sense that e −αT ≪ 1 and e −βT ≪ 1. Using (B3) and (B7), we can write the explicit analytic from as For self-consistency we need to check that y c (t) given by (32) is small, compared with x c (t), if we are to justify the linearization procedure which led to (29) and (30). From (32) it is straightforward to check that y c (t) has a maximum at t = t ′ , where t ′ = T , if X < X min and T − t ′ = α −1 ln(X/X min ), if X > X min . As discussed in Sec. III, X min is the minimum value of x at which state selection can be said to have been completed, and therefore we ask that X ≥ X min . The maximum value of y c (t), increases rapidly with X and is already extremely large when X ∼ 5X min . It might therefore be tempting to argue that we need to take X to be much less than this value if the linearization procedure is to be valid. Since (33) has its least value when X = X min , this approximation would seem to be best for this choice of X, and in fact this was the choice made in our earlier work [1]. In fact, there are a variety of reasons why choosing X to be X min is not suitable. First, the argument above, namely that if X is too large then the linearization procedure is invalid, is in fact incorrect. We will find that the width of the distribution is so narrow in the y direction that the range of integration for y f need only be tiny -in fact only going out to values of y f such that y 2 f exp γX 2 /α ∼ 1. From (33) we see that under these conditions y c (t) is indeed small. We will come back to this point again in Secs. VI and VII.
The second reason is in fact more profound. Intuitively, we expect that for large enough values of X the probability flux through a given valley should be independent of X. This is due to the fact that in the regime of weak noise, once state selection has occurred there is no flux leakage out of a valley. This implies that we should not have to make any choice for X because our results should be independent of X as long as X is large enough. In fact, this is precisely what we find, namely that the probability flux through the x-valley tends to an asymptotic value as X tends to infinity. The remarkable cancellations that occur to produce an X independent flux is a good indicator of the correctness of our calculational scheme. Again, we will come back to this point in Secs. VI and VII.
Finally, throughout our calculation we use the approximation that T is large in the sense that e −αT ≪ 1 and e −βT ≪ 1. This is apparently a problem, since we need the form of the distribution for all T in order to perform the time integral of the probability flux in the calculation of the total flux through the valley. The way out of this impasse is to take X sufficiently large that the probability current is essentially zero at small times, and only starts making appreciable contributions to the T integral when T is such that e −αT ≪ 1 and e −βT ≪ 1. It is not clear how large X will have to be, but as discussed above, the probability flux through the x-valley tends to an asymptotic value as X grows, so we are assured of always finding a value of X for which the approximation that T is large is valid.
At first sight, the large T approximation might appear to be problematic because it ignores the shorter state selection time scales α −1 and β −1 . However, this is not the case. The reason is that the integrals over time that are performed in the calculation of the action and Jacobian prefactor include these earlier times. In other words, the entire classical path, in particular the spatial region near the origin, is included in the calculation.
After these rather technical asides, it is worthwhile summarizing what we have deduced about the equation for the optimal path which leaves (0, 0) at t = 0 and arrives at (X, y f ) at time t = T . It is given by (31) and (32) under the assumptions that (i) T is such that e −αT ≪ 1 and e −βT ≪ 1, (ii) X is larger than a few X min , and (iii) t m and (T −t m ) are large (we will later find that we require that e −2αtm ≪ 1, e −2βtm ≪ 1, e −2α(T −tm) ≪ 1 and e −2β(T −tm) ≪ 1).
We are now in a position to calculate the leading order contribution, P 0 (X, y f , T ), to the conditional probability distribution (21), by finding the action of the optimal path given by (31) and (32). We begin by noting that, from (26), the classical paths are solutions ofẍ = −∂U/∂x andÿ = −∂U/∂y. Multiplying the first equation byẋ, the second one byẏ, adding them and integrating the result, gives 1 2 (ẋ 2 c +ẏ 2 c ) + U = E, a constant. Thus another form for the action of classical paths is Details of the calculation are given in Appendix B. From (B20) we find where S y = 1 2α The two terms in (35) have different characters: as remarked earlier, we will eventually integrate over y f to obtain the probability flux through the valley, but the first term will remain, giving a leading order contribution of exp −(αX 2 /4D)(cothαT − 1) to (21). When performing similar calculations to find the escape rate from one metastable state to another, the leading order result e −∆V /D gives a reasonable estimate, and the need to go to next order and to calculate the prefactor only arises if increased accuracy is required. The situation is very different in the present case of the decay from an unstable state. As mentioned already, we expect that the main contribution to the probability flux through the valley will occur when T ∼ (2α) −1 ln(αX 2 /D). The conditional probability distribution (21) is peaked around a value of T of this order because of a balance between the leading (classical) term and the prefactor (a fluctuational) term. For this reason, the calculation of the prefactor in this problem, unlike in escape problems, is vital if the essential structure of the state selection probabilities is to be captured. We therefore turn to the calculation of this quantity.

V. CALCULATION OF THE PREFACTOR
In this section we will calculate the next to leading order contribution in (21). From (22) we see that it consists of two distinct contributions: from the Jacobian J[ r] evaluated at the optimal path, and from the fluctuations around the optimal path which give rise to the determinant of the second functional derivative of the action (23). We begin by evaluation of this determinant.
Starting from (26), and using the notation (23), we obtain where, in line with our previous assumption, we have neglected terms which are O(y 2 c ) down on the terms shown in (37). We have omitted any time dependence in (37) for the sake of clarity, but it should be remembered that the classical solutions x c and y c are functions of t and the matrix is multiplied by an overall factor of δ(t − t ′ ). The diagonal entries are as in the y c = 0 case [1] -only the off-diagonal entries are different. However, it is easy enough to see that the eigenvalues of (37) are even in y c , and since the matrix with y c = 0 has no zero eigenvalues, the off-diagonal entries only provide O(y 2 c ) corrections to the eigenvalues of the matrix with y c = 0. Therefore neglecting these terms as before, we conclude that within the approximation we have adopted in this paper, we may set the off-diagonal entries in (37) to zero. From (23) and (37) we now find where L x and L y represent fluctuations in x and y respectively: To evaluate the determinants of the operators (39) and (40) we make use of a well-known formula for such determinants [28]. Let L be an operator of the type −d 2 /dt 2 + φ(t) and suppose that the eigenfunctions of L are required to vanish at the boundaries at t = 0 and t = T , as in the case of interest to us here. Now let h 1 (t) and h 2 (t) be two independent solutions of the homogeneous equation Lh = 0. Then The constant of proportionality in (41) may be omitted for the same reason the overall constant was omitted in (20): the probability of ending up in the positive x-valley will be normalized by the sum of the probabilities of ending up in the x-and y-valleys. Furthermore, if h 2 is the solution which vanishes at t = 0, then the formula for DetL involves only this solution: The solution of L x h = 0 which satisfies h(0) = 0 is h 2 (t) = F sinh αt, where F is a constant. Therefore To find DetL y we need only note that the solutions of the homogeneous equation L y h = 0 are also the solutions of the classical equation (30). But the solution which vanishes at t = 0 has already been found in Appendix B: it is given by (B3) and (B7), where G is a constant. Sinceḣ 2 (0) = βG, use of (42) yields, It only remains to calculate the Jacobian at the optimal path. From (18) this is given by Neglecting the O(y 2 c ) term as before, one finds that J( r c ) = J x J y where since we are assuming that e −αT ≪ 1. Since both the determinants and the Jacobians factorize into contributions associated with the x coordinate and contributions associated with the y coordinate, we may summarize the results of this section so far as It is straightforward to check that these results agree with those of Ref. [1] under the assumptions that e −αT ≪ 1 and e −βT ≪ 1. While neglecting corrections to the leading behavior of the prefactor which are of the type e −2αT or e −2βT is in line with this large T approximation, we will see later that the final integration over T will also provide a posteriori justification for the neglect of these terms. From (21), (35), (49) and (50), we can now write down the conditional probability for the system to be at r f = (X, y f ) at time T , given that it was at the origin at time t = 0, as As we discussed in detail in Secs. II and III, when using the modified potential (8) the question is not what are the paths of least action from the origin to the minima of the potential, since those minima cease to exist when the limit δ, ǫ → 0 is taken. Instead, the question is what is the relative flux through one valley compared to the other. To calculate the flux, it is not the conditional probability distribution (51) that we need to know, but the probability current J = (J x , J y ). This current is related to the conditional probability distribution through the Fokker-Planck equation [29] ∂P ∂t where In order to find J from P given by (51), we first note from (21) that Therefore we only need to differentiate V and P (0) with respect to X and y f in order to find J . Carrying this out we obtain In the next section we will calculate the flux through the positive x-valley. This will involve only the component of J normal to the line x = X. Therefore, the y-component of the current, J y , will not contribute and need not be considered any further; the entire contribution to the flux will come from integrating J x , given by (55), over y f . By virtue of the y 2 f term in the exponential in (51) this will be a Gaussian integral. While we would naturally neglect the O(y 2 f ) term in (55) in line with previous approximations in this section, we now see that it would in any case give a contribution of order D once the y f integral is performed. This is a further justification for ignoring such terms. Finally, we have been neglecting e −2αT type corrections in the prefactor, which means that we should replace cothαT by 1. These approximations lead to the result J x = αXP and so we find from (51) This is the key result of this section. We need now only use it to calculate the total flux through the positive x-valley at x = X and compare it with the analogous quantity through the positive y-valley. This will give us the relative probability of the x-mode being selected.

VI. RESULTS
The idea underlying the method we use to calculate the relative probability of one of the states being selected is most easily understood by first describing the analogous procedure used when carrying out Monte Carlo simulations. As was briefly alluded to in Sec. III, for each run the particle starts at the origin and subsequently follows the Langevin dynamics of the reduced problem until it reaches x = ±X in one of the x-valleys or y = ±Y in one of the y-valleys. We used X = 5X min and Y = 5Y min to be certain that state selection had occurred. At this point that particular run ceases. If the particle has ended up in a y-valley, for instance, 1 is added to the total number of runs which have selected the y-mode. Another run is initiated and the result of that is added to the totals. After a large number of runs the proportions selecting the x-mode and y-mode are used to calculate N x and N y , the relative probabilities of these states being selected.
Let us focus on the runs ending up in the positive x-valley, exactly as we have been doing in the analytic treatment in Secs. IV and V. For each of these runs the final value of y -called y f in the analytic treatment above -will vary. We would expect most of the runs to end near to the x-axis, with off-axis end-points becoming less and less common as we move away from the axis. This is reflected in the y f dependence of J x given by (57). Just as in the Monte Carlo simulation, where we add up all contributions with differing final y-coordinate at x = X, so we need to integrate (57) over all y f at x = X. Similarly, just as in the Monte Carlo simulation, where we add up all of the contributions, no matter how long they took to get to the end-point, so in the analytic treatment we have to integrate over all T to obtain the total flux. Therefore we need to calculate in order to calculate the state selection probabilities. While the above justification for (58) as the quantity which we need to calculate seems intuitively plausible, it is worthwhile formally proving this. We first need to specify what "ending up in the positive x-valley" means in terms of a mathematical expression. Since, in the limit T → ∞, trajectories will have entered one of the four valleys, we define the probability that it has entered the positive x-valley as Here the positive x-valley is defined as the entire potential surface to the right of the line x = X. On the other hand from the continuity equation (52) lim T →∞ Integrating (60) over the "volume" {x ≥ X, −∞ < y f < ∞}, using that fact that at T = 0 the end-point is at the origin and making use of the divergence theorem, gives The "surface" integral on the right-hand side of (61) only gives a contribution at x = X, where d S is in the direction of the outward normal, i.e. in the negative x direction. Thus J .d S = −J x dy f and so the required probability (59) is equal to (58). Substituting (57) into (58) and performing the y f integration yields As remarked in Sec. V, any O(y 2 f ) corrections in the prefactor give O(D) corrections to (62), as do O(y 4 f ) corrections to the action, which justifies their omission. The final integral over y f also provides justification for the linear approximation. An examination of J x (X, y f , T ) shows that this function is effectively non-zero only for small values of |y f |, which get still smaller as X increases, so that unless |y f | < ∼ exp −γX 2 /2α , the flux is essentially zero. Thus either the conditional probability distribution or the conditional probability current resemble a thin wafer centered on the x-axis when plotted at fixed T . Moreover, this wafer gets very much thinner with increasing X. This means that the limits on the y f integral in (58) are effectively ± exp −γX 2 /2α , which is tiny for large X. It is more instructive to change variables from y f to y * f = y f exp γX 2 /2α , so that the integral now has the range (−1, 1). But, as discussed in Sec. IV, this means that (33) now reads which has a magnitude less than one for −1 < y * f < 1. Thus the linearization procedure is justified. In some sense, the variable y * f is more appropriate in this problem than y f itself. An examination of the integrand in (62) shows that is has a maximum when T = T * ≡ (2α) −1 ln(α 2 X 2 /[(α+ β)D]), which is slightly different to the value mentioned in Sec. IV. This is because that value came from the estimates of Sec. III and Appendix A, which used the value for the maximum which is found in the analogous one-dimensional stochastic process in x(t). The more refined calculation we have been discussing here has an additional e −βT contribution coming from fluctuations in y(t) and from J y . It is the occurrence of this term in the integrand of (62) which changes the value at which the maximum occurs by this very slight amount. From (62) we can also check the conjecture made in Sec. IV, that if X is greater than a few X min , then even the tails of the integrand will be at values of T which are large enough for the approximations made in this paper to be valid. A numerical investigation of the integrand shows that it has typically already fallen by two orders of magnitude from its maximum value for T satisfying e −αT , e −βT ≪ 1, when X is larger than a few X min .
To perform the integration in (62) we (i) replace (cothαT − 1) by its large T form, 2e −2αT (justified below), and (ii) change variables to ξ = (αX 2 /2D)e −2αT . Then the integral becomes where Ξ = (αX 2 )/2D. The integral is an incomplete gamma function: it equals Γ([α + β]/2α) plus exponentially small corrections in D coming from the large, but finite, upper limit Ξ [27]. Neglecting these exponentially small terms we therefore obtain for the total flux through the positive x-valley at X By using the approximation (cothαT − 1) ≈ 2e −2αT , terms of order D −1 e −2(n+1)αT , (n = 1, . . .) were neglected, but these are of the order of D n ξ n+1 and give contributions which are O(D) down on the leading term (65). Similarly, any corrections to the prefactor of order e −αT or e −βT would also give contributions which are down compared to (65). Once again, we see that these final integrals justify approximations which we made earlier.
The total flux given by (65) depends on X, which is still undetermined, apart from the requirement that it should be greater than a few X min . However, if we investigate the X dependence of the expression in (65), we find that as X increases from X min , it first decreases until X equals two or three times X min , where it reaches a constant value and subsequently remains at this value as X → ∞. Of course, this constancy of the flux is exactly what we would expect: once the system has selected the positive x-valley, it remains in that state and so the flux should be conserved within the valley, i.e. be the same for all X. We will denote this constant value simply as F x ; it is given by the X → ∞ limit of (65): The integral in (66) may be expressed in terms of gamma functions: it equals αβΓ(β/α). Therefore The result (67) only depends on the four parameters α, β, γ and D, as we would expect. An exactly analogous calculation of the total flux through the y-valley at Y gives the same form but with α and β interchanged: Throughout the analysis we have ignored constant factors which multiply both the results coming from paths ending up in the x-valley and those ending up in the y-valley. The reason given was that eventually we would normalize the probabilities of state selection and therefore common factors were irrelevant. We are now at the point where we can impose this normalization, but before doing so, let us put in the constant multiplying factor explicitly, and in doing so absorb the additional factors of 2D (2π) 1/2 which appear in (67) and (68) into it. We therefore the final forms for the total fluxes along the x-and y-valleys as where K is the overall constant. The relative probability of flowing into an x-valley has already been introduced in Sec. III. In terms of the total fluxes it is given by with N y = 1 − N x . From Eqns. (69), (70) and (71) we obtain: whereD Equations (72) and (73) are the desired results. Note that the only place where the coupling constant γ enters is multiplying the noise strength D. It is simple to see that this will always be the case since, if x and y are measured in units of X min and Y min respectively in the reduced problem, the only place where γ appears in the Langevin equation is multiplying D. Thus the effect of the interaction, specified by γ, is to renormalize the noise. This means that the probability of either the x-mode or of the y-mode being selected depends only on the three quantities α, β and γD. It is clear that N x and N y given by (72), have the right behavior in various limits. First of all, if α = β they are equal. If α and β are not approximately the same, the relative magnitudes of N x and N y are largely governed by their dependence onD: if α > β, then σ > 0, and N x < ∼ 1 while N y < ∼D σ ≪ 1. A more detailed comparison with the results of Monte Carlo simulations is given in Fig. 6. The final result (72) is in good, but not perfect, agreement with simulations. We believe that the only restriction on our method is that D should be small, so that the method of steepest descent is appropriate. For D = 0.01, the neglect of O(D) corrections to (69) and (70) should mean that our result is within approximately 1% of the simulation results. It is not clear from an examination of Fig. 6 whether the slight discrepancy is due to the O(D) corrections or to additional contributions which have not been accounted for. Since the agreement is still not perfect for smaller D, we should take the latter explanation seriously and carry out a search to see if we can find such additional terms. A possible source is the following. We have conjectured that the solution of (24) and (25) which has the least action can be obtained first by finding the x c (t) solution with y c = 0, and then linearizing about it. However, there are undoubtedly higher action solutions of (24) and (25) which have a non-trivial structure in x and y, that is, solutions which cannot be obtained as a simple expansion about the y = 0 solution. In other applications of the steepest descent method, we would neglect these solutions as they would give contributions which were exponentially smaller than the contributions from the least action solutions, and thus be completely negligible for small D. However, in this problem, after integrating over T , the previously exponentially small contribution becomes a power law (see Eqn. (64)), or more precisely, D divided by the a quantity related to the action of the solution, all raised to a power. Although we do not know what power the higher action solutions will be raised to without further analysis, it is clear that solutions to (24) and (25) with, for instance, twice the action of the solution considered in this paper, have at least the potential of changing our result by a not insignificant amount. Thus it appears that the neglect of higher action solutions, which is common in most problems involving optimal paths, may not be so well-founded in the situation which we are considering in this paper.

VII. CONCLUSIONS
In this paper we have presented the motivation, justification and calculational details of a method for determining the probabilities of various metastable states being selected for occupation, when an unstable state decays. Although, as we hope we have made clear, this is a widespread phenomenon, theoretical progress seems to have been hampered by a lack of sufficiently powerful tools with which to attack the problem. As a consequence, most studies were carried out some time ago on one-dimensional systems (or quasi-one-dimensional systems). Yet, these are actually the systems of least interest. The really interesting aspects of the problem arise in two or more dimensions, where the system can perform Brownian motion near an unstable state, and then select a final state through a mixture of randomness (noise) and dynamics (which will necessarily be non-linear in all non-trivial cases). It is this interplay of randomness and determinism, both essential to the process, which makes the problem so difficult.
Our main purpose in this paper has been to show how the use of path-integral techniques, and especially the notion of optimal paths, can be used to analyze this technically hard problem. In order to stress that relatively naive approaches are not able to capture much of the subtlety of this problem, we carried through what seemed to us the most simple-minded approach to the problem in two dimensions, and showed that it was not able to give satisfactory results. Although the extent of the agreement with numerical work is always open to interpretation, there were two features which signified that the naive method was defective. Firstly, there was no systematic way to proceed with the calculation; one was faced with ad hoc choices which had to be made at several junctures in the calculation. Secondly, the crude approximation (10) to the already unsystematic result (12), actually gives better agreement with simulations than (12) itself, which indicates the dubious nature of these results. It should also be borne in mind that the prediction that N x = 1/2 when α = β is assured by symmetry and that in almost all schemes N x → 1 as α gets very much larger than β, simply because the noise will have very little effect in this case. If we add to this the expectation that N x (α) will be a smooth function, we see that even the most simple minded approximation cannot be expected to be too far out. These observations indicate that we should focus more on the systematics (or lack thereof) of the calculation, rather than getting perfect agreement with simulation data.
The nature of the optimal paths used in the treatment presented here is different from those in the more familiar problem of escape from one potential well to the other. In that situation it takes an exponentially large time (exp ∆V /D, where ∆V is the height of the barrier), to make the transition from one potential well to another, and this is reflected in the fact that the optimal paths are essentially of infinitely long duration (T → ∞). In the problem studied in this paper we are focusing on phenomena on much shorter time scales: the decay of a state with no barrier to a metastable state. Any subsequent process (presumably a noise-activated escape to a different well, of the kind we have just discussed) is of no interest to us here. Actually, as explained in Sec. II, we do not even ask that the end point of the optimal path is a metastable state, but only that it lies somewhere in the valley which leads to the metastable state in the full problem. The optimal paths in these cases are of finite duration. To achieve this, the initial velocity of the particle in the mechanical analogy has to be non-zero at the origin. As T increases this initial velocity decreases, and eventually tends to zero if T tends to infinity.
Many of the other details of the calculation involving optimal paths turn out to be remarkably subtle. For example, the use of the reduced potential simplifies the problem, but eliminating the actual minima means that we have to specify an arbitrary final value of x (if we are considering state-selection into the positive x-valley), which we denote by X. However, the probability of the positive x-valley being selected should not depend on X, as long as X has a value which lies in the valley. Reassuringly, this is what we find in our calculation, but it requires us to take into account the y-dependence of the path for this "common-sense" condition to be observed. In our earlier treatment [1], we sought to give only the outline of our method, and made the assumption that the total probability flux through a particular valley was proportional to (a) its value on the axes, and (b) its maximum value, which occurs at T = T * -the time for which P (X, 0, 0| 0, 0) is a maximum at a fixed X. Within this scheme there was no way of determining X and Y , and so we fixed them to our simulation results. The best values were X > ∼ X min and Y > ∼ Y min , in line with our expectations. However, we have found here that it is possible to get a form for the y-component of the optimal path only if we assume that X is sufficiently large -which comes from the requirement that T must be sufficiently large -and so setting X ≈ X min and Y ≈ Y min is not possible in the approach we have adopted here. Needless to say, we expect that if we were able to find the full solution for the optimum path, then we would not need to constrain X or Y to be large: the expressions for the fluxes F x and F y would be valid for all X ≥ X min and for all Y ≥ Y min . Presumably the fact that F x decreases from its value at X = X min , and only becomes constant when X equals a few X min , is not due to a violation of flux conservation (this would necessitate some flux leakage, and it is difficult to see where it would go to), but instead to the inapplicability of our calculational scheme for small X, which simply gives the incorrect expression for F x . It should be stressed that this is not a problem for our method: we simply take X and Y large enough so that F x and F y respectively, do tend to asymptotic values. These are the required values and are those which we use in our determination of N x and N y .
An even more subtle aspect of the calculation concerns the applicability of the linear approximation in determining the y-component of the optimal path. Although at first sight y c (t) is not small for large X, it turns out that the range of final values of y c (t), denoted by y f , decreases so fast as X increases, that the small y c (t) approximation is still valid. Put more intuitively: if T is large, the path has to take a wide detour to large y if it is to take sufficiently large time to reach (X, y f ) (recall that T is fixed for these paths). But "large y" means large relative to y c (T ) = y f , and, as discussed earlier, the current J x (X, y f , T ) is effectively non-zero only for tiny values of |y f |. So while it is true that on the scale of y f the excursions of the optimal path spread deep into the plane, in fact the range of y f is so minute, that the linearization assumption remains valid. The result that the probability is so concentrated about the axes is used in other places in the calculation. For example, in Sec. VI we argued that the "volume" integral in (60) could be replaced by the "surface" integral in (61), and then we proceeded to ignore some of the contributions from the "surfaces" in the later integral. We specified the volume in the case of the positive x-valley as the entire plane to the right of the line x = X. Actually, this is not quite correct, since we need to define a similar volume for the y-valleys, and these will overlap in the sectors where x and y are both large. We can be reasonably vague about how precisely to define these volumes and surfaces simply because J x (and J y ) fall away so fast that they have an utterly negligible contribution however we define them. In any reasonable definition, the only contribution to the surface integral (61) will be on the line x = X, very close to the x-axis.
We believe that one of the best indicators of the correctness of our approach is the cancellation of the various factors of X coming from several different sources: (a) the fluctuations about the y c (t)-solution, (b) the action of the y c (t)-solution (after integration over y f ), (c) the action of the x c (t)-solution (after integration over T ), and (d) the definition of the current in terms of the probability distribution. What remains after these cancellations is essentially independent of X, as would intuitively be expected. As we have repeatedly stressed, many of these contributions come from assuming that T is large. This seems somewhat paradoxical, since we would expect that state selection is determined at earlier times. This can be understood if one realizes that the large X calculation determines the normalization of the flux in (65), or equivalently, the quantity ζ(α, β) in (72). To obtain a finite X → ∞ limit, it is essential to include a non-trivial y c (t) optimal path in the analysis. If the optimal path is taken to be y c (t) = 0, then only theD σ factor in (72) is determined. This was essentially the approach we adopted in our earlier paper [1]. While we believe that the treatment given here is a great improvement on that reported in Ref. [1], there is no doubt room for further improvement. For instance, the linearization approximation seems very reasonable on physical grounds, but it would be useful to put it on a sounder mathematical footing. A deeper understanding of the origin of the small scale set by exp (−γX 2 /α) would also be valuable.
The competition between new modes, when a previously stable mode becomes unstable, is responsible for much of the emergent order found in systems far from equilibrium. The continual branching to more complex structures that this entails [30], is not only governed by the dynamics of the process but also by the random fluctuations, or noise, generated by the large number of other degrees of freedom of the system not explicitly included in the description. The model which we have investigated in this paper is, as has already been emphasized, the simplest showing enough of the complicated features of this process of state selection to serve as an illustrative example of our method. In forthcoming papers we will show how the method can be applied to more complex examples, such as state selection in lasers and the study of population dynamics in a fluctuating environment. We believe that the ideas and techniques developed here will not only be applicable to these situations, but to many others where multiple states compete for occupation.

ACKNOWLEDGMENTS
AM wishes to thank the Department of Physics at the University of Illinois at Urbana-Champaign, where part of this work was carried out, for hospitality, and EPSRC for financial support under Grant No. K/79307. MT gratefully acknowledges support from the MRSEC program of the NSF under award number DMR 9808595.

APPENDIX A:
In this appendix, we explore the consequences of the elementary theory put forward at the end of Sec. III. We will obtain simple expressions for the probability that the particle ends up in an x-valley, which are compared to Monte Carlo simulations in Sec. III.
We begin with the first ingredient in this simple theory. Since γ = 0 we have only to solve the Langevin equationṡ x = αx + η x (t) andẏ = βy + η y (t). Such linear problems can always be exactly solved. Since x and y are linearly related to η x and η y , they are also Gaussian random variables. It is easy to show that x(t) = y(t) = 0 and therefore the probability that the particle is at (x, y) at time t, given that it started at the origin at time t = 0 is P (x, y, t|0, 0, 0) = N exp − x 2 2 x 2 (t) − y 2 2 y 2 (t) , where N is a normalization constant. For later use we will also need To quantify the second ingredient, first suppose that the coordinates of the particle (x, y) are in the first quadrant. If the angle it makes with the positive x-axis, tan −1 (y/x) is less than tan −1 (Y min /X min ), then according to our criterion, it goes into the x-valley. Since the tan function is monotonic in the interval (0, π/2), and generalizing to other quadrants by using |x| and |y| rather than x and y, the criterion becomes: These ingredients can now be put together. The question reduces to how often the particle is in the sector specified by (A3) at a given time t and how often it is in the sector specified by (A4) at this time. This probability is found by integrating (A1) over all x and all y with the constraint that it satisfies either (A3) or (A4) depending on whether we want the probability of ending up in the x-valley or y-valley. Restricting ourselves to the positive quadrant, which we can obviously do on symmetry grounds, the probability of ending up in an x-valley, N x and in a y-valley, N y , are − y 2 2 y 2 (t) and (A5) where the factor 4 comes from including the contributions from the other quadrants and λ = Y min /X min . In the N x integral, change variables from y to z where y = λxz at fixed x. The x integral can now be easily performed, giving where κ 2 (t) = λ 2 x 2 (t) y 2 (t) = Y 2 min x 2 (t) X 2 min y 2 (t) .

(A8)
In a similar way From (A2) we see that if α and β are not very different from each other, then κ(t) will be neither very large nor very small, and we can ignore the integrals in (A7) and (A9) as a first approximation. Then determining N by asking that N x + N y = 1 and ignoring the factor 1 compared to the exponentials in (A2) yields N x = e 2αt e 2αt + e 2βt ; N y = e 2βt e 2αt + e 2βt .
There are several simple arguments which give formulas resembling this. It is not difficult to improve on this by including the integrals in (A7) and (A9). Once again determining N by normalization gives where we have used the result tan −1 (w) + tan −1 (1/w) = π/2. If x 2 (t) ≫ y 2 (t) or y 2 (t) ≫ x 2 (t) , which can happen even when α and β are not too different, we may replace (A11) by a simpler form which resembles (A10), but with constants multiplying the exponentials. The simple result (A10) does not depend on the parameter γ, which represents the non-linear part of the potential, or on the noise strength, D. This is also true of (A11) if we ignore the factor 1 compared to the exponentials in (A2). To understand this, let us suppose that α > β. From (A10) or (A11), we see that when t is small although there is a slight bias in favor of the x-valley over the y-valley, it is not that marked. However, as t increases the asymmetry becomes more marked, until at large times the particle has only a very small chance of being in the y-valley. We need to simulate the role that the saddle points play in state selection by choosing the time to have the value t c , which is the time at which it is most likely that the x and y coordinates of the particle will have magnitudes X min and Y min respectively. This is implemented by asking that x 2 (t c ) = X 2 min or y 2 (t c ) = Y 2 min . Note that we have used the word "or", since only one relation is required to determine t c , and the two relations given give incompatible values for t c unless α = β. Once again, as is common in many aspects of the naive approach to state selection which we have been summarizing in this appendix, there is a considerable degree of arbitrariness in the choice of t c . We will assume that for α > β, when the particle is most likely to go into the x-valley, t c is determined from the condition on x 2 (t c ) and vice-versa. Then from (A2), which has the characteristic form − ln D as D → 0. If we substitute (A12) into (A10) we obtain the probability of ending up in an x-valley. In terms of the quantitieŝ D ≡ γD αβ and ρ = (α − β) α , we obtain If we substitute (A12) into (A11) we obtain the more general form Equations (A14) and (A15) are the principal results of this appendix. The two main points we wish to make is (i) their derivation is somewhat arbitrary -there are several other similar assumptions we could make which would give us slightly different formulas, and (ii) although, as is discussed in Sec. III, they show the right qualitative features, they are are not in good agreement with Monte Carlo simulations. Thus a more systematic and sophisticated theory is required. This is given in Secs. IV, V and VI of the paper.

APPENDIX B:
The purpose of this appendix is to explore the solution of (30), where x(t) is given by (31) and where y(0) = 0 and y(T ) = y f . Throughout we will assume that T is large in the sense that e −αT ≪ 1 and e −βT ≪ 1, which will allow us to obtain a rather simple form for the solution.
Clearly the zeroth-order solution which satisfies y(0) = 0 is y 0 (t) = A sinh βt, where A is an arbitrary constant. It is straightforward to determine y 1 (t), but since a short analysis shows that y 1 (t)ζ ≪ y 0 (t) for all t < ∼ T /2, we conclude that, it is sufficient to take the small t of the solution of (30) -which we denote by y < -to be