Nonequilibrium orientational patterns in two-component Langmuir monolayers

A model of a phase-separating two-component Langmuir monolayer in the presence of a photo-induced reaction interconvering two components is formulated. An interplay between phase separation, orientational ordering and treaction is found to lead to a variety of nonequilibrium self-organized patterns, both stationary and traveling. Examples of the patterns, observed in numerical simulations, include flowing droplets, traveling stripes, wave sources and vortex defects.


I. INTRODUCTION
Spurred by experimental and technological developments in supramolecular chemistry and biophysics, there is a raising interest to study nonequilibrium structures related to self-organization phenomena in weakly condensed systems. In such systems, attractive potential interactions between constituent molecules are essential. Equilibrium structures in soft matter correspond to minima of free energy, and are caused by the competition between short-range attractive and long-range repulsive interactions [1]. In contrast to this, nonequilibrium patterns require permanent supply of energy and correspond to persistent (stationary or time-dependent) kinetic states of a system [2]. Typically, they result from an interplay between reactions, diffusion and phase transitions. Examples of nonequilibrium pattern formation in reactive soft matter include stationary Turing-like patterns in phase-separating binary polymer mixtures with chemical reactions [3][4][5] and in monomolecular adsorbates on metal surfaces [6,7]. Theoretical analysis for twocomponent reactive adsorbates [8,9] and for three-component reactive polymer systems [10,11] has further shown that not only stationary patterns, but also traveling and standing waves are possible.
Langmuir films are monomolecular layers of amphiphilic molecules on an air-water interface. Such films are a classical example of soft matter and, generally, it should be expected that, in the presence of chemical reactions and energy flows, they would form nonequilibrium patterns. Though equilibrium properties of Langmuir monolayers are thoroughly investigated (see [12]), nonequilibrium phenomena in these systems still remain less explored. In an interesting series of experiments, Tabe and Yokoyama have studied illuminated Langmuir monolayers of amphiphilic derivatives of azobenzene by means of the Brewster-angle microscopy [13][14][15][16]. In these experiments, transitions between trans and cis conformations of individual molecules were photo-induced by polarized light of a selected wavelength. Since physical properties of the two conformations are different, the trans and cis isomers essentially represented two different species. In addition to photo-induced periodic stationary patterns, these experiments have shown, for the first time, spontaneously emerging patterns of propagating waves of molecular reorientation under appropriate illumination conditions.
In the Letter [17], we proposed a model of reactive two-component Langmuir monolayers with orientational ordering. This model was taking into account phase separation in the two-component system, reaction interconverting both species, diffusion of reactants and processes of orientational ordering. Interactions between the components of the monolayer resulted not only from the positional, but also from the orientational order of the hydrophobic tails of constituent molecules, determined by their tilt. We have shown that this model already reproduces nonequilibrium traveling-wave patterns which arise as a consequence of a Hopf bifurcation with a finite wavenumber. To simplify the analysis, it was assumed that the azimuthal orientation of molecules remained fixed and uniform, so that only their tilts could vary. Moreover, some orientational-order contributions to the free energy of the monolayer were neglected. The inclusion of azimuthal variations is, however, important for a comparison with the experimental data yielded by the Brewster-angle microscopy that is sensitive to the local azimuthal orientational ordering.
The aim of the present paper is to formulate and to study a more general model that contains both orientational variables, and includes bend and splay distortion terms. Our analytical and numerical investigations show that this model has a significantly different phase diagram and new kinds of nonequilibrium patterns are possible here. After introduction of the model in Sec. II, we investigate in Sec. III the behaviour of the system in the equilibrium case, when the reaction is absent. The bifurcation analysis of the uniform steady state of the nonequilibrium system under illumination is presented in Sec. IV. Numerical simulations, revealing the formation of such nonequilibrium patterns as traveling droplets and stripes, stationary splay-defects, and complex azimuth reorientation kinetics, are reported in Sec. V. The paper ends with the conclusion and the discussion of the obtained results.

II. THE MODEL
We consider a monolayer formed by two diffusive immiscible components, A and B, that have strongly different shapes. Modeling the situation in the photo-isomerization experiments, molecules A are supposed to have an elongated shape with a long tail (the trans isomer) and molecules B to have a crumpled conformation (the cis isomer). According to this assumption, only molecules A are subject to orientational order, whereas molecules B play essentially the role of passive dilution with respect to such ordering. Furthermore, a photo-induced reaction interconverting A and B molecules is considered. The total concentration of components A and B in the monolayer is assumed constant. Therefore, the local composition of the monolayer is characterized by the concentration order parameter c representing the local fraction of molecules A (so that 1 − c gives the local fraction of molecules B). The local orientational order is described by the order parameter a, that corresponds to the projection of the local mesoscopic average of the unit molecular director n of the elongated molecules A onto the monolayer plane (see Fig. 1). The vector a is defined by its modulus a = | a| = sin η, where η is the tilt angle, and its azimuth angle ϕ. After we have introduced the two order parameters c and a needed to describe the system, the kinetic equations for their evolution should be formulated. Following the mesoscopic approach [17], the evolution equation for c reads as Here D is the diffusion coefficient, T is temperature, and µ is the chemical potential defined as µ = δF /δc, where F is the free energy functional that will be specified below. In the presence of nonpolarized light, the reactive term R(c) is given by [18][19][20], where the first term corresponds to the conversion of molecules B into molecules A and the second term describes the reverse conversion process. In this equation, I is the light intensity, λ is the wavelength of light, and f (λ) and g(λ) are related to the surface molar extinction coefficients and quantum yields for the conformations of A and B, respectively. For moderate and high illumination intensities, photo-induced conversion dominates over the thermal conversion processes, so that the thermal rate constants k 0 ±1 can be neglected. In this case, the ratio of conversion rates for forward and backward processes is independent of temperature or light intensity, and determined only by the light wavelength λ. For simplicity, we shall assume in our subsequent analysis that the wavelength λ is chosen in such a way that f (λ) = g(λ) . Under this condition, equation (2) takes the form where k = If (λ) = Ig(λ) is the reaction constant proportional to the intensity of the applied light. The kinetic equation for the local orientation vector a is chosen as The first term on the right side corresponds to relaxation with a typical relaxation time Γ −1 . The second term takes into account that the reaction, interconverting one molecular conformation to another, also affects the local average orientation of molecules. The choice of the function q(c) should be based on the kinetic analysis with respect to the orientation of molecules involved in the reactive process. In our model, the reaction B −→ A, transforming crumpled molecules into the elongated form, is assumed to be strongly energetically activated by light. This means that new molecules A would be created with high initial energy and can adopt any orientation a, independent of the orientation of surrounding molecules. Therefore, the orientation of newly created molecules A is random and the statistical average of a is zero. Because any conversion from an elongated molecule A to the crumpled molecule B decreases the local order and the reverse conversion process does not increase it, the overall reaction effect is that it tends to destroy orientational ordering. Under such assumptions, local evolution for the orientation momentum c a is described by a simple equation ∂ t (c a) = −kc a, that corresponds to the loss of c a when molecules A transform into B, having no defined orientation. Splitting ∂ t (c a) and substituting the local variation of c due to reaction, ∂ t c = k(1 − 2c), leads to ∂ t a = −k a(1 − c)/c. Therefore, we have q(c) = (1 − c)/c. This function q(c) is used below.
The energy functional F is constructed in terms of the order parameters c and a of the system. We decompose the energy functional into two parts: one that accounts for the affinity between the two isomers in the monolayer, and the other due to the tail-tail orientational interactions. The first of these contributions, F c , depends exclusively on c, whereas the tail-tail interaction, F or , is sensitive to both c and the orientation of the elongated molecules. Accordingly, F = F c (c) + F or (c, a). Note that entropic contributions are not considered, since they have been directly included as the diffusive term in the kinetic equation (1). The expression for F c reads as [17], and describes phase separation due to lateral interactions between molecules. In the limit of short-range lateral interaction potentials,χ 0 is determined by its strength, andχ 2 can be estimated asχ 2 ≈ 1 2χ 0 r 2 0 , where r 0 is the characteristic radius of the interaction [21]. Near the critical point of the equilibrium phase separation, the combination of Eqs. (1) and (5) is equivalent to the usual Cahn-Hilliard equation with the Landau free energy functional [7].
The part of the free energy functional associated with the distortion of the orientation of tilted molecules in the monolayer can be written as The first two terms correspond to the Landau expansion up to the quartic term for the modulus of the orientation vector [12]. This expansion is only justified for sufficiently small η, since in this case a = sin η ≈ η is small as well. Generally, all coefficients in the Landau decomposition should depend on the local concentration c. We shall, however, consider only weakly nonuniform states, where local deviations of the concentration c from the uniform stationary state c = 1/2 are small. Therefore, we neglect the dependence of the coefficientβ on the variable c. However, the coefficientp of the quadratic term in the Landau free energy is already small near the instability and its dependence on c must be retained. For condensed phases, lowering the lateral pressure of a Langmuir monolayer leads to an increase of its equilibrium tilt [12]. Since we assumed that molecules B play a role of passive dilution for the tilted molecules A, decreasing c is roughly equivalent to decreasing the lateral pressure. Thus, we choose a linear dependencẽ p(c) =π 0 + 2α(0.5 − c), whereα is a positive coefficient, andπ 0 is a decreasing function of the lateral pressure which determines the equilibrium tilt in the monolayer with c = c. The other contributions in Eq. (6) stand for the bend and splay distortion terms [22]. More specifically, the third and fourth terms correspond to the classical Frank elasticity terms, that account for the splay and bend distortions, respectively. Normally, one takes the single Frank constant approximation (K =K s =K b ). The fifth term is the lowest-order splay contribution, ( ∇ · a), that appears coupled to the composition order parameter c. Although some authors [23,24] prefer to couple the linear splay term to the tilt angle, we follow the suggestion by Selinger et al. [25] for a two-component monolayer undergoing phase separation, as it is in our case. A similar approach was taken by Tabe et al. [14,15] who coupled the linear splay term with a certain order parameter related to the molecular density that varies across the monolayer. In general, for sufficiently strong couplingΛ, the linear splay term destabilizes uniformly oriented phases, leading to equilibrium nonhomogeneous splayed states (see Sec. III). Finally, notice that a term ( ∇ × a) linear with respect to bend distortion is not considered, because it is not permitted by symmetry in nonchiral Langmuir monolayers. Summarizing, the model presented here can be viewed as a Cahn-Hilliard equation for the composition variable c, coupled to a relaxational equation for the orientation order parameter a. The contributions to the free energy have been derived considering that c is close to its stationary uniform solution and for sufficiently small tilt angles η.
The analysis can be simplified by appropriately adimensionalizing energy, time and space. Energy is measured in units of k B T , time in units of the relaxational time (Γk B T ) −1 , and spatial coordinates are rescaled with the relaxational length D/(Γk B T ). The model is then characterized by the dimensionless parameters χ 0 =χ 0 /k B T , With this choice, the final equations for the evolution of c and the two components of vector a read as where . This is the mathematical model which will be investigated below.
In order to obtain numerical results from the model, we numerically integrate Eqs. (7) on a 100 × 100 square grid, using an explicit Euler scheme with a mesh size ∆x and a time step ∆t that assure a good numerical accuracy. Periodic boundary conditions are chosen to model the behaviour in a large system far from the boundaries. As initial conditions, small random perturbations around the stationary uniform states of the system c and a, and a random distribution of azimuth angles ϕ ∈ (0, 2π) are taken. To display simulation results, snapshots of the patterns after the transients are given in the figures. Each figure consists of two panels: the left panel shows in grey scale the spatial distribution of the composition variable c, with larger values corresponding to the darker color, and the right panel is used for visualization of the orientational field a. The local directions of this field are visualized by using small arrows. Note that, for technical reasons, such arrows could not be used to indicate the states of all grid points in the simulations and therefore the visualization of the azimuthal orientation is rough. The grey color in the right panels is used to display the local tilt a, and again, darker regions correspond to the larger values of this variable. In some cases, videos of time-dependent patterns are additionally provided.

III. EQUILIBRIUM PATTERNS
Before addressing the nonequilibrium cases, we show some examples of equilibrium pattern formation in the considered system. The equilibrium conditions correspond to absence of illumination and are realized if all reactive terms in Eqs. (7) are omitted. One of the limitations of our previous simpler model [17] was that it did not exhibit any equilibrium pattern formation, despite the experimental evidence of spontaneous generation of striped patterns in non-illuminated monolayers [13][14][15]. The linear splay term included in Eq. (6) can already lead to the formation of equilibrium orientational structures.
In the limit of α → 0 (i.e., in the absence of the quadratic Landau term for the tilt variation), the nonreactive version of the present model is similar to the description proposed by Selinger et al. [25] for nonchiral Langmuir monolayers. In their study, three different nonuniform phases (sinusoidal stripes, soliton stripes and square lattice of vortices) were found for smectic films undergoing chiral symmetry breaking under variation of a control parameter (corresponding to temperature), and similar phases for multicomponent Langmuir monolayers were predicted. As illustrated in Figs. 2, 3, 4 and 5, our model reproduces these regimes under appropriate choices of the parameter χ 0 specifying the relative intensity of energetic lateral interactions.
When the characteristic energy of lateral interaction is much weaker than the orientation energy, a square array of alternating "inward" and "outward" splay defects (vortices) is formed (Fig. 2). The "inward" defects represent regions rich in elongated molecules A which are oriented towards the center of the defect. In analogy to the "escape to the third dimension" found in defects in 3D nematics [26], the tilt in the center of these defects is almost zero, reducing the Frank elastic energy near that point. The "outward" splay defects are poor in elongated molecules A, but also exhibit vanishing tilt in their centers. Increasing the lateral interactions with respect to the splay coupling, the system organizes into a stripe phase with smooth variations of c and a (see Fig. 3). The profiles of variation of c, a and ϕ across a stipe are displayed in Fig. 4. The equilibrium stripe patterns are similar to those observed in the experiments in absence of illumination [14,15,27]. The modulation of the tilt (absent in the model by Selinger et al. because tilt variations were not included there) has a spatial frequency twice that of the azimuthal and concentration modulations. Moreover, in agreement with the experiments the amplitude of the tilt angle modulation depends on the stripe size, decreasing as the stripe widens (this is observed, for example, by decreasing the parameter Λ).
Strong lateral attractive interactions (i.e., large interaction strengths χ 0 ) lead to the formation of a striped phase with sharp wall domains ("soliton stripes") in the modulation of composition, tilt and azimuth, which has also been reported by Selinger et al. [25]. In Fig. 5 profiles of variation of c, a and ϕ across such a stripe are plotted. Note again the double frequency modulation of the tilt with respect to the composition and azimuth variation.   We have checked that the inclusion of the quadratic Landau term for the tilt variable, neglected in the simulations described in this section, does not lead to significative differences in the properties of equilibrium patterns.  (7) are c = c = 1/2 and a = a = (π 0 − κ) /β, provided that π 0 > κ. The azimuthal orientation is arbitrary in such a state. For convenience, we choose a x = a and a y = 0. The linear stability of these uniform solutions is performed by adding small plane-wave perturbations δc, δa x and δa y proportional to exp(iq x x + iq y y + γ( q)t), and linearizing Eqs. (7). This leads to the linearization matrix where q = (q x , q y ) is the wavevector of the considered mode. The first line in the matrix corresponds to the composition variable δc and the next two lines stand for the orientation components δa x and δa y .
The (complex) rates γ( q) of various modes are determined by the eigenvalues of the linearization matrix L. The unstable modes are identified like those with Re(γ( q)) > 0. If the imaginary part of γ( q) is not zero for the first unstable mode, we have a wave instability (a Hopf bifurcation with a finite wavenumber), resulting in traveling or standing waves (cf. [28]). On the other hand, if Im(γ( q)) = 0 for the first unstable mode, a Turing instability leading to nonequilibrium stationary periodic patterns is realized. The values of γ( q) are yielded as the roots of the characteristic equation associated with the matrix L. This characteristic equation is, however, cubic and therefore its analytical solution is possible only in some special cases.
If splay coupling is absent (Λ = 0), the stability analysis is simplified and becomes equivalent to that of the previously studied reduced model [17]. The phase diagram for the non-splay case in the plane (π 0 ,κ) is presented in Fig. 6a. We see five different regions, whose boundaries and marginal wave numbers can be obtained analytically [17]. Region I corresponds to the wave instability regime, where traveling or standing waves are expected. In region II we have the Turing instability region, where stationary droplet-like structures with periodic variation of both local concentration and tilt are found. Region III corresponds to the uniform tilted phase. When π 0 < κ, only nontilted phases (a = 0) exist: a nonuniform phase due to a Turing instability in region IV a and stable uniform phase in region IV b. In order to understand the effects of the linear splay term, that was neglected above, we perform the stability analysis of the full model equations for a fixed wave plane q y = 0. In this case, the variable a y is decoupled from the other two degrees of freedom (actually, its dynamics is exclusively governed by the elastic damping term, see Eq. (8)), so that the stability analysis is reduced to a 2×2 problem. The following results, summarized in Fig. 6b, are obtained: For nontilted phases (π 0 < κ), the effect of the linear splay term is to move to larger reaction rates κ the boundary between the Turing-instability region (IV a) and the uniform nontilted phase (IV b). For tilted phases (π 0 > κ), the changes are more profound. All unstable modes in this area of the parameter space have now a nonvanishing imaginary part Im(γ( q)). Therefore, the region II with Turing instability and stationary tilted patterns completely disappears and becomes replaced by the region I with traveling waves. Additionally, the region III occupied by the tilted uniform phase is reduced.
Based on this (limited) stability analysis, one can expect that the traveling patterns would be found more often when the splay interactions are taken into account. Furthermore, it can be expected that the traveling patterns in the parameter region occupied by stationary Turing-like structures in absence of splay interactions (region II in Fig.  6a), would be different for weak splay interactions from the traveling patterns in the parameter region where traveling waves are observed even without the splay effects (region I in Fig. 6a). As Λ → 0, traveling waves in the former region II should slow down and become frozen at Λ = 0. Indeed, the analytical stability investigation shows that, in this region, the velocity of the most unstable mode is proportional to κΛa. On the other hand, the velocity of patterns in the former region I remains finite in the limit Λ → 0.
Numerical simulations of the model, which represent the main part of the reported study, agree with the predictions of the stability anaylsis.

V. NONEQUILIBRIUM PHENOMENA: NUMERICAL RESULTS
To facilitate the comparison with previous simulations of the model without azimuthal variation [17], we choose here the same numerical values of the common parameters χ 2 = 0.0052, χ 0 = 2.27, α = 1.5 and β = 2. The above mentioned changes in the phase diagram due to the inclusion of the linear splay term indicate that this parameter region is worth examining with the model presented here in order to obtain spatio-temporal behaviors involving now composition, tilt and azimuth modulations. We examine numerically such parameter region, and the results are summarized in the three following subsections according to the value fixed for the splay coupling constant Λ. All simulations are perfomed for a system with a linear size L = 5.

A. No splay coupling
When splay coupling is absent (Λ = 0), azimuthal orientation of elongated molecules is not influenced by variations of the local composition c or the local tilt a. Indeed, by looking at the second of equations (7) we notice that then all terms on the right hand side of this equation, except for the elastic term K∇ 2 a, are proportional to a and therefore cannot change the direction of this vector. If the azimuthal orientation is initially uniform (ϕ = const), this state is maintained at all times. For such a state, Eqs. (7) reduce to the model which was already investigated in Ref. [17] and new numerical simulations are not needed. If the initial azimuthal orientation is not uniform (ϕ = const), the subsequent evolution of the orientation field is determined only by the elastic term.
In numerical simulations, we choose the parameters inside the region I with traveling waves in the phase diagram shown in Fig. 6a, by fixing π 0 = 1.5, κ = 1 and taking a small value K = 0.001 of the elastic interaction constant. Random distribution of azimuthal orientations is chosen as the initial condition. The simulation results are presented in Fig. 7 (see also video Fig7.mpg). Because of elastic interactions, the molecules tend to have parallel orientation, and this leads, after some time, to a pattern characterized by a number of orientational defects that remains stationary henceforth. In the center of a defect, the molecules are non-tilted; the azimuthal direction changes by 2π after passing around a defect. The waves, similar to those described in Ref. [17], travel on the background of the stationary orientational pattern, and their motion is not generally influenced by the azimuth angle distribution. The waves, however, break when they pass through a defect. Inside the region II in the phase diagram in Fig. 6a, stationary Turing-like patterns representing arrays of droplets are observed.

B. Weak splay coupling
The most important change in the phase diagram due to the inclusion of the linear splay term is that the region II occupied by the Turing-like patterns in absence of splay coupling (Fig. 6a) is transformed into region I where traveling patterns should be observed (Fig. 6b). However, as we shall see, the properties of traveling patterns in the parameter domain are different, for weak coupling, from those of the traveling waves in the old region I. Figure 8 and video Fig8.mpg show the traveling pattern yielded by numerical simulation with the parameters κ = 0.5, π 0 = 1.25, and Λ = K = 0.01 that correspond to the former region II. The pattern looks like a flow of droplets, with the local direction of flow determined by the azimuthal orientation of the elongated molecules in the monolayer. Such droplet-like structures of large c and small a emerge and move following the local molecular orientation path. Occasionally, rupture of the droplets, when they happen to approach divergent points of the molecular orientation field, is observed. Some of the droplets are pinned by the orientational defects and exhibit only rotation, but not a translational motion. Remarkably, the pattern of azimuthal orientation remains frozen under weak splay coupling, as it was the case in its absence.
In agreement with the linear stability analysis, the velocity of the emerging droplet structures is an increasing function of the splay coupling constant Λ and the reaction rate constant κ. Moreover, we have found that the size of such traveling structures is affected by the strength of elastic interactions: the larger K the bigger are the droplets. For large K, however, the uniform state becomes stable, and no spatial organization is observed. Figure 9 and video Fig9.mpg show traveling waves for a larger value of the reaction rate constant (κ = 1), such that we are now inside the old region I in Fig. 6a. The morphology of the pattern is now different and it resembles the pattern of traveling waves in absence of splay coupling (Fig. 7). The pattern is formed by stripes that move along the direction determined by local azimuthal orientation of the elongated molecules. Two orientational defects are seen in the upper right corner in Fig. 9,and they also correspond to defects of the traveling stripe pattern of the composition field. The waves are rotating around these orientational defects. Again, the pattern of azimuthal orientation becomes frozen after a rapid initial transient and it is not significantly affected by the traveling waves. Note from the simulation videos that the stripes move much faster than the droplets in Fig. 8.
Next, we examine more closely the profiles of the composition c and the tilt a in different traveling patterns. Figure  10a shows such profiles for a single traveling droplet from Fig. 8, also displayed in the inset in the left upper corner of this figure. The profiles in the cross sections which are parallel and perpendicular to the motion direction are presented here. The droplet corresponds to a local increase in the concentration of elongated molecules A and a local decrease in the tilt of such molecules. It can be noticed that the tilt is also slightly increased along a ring surrounding the droplet, and that the droplet is not axially symmetric. Comparing the profiles for the perpendicular and parallel cross sections in Fig. 10a, we find that the droplet is slightly elongated in the direction parallel to the azimuth molecular orientation and, moreover, the tilt peak in the rear part of the moving droplet is higher than that at the front. This asymmetry determines the propagation direction of the droplet. Furthermore, we have checked that the asymmetry gets stronger when the splay coupling coefficient Λ is increased.
The profiles of a traveling stripe from Fig. 9 are displayed in Fig. 10b. We see that the profiles exhibit more smooth variation in this case, and the variation of c and a is closer to harmonical. Note that the spatial profile of the tilt a is retarded with respect to that of the local composition c and this again determines the motion direction. . The respective patterns and the orientations of lines used to make the cross sections are displayed in the insets. For the droplet, both the profiles along the lines parallel ("para") and perpendicular ("perp") to the motion direction are presented. The patterns move from left to right.
In the simulations described above random initial conditions were chosen. We have also studied some special initial conditions for the orientation distribution that give rise to domain sinks or sources. In the simulation shown in Fig. 11 and video Fig11.mpg, the initial condition with two orientational domains (ϕ = 0 for x < L/2 and ϕ = π for x > L/2) has been chosen. The boundary between the domains plays here a role of source emitting waves. For other initial conditions, corresponding to outward vortices of molecular orientation, we have observed droplets that are generated in the center of a vortex and travel in the radial direction out of it. Generally, linear and point defects with positive (negative) splay ( ∇ · a) act as sources (sinks) for the traveling structures that involve tilt and composition variations.

C. Strong splay coupling
When the coefficient Λ, specifying the intensity of splay coupling, is further increased, the azimuthal orientation of molecules becomes influenced by the traveling or stationary patterns. For strong splay coupling, stationary nonequilibrium patterns are usually observed. An example of such a pattern is shown in Fig. 12. Starting from random initial conditions, the system first develops a pattern of traveling stripes. Subsequently, the stripes undergo breakdown and a frozen array of orientation defects is produced. Inside each defect, the concentration of the elongated molecules is increased, and in general, these molecules are oriented towards the center of a defect. This stationary pattern is favoured in the case of strong coupling, since it minimizes the splay contribution towards the free energy of the system. To avoid the formation of immobile splay defects, a region in the parameter space, characterized by the wave instability but lying closer to the boundary of a transition to the nontilted state, can be considered. Because the elongated molecules are only slightly tilted in this region, their reorientation is energetically easier than far from the boundary a = 0. If a relatively large value of the reaction rate constant κ is additionally chosen, so that the system is driven further away from thermal equilibrium, traveling waves that are accompanied by azimuthal reorientation of molecules are formed.
An example of the developing wave pattern is shown in Fig. 13 and video Fig13.mpg. The temporal evolution of this pattern along the central vertical cross section is additionally displayed in Fig. 14. After a short transient starting from random initial conditions, the system soon developes a regular pattern of stripes traveling at a constant velocity.
The profiles of the composition c, the tilt a and the azimuth angle ϕ along the line orthogonal to the propagation direction are shown in Fig. 15. An important difference with the profiles of equilibrium stationary stripes is that, for traveling stripes, the azimuth angle does not undergo 2π rotations from one stripe to another (seen in Fig. 4). Instead, a periodic angular modulation within the interval from zero to π is observed here. Moreover, it can be noticed that the composition c and the tilt a vary completely in phase for such traveling stripes, i.e. the maximum of the concentration field coincides with a minimum of the tilt distribution. This differs from the behaviour of such variables in traveling stripes at weak splay coupling (cf. Fig. 10b). The propagation direction is now fixed by a small shift observed in the azimuth angle variation with respect to the profiles of variables c and a.  Fig. 13. Here, the grey-scale coding is such that ϕ = 0 is shown as white and ϕ = π is shown as black. The azimuth angle is defined in such a way that ϕ = 0 corresponds to the propagation direction of the stripes.

VI. CONCLUSIONS
We have formulated a theoretical model to describe pattern formation in illuminated two-component Langmuir monolayers. In this model, nonequilibrium pattern formation results from an interplay between phase separation, optically induced transitions between the two immiscible conformations, and coupling of phase separation to orientational ordering of molecules in the monolayer. In comparison with the previous model [17], the current description is more complete because variations of the azimuthal variation of the elongated molecules in the monolayer are allowed. We have also incorporated into the model the splay coupling between azimuthal orientation and variations of local concentration. The situation corresponding to illumination with non-polarized light has been considered.
In contrast to the earlier simple model [17], the investigated system is characterized, in absence of illumination, by the formation of equilibrium stationary patterns representing arrays of orientational defects (vortices), stripes or soliton-like structures with sharp domain walls. All these patterns are caused by splay interactions and have a purely energetic origin, so that they can be also interpreted by considering minima of free energy (cf. [14,15,25]).
Linear stability investigation of the uniform state of the model under illumination conditions has been performed. It reveals that, generally, splay interaction favour the appearance of traveling patterns. Such patterns are observed even in the parameter regions occupied by stationary structures in the earlier model [17].
Numerical simulations of the proposed model have been undertaken. They show a rich spectrum of spatio-temporal structures. For weak splay coupling and relatively low reaction rates, a pattern of droplets slowly flowing along the local directions determined by azimuthal orientation of molecules has been observed. If the reaction is stronger, the droplets are replaced by a pattern of traveling stripes following the local azimuthal orientation. Interactions between droplets and stripes and the orientational defects have been seen. It is found that linear and point orientational defects play a role of sources or sinks of traveling structures in this system.
When splay coupling is weak, the spatial pattern of azimuthal orientation is not affected by local concentration variations. The temporal evolution of the orientation pattern in this case is governed by the elastic interactions between molecules. Starting from a random initial distribution, they produce after a short transient a stationary pattern of azimuthal orientation that remains frozen afterwards. In contrast to this, strong splay coupling with the concentration field leads to the appearance of stationary and spatiotemporal patterns where variations of azimuthal orientation are correlated with the the changes in the local composition of the monolayer. Often, stationary arrays of splay defects are observed in this case. However, when the system is close to the orientational phase transition (from the tilted to the non-tilted phase), patterns of traveling waves accompanied by strong variation of azimuthal orientation are found in the model. This investigation, as well as the earlier Letter [17], are motivated by the experimental discovery of traveling structures in illuminated Langmuir monolayers by Tabe and Yokoyama [13,16]. The extended model, which we have now explored, is able to reproduce some essential features of the experiments. In absence of illumination, equilibrium orientation patterns are yielded by the model. Under illumination for sufficiently strong splay coupling, traveling waves accompanied by changes in azimuthal orientation of molecules are also found. The experiments [13] were performed by using polarized light, and light polarization significantly influenced the properties of the observed patterns. In contrast to this, our model does not include the effects of light polarization and thus corresponds to illumination with non-polarized light. Though the inclusion of anisotropy effects due to light polarization is straightforward, the resulting model is more complicated and not yet discussed in the present paper.
Langmuir monolayers represent a classical example of soft matter and are closely related to biomembranes, playing a fundamental role in cell biology. Therefore, investigations of nonequilibrium pattern formation in such monolayers under reactive conditions can help in understanding of general mechanisms of nonequilibrium self-organization in soft matter. Our study provides evidence that traveling orientational wave patterns may represent a generic property of Langmuir monolayers that are subjected to chemical reaction and composition changes. We have found that various patterns, representing traveling droplets or stripes, wave sources and sinks, and orientational defects interacting with the traveling structures, are possible in such systems.