{'fraction_non_alphanumeric': 0.03350227420402859, 'fraction_numerical': 0.012267706302794022, 'mean_word_length': 4.6210372534696855, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Doping front instabilities in organic semiconductors: a means for optimizing optoelectronic devices

Recently, it was demonstrated that electrochemical doping fronts in organic semiconductors ex- hibit a new fundamental instability growing from multidimensional perturbations [Phys. Rev. Lett. 107, 016103 (2011)]. In the instability development, linear growth of tiny perturbations goes over into a nonlinear stage of strongly distorted doping fronts. Here we develop the nonlinear theory of the doping front instability and predict the key parameters of a distorted doping front, such as its velocity, in close agreement with the experimental data. We show that the instability makes the electrochemical doping process considerably faster. We obtain the self-similar properties of the front shape corresponding to the maximal propagation velocity, which allows for a wide range of controlling the doping process in the experiments. The theory developed provides the guide for optimizing the performance of organic optoelectronic devices.


I. INTRODUCTION
The discovery and development of π-conjugated polymers -plastics with properties of semiconductors, often termed organic semiconductors (OSCs) -has revolutionized the field of electronics, and was awarded with the Nobel Prize in Chemistry in 2000 [1][2][3] . The OSCs, in contrast to crystalline inorganic semiconductors, are simple to process and relatively inexpensive in production, while their soft and conformable character enables large number of novel technological applications ranging from optoelectronics and photovoltaics to sensors and artificial muscles [4][5][6][7][8][9] . The unique combination of electronic and ionic conductivity in OSCs introduces the exciting functionality that triggers experimental and technological activities in the nascent field of so-called organic (i.e. based on OSCs) electronics. At the same time, as formulated in Ref. 8 , "the understanding of the fundamental processes that take place [in OSCs] is still in its infancy", especially if compared to the theoretical fundament of the field of crystalline inorganic semiconductors 10,11 . The inadequate understanding of fundamental OSC properties limits the technologies to "cut-and-try" methods and demands development of the theory on the subject. In line with this demand, the goal of the present work is to attain better fundamental understanding of the doping process in OSCs, which may provide a guide to optimizing performance of optoelectronic devices such as, for example, light-emitting electrochemical cells.
In the process of electrochemical doping the OSCs can be transformed from an almost insulating state to a metallic-like state with the conductivity increased by several orders of magnitude 4,12,13 . Apart from conductivity, this doping transformation also allows for tuning other key properties of OSCs, such as color, volume, surface energy and photoluminescence capability, in a controllable way. This possibility can be profitably employed in numerous applications, e.g., in optoelectronics [3][4][5][6][7][8] . Operationally, the electrochemical doping can be performed in-situ by applying a potential difference between two metal electrodes in contact with an OSC, blended with a solid-state electrolyte, as illustrated in Fig. 1 (a). Because of the applied electric field, positive and negative ions redistribute in the active material and allow the injection of electrons (n-doping) and holes (p-doping) into the OSC from the cathode and anode, respectively [14][15][16][17][18] . As the electronic charges populate the OSC, its conductivity increases by several orders of magnitude in comparison with the initial conductivity supported by ions. It has been demonstrated experimentally that the doping process develops in the form of two counter propagating fronts illustrated in Fig. 1 (a), which separate the regions of p-and n-doped material from the undoped one [17][18][19] . When the two fronts meet, a p-n junction is formed in which the recombination of holes and electrons results in the visible light 4,[17][18][19][20][21] . Thus, the velocity of doping front propagation determines the turn-on time of the light-emitting electrochemical cells.
As pointed out above, the present understanding of the doping process in OSCs is far from being adequate, with only few theoretical papers focusing on the subject 19,[22][23][24][25] . In particular, the analytical formula for the velocity of a planar doping front has been found only recently in Ref. 25 , which demonstrated good agreement with the experimental observations. The analytical theory has shown that the doping rate is controlled by the mobility of heavy ions that results in a slow turn-on time of light-emitting electrochemical cells. At present, one of the most difficult problems to solve is to speed-up the doping process and to reduce the turn-on time. The speed-up could be achieved with the help of a new basic phenomenon -the multidimensional doping front instability -discovered recently in OSCs in Ref. 19 . This new fundamental instability distorts the shape of doping fronts, increases locally the electric field ahead of the fronts and hence increases their propagation velocity. The recent paper 19 has developed a linear theory and provided the experimental evidence of the instability; it has also proposed a possible way to amplify the instability by means of minor humps introduced at the electrode surfaces. Figure 1  Ref. 19 . Still, the dynamics of the strong instability at the nonlinear stage remained unclear and demanded deep theoretical studies.
In the present paper we develop the nonlinear theory of the doping front instability, focusing on the velocity and shape of the corrugated fronts. For the 2D geometry of the OSC films, the theory predicts increase of the doping front velocity by a factor of two in comparison with the planar front velocity, which is in a good agreement with the previous experimental data. We obtain self-similar properties of the front shape corresponding to the maximal velocity, which promises a wide range of controlling the doping process in the devices, e.g., through the design of the electrodes. We show that moderately corrugated doping fronts are preferential for achieving higher velocity as compared with the strong corrugations obtained in Ref. 19 and illustrated in Fig. 1. The developed theory provides a guide for optimization of organic optoelectronic devices utilizing electrochemical doping.
In addition, it is expected that doping fronts in OSCs will become a new popular object of nonlinear science alongside with other fronts and interfaces in fluids, gases and plasmas [26][27][28] . In particular, in the present work we discuss similarity and difference between the doping front instability in OSCs and the hydrodynamic instability of a flame, i.e. a slow combustion front [29][30][31] . The electrochemical doping process provides a fascinating opportunity to observe directly and investigate a number of interesting nonlinear phenomena: front propagation, multidimensional instabilities and pattern formation, with surprisingly different outcome for p-and n-fronts. In addition, the present analysis suggests different regimes of the doping front instability for the 2D geometry of OSC films and for the 3D geometry of the bulk optoelectronic devices 18 : thus, the nonlinear dynamics of the instability depends also on the system dimension.
The paper is organized as follows: In Sec. II we remind basic results on the linear stage of the doping front instability needed for the subsequent nonlinear analysis. In Sec. III we derive the nonlinear equation for the p-doping front bent by the instability. In Section IV we analyze stationary periodic solutions to the nonlinear equation. In Sec. V we compare the theoretical results to the experimental data and discuss technical implications of the theoretical predictions. The paper is concluded by a brief summary in Sec. VI.

II. LINEAR STAGE OF THE DOPING FRONT INSTABILITY
In this section we remind basic results of Ref. 19 on the doping front instability at the linear stage, when perturbation amplitude is infinitesimal. Details of the linear analysis are needed to understand the nonlinear theory, which we develop below in Sec. III and IV. The instability arises because speed of a doping front depends on the electrical field ahead of the front as it was derived in 25 . We present the p/n-doping front velocity in the vector form as Here E = −∇φ is the electric field intensity just ahead of the front, φ the electric field potential, n 0 the initial ion concentration in the active material, n h,e the final concentration of the holes/electrons behind the front determined by the thermodynamic properties of a particular OSC, µ + and µ − the mobilities of the positive and negative ions, respectively. The "±" sign in Eq. (1) indicates that the p-type (+) and n-type (−) doping fronts propagate in opposite directions. We start our analysis with a single p-doping front propagating stationary with velocity U p in a constant uniform external electric field E =â x E 0 pointing along x-axis, whereâ x is the unit vector in x-direction. The analysis of the n-front instability may be performed in a similar way. In the case of a planar stationary front, the general 2D equation for the p-front position x = X p (t, y) reduces to X p (t) = U p t. Next, we consider small perturbations of the stationary p-doping front whereX p ∝ exp(σt + iky) stands for one Fourier mode at the linear instability stage, k = 2π/λ is the perturbation wave number characterizing infinitesimal bending of the doping front (λ is the perturbation wavelength), σ is the instability growth rate, see Fig. 2. The purpose of the linear stability analysis is to find the dispersion relation σ(k). In the case of Re(σ) > 0 valid for a certain domain of perturbation wave numbers we conclude that the planar front is unstable with respect to multidimensional perturbations. Small bending of the doping front induces also perturbations of the electric field in the undoped regionφ(t, y, x) ∝ exp(σt+iky). The doped region may be treated as equipotential due to high conductivity. Thickness of the doping front may be characterized by a length scale L p , of about 10 −4 −10 −3 mm, determined by ionic diffusion 23,25 . Since the characteristic size of the experimentally observed p-front perturbations (approximately 0.2mm) is much larger than L p , we may treat the front as infinitesimally thin. Most of the time the doping fronts are sufficiently far away from each other, k(X n − X p ) ≫ 1, therefore, the instability of one front is not affected by the other. The linearized form of Eq. (1) may be written in components as Solving the Laplace equation ∇ 2 φ = 0 for perturbations of the electric potential in the udoped region we find the solutioñ φ ∝ exp(σt + iky − |k|x).
Eqs. (3), (4), (6) may be reduced to a single linear equation describing the perturbation growth with time Thus, we obtain an instability for the p-type doping front, where the initially small perturbations grow exponentially asX p ∝ exp(σt), with the growth rate The same result holds for the n-type doping when U p is replaced with U n . The obtained dispersion relation σ ∝ |k| is mathematically similar to the Darrieus-Landau (DL) instability encountered in combustion, astrophysics and laser fusion [29][30][31][32] . The role of the electric field in the new instability may be also compared to that of the gas velocity field in flames. According to the above analysis, the doping front instability develops because of a positive feedback of two effects: First, the local doping speed depends on the electric field, Eq. (1); and, second, the electric field increases close to the humps at a curved surface of the conducting material 33 . The latter is reminiscent of the St. Elmo's fire phenomenon, which visualizes strong increase of electric field at sharp conducting surfaces by numerous discharges in the surrounding air. In the present configuration of a curved doping front, any small hump at the front produces local increase of the electric field in the undoped material close to the hump, which, in turn, facilitates further growth of the hump, thus leading to the instability.

III. THE NONLINEAR EQUATION FOR THE P-DOPING FRONT BENT BY THE INSTABILITY
A smooth shape of the p-doping front observed in the experiments 19 allows reducing the full diffusion-mobility model 25 to a single nonlinear equation. Similar to the previous section, we consider a single p-doping front propagating in the initially uniform external electric field E =â x E 0 . In developing the nonlinear theory of the doping front instability for perturbations of considerable amplitudes, we have to take into account interaction of modes with different wave numbers. Then Eq. (7) has to be presented in the form where the DL operatorĴ implies product by |k| in Fourier space, i.e.ĴX in the 2D geometry of an OSC film. Equation (9) will be obtained below as the principal linear terms in the nonlinear equation for the curved p-doping front.
We study dynamics of the p-doping front as a geometrical surface using the identity which leads to the normal front velocity 30,31 Here we have introduced the unit normal vector to the surfaceâ On the other hand, the normal velocity of the curved pdoping front is determined by the electric field just ahead of the front, see Eq. (1), which leads tô Taking into account the relation between the planar front velocity U p and the external electric field E 0 , we write Eq. (15) in components as We also employ the solution to the Laplace equation for one Fourier mode of the perturbed electric fieldφ ∝ exp(iky − |k|x) and find the general identity where the Hilbert operatorĤ means multiplication by |k|/ik in Fourier space. Then Eqs. (16), (17) may be reduced to Within the limit of weak second-order nonlinearity we simplify Eq. (19) as The first two terms in Eq. (20) describe the linear instability growth similar to Eq. (9). The third term corresponds to the nonlinear Huygens stabilization of the instability due to the front propagation, which is wellknown, e.g., in combustion science [29][30][31] . Because of the Huygens stabilization convex parts of the front (humps) become smoother, while concave parts converge to sharp cusps, as illustrated in Fig. 3. Propagation speed of the cusps increases in comparison with the planar front velocity U p , and the cusps tend to catch up with the humps. When linear growth of the humps because of the instability is balanced by the nonlinear increase of the cusp velocity, we should expect formation of a curved stationary or quasi-stationary p-doping front propagating with an increased velocity in comparison with U p . Still, one should employ the similarity between combustion fronts and the doping fronts in OSCs with caution, since these two phenomena demonstrate also a number of important different features apart from similar ones. In particular, the normal velocity of a thin combustion front may be typically treated as a constant value independent of the large-scale gas flow of the fuel mixture [29][30][31] . On the contrary, the normal velocity of a doping front is determined by the electric field in the undoped material just ahead of the front, see Eq. (1). Finally, the last term in Eq. (20) has not been encountered before in combustion science or any other study of front dynamics. This term is intrinsic for the doping fronts in OSCs, and, as we will see below, it provides additional nonlinear destabilization of the front since it has the same origin as the linear destabilizing termĴX p . The nonlinear destabilization may be quite strong in the 3D geometry so that the whole approach of weak nonlinearity fails and one has to go back to the fully nonlinear equation Eq. (19), which is much more difficult for solution than Eq. (20). Still, in the 2D geometry of the OSCs films studied experimentally in Ref. 19 the approach of weak nonlinearity provides reasonable results able to explain the experimental data both qualitatively and quantitatively.
In the studies of the unstable p-doping fronts we also have to take into account linear stabilization of the instability at short wavelengths below the cut-off λ p ∝ L p . The short wavelength stabilization is related to the small front thickness and, therefore, was not accounted in Eq. (7). Though negligible at the experimentally observed length scales, the stabilization prevents the ultraviolet divergence of the perturbation spectrum in the doping front dynamics by use of the modified dispersion relation The above dispersion relation may be also presented in the form of a linear integral-differential equation similar to Eq. (9) The linear stabilization terms have the physical meaning of the local front velocity depending on the front curvature [29][30][31] ; these terms may be also incorporated as curvature corrections to the doping front velocity in Eq.
(1). Then, taking into account the linear stabilization term in the nonlinear equation Eq. (20), we come to the final equation for the p-doping front We point out thatĴf is the key term describing the front instability at the linear stage; thus the product (∆U/U p )Ĵf with ∆U > 0 means additional nonlinear destabilization of the front. It should be stressed that such a term has not been encountered before in the combustion studies of the DL instability, and it modifies the nonlinear dynamics of the p-doping fronts quite strongly as compared to the combustion fronts. By employing the method of pole decomposition 36-40 , we look for λ-periodic solutions to Eq. (24) in the form where the complex values y α represent the poles, unknown so far. Substituting the pole solution Eq. (25) to the stationary equation (24) we derive the set of equations for the poles y α . In these calculations, action of the DL operatorĴ on Eq. (25) requires additional comments.
The definition of the operatorĴ may be re-written as followŝ where the integral in Eq. (26) is similar to the Fourier transform. We integrate Eq. (26) by parts using the following identities 41 2 ∞ n=1 e −nz cos ny = sinh z cosh z − cos y − 1, z > 0. (28) As a result, we can specify action of the DL operator on the pole solution Eq. (25) aŝ with the designation ε α = sign [Im(y α )] introduced for brevity. Then, we substitute Eq. (25) to Eq. (24) and perform straightforward calculations to find the set of equations for the poles where N pairs of poles are involved in the solution. In the derivation of Eq. (30) we also employed the following trigonometric formula Apart from the equations for the poles, the calculations yield an equation for the velocity increase where N is an integer standing for the number of the pole pairs. Solution to Eq. (32) is illustrated in Fig. 4 (a) for N = 1: the solid lines present the solutions with velocity increase limited for any λ p /λ; the dashed lines show the solutions with unlimited velocity increase. The latter, presumably, does not have a physical meaning and hence will be omitted in the following. The solution presented by the solid line is determined by a simple formula which follows from Eq. (32). We can also obtain the same result for the velocity increase by taking the spatial averaging of Eq. (24) as The velocity increase is presented in Fig. 4 (b) versus the scaled period of the solution λ/λ p for different numbers of the pole pairs. As we can see from Fig. 4 (b), a curved p-doping front is possible only for a sufficiently large period λ exceeding the instability cut-off λ ≥ λ p ; the branch starting at λ = λ p corresponds to one pair of poles, N = 1. The solution with N pairs of poles becomes possible for λ > N λ p and it provides the maximal front velocity within the interval of N < λ/λ p < N + 1. As we increase the period λ in comparison with λ p , multiple stationary periodic solutions of the same period become possible with different numbers of the pole pairs. In Fig.  4 (b), the solid lines indicate the solutions corresponding to the maximal velocity increase for a fixed period λ; other solutions are presented by the dashed lines. The solutions with maximal velocity involve also the largest possible number of poles for a fixed period. Presumably, the solutions with the maximal propagation velocity and the maximal number of poles are stable, while other solutions (shown by the dashed lines) are unstable. In the context of combustion fronts, such a result has been demonstrated in Ref. 38 by investigating stability of the pole solutions; we expect the same conclusion to hold for the doping fronts in OSCs as well, though the detailed stability analysis requires much work and is left for the future. Figure 4 (b) demonstrates also that the maximal velocity increase tends to the limiting value max∆U/U p = 1 (35) as the period λ becomes much larger than the cut-off λ p , i.e. for λ/λ p → ∞. The limit of large scales in comparison with the cut-off is of the most importance for the analysis of the experimental results 19 . Thus, we expect strongly curved p-doping fronts in OSC films to propagate twice as fast in comparison with the planar fronts. Next, by using Eq. (30), we investigate the front shape predicted by the stationary solutions. Translational invariance of the solution Eq. (25) implies that we can chose arbitrary Re(y α ); hence we set Re(y α ) = 0, so that the poles y α become pure imaginary y α = ±iB α . Taking into account complex conjugate pairs, we rewrite Eq. (25) in a more convenient form and the set of equations for the poles, Eq. (30), reduces to  Fig. 5 (a)). We observe in the figure that the pole pair with the minimal number α = 1 plays the dominant role in building the amplitude of the final hump, while the poles with higher numbers are responsible mostly for fine details. In that sense, the role of poles in the solution (36) may be compared qualitatively to the role of harmonics in Fourier series. Substituting the solutions for the poles into Eq. (36) we find the shapes of the curved p-doping fronts. In agreement with the experimental observations, Fig. 1, the shapes are built of smooth humps facing the undoped material and sharp cusps pointing to the doped one, e.g. see Fig. 6 (a) with the solution plotted for λ/λ p = 2. For illustrative purposes, Fig. 6 (a) presents three periods of the solution. We remind that in this section we focus on the stationary solutions with humps and cusps reproduced periodically. Still, these humps provide the key to understanding irregular front patterns as well like those observed experimentally and shown in Fig. 1. Figure 6 (b) presents the stationary humps corresponding to the maximal front velocity and the maximal number of poles found for different values of the scaled period λ/λ p = N = 1 − 10. In Fig. 6 (b) we observe only minor variations of the hump shape with increasing period. The variations concern mostly the cusps, which become sharper as we increase the hump size λ in comparison with the instability cut-off λ p . We remind that the limit of large length scales, λ/λ p → ∞, is the most interesting from the experimental point of view. As we can see in Fig. 6 (b), the shape of a hump becomes self-similar in the limit of large length scales, and this shape may be approximated with a very good accuracy by a solution with a limited number of pole pairs. In particular, Fig. 6 (b) demonstrates negligible difference between the solutions obtained for λ/λ p = N = 5 and 10.
At the end of this section we point out that the pole solutions may be written in an explicit analytical form for small numbers of pole pairs N = 1, 2. For one pair of poles (N = 1) we find from Eq. (37) which leads to the front shape In the same way, for two pairs of poles (N = 2) we obtain where the designation Λ = [3λ p /(λ + 2λ p )] 2 has been introduced for brevity.

V. COMPARISON TO THE EXPERIMENTS AND TECHNICAL IMPLICATIONS
In this Section we compare the theoretical predictions to the experimental results of Ref. 19 and demonstrate how theoretical understanding may be employed for better control and optimization of the electrochemical doping process in optoelectronic devices. In particular, Figure 7 compares positions of the doping fronts as predicted by the theory and measured experimentally. Here, the solid lines present front positions assuming the planar front shape (labeled "theory 1") and accounting for the instability and the curved shape of the p-doping front (labeled "theory 2"). The markers show average positions of the strongly corrugated fronts obtained experimentally in Ref. 19 and illustrated in Fig. 1; the error-bars indicate average spreading of the front brush. Details of the theoretical calculations for the geometry of the planar fronts may be found in Ref. 25 . Here we remind two important features of the process taken into account in the calculations: (i) according to the experiments 19 , the n-doping front starts later and propagates noticeably slower than the p-doping front; (ii) the fronts accelerate towards each other because of much higher conductivity of the doped material in comparison with the undoped one. Due to the relatively high conductivity of the doped material, the potential difference applied to the electrodes drops effectively at the gap between the fronts. As the distance between the p-and n-fronts decreases, the electric field in the gap increases, which leads to increase of the front velocities according to Eq. (1). Theoretically, the electric field and the planar front velocities diverge as the fronts meet; the front propagation comes to a complete halt at that point. The theory of Sec. IV predicts that, in the experimental conditions, the doping front instability and the corrugated shape may increase the propagation velocity of the p-doping front by a factor of 2. We employed this theoretical result for calculating the average position of the corrugated p-doping front labeled as "theory 2" in Fig. 7. Since the theory is valid only for p-doping fronts, then we use the model of a planar n-doping front as in the previous calculations. Still, the n-doping front sweeps only about 10% of the active material in the described experiments, and, therefore, a particular model for the nfront dynamics produces only slight modifications of the results.  Figure 7 shows that the well-developed front instability and the curved shape of the p-doping front may speedup the doping process considerably thus improving performance of the optoelectronic devices. The theoretical predictions for the p-front dynamics are in a good agreement with the experimental data during most part of the doping process (for the time instants up to ∼ 65 s). The front velocity found experimentally is only a little lower than the theoretical expectations for the maximal front velocity within that interval. Some deviations in the front position and velocity between the theory and the experiments observed in Fig. 7 are quite natural, since the periodic solutions studied theoretically in Sec. IV are not identical to the irregular pattern of Fig. 1 obtained experimentally. As we discuss below, the difference in the front shape should lead to larger velocity of the periodic curved fronts of Sec. IV in comparison with the front pattern of Fig. 1. The deviations between the theory and the experiments become considerable at the end of the process, for the time instants about (70−80) s. On the other hand, the very approach of the average front position is misleading at the end of the doping process, since the gap between the p-and n-fronts becomes comparable to the size of the humps. Figure 8 compares shape of the humps at the curved p-doping front as predicted by the theory for a stationary periodic solution and obtained in the experiments 19 . In order to perform the comparison, we shift and re-scale the humps labeled in Fig. 8 (a) keeping equal scales in both directions. We point out that comparison is possible only for the top part of the humps, since the bottom parts are strongly affected by particular irregular conditions of front initiation. Figure 8 demonstrates a reasonable agreement between the theory and the experiments. Still, we can see that the periodic stationary humps studied in Sec. IV are somewhat smoother and the cusps in between them are more shallow than those observed in the experiments. In part, the difference may be attributed to the limited accuracy of the assumption of weak nonlinearity used for the derivation of Eq. (23) from the fully nonlinear equation Eq. (19). However, the fully nonlinear equation Eq. (19) is much more difficult for the analysis than Eq. (23); for this reason, investigation of Eq. (19) is left for future studies. Large amplitude of the humps in Fig. 1 may also originate in a particular way of initiating the instability in the experiments 19 . As demonstrated in Fig. 1, the strong instability of the p-doping front has been initiated by introducing minor corrugations at the anode surface. Because of the local increase of the electric field at the corrugations, the doping process started faster at the corrugation tips. For a while, the regions of p-doping developed separately in the form of isolated humps close to the anode corrugations without forming a single continuous p-front. When the front was finally formed, the employed way of initiating the instability produced extremely deep cusps -chasms -in between the humps, which remained quite deep till the end of the doping process thus reducing the average front position and the average front speed. Here we would like to stress an important difference between the doping fronts in OSCs and combustion fronts one more time; the similarity between these two phenomena is limited and may be often misleading. In the case of combustion, normal velocity of the front propagation is approximately constant, so that deep chasms at a combustion front imply strong increase of the front speed, and the Huygens stabilization mechanism leads to fast leveling of the chasms as shown in Fig. 3. In contrast to combustion, such effects do not happen at the doping fronts in OSCs. The doping front velocity is determined by the electric field just ahead of the front, which is strongly reduced in the chasms. In particular, Fig. 9 presents the electric field (a) and potential (b) calculated for the corrugated doping front of Fig. 8 (a). The electric field in Fig. 9 varies by 3 − 4 orders of magnitude from ∼ 12kV/m at the tips of the humps to ∼ (10 − 100)V/m at the bottom of the chasms.
For this reason, the Huygens stabilization mechanism is much weaker for the doping fronts, and deep chasms reduce the average front velocity instead of increasing it.
As a result, moderately corrugated doping fronts studied in Sec. IV provide higher average propagation velocity than the strongly corrugated fronts of Fig. 1. The electric field (a) and potential (b) in the undoped material ahead of a moderately corrugated doping front studied in Sec. IV is shown in Fig. 10; in that case the electric field between the humps and the cusps varies much less than in Fig. 9, approximately by an order of magnitude, from ∼ 5kV/m to ∼ 0.5kV/m. Thus, in order to obtain further increase of the doping rate in comparison with the experiments 19 , one should avoid extremely deep chasms when initiating the doping front instability. This can be achieved, for example, by introducing smoother corrugations at the electrodes imitating the curved periodic doping fronts obtained in Sec. IV.
Another important feature of the process is the selfsimilar shape of the humps at the corrugated p-doping front, which allows additional control of the front corrugations and, hence, of the final light-emitting p-n junction curve. As we have shown in Sec. IV, in the limit of large length scales λ/λ p ≫ 1 inherent to the experiments 19 , the maximal velocity increase tends to the limiting value ∆U/U p → 1 independent of a particular length scale λ, and the humps attain a selfsimilar shape presented in Fig. 6. We point out that the sample sizes in the experiments 19 , namely, breadth of 1 mm and length of 1 cm, were many orders of magnitude larger than characteristic width of the doping front, L p ∼ (10 −4 −10 −3 )mm, determined by ionic diffusion 23,25 . For this reason, humps of both large and small size in comparison with the sample sizes may satisfy the limit λ ≫ λ p ∝ L p and provide maximal velocity of front propagation. As an illustration, Fig. 11 shows two p-doping fronts with two different periods of the humps, λ = 0.3 mm and 1.2 mm, which propagate with the same maximal velocity. Visually, front (b) demonstrates strong corrugations, while front (a) is almost planar in average when compared to the sample sizes. Respectively, in the case (b) one should expect a strongly curved lightemitting p-n junction, while in the case (a) the p-n junction will be almost planar with slight wrinkles, though both cases demonstrate the same maximal doping front velocity increased because of the instability.
Finally, we discuss implications of the present analysis for the 3D geometry of the doping fronts. Though the experiments of Ref. 19 concerned only 2D OSC films, the 3D bulk geometry of the doping process is also technically attainable for transparent electrodes made, e.g., of chemically derived graphene 18 . In the 3D geometry, the stationary weakly nonlinear equation (24) becomes where ∇-operator concerns two transverse coordinates, y and z. By using the methods of combustion theory 31,40,42 , we may look for a particular solution to Eq. (42) in the form The solution (43) describes only rectangular 3D patterns at the doping fronts as illustrated in Fig. 12. Still, comparison of triangular, rectangular, hexagonal and axisymmetric patterns for a similar combustion problem demonstrated that increase of the front velocity depends slightly on a particular pattern 43 . The combustion studies suggest also the evaluation max(∆U 3D ) ≈ 2 max(∆U 2D ), which would imply velocity increase by a factor of 3 for curved p-doping fronts in the 3D bulk OSC devices. However, careful investigation of Eq. (42) indicates that velocity increase of a curved p-doping front with 3D patterns is even stronger than that suggested by the above estimate. Indeed, by substituting the solution Eq. (43) into the stationary equation (42), we come to the following equation for the maximal velocity increase in the 3D geometry which looks similar to Eq. (32), but, surprisingly, has no physical solution. The last fact implies some important consequences for curved doping fronts in optoelectronic devices with the 3D bulk geometry: (i) the approach of weak nonlinearity of Eq. (23) does not work in the 3D geometry, and one has to employ the fully nonlinear equation (19); (ii) presumably, the doping front instability is so strong in the bulk geometry, that it cannot be not stabilized at the nonlinear stage to any stationary/quasi-stationary solution. Instead of the stabilization, one may expect acceleration of the 3D humps at the doping fronts resembling the breakdown in lightning and development of streamers 44,45 . Such an acceleration of the doping fronts with 3D patterns may be also compared to acceleration of curved combustion fronts in deflagration-to-detonation transition 46,47 . Similar to the present work, combustion studies have encountered an example of curved flames in channels with the propagation regime determined by the geometry dimension: in open tubes, the combustion fronts propagate stationary in average in the 2D geometry, but accelerate in the 3D case, see 48,49 . Thus, one should expect new interesting effects for doping front dynamics in the 3D bulk optoelectronic devices, which are qualitatively different from the case of 2D OSC films, and which may further improve the device kinetics.

VI. SUMMARY
In the present work we have developed a nonlinear theory of the doping front instability discovered recently in Ref. 19 . Our nonlinear studies predict the possibility of quasi-stationary curved p-doping fronts in OSC films, such as the ones employed, e.g., in light-emitting cells. The curved doping fronts may propagate with a velocity that is increased by a factor of 2 in comparison with planar fronts, thus allowing for an optimization of the optoelectronic devices. The theoretical predictions for the velocity increase are in good agreement with the experimental measurements, see e.g. Ref. 19 . According to the theory, further increase of the front velocity may be achieved by fine tuning the initiation of the instability using minor smooth humps introduced at the electrodes.
The theory also predicts the front shape in the form of smooth humps facing the undoped material and sharp cusps facing the doped substance, in agreement with experimental observations. At the same time, in contrast to analogous studies of combustion front dynamics, extra deep cusps at the doping fronts do not produce additional increase of the front velocity; instead, they lead to certain velocity reduction in comparison with the moderately corrugated doping fronts. In the experimentally important parameter domain, humps at the p-doping fronts demonstrate self-similar properties, which allows for better control of the front shape and the doping rate. In particular, we suggest that a maximal doping rate may be achieved both for strongly corrugated fronts, and for fronts which look planar on average on large scales but demonstrate corrugations on small scales. This interesting property of the doping front instability allows control of the shape of the light-emitting p-n junction.
Finally, the present analysis suggests that dynamics of the curved p-doping fronts in optoelectronic devices with a 3D bulk geometry may be qualitatively different from that in the 2D OSC films studied in this paper. We expect that p-doping fronts with 3D patterns will not be stabilized to a stationary (or quasi-stationary) nonlinear pattern even in the case of a single front propagating in an initially uniform external electric field, but they will rather accelerate, similar to breakdown in lightning. Such an acceleration may reduce the turn-on time in the 3D optoelectronic devices considerably in comparison with the light-emitting cells based on 2D OSC films.