Kinetic Monte Carlo Simulations for Birefringence Relaxation of Photo-Switchable Molecules on a Surface

Recent experiments have demonstrated that in a dense monolayer of photo-switchable dye Methyl-Red molecules the relaxation of an initial birefringence follows a power-law decay, typical for glass-like dynamics. The slow relaxation can efficiently be controlled and accelerated by illuminating the monolayer with circularly polarized light, which induces trans-cis isomerization cycles. To elucidate the microscopic mechanism, we develop a two-dimensional molecular model in which the trans and cis isomers are represented by straight and bent needles, respectively. As in the experimental system, the needles are allowed to rotate and to form overlaps but they cannot translate. The out-of-equilibrium rotational dynamics of the needles is generated using kinetic Monte Carlo simulations. We demonstrate that, in a regime of high density and low temperature, the power-law relaxation can be traced to the formation of spatio-temporal correlations in the rotational dy- namics, i.e., dynamic heterogeneity. We also show that the nearly isotropic cis isomers can prevent dynamic heterogeneity from forming in the monolayer and that the relaxation then becomes exponential.


I. INTRODUCTION
The possibility to control organic and inorganic materials at the molecular nanoscale level is crucial for a large variety of technological application and for a deeper understanding of matter 1-8 . Among possible tools for molecular control, light is one of the most promising. Some of the appealing applications include: illuminating the metallic tip of a scanning force microscopy to precisely control the position of single molecules 2 , using nanowires to build miniaturized photonic devices 4 , and inscribing nano-sized geometrical patterns on a surface by photolithography 9 .
Photochromic molecular switches, molecules that undergo configurational changes between two (or even more) isomeric states when irradiated by light 10 , offer yet another appealing way to control material properties with light. For example, they are used to fabricate functional surfaces with tunable chirality, wettability, conductivity etc. 11,12 . Among other things 13 , photoswitching molecules illuminated by light can reorient a nematic liquid crystal or directly control both the formation and relaxation of orientational order in a monolayer 1, [14][15][16][17] . Thus, many material properties such as mass transport, mobility, and viscosity can be efficiently tuned while producing hardly any heat. Although these phenomena have attracted much interest for many years, details of the microscopic dynamics still need to be clarified.
In this article we consider a self-assembled monolayer of light-switchable molecules tethered to a surface. Instead of the method mentioned above, we perform kinetic Monte Carlo simulations for a molecular model, where we approximate the two isomeric states, called trans and cis, by a straight and a bent needle, respectively. The simplicity of the model allows us to study the long-time collective dynamics of a statistical ensemble consisting of 10,000 molecules 28 , much more than atomistic molecular dynamics simulations can handle. In a previous work we have investigated the phase ordering of bent needles with varying shape 29 .
Our work is motivated by a recent experimental study of Fang et al. 1 on the glasslike orientational dynamics of a self-assembled monolayer of photo-switching molecules. After aligning the molecules with light, the authors observed the decay of orientational order (or birefringence) under either thermal erasure or erasure with circularly polarized (CP) light. In both cases they find that the relaxation of birefringence follows a power law, which is typical for glasslike dynamics.
Within our relatively simple model we can reproduce this feature in a system containing straight needles alone (trans molecules), if the density is sufficiently high and temperature is low. We demonstrate that the needles, when randomizing their orientations, develop dynamic heterogeneities in space and time 30,31 , which ultimately cause the power-law decay. The presence of cis molecules, which have a rather isotropic shape, can prevent the formation of such spatio-temporal variations in the local structure and the birefringence relaxation then becomes exponential. In the following, we clarify under which conditions our model nevertheless reproduces the experimental observation of a power-law decay by tuning isomerization probabilities.
The plan of this paper is as follows. In Sec. II we thoroughly review the experimental motivation for our work and explain the model and simulation method in Sec. III. The system of pure trans molecules is studied in detail in Sec. IV and in Sec. V we discuss the birefringence relaxation in a system with light-switchable molecules. We close with a summary of the results and a conclusion in Sec. VI.

II. EXPERIMENTAL MOTIVATION
The work presented in this paper is strongly inspired by recent experiments of Fang et al. 1 that have attracted considerable attention. In this section we first shortly summarize their results.
A self-assembled monolayer (SAM) with glasslike dynamics is realized by covering at high in-plane density a glass surface with dye Methyl-Red (dMR) molecules and tethering them with covalent bonds at random positions. The molecules are then free to rotate but their translational freedom is constrained to ≈ 1nm displacements. In the ground state dMR molecules assume a rod-like trans configuration with anisotropic optical properties. The light-induced cis configuration has a bent-core shape and is nearly isotropic. Because of their photo-switchable azo-core, isomerization between the two configurations can efficiently be induced by illumination with light at a 514 nm wavelength. The absorption spectra of the cis and trans isomers partially overlap and illumination with light of this wavelength induces both cis-to-trans and trans-to-cis transitions.
In Ref.
[1] the initial orientation of the molecules in the SAM was random. Under illumination with linearly polarized light (writing process), molecules preferentially aligned along the light polarization vector more likely switch toward the cis configuration, while the remaining trans molecules form nematic order perpendicular to the light polarization (hole-burning) [32][33][34] . The monolayer thus exhibits some birefringence Q(t). Since light illumination also induces the reverse, cis-to-trans transition, in steady state a mixture of both isomer exists at a relative concentration that mostly depend on the light intensity and the monolayer density.
Reference [1] investigated the relaxation of the birefringence under two different illumination conditions following the writing process: 1) The SAM is left in the dark, resulting in the relaxation being driven by thermal fluctuations (thermal erasure). At room temperature all molecules in the cis configuration relax back to the trans form after a characteristic time 35 . They assume random orientations and the nematic order is lost.
2) The SAM, immediately after the end of the writing process, is instead illuminated with circularly polarized light at 514 nm, which makes the molecules cycle between the trans and cis configurations and speeds up the relaxation dynamics of the birefringence (CP erasure). 1. The birefringence inscribed in the monolayer by illuminating it with linearly polarized light can be erased either by thermal fluctuations at room temperature (thermal erasure) or by illumination with circularly polarized light (CP erasure). In both cases the relaxation is well described by the long-time power law decay of Eq. (1). Parameters τt are shown as a diamond marker. Other parameters are η = 0.25 for the thermal erasure curve and η = 0.63 for the CP erasure curve.
A schematic of the typical experimental results is given in Fig. 1, where the temporal evolution of the birefringence Q(t) is shown for both thermal and CP erasure processes. Fang et al. 1 showed that the relaxation of the birefringence is non-exponential, and is instead accurately described by a functional form of the type (1) An asymptotic power-law relaxation starts at t ≈ τ t and proceeds as Q(t) ≈ (t/τ t ) −η for t > τ t . Non-exponential relaxations are characteristic of a glassy state, i.e., of systems sampling a rugged free energy landscape [36][37][38][39] . The underlying dynamics does not possess a single characteristic time scale, but rather a distribution of them. It has been suggested 1,16 that the distribution of energetic barriers that leads to the powerlaw relaxation in the SAM originates from i) the high packing density and ii) the molecular rotations proceeding in discrete jumps as molecules pass each other by stretching or squeezing the covalent bonds.
In the following we describe a molecular model, which incorporates the essential features of the experiment outlined in this section and consider its dynamical behavior under conditions that mimic those studied experimentally 1 .

III. MODEL AND SIMULATION METHODS
In this section we detail the molecular model and the kinetic Monte Carlo simulations that we implemented to investigate the relaxation of the birefringence in the experimental system.

A. Molecular model
As explained in Sec. II, when the monolayer of tethered molecules studied in Ref.
[1] is illuminated by light, molecules cycle between an anisotropic, rod-like trans configuration and a nearly isotropic, bent-like cis configuration [see Fig. 2(a)]. We model the trans isomer as an infinitely thin, straight needle of unit length L = 1 and the cis isomer as a bent, infinitely thin needle also of total length L [see Fig. 2(b)]. In the cis configuration we fix the angle between the central and the tail segments to γ = π/3, while each tail segment has a length of 0.35L. The total number of trans (N t ) and cis (N c ) molecules is fixed, i.e., N t + N c = N . Since the transition moment of the dMR molecules, to which the light polarization couples, is nearly parallel to the monolayer surface, we consider the system to be purely two-dimensional. To mimic the effect of the covalent tethering, molecules are allowed to rotate within the plane but cannot translate.
The isomerization of an azobenzene is a very complicated process when considered in full atomistic detail. The transition between the two isomeric states occurs via a number of intermediate excited states while sev-eral conformational degrees of freedom of the molecule change 40,41 . The situation is even more complicated when the isomerization occurs in a crowded environment because neighboring molecules then most likely interfere 42 .
In our model the isomerization process is drastically simplified: it consists of a simple switch from the straight to the bent needle and vice-versa. Also, we assume that the conformational change happens instantaneously, a reasonable assumption given the relative time scales involved. The relaxation of the birefringence happens at least on the second scale (see Fig. 1), while the isomerization occurs on the picosecond scale even in relatively dense organic solvents 42 .
In order to incorporate isomerization as described in Sec. II in the kinetic Monte Carlo simulation, we define the following set of rules: 1) thermal erasure: only the cis → trans spontaneous transition is allowed and isomerization occurs with probability P th (c → t ).
2) CP erasure: under illumination with circularly polarized light the isomerization rate is proportional to the light intensity, the quantum yield of the transition process, and the absorption cross sections of the isomers 23,24 . To limit the number of free parameters, we incorporate all these factors in the two isomerization probabilities P CP (t → c) and P CP (c → t ) for trans-to-cis and cis-totrans isomerization, respectively. When birefrigence has relaxed towards zero, a steady state is reached, where the numbers of cis and trans isomers fulfill N c P CP (c → t ) = N t P CP (t → c). Hence, the ratio is an essential parameter of our model. Previous measurements show that the absorbance of the two isomers is very similar at 514 nm, therefore R should be in the order of unity 35 . Because light-induced isomerization cycles occur at a much faster rate than the spontaneous cis-trans relaxation, the latter is neglected during CP erasure. Importantly, we assume that after an isomerization event, the orientation of the molecule is chosen at random 32,33 .
As anticipated in Sec. II, because of the high in-plane packing density of the SAM, molecules must overlap during the relaxation process in order to become randomly oriented. Following again a minimal approach, and because we restricted molecular motion to a plane, we allow the molecules to overlap by introducing a simple interaction potential U (O) = U step O, which is proportional to the number of overlaps O and the energetic cost U step of each overlap, which sets the unit of energy. We clarify below how U (O) affects the rotational dynamics of our model molecules.

B. Kinetic Monte Carlo simulation
The system dynamics is generated using a kinetic Monte Carlo algorithm 28 . Rotational dynamics is implemented by picking a molecule at random and rotating it by an angle δθ chosen with equal probability from the interval [−Θ, Θ]. The maximum rotational step size Θ is connected to the molecular self-diffusion constant D θ via a Monte Carlo time step dt 43 The Monte Carlo time step dt is set such that a single Monte Carlo trial move is accepted with a rate close to one, which avoids non-local moves and guarantees a reliable dynamics 28,43,44 . In Sec. IV we fix dt = 10 −4 while in Sec. V we fix dt = 10 −5 to increase the time resolution. Both values ensure an acceptance probability close to one. For the self-diffusion constant of our simple molecular model, we rely on the result for a very long cylinder 45,46 where L is the cylinder length, σ the cylinder aspect ratio and D 0 = k B T /(η S L) with k B the Boltzmann constant, T the temperature, and η S the shear viscosity of the fluid. Since our needles are infinitely thin, we choose σ = 1000 and for simplicity, the rotational diffusion constants of trans and cis molecule is assumed to be the same. Both rotational motion and isomerization take place under the influence of the interaction potential U (O). At each Monte Carlo step a molecule is picked at random and rotated and isomerized using the set of rules defined in the previous section. After a trial move of a single molecule, the number of overlaps and the energy U (O new ) of the new configuration are evaluated and compared to the old configuration with U (O old ). Following the standard Metropolis scheme, the move is accepted with probability In the following we express k B T in units of U step .
A complete sweep consists of N trial moves and the running Monte Carlo time is measured in units of dt. Thus, t MC = N s , where N s is the number of Monte Carlo sweeps. The number density of the model molecules is defined as ρ = N/l 2 , where l, in units of L, is the side length of the square simulation box. For all the results presented in the following, we use a total of N = 10000 molecules under periodic boundary conditions.
As initial condition, we use a configuration, in which both isomers are equally present. The trans molecules exhibit orientational order while cis molecules are randomly oriented. Starting from this configuration, we then follow the relaxation of the birefringence towards equilibrium.

C. Birefringence Relaxation
The degree of alignment within the system is evaluated by the nematic order parameter S(t). Since trans isomers have much higher shape anisotropy than cis isomers 47 , we evaluate the nematic order parameter only for the trans isomers. We calculate the nematic order parameter as the positive eigenvalue of the tensor order parameter where u i α (t) is the α-th Cartesian coordinate of the unit vector pointing along the central segment of the i-th molecule in the trans configuration at time t and . . . denotes non-equilibrium averaging, for which we used at least 10 different runs. For a perfectly aligned system S = 1. The monolayer birefringence ∆n(t) is proportional to both the degree of molecular alignment and the number of rod-like molecules 48 . Therefore, to monitor birefringence relaxation, we keep track of The data points in Fig. 3 give a typical temporal relaxation of Q(t) from our kinetic Monte Carlo scheme. We now illustrate how we discriminate between three possible functional forms of the relaxation dynamics of Q(t). Recall that simple relaxation processes are expected to show an exponential form φ(t) = e −t/τ (a Maxwell-Debye relaxation), which is characterized by a well-defined time τ that fully determines the kinetics of the system. In some systems, however 36 , the relaxation significantly deviates from an exponential form and is described either by a stretched exponential decay φ(t) = e −(t/τ ) β with 0 < β < 1 or by an asymptotic power law as in Eq. (1).
The following procedure is used to discriminate between an exponential, a stretched-exponential, and an asymptotic power-law relaxation of Q(t). For each choice of model parameters, we run kinetic Monte Carlo simulations until Q(t) has reached a clear steady-state value Q eq (illustrated as an horizontal line in Fig. 3). The equilibration time t eq is the first time for which Q(t eq ) = Q eq . Q(t) is then fitted by least-square minimization over the range t ∈ [0, t eq ] with the three functional forms given in Fig. 3.
For each fit curve a goodness-of-fit test 49 is performed. The most reliable fit function is chosen as the one with the value of the reduced χ 2 closest to 1. In the example given in Fig. 3, we obtain χ 2 = 9.65 for the exponential function, χ 2 = 2.32 for the stretched-exponential, and χ 2 = 0.54 for the asymptotic power law. The power law thus clearly provides the best fit of the simulation data. The results of the fitting procedure in different regions of parameter space and under different initial conditions will be discussed in the following two sections.

IV. RESULTS: PURE TRANS SYSTEM
In this section we present the results of the kinetic Monte Carlo simulations for a system that only contains trans molecules.
First, we characterize the equilibrium properties of the system with temperature T and density ρ.  The logarithmic derivative of ∆θ 2 (t) to extract the local exponent ν in ∆θ 2 (t) ∝ t ν . The extent of the subdiffusive regime (ν < 1) increases with decreasing temperature. The inset shows the short-(blue circles) and long-time (red squares) diffusion constants as a function of kBT .
shows the equilibrium value of the nematic order parameter, S eq , determined after equilibration. The main plot refers to molecules with hard-core interactions that cannot overlap at all. The steady-state value of the nematic order parameter shows a steep increase with ρ. The inset in Fig. 4 instead plots S eq versus k B T at ρ = 20. If the temperature is sufficiently high, molecules can pass over their neighbors by creating overlaps and thereby drastically reduce S eq . So, allowing the molecules to overlap, results in an isotropic state at sufficiently high temperatures even at high in-plane packing density.
We now characterize the relaxation dynamics of the birefringence at different temperatures. We fix the initial degree of the nematic order to S(0) = 0.6 (the same qualitative results are obtained using values from S(0) ≈ 0.5 to S(0) ≈ 1.0) and follow the temporal evolution of S(t) while it relaxes back to its equilibrium value S eq . First, we analyze the rotational diffusion of the molecules, which results in the decay of Q(t), by looking at the rotational mean square displacement where θ i (t) is the orientation angle of the i-th molecule at time t. In Fig. 5(a), we plot ∆θ 2 (t) versus t MC for different temperatures in a system with fixed density ρ = 20. In addition, Fig. 5(b) shows the logarithmic derivative of the rotational mean square displacement, which gives the local exponent ν in ∆θ 2 (t) ∝ t ν . Anomalous diffusion has ν = 1. Initially, the molecules diffuse with a diffusion constant D s that increases linearly with temperature, as expected from the Einstein relation [see inset of Fig. 5(b)]. When the needles start to overlap, a subdiffusive regime emerges and its extent increases with decreasing temperature [see Fig. 5(b)]. Ultimately, normal diffusion is recovered even in systems where nematic order is well developed (for instance at k B T = 2.5 one finds S eq ≈ 0.6 in Fig. 4). In order to diffuse, molecules have to pass each other by creating overlaps. At low temperatures, crossing energy barriers makes rotational diffusion an activated process. Hence, the long-time diffusion constant D l is no longer linear in temperature [see inset of Fig. 5(b)]. Rotational diffusion causes initial nematic order to fully decay for k B T ≥ 4. Figure 6(a) shows the relaxation of Q(t) = S(t)/S(0) at fixed ρ = 20 and different k B T , while in Figure 6(b) we keep temperature at k B T = 4.5 and vary density ρ. In a high-temperature or low-density regime the relaxation is very well fitted by an exponential function [k B T ≥ 7.0 in Fig. 6(a) and ρ ≤ 14.0 in Fig. 6(b)], while in a high-density, low-temperature regime [k B T ≤ 5.0 in Fig. 6(a) and ρ ≥ 18.0 in Fig. 6(b)] the power-law decay defined in Eq. (1) provides an excellent fit of the simulation data. In this regime both the exponential and stretched-exponential functions give significant residual errors. Our Monte Carlo results are best fitted by a stretched exponential function with β ≈ 0.8 in an intermediate regime [k B T = 5.5, 6 in Fig. 6(a) and ρ = 16.0 in Fig. 6(b)]. To summarize, the interaction with neighboring molecules, which is more relevant at low temperatures and high densities, causes a transition from an exponential decay of Q(t) to a power-law relaxation.
As anticipated in Sec. III, non-exponential relaxation originates from the presence of a wide distribution of relaxation times in the system. In the concept of dynamic heterogeneity 30,31 , this distribution is traced back to spatial and temporal variations in the local structure of the system, which then determines its dynamic evolution. Molecules diffusing slower or faster than the average become spatially correlated, giving rise to regions with slow and fast dynamics. Hence, averaging the dy- namics over this heterogeneous environment leads to an overall non-exponential relaxation. A typical quantity to monitor this dynamic heterogeneity is a four-point correlation function, which we introduce here for the angular displacement following Ref. [30]. We define the mobility c i (t), in order to quantify how mobile the molecule i is, where ∆θ i (t) 2 is defined as in Eq. (8). The variable c i (t) has zero mean and is positive (negative) if the ith molecule moves less (more) than the average. The

four-point correlation function is then
where | r ij | is the distance between molecule i and j and δ is the Dirac delta function. The correlation function G 4 (r, t) measures the spatiotemporal correlations in the dynamics of the molecules over a distance r at time t. The presence of correlated domains in space and the degree of dynamic heterogeneity is monitored by the dynamical four-point susceptibility for the angular displacement 31,50 , It increases with increasing size of correlated domains.  Figure 7(a) shows the corresponding four-point correlation function G 4 (r, t) for k B T = 4.0 at t MC = 3.0 × 10 4 , where the correlations in time and space are largest compared to other times. In order to monitor the complete temporal evolution of the dynamical heterogeneities, we present the susceptibility χ 4 (t) in the main graph of Fig. 7(a) for different temperatures and at density ρ = 20. The dynamical susceptibility shows a clear peak, which coincides with the time window, during which the power-law relaxation of Q(t) is observed in Fig. 6(a). The color plot of Fig. 7(b) is obtained at the maximum of χ 4 (t) at k B T = 4. The maximum increases with decreasing temperature and therefore demonstrates that regions of correlated motion become more relevant at low temperatures.
Spatio-temporal variations in the local environment of the pinned molecules are also responsible for the anomalous diffusion, which we demonstrate in Fig. 5. In Fig. 8(a) we show individual molecular trajectories for k B T = 4.0 and ρ = 20.0. One already senses that some molecules rotate significantly faster than the average (shown as a thick black line), while all trajectories indicate that rotational diffusion proceeds by sudden jumps. This feature, together with the spatio-temporal correlations discussed before, are common signatures of a glass-like dynamics [51][52][53][54] . In particular, it has been suggested recently that a non-Gaussian distribution of molecular displacements is a universal feature of finitedimensional glass-like dynamics 55 . To investigate this point, we evaluated the non-Gaussian parameter α 2 of the distribution of molecular displacement P (∆θ). In two dimensions it is calculated using the second and fourth moment of P (∆θ), where for a Gaussian distribution one has α 2 = 0. In Fig.  8(b), we show α 2 for different temperatures k B T and at density ρ = 20. The distribution P (∆θ) becomes highly non-Gaussian as also demonstrated in the inset, which shows P (∆θ) for k B T = 4.0 and at the time, when α 2 (t) maximal. Fast rotating molecules are responsible for the non-Gaussian tails of the distribution. The peak in α 2 increases and shifts to later times upon cooling, which again coincides with the developing non-exponential relaxation of the birefringence. Interestingly, the characteristic times in the relaxation laws for Q(t), shown as diamond markers in Fig. 8(a), agree nicely with the locations of the maxima in α 2 (t).

V. RESULTS: RELAXATION WITH ISOMERIZATION
In this section we explore the relaxation of birefringence Q(t) in a system, in which both trans and cis molecules are present. In particular, we investigate how the nearly isotropic cis molecules influence the relaxational dynamics of Q(t) under two different isomerization scenarios: thermal erasure and CP erasure.
The model has a rather rich parameter space and its full exploration is beyond the scope of this paper. In the following, we fix the initial conditions to be close to the experimental values and explore the relaxational dynamics of Q(t) for different isomerization rates. We prepare the initial state with an equal number of trans and cis molecules, N t (0) = N c (0), and order parameter S(0) = 0.6; both values are close to the estimates carried out in Ref. [1]. The cis molecules are randomly oriented at t = 0 and do not contribute to S(t). In this section we choose the time step dt = 10 −5 to increase the temporal resolution, which results in an equilibration time of up to 9 weeks on an Intel Xeon X5550 machine with a 2.66 GHz CPU. In order to compare the results of this section directly with the results of Sec. IV, we still give the Monte-Carlo time t MC in units of dt = 10 −4 .

A. Thermal erasure
Here, we discuss the relaxation of the birefringence, when the monolayer is not illuminated. Since the trans isomer is the ground-state of the dMR, all the molecules in the cis configuration will isomerize back to the trans state after some characteristic time. As discussed in detail in Sec. III A, the cis to trans isomerization rate of isolated molecules is the isomerization probability P th (c → t ).
In view of the results discussed in Sec. IV, we expect the relaxation dynamics to be exponential in the hightemperature and low-density regime, regardless of the isomerization probability. Therefore, we fix both temperature and density to k B T = 4.0 and ρ = 20.0, respectively, where the pure trans system shows a clear power-law decay of the birefringence, and monitor how this decay is influenced by the isomerization rate. Figure 9(b) shows the temporal evolution of the number of trans isomers for different P th (c → t ). We expect an exponential relaxation and, indeed, the Monte Carlo data (circles) are well fitted with where the fit parameter, the relaxation time τ c , is proportional to the inverse isomerization probability, τ c = 1.54 P −1 th (c → t) [see inset of Fig. 9(b)]. We find τ c > P −1 th (c → t) because in a crowded environment some of the attempted isomerization events are rejected as they generate more overlaps between the molecules. In Figs. 9(a) and (c) we show the respective temporal evolutions of the nematic order parameter S(t)/S(0) and the birefringence Q(t) ∝ N t (t)S(t), which originate from the alignment of the trans isomers. We also evaluate the dynamical susceptibility χ 4 (t) defined in Eq. (11), but only on the subset of molecules that are in the trans configuration at the beginning of the simulation at t MC = 0. The results are plotted in the inset of Fig. 9(c), where we also include the dynamical susceptibility for the pure trans system at ρ = 20 at k B T = 4.0 (black circles) as a reference. This will allow us to to quantify how the presence of the cis isomers influences the development of dynamic heterogeneities. The temporal relaxation of both S(t) and Q(t) strongly depends on the isomerization rate. Analyzing Figs. 10(a), (b) and (c), we identify four different regimes: The full lines are the best fits to an exponential (magenta), a stretched exponential (green), or a power law (red). Diamonds indicate the characteristic times of the power law relaxation. The inset shows the dynamical susceptibility χ4(t) for trans molecules already present at t = 0. As a reference the black circles give χ4(t) for the pure trans system at the same ρ and kBT . 1) For sufficiently small relaxation rate [P th (c → t ) = 10 −5 in Fig. 9], the birefringence Q(t) relaxes exponentially. On the time scale of the declining S, the number of cis molecules stays constant [compare plots (a) and (b)]. They thus create a more uniform environment, as indicated by the nearly vanishing dynamical susceptibility χ 4 , in which S(t) and Q(t) relax exponentially and much faster than the pure trans system for the same temperature and density in Sec. IV. Note that at late times Q(t) first falls below its equilibrium value (dashed horizontal line) and then increases again. The reason is that the equilibrium value of S depends on the number of trans molecules [inset of Fig. 9(a)] because the presence of the nearly isotropic cis molecules decreases S eq . Since N t increases for t MC > 10 4 , Q(t) also increases.
2) For P th (c → t ) = 10 −4 and 2.5·10 −4 the best fitting function for the temporal evolution of Q(t) is provided by a stretched exponential with β ≈ 0.8 indicating the transition to the power-law decay.
3) At intermediate isomerization rates P th (c → t ) = 3.5 ·10 −4 and 5.0 ·10 −4 the power-law relaxation provides the best fit of the simulation data. Here, the isomerization of the randomly oriented cis molecules into the trans state happens on the same time scale as the relaxation of the nematic order parameter. Thus, Q(t) does not decay in the static disordered distribution of cis molecules but in a dynamic and heterogeneous environment, as demonstrated by the clear peak in χ 4 (t). As a result, the relaxation of Q(t) follows a power law. Its characteristic time τ t does not change significantly, but the isomerization probability seems to control the power-law exponent. We find η = 0.281 at P th (c → t ) = 3.5 · 10 −4 and η = 0.225 at P th (c → t ) = 5.0 · 10 −4 .
4) The situation changes again if the isomerization rate is very large [P th (c → t ) = 10 −2 in Fig. 9]. All the cis isomers rapidly isomerize into trans molecules with random orientation and S(t) drops to S(0)/2. This is compensated by the resulting increase in N t (t), which ultimately generates a bump in Q(t) at t MC ≈ 10 2 . Once isomerization is completed at t MC ≈ 10 3 , the system is composed of only trans molecules the aligment of which relaxes via rotational motion. This is the same situation as discussed in Sec. IV. The environment is heterogeneous as indicated by the peak in the dynamical susceptibility, which is nearly as large as in the pure trans system of Sec. IV. Without the initial bump, Q(t) is fitted well by a power-law with a larger characteristic time and a larger power-law exponent η ≈ 0.5 as compared to case 3) but similar to the pure trans system.
Summarizing the results in Fig. 9, we find that cis isomers pinned at random positions not only accelerate the birefringence relaxation compared to the pure trans sys- tem but also prevent the development of the dynamical heterogeneity that is responsible for a non-exponential decay. We attribute this behavior to the nearly isotropic shape of cis isomers, which, regardless of their orientations, create a similar environment for the exisiting trans isomers. In contrast, newly formed trans isomers can adjust their orientations to their neighbors and thereby allow for the formation of dynamical heterogeneities. In the experiments of Ref. [35] the lifetime of cis isomers during thermal erasure is estimated as 0.7s, while the characteristic time of the power-law decay is measured as τ th ≈ 2s. These values are achievable by isomerization rates between our cases 3) and 4).

B. CP erasure
Illumination of the SAM with CP light induces transcis isomerization cycles at a rate much faster than the spontaneous relaxation due to thermal erasure. Therefore, we neglect the thermally induced cis-to-trans transition when modeling the CP erasure process. As discussed in Sec. III A, in order to take into account the different light absorption of the two isomers, we intro-duce the two respective probabilities, P CP (t → c) and P CP (c → t ), for trans-to-cis and cis-to-trans isomerization. Since light is circularly polarized, the isomerization probabilities do not depend on the molecular orientation.
In order to limit the computational cost, we do not explore the full range of isomerization probabilities. Instead, we choose their values such that the characteristic times of birefringence relaxation for thermal and CP erasure matches the experimental observations, where they differ by approximately two decades in time (see the diamond markers in Fig. 1). Because illumination of the monolayer produces a very negligible amount of heat 1 , we fix both density and temperature at ρ = 20 and k B T = 4, exactly as during thermal erasure discussed in Sec. V A.
In Fig. 10 we show the relaxation of Q(t) and N t (t) towards their steady-state values starting with N t = 0.5 at t = 0. In Fig. 10(a) we set P CP (t → c) = 5.0 · 10 −1 while in Fig. 10(b) P CP (t → c) = 10 −1 . The backward isomerization rates, P CP (c → t ), are chosen by the ratio R = P CP (t → c)/P CP (c → t ), which also determines the number of isomers in steady state: R = N c (t → ∞)/N t (t → ∞), as discussed in Sec. III A.
The implementation of our kinetic Monte-Carlo simulations suggests that the number of trans isomers evolves according to the following kinetic equation, where τ ct and τ tc are characteristic relaxation times and N c (t) = N − N t (t). It is solved by with a = 1 1 + τ ct /τ tc , a + b = N t (0) N , and 1 τ = 1 τ tc + 1 τ ct (16) Since a is the steady-state value at t → ∞, we also have a = 1/(1 + R), which is confirmed in Fig. 11(a). Thus, the ratio of isomerization probabilities also determines the ratio of the two relaxation times, R = τ ct /τ tc . The fits of the simulated N t (t)/N in Figs. 10(a) and (b) excellently confirm the kinetic model. However, as in the case of thermal erasure, we find τ ct > P −1 CP (c → t) [see Fig. 11(b)], because in a crowded environment some of the attempted isomerization events are rejected.
In Fig. 10 we show the best fitting functions for Q(t) as continuous lines. In both Figs. 10(a) and (b) the relaxation of the birefringence is exponential for R < 0.7, stretched-exponential for 0.7 ≤ R < 0.3, and follows a power-law for R ≤ 0.3. The characteristic times of the relaxation are shown as diamond markers and the exponents are given in Table I. As expected, a larger isomerization rate P CP (t → c) shifts the birefringence relaxation to smaller times [compare Figs. 10 (a) and (b)] because aligned trans molecules are faster transformed to the cis state. Unlike the case of thermal erasure, the steady value of N t (t)/N is reached well before stretched-exponential or power-law relaxation sets in. The relevant characteristic times τ t are indicated   Fig. 10. Table  on the left gives the exponent for the stretched-exponential fit (green lines) and table on the right gives the exponent for the power-law fit (red lines).
by diamond markers in Fig. 10. Interestingly, for constant P CP (t → c) the ratio of isomerization rates R controls the functional form of the relaxation. In the powerlaw regime, the characteristic times τ t do not change significantly, while the power-law exponent η heavily depends on R (see Table I).
At a first glance, the behavior in Fig. 10 seems surprising. For increasing R the isomerization rate P CP (c → t ) = P CP (t → c)/R decreases and consequently the molecular orientations after isomerization become randomized less frequently. So, we expect the birefringence relaxation to become slower in contrast to the results presented in Fig. 10. However, we know already from thermal erasure that cis isomers create a uniform environment, where the orientation of trans molecules relaxes faster and exponentially. This is the case for R > 1, where the cis isomers are in the majority in steady state. Decreasing R increases the number of rod-like trans molecules, N t . As discussed already for thermal erasure, they hinder the orientational relaxation of their neighbors more efficiently than cis molecules. But they also create a more heterogeneous environment, where the birefringence relaxation first follows a stretched exponential and then for further decreasing R becomes a power law. However, we were not able to quantify the dynamic heterogeneity using the dynamical susceptibility of Eq. (11) for a subset of molecules as in Sec. V A since they continuously cycle between their two configurations.

C. Comparison between thermal and CP erasure
In Fig. 12 we show the power-law relaxation of Q(t) for both thermal and CP erasure using typical parameters. Comparing it with Fig. 1 gives an idea of the degree of agreement between our model and the experimental results. Birefringence relaxation is efficiently accelerated by the isomerization cycles induced by illumination with CP light. A similar speedup was demonstrated in Refs. [56] and [57], where isomerization of an azo-dye embedded in a molecular matrix significantly in- creased translational diffusion of surrounding molecules. By tuning the isomerization probabilities for CP erasure, we achieved a difference between the characteristic times τ t in the power-law decay of approximately two orders of magnitude, in good agreement with the experimental results.
The larger isomerization rates during CP erasure and the presence of cis isomers in steady state give a faster power-law relaxation with a larger exponent η, again in qualitative agreement with the experimental results. Our model also accounts for the smaller steady-state value of Q(t) after illumination with CP light as reported in Ref.
[1]. This is due to the presence of the nearly isotropic cis molecules that destabilize orientational order.

VI. SUMMARY AND CONCLUSION
Experiments on SAMs with tethered light-switchable dye molecules show power-law relaxation of initial birefringence during both thermal erasure and the faster CP erasure. Despite its simplicity, the molecular model discussed in this paper is able to reproduce the experimental results and to identify dynamic heterogeneity as the main cause for the power-law decay.
First, we studied a system of pure trans molecules. Here, the non-exponential, glass-like relaxation of the orientational order inscribed in the monolayer emerges naturally at high densities and low temperatures due to the presence of dynamic heterogeneity. Rotational motion develops a transient sub-diffusive regime upon cooling. At the same time, molecules with dynamics faster and slower than the average become spatially correlated. The spatial average over these regions with different ori-entational mobilities results in the power-law decay of birefringence.
In a second step we included the possibility that the molecules can assume two different isomeric forms. During thermal erasure, the nearly isotropic cis isomers create a uniform environment because they do not align locally with their neighboring molecules. Hence, the orientation of trans molecules relaxes exponentially. The experimental power-law relaxation of birefringence is only recovered if a sufficient number of trans isomers is present. They slow down orientational relaxation but also initiate the formation of dynamic heterogeneities as in the pure trans system.
During CP erasure light adsorption induces a fast isomerization cycle between cis and trans isomers and thereby the overall orientational relaxation becomes faster in agreement with experimental results. The functional form of the birefringence relaxation is controlled by the ratio of the two isomerization probabilities, which determine the number of trans and cis molecules in steady state. As in thermal erasure, a larger number of cis isomers speeds up the exponential birefringence relaxation, whereas trans isomers in the majority hinder relaxation and ultimately give rise to a power-law decay. Finally, the presence and nearly isotropic cis isomers also explains the smaller steady-state value, which birefringence reaches during CP erasure. All these findings strongly suggest the possibility to change the monolayer dynamics by tuning the absorption properties of the molecules and their geometrical shape.
To reproduce the power-law decay of birefringence in our model, we have to fine-tune the parameters within a range of values that are experimentally reasonable. It remains to be demonstrated if this is due to our simplified model or a general feature of light-switchable molecules. Future work should address this question by exploring the effect of more complex molecular geometries with more realistic molecular interactions. Another promising direction for exploring how light can be used to control material properties are light-switchable sufactants [58][59][60][61] . They accumulate at fluid interfaces. By switching locally between the two isomeric states, the surface tension changes and its gradient drives Marangoni flow. This moves emulsion droplets along a surface or in bulk with interesting non-linear dynamics 62,63 .