Defect dynamics in active smectics induced by confining geometry and topology

The persistent dynamics in systems out of equilibrium, particularly those characterized by annihilation and creation of topological defects, is known to involve complicated spatiotemporal processes and is deemed difficult to control. Here the complex dynamics of defects in active smectic layers exposed to strong confinements is explored, through self-propulsion of active particles and a variety of confining geometries with different topology, ranging from circular, flower-shaped epicycloid, to hypocycloid cavities, channels, and rings. We identify a wealth of dynamical behaviors during the evolution of complex spatiotemporal defect patterns as induced by the confining shape and topology, particularly a perpetual creation-annihilation dynamical state at intermediate activity with large fluctuations of topological defects and a controllable transition from oscillatory to damped time correlation of defect number density via mechanisms governed by boundary cusps. Our results are obtained by using an active phase field crystal approach. Possible experimental realizations are also discussed. It is intriguing but challenging to control the complex dynamics of topological defects that emerge in nonequilibrium spatiotemporal systems. Through an effective modelling this work explores a viable way to controllably vary the properties of persistent defect dynamics in self-propelled active smectics via the competition of activity and boundary confinements with different topology and geometry.


Introduction
Defects in ordered or pattern-forming systems are of great interest both from a fundamental physics point of view highlighting the role of topology in condensed matter [1][2][3][4][5] and for applications since they largely control the material properties.For the example of liquid crystals, most of the studies have focused on topological defects in the orientational ordered nematic phase of lyotropic or thermotropic liquid crystals and by now it has been well understood how to trigger them by external influences [6][7][8][9][10] and confinements [11][12][13][14][15][16].In the layered smectic phase, defect characterization is more complex due to the additional positional ordering [17][18][19][20] but can be varied by confinement as well [21,22].
In recent years, active particles that are self-propelled intrinsically and relevant for biological swarms, motor proteins and biofilaments have been of tremendous interest [23][24][25][26].These particles self-organize into fascinating "active" liquid crystalline phases which are qualitatively different from their passive counterparts in or near equilibrium and are governed by nonrelaxational dynamical processes that are the characteristics of far-from-equilibrium pattern-forming or chaotic systems [4].One important example of active matter systems involves active nematics showing orientational order.The dynamics of topological defects in active nematics have been analyzed to a large extent and a plethora of new phenomena were discovered (see, e.g., Refs.[27][28][29][30][31][32][33][34][35][36]).In particular, effects of confinement were explored for channels and capillaries [37][38][39][40][41][42][43][44], resulting in e.g., nontrivial dancing motion of defects [45].However, confinements different from channels were only rarely addressed [46].In parallel, effects of topology in the confinement have been studied for active nematics [47,48].Also active layered smectic-like states which exhibit an additional positional order have been examined [49,50].These systems are modeled by aligning [51] and nonreciprocal [52] interactions or nonlinear feedback [53].Due to the nontrivial coupling between orientational and positional degrees of freedom with the active driving force, active layered smectic systems and their defect dynamics are much more complex than their nematic counterparts.To the best of our knowledge, studies of defects in active smectics, including their dynamics and controllability by external constraints, are still sparse.Moreover, the effect of sharp cusps in the confinement shape has never been addressed for any active liquid crystalline or active pattern systems.
Here, we contribute to fill this gap and explore defect dynamics of active smectic systems by using an active phase field crystal (PFC) model [54][55][56][57][58][59][60][61] in a parameter setting where a traveling stripe phase is stable.We expose this state of self-propelled smectics to strong confinements where the smectic layer width is getting comparable to or not far from the confining length scales.A rich variety of boundary geometries with different topologies, including cavities and open or closed channels that are of various convex, concave, and cusped shapes, are considered.Our purpose here is to attain a systematic understanding of the complex dynamics of defects in the confined active smectics, particularly the effects of both geometry and topology.The confinements examined can be classified into two types of topology, closed cavity vs channel/ring (including closed rings and open channels with periodic boundary condition along open ends).Each of them includes different types of confining geometries with different number of sharp, singular cusps, either inward cusps (in epicycloid cavities or rings or cycloid open channel) or outward cusps (in hypocycloid cavities), or both combined (as in hypocycloid rings).Each type of topology also involves smooth boundaries with the lack of cusps (such as circular cavity, S-shaped open channel, and annulus).
In this study a wealth of nonequilibrium defected states are found in the evolution towards complex spatiotemporal patterns arising from a competition of activity and confinement.The mechanisms governing the dynamics of defects at large enough activity strongly differ from those of passive patterns and active nematics.The persistent defect dynamics identified here, especially the highly fluctuating defect state, goes far beyond the traditional classification familiar in passive systems, and is shown to be induced by the confining shape and topology of the cavity or channel and the degree of particle self-propulsion.This dynamical regime of high defect fluctuations occurs at intermediate activity as characterized by the perpetual process of defect creation-propagation-annihilation, showing as intermittently varying time stages with bursts of defects emerging in most of boundary confinements other than the smooth-boundary channels without cusp (i.e., S-shaped channel and annulus).In particular, the presence of cuspated boundaries induces or annihilates defects, which can be utilized to control the dynamics of defect creation and the quantitative behavior (oscillatory vs damped) of time correlation of defect density.Our predictions can be verified for confined dense vibrated granular rods [62][63][64] or self-propelled colloidal Janus particles [65,66] exposed to strong confinements [21].

Model
We describe the evolution of active smectics under confinement based on a continuum density-field theory, i.e., the active PFC model which can be derived from dynamical density functional theory [54,55] and also from a particlebased microscopic description [60].It reads where ψ is the particle density variation field, the polarization P represents the local orientation vector field, v 0 measures the strength of particle selfpropulsion, and D r is the rotational diffusion constant.The above dynamical equations have been rescaled, with a diffusive timescale and a length scale set via the pattern periodicity.The model considers conserved dynamics for ψ, as seen in Eq. ( 1), and thus the average density ψ 0 remains unchanged during the system evolution.We set F = F aPFC + F anch , where F aPFC is the rescaled free energy functional of active PFC [54,55] with < 0, the characteristic wave number q 0 = 1 after rescaling, and C 1 > 0 tending to suppress any spontaneous ordering of orientational alignment.and the average density ψ 0 are chosen to give rise to the resting or traveling active smectic phase [54].
We represent the effect of boundary confinement via an anchoring energy to effectively satisfy both Dirichlet and Neumann boundary conditions ψ = ψ b , P = P b , and n • ∇ψ = 0 (with n the local unit normal) at any implicitly defined domain boundary, with where V b0 gives the anchoring strength, r s (r) is the signed distance function to the domain boundary (with r s < 0 inside the domain and > 0 outside), and ∆ sets the thickness of boundary interface.This approach we develop here combines an approximation of domain interface energy for imposing the boundary conditions and a setup of boundary via V b (r) = V b0 φ(r) with an auxiliary phase field function φ(r) used in the diffuse domain method [67] to control the confinement geometry implicitly.More details, including the specific analytical forms for different geometries of the cavities or channels simulated, are given in the Methods section.Among them, the geometry of epicycloid or hypercycloid with integer n cusps is described as a closed plane curve formed by the rolling of a small circle of radius b on the outside or inside of another larger fixed circle of radius a = nb, respectively, with their parametric equations given in Eqs.(26) and (27) of the Methods section.Equation (4) produces the condition of planar anchoring as found in experiments.Its last term is analogous to the Rapini-Papoular form of surface potential [3,14].In our simulations (starting from random initial conditions), we set ( , ψ 0 , D r , C 1 ) = (−0.98,0, 0.5, 0.2) in the strong segregation regime of stripe phase, and (V b0 , ∆, ψ b , P b ) = (1, 0.1, 0, 0) for the cavity or channel boundary setup.We have tested stronger anchoring strength with larger value of V b0 and found that larger self-propulsion strength v 0 would then be needed to overcome the stronger boundary confinement particularly for defect nucleation, while the corresponding results obtained are qualitatively similar.

Defects dynamics in closed cavities
The emergence of topological defects, including disclinations, dislocations, and grain boundaries, is observed in our simulation systems (see Fig. 1 for some typical smectic defects, similar to those found in passive smectic or stripe-pattern systems [1][2][3][4][5]).Their complex behaviors of dynamical evolution is further complicated by the effects of active self-driving and boundary confinement.In the confined cavities of different geometries, our results presented in Fig. 2 show that the evolution of these defects in active smectics is governed by three intrinsically different dynamical regimes.At weak enough self-propulsion strength v 0 (first column of Figs.2(a)-2(c)), the stripes (smectic layers) remain perpendicular to the boundary, showing planar anchoring with tangential alignment of constituent particle orientations.The defects emerging from the early stage of system evolution become mostly pinned, with extremely slow local dynamics, similar to the glassy state observed in strongly segregated passive stripe patterns showing no long range orientational order as a result of defect pinning by the pattern-periodicity induced potential barrier [68].
In the other limit of large enough v 0 , the particle self-propulsion completely overcomes the defect pinning barrier, facilitating the fast annihilation of defects as accelerated by the effect of self-driving on the smectic layer alignment.As shown in the last column of Figs.2(a)-2(c), large particle activity also enables the overcoming of boundary anchoring constraint, leading to the violation of local planar anchoring, and depending on the boundary geometry, even to (partially) homeotropic anchoring.At late stage such a strong self-driving induces the rotation of smectics inside the cavity, either clockwise or counterclockwise, and hence the persistent self-rotation of the remaining multi-core spiral defects trapped at the cavity center.The transition between these different dynamical regimes occurs within a narrow range of the self-propulsion strength v 0 (middle column of Figs.2(a)-2(c)).To understand this transition between localized and traveling smectic patterns, we focus on the bulk state and neglect the boundary anchoring energy.This allows us to rewrite Eqs. ( 1) and ( 2) via defining a local polarization divergence field S = ∇ • P, i.e., Working in a comoving frame with r → r = r − v m t and assuming ψ = ψ(r − v m t) and S = S(r − v m t) in the nonequilibrium steady state of pattern dynamics, where v m is the migration velocity of smectic layers, we can obtain the equations governing the Fourier components ψq and Ŝq of the ψ and S fields as In one-mode approximation the density field for a perfect stripe phase is given by ψ , where q s is the selected wave vector of the pattern and A = A 0 e iφ 0 with a constant phase φ 0 .
Here we have assumed that the model parameters (including , v 0 , C 1 , and D r ) are within the range of the traveling state of active smectic pattern examined in this work.For mode q s Eqs. ( 8) and ( 9) then become Separating real and imaginary parts of Eq. ( 10) leads to the following solutions for velocity v m and amplitude A 0 : When |v 0 | ≤ v 0c , we have v m = 0 corresponding to localized stripe patterns, with while at large enough activity |v 0 | > v 0c , the pattern travels with nonzero velocity v m as determined by where v 0c is the critical threshold of activity for the transition, as given by For a nonpotential, nonrelaxational system like the active system studied here, the selected wave number q s of an ordered pattern in the long-time steady state cannot be determined from free energy minimization and is difficult to identify through analytical calculations [4].In the active PFC model the value of q s is expected to be near q 0 as also found in our numerical simulations.
Substituting the parameter values used in our simulations (i.e., C 1 = 0.2 and D r = 0.5) and approximating q s ∼ q 0 = 1, we get v 0c = 0.3 from the above analytic result, the same as what has been found in simulations and presented in Figs.2(a)-2(c).
The direction of pattern traveling velocity v m and the orientation of the selected wave vector q s (which tend to locally align with each other) highly depend on the initial and boundary conditions, and could vary in different parts of the system due to the effect of boundary confinement, as seen in our numerical simulations (e.g., Fig. 2 and Supplementary Movies 1-7).
Defect dynamics in these different regimes (localized and traveling smectic patterns and their transition) reveal the competition between rigid boundary confinement restricting the local smectic orientation and the tendency of bulk alignment of stripes [69,70].An interesting type of dynamics occurs in the transition regime when such two incompatible boundary and bulk effects are of comparable strength, giving rise to a highly fluctuating state (middle column of Figs.2(a)-2(c)).Although the planar boundary anchoring is maintained at very early stage, the deviation occurs at later times as caused by self-propulsion, leading to local distortion of stripes inside the cavity as a result of confinement-alignment competition.Importantly, in addition to the annihilation of defects (including dislocations and disclinations, majorities of which occur at cavity boundaries), new defects can nucleate from the boundary, propagating into the bulk, evolving and generating a subsequence of more new defects like a chain effect, as seen in the circled regions of Fig. 2 and Supplementary Movies 1-7 for epicycloid and hypocycloid cavities.This results in the repeated succession of tranquil and active time stages in terms of defect density and dynamics, with some examples for circular and 6-cusp epicycloid and hypocycloid cavities given in Fig. 2(d).In contrast to the fully bulk state without any boundary confinement (thus with the absence of defect generation) which shows a monotonic time decay of defect number, the cavity confinement induces an intermittency-type behavior with seemingly irregular bursts of number of defects.The boundary cusps appear to enhance the creation of new defects, yielding higher defect density peaks, as compared to the smooth boundary of circular cavity.It is noted that the overall system activity is governed by v 0 which is kept unchanged for different confinement geometries and cusp number, although there would be local variations of effective activity as a result of complex nonlinear defect dynamics particularly defect creation and annihilation.
The property of this transition zone with dynamical fluctuations of defects can be further quantified through the normalized time autocorrelation function of the defect number density n d , i.e., where the averages are conducted over a long time series in the steady state (e.g., t = 10 5 − 10 6 in our calculations) for each simulation run, assuming ergodicity of the corresponding probability measure [4].Some results of C n (τ ) are presented in Fig. 2(e), showing a decay behavior for circular and 6-cusp hypocycloid cavities.Interestingly, for epicycloid cavity with n = 6 cusps a weak oscillation around a negative minimum correlation (near time scale τ m ∼ 5000) appears, implying a correlated behavior between the burst (active) and low-number (tranquil) regimes of defect density and dynamics.For a given geometric type of confinement, the behavior of defect autocorrelation can be qualitatively changed through different number of boundary cusps.As shown in Fig. 3, for epicycloid cavities decreasing the cusp number n from n = 6 to 3 leads to the variation of C n (τ ) from local negative minimum to positive maximum of correlation.A more dramatic change occurs when lowering the cusp number of hypocycloid cavities.When n = 5 and 6 a damped time correlation of defect density is observed, while the n = 4 (i.e., astroid) cavity is featured by an oscillatory behavior (within the statistical error) of time correlation, indicating a cyclic state of defect density variation with periodic creation and annihilation of defects over a characteristic time period τ T ∼ 4300.In this cyclic defect dynamics, although generally the spatial locations of boundary defect nucleation and annihilation seem uncorrelated, statistically the periodicity in autocorrelation C n (τ ) can be attributed to the propagation of defects between different sides of boundary within the defect creation-annihilation time interval.
Figure 4 gives the corresponding power spectrum of the defect number density for some closed cavities with different number of cusps.It has been averaged over independent simulation runs as calculated via where nd (ω) is the temporal Fourier transform of defect density n d (t) and L t is the length of n d time series which needs to be long enough particularly for non-periodic variation of n d (in our calculations a time range of t = 10  behavior of autocorrelation C n (τ ) then corresponds to a single high peak in the power spectrum plot, as shown in Fig. 4 for the n = 4 hypocycloid cavity with a peak located around ω ∼ 2π/τ T = 0.00146, consistent with its C n (τ ) result given in Fig. 3.
The complex behaviors of defect dynamics revealed above for active smectic systems are governed by intrinsically different mechanisms compared to those examined in previous works of defect dynamics in stripe or convection-roll patterns of passive systems [4].Most of those passive model systems, either potential or nonpotential, are governed by nonconserved dynamics, such as the original or generalized Swift-Hohenberg equations with or without the coupling to hydrodynamic mean flow and vertical vorticity [68,[71][72][73][74][75].For conserved dynamics as studied here, the treatments are more complicated and some related developments for the study of dislocation motion in passive crystalline systems were available only recently [76,77].The considered approach was based on amplitude equation expansion in the weakly nonlinear regime and applied to the pattern near onset.This is different from the active system studied here which is nonpotential and nonrelaxational and is far from onset with strong segregation.The corresponding amplitude equation formulation and the related analysis for strongly segregated patterns would be much more involved and need the incorporation of nonadiabatic effects for lattice pinning [68,78,79].
The defect dynamics here also differ from those known in active nematics, where they strongly depend on the type of the defects and the interactions between them [26][27][28][29]34].In active smectics the defect type and the interactions between individual defects are only of secondary effect for the parameter range and the corresponding dynamical regimes examined in this work.At low activity (i.e., v 0 ≤ v 0c with the threshold v 0c determined by Eq. ( 15)), since the system examined here is in the strong segregation regime (with = −0.98), the defects are pinned to the underlying periodic structure of the stripe pattern (similar to that of a passive potential system exhibiting a stripe phase in Ref. [68]); thus the driving forces caused by the intrinsic interactions between defects (between either different dislocations or disclinations) are much smaller than the pinning force and do not overcome the pinning barrier.When v 0 > v 0c at large enough activity which is the main focus here, the active self-driving force overcomes the pattern pinning effect and when combined with the effect of boundary confinement, dominates the evolution and motion of all the defects as it clearly exceeds the pinning force and hence far exceeds the defect-defect interactions which now play a secondary role.This can be observed in numerical simulations (see Supplementary Movies 1-7), where at the leading order the motion of defects mostly follows the traveling of local stripe layers as driven by the self-propulsion, accompanied by further defect generation/splitting or annihilation through the coupling to rigid boundary confinement and the effect of large activity.Some sample trajectories of defect motion are illustrated in Fig. 5 for an n = 4 hypocycloid (astroid) cavity at v 0 = 0.31.The dynamics starts from a single defect created at the upper-right boundary of the cavity; then a chain of new defects is generated sequentially during the traveling into the bulk (see Fig. 5(b) and the corresponding Supplementary Movie 5).This procedure of new defect generation is caused by the active driving and local distortions of the smectic layers which depend on the specific geometry and topology of the boundary confinement (with more related studies given in the next section), and their motion is subjected to the flowing of local stripes as observed in Supplementary Movie 5.
The self-driving force for the time evolution of density field and the resulting defect dynamics is determined by the local polarization divergence term v 0 ∇ • P as can be seen from Eq. ( 1) or (6).As an example, in Fig. 5 we show the spatial distribution of the polarization field P and its divergence field S = ∇ • P (which is in turn coupled to the local variation of density field through v 0 ∇ 2 ψ as given in Eq. ( 7)) around a defect core in a four-cusp astroid cavity.As illustrated in Fig. 5(c), the vector field P represents the distribution of local orientational order of active-driving directions, with the net self-propulsion determined by the asymmetric distribution of the P field surrounding the local density peak of ψ and the corresponding net orientation of P [54,55].The information of self-driving can then be obtained from the local spatial gradients of P, with ∇ • P being of similar pattern as the density field ψ (see Fig. 5(d)).The only difference is a nonzero phase shift between them giving the effect of self-propulsion.This can be also seen from Eq. ( 9) or (11) showing their Fourier components ( ψq and Ŝq ) being proportional to each other but with a phase difference when the net migration velocity v m = 0.

Topology-dependent defect dynamics and cusp-induced mechanism of defect generation
To further understand the mechanism of defect creation and correlation and hence the accessibility for varying defect dynamics, we examine the process of defect flow for a different topology including two types of open channels, the S-shaped channel with smooth boundary and the cycloid channel with a single cusp (see Fig. 6, noting the periodic boundary condition along the vertical open ends).No new defects are nucleated from the smooth boundary of S-shaped channel and the defect number decreases with time, with few defects left at late stage, as seen in Figs.6(a) and 6(b).In contrast, during the flow of stripes in the cycloid channel (Figs.6(c)-6(e)), the cusp singularity enhances the local distortion of the smectic layers, thus enabling the formation of defects at boundary (not necessarily at the cusp location).Further distortion during flow can facilitate the defect motion into the bulk and induce a chain of new defects (see Fig. 6(e)).This burst of defects will be diminished due to their annihilation when traveling to the channel boundary, with few or none remained.After then the similar nucleation-propagation-annihilation process repeats, resulting in the periodic variation of defect number density as shown in Fig. 6(a).
For the confined cavities studied in the last section, the enclosure constraint without open ends, competing with self-propelled alignment and domain flow in the bulk, leads to a high degree of local pattern distortion which enables the defect nucleation at boundary even for no-cusp, circular cavity of large enough aspect ratio between the lateral cavity dimension and the smectic layer spacing ( 25 as estimated from simulations; see Fig. 2).This is different from the above results for open channels (and the results below for rings or annulus) where no defect nucleation would occur at smooth boundaries without cusps, due to different confinement topology.The mechanism originated from cusp singularity of confinement would play a key role on enhancing defect generation, and importantly, on controlling the time-correlated property of defect variations, including the oscillatory behavior of defect correlation for cavities with small number of cusps as observed in simulations.When the cusp number increases the oscillation of time correlation function would be damped, and thus the degree of defect periodic variation be reduced, as a result of the interference between the effects induced by different individual cusps.This interference effect can account for the transition from oscillatory to non-oscillatory decay of the correlation shown in Fig. 3.In the other limit of zero cusp without the oscillatory mechanism, such as the circular cavity, faster decay of correlation is found (Fig. 2(e)).These results further indicate that the mechanism generated by cusp singularity provides an effective route for controlling the collective property of defect dynamics and correlation.Applying this mechanism, one can expect that further confining via both inner and outer boundaries, i.e., a closed channel with a void, would result in the damping of correlation due to greater constraint and larger degree of interference between boundaries, although with similar processes of defect generation, traveling, and annihilation.This is seen in Fig. 7 where the single-cusp epicycloid ring can be viewed as a curved analog of the cycloid channel given in Fig. 6.For small system size (e.g., 256 × 256 grid points, same as Fig. 6) with narrow channel width, much weaker oscillation and a damped behavior of time correlation occurs, as compared to the time-periodic behavior for the open cycloid channel (see Fig. 7(b)).Interestingly, wider channel width leads to longer period of time correlation of defect number variations, as seen from the C n (τ ) result presented in Fig. 7. (It is possible that defects could be trapped inside a wide enough epicycloid channel, without any new boundary defects nucleated, as found in roughly half of the independent simulation runs conducted for system size of 512 × 512 grid points.Those cases are not used in the calculation of C n (τ ) in Fig. 7.) In comparison, for annulus with smooth boundary, no defect multiplication occurs, analogous to the case of smooth Sshaped channel, and the bulk defects could self-circle persistently as a result of self-propulsion (see Supplementary Movie 9 for an example of annulus at v 0 = 0.31).
Higher degree of correlation damping is expected for closed channels/rings with more cusps interference, as verified in Fig. 8 for hypocycloid rings with n varying from 3 to 7. Weak oscillation of C n (τ ) is found for n = 3 and 4; between them the n = 3 hypercycloid ring exhibits faster decay of correlation over time, which can be attributed to its narrower channel width.Larger number of cusps (e.g., n = 6, 7) results in non-oscillatory, damping behavior of correlation due to increasing interference between cusps, with faster damping for larger n (see Fig. 8 if comparing n = 7 to n = 6), as expected.On the other hand, the confinement of closed channel leads to the increase of the correlation time as compared to the closed cavity, as seen from the C n (τ ) plots for n = 6 hypocycloid cavity vs ring in Figs. 3 and 8.
The above results for cusped channels or rings correspond to the intermediate regime of activity with v 0 larger than but near v 0c showing high defect number fluctuation and defect creation-annihilation.In another dynamical regime of larger activity, the effect of self-driving dominates over that of rigid boundary confinement, leading to the absence of defect multiplication as in the case of closed cavities.In the closed channels of annular or ring geometries the defect dynamics manifests itself either via local persistent variations involving bulk grain boundaries and spirals confined by the channel walls, or interestingly, in the form of persistent self-circling current of bulk defects (similar to that found in Supplementary Movie 9 for annulus at smaller v 0 ).
We finally remark that defect pinning or local trapping can be identified under two different conditions.(i) For small self-driving strengths with v 0 < v 0c , defects are pinned near the cusps as well as in the bulk (see the first column of Figs.2(a)-2(c)), with similar results obtained for different confining geometries and topologies including closed cavities and open and closed channels as well as different cusp shapes such as the polygonal confining geometry as found in our simulations (see Supplementary Note 1 for some sample results).(ii) In the intermediate range of v 0 with states of highly fluctuating defect dynamics, it is possible that some defects could be temporarily trapped close to the cusp as caused by strong boundary confinement particularly in narrow closed channels, as shown in an example of an n = 6 hypocycloid ring in Fig. 8 (see the white-circled corner) and Supplementary Movie 8.Although these look similar to that found for self-propelled rods with smectic defect structures trapped in a wedge by simulations and experiments [80,81], the underlying mechanisms are different.In those previous studies the trapping or escape of active rods and smectic defects was for an isolated cusp/wedge in an open environment of particles.While there is an active smectic structure pinned by the open cusp, the outside region is pretty dilute in active particles and as a consequence the cusp angle and shape play a crucial role [80,81].In contrast, here the active smectics fill the whole cavity and are overall confined.Therefore defect pinning or trapping either occurs as a result of too small activity v 0 to overcome the large pinning barrier in the full strongly segregated pattern, the occurrence of which is independent of cusp geometry, or at the intermediate value of v 0 , could temporarily appear due to local strong confinement of narrow channels that hinders the propagation of boundary-nucleated defects and any subsequent defect multiplication.Finally we also mention that the classical sources of defect generation and multiplication by crystal deformation such as the Frank-Read source [82] would not directly apply here since our effects are controlled by a combination of activity (with the self-propulsion of smectic layers) and confined boundary conditions.

Conclusions
We have examined the dynamics of topological defects in active smectic systems subjected to three types of boundary confinements, i.e., closed cavities, open and closed channels with various geometries.Our simulations based on active PFC modeling indicate a viable way to effectively vary or control the complex dynamics and collective time-variation properties of defects through both particle self-driving and the geometry and topology of strong confinement, with the underlying mechanisms intrinsically distinct from those of passive patterns and active nematic systems.These confined nonequilibrium smectic systems are featured by three distinct regimes of active and persistent defect dynamics, including defect pinning in a glassy state with ultraslow evolution, the fast self-rotating of spiral defects in cavities or local defect variations or persistent self-circling defect currents in closed channels, and interestingly, a dynamical state governed by far-from-equilibrium, nonrelaxational processes with large defect fluctuations in the transition between localized and traveling smectic patterns.For the latter, a key factor is the intermittent but perpetual creation of new defects as enabled by the confinement boundary and enhanced by cuspate boundary geometry.A transition from random to time-periodic process of defects creation and annihilation can be made possible through the control of boundary cusp singularity as the mechanism of confinement-induced defect generation.These predictions can be examined and achieved in experiments on e.g., dense self-propelled rods [83] which form an active smectic phase.Examples range from vibrated granular rods which can be exposed to circular [62], epicycloid-like flower-shaped [63] or annular [64] confinements, to active colloidal rods [65,66] in channels and cavities.

Confinement geometries
As described above in Eqs. ( 4) and (5), to implement the planar anchoring condition of boundary confinement we develop an approach by combining a formulation of surface/interface free energy for imposing the rigid boundary conditions (i.e., ψ = ψ b , n • ∇ψ = 0, and P = P b , with ψ b and P b the boundary values of density field ψ and orientation field P), with the control of confinement geometry via the spatially dependent interface energy amplitude V b (r) = V b0 φ(r), where V b0 is the anchoring strength.Here φ(r) is an auxiliary phase field function used in the diffuse domain method [67] and is approximated as where ∆ is the interface width of the boundary and r s (r) is the signed distance function from any location r to the domain boundary ∂Ω.r s (r) = −d(r) < 0 inside the domain Ω and r s (r) = d(r) > 0 outside Ω, with d(r) the distance function to ∂Ω.Various methods or algorithms for calculating signed distance functions have been available.For 2D domain boundaries studied here, we can directly obtain analytic forms of r s for some simple geometries (see below).In the cases of complex boundary geometries, a straightforward way of approximation, as used in our simulations, is to numerically compute the distances d i = |r − r bi | to the points r bi ∈ ∂Ω (i = 1, 2, ...) discretized on the boundary curve and find the shortest distance among them, i.e., which also determines the boundary point r b corresponding to each r.The local unit normal to the boundary is then The following types of confinement geometries have been examined in our simulations: (i) Circular cavity and annulus: For a 2D circular cavity of radius r 0 , with cavity center located at r c = (x c , y c ), we have |r| = r = [(x−x c ) 2 +(y−y c ) 2 ] 1/2 , n = r, r s = r − r 0 , and For an annulus with inner and outer radius of r in and r out respectively, n = r and (ii) S-shaped open channel : Assume the channel is aligned vertically, with its center at r c = (x c , y c ) and the average locations of right and left boundaries at x − x c = ±x 0 .The boundary curves are of the form where S 0 is the amplitude and 2π/q s is the periodicity of the S-shaped modulation.The boundary normal is given by n = (−1, dx/dy) The area of this open channel in a system of vertical length l y is equal to 2x 0 l y (when l y is set as an integer number of modulation periodicity), which is used in the calculation of defect number density.The corresponding phase field function is approximated by Note that although r s = |x − x c | − |x b | for this S-shaped channel, the above equation is still a good approximation for φ(r) when ∆ is small enough (i.e., for sharp boundary interface).
(iii) Epicycloid cavities with n cusps: The corresponding parametric equations are given by where the parameter θ (not the polar angle) ranges from 0 to 2π, and b = a/n for a n-cusped epicycloid.The area of the enclosed cavity is (n + 1)(n + 2)πb 2 for integer n.The related phase field function φ(r) is calculated via Eq.( 18) with r s (r) and the unit normal n identified numerically as described above.Some sample results obtained from our simulations are presented in Figs. 2  and 3 and in Supplementary Movies 1-4 for different cusp number n, including n = 3 (of shape similar to trefoil), 4 (similar to quatrefoil), 5 (ranunculoid), and 6.
(iv) Hypocycloid cavities and rings: For a hypocycloid cavity with n cusps, the parametric equations are written as where b = a/n.The cavity area is equal to (n − 1)(n − 2)πb 2 .At each position inside or outside of cavity, values of phase field function φ, r s , and n are computed numerically by following the above procedure.Some sample simulation results for n = 4 (astroid), 5, and 6 hypocycloid cavities are given in Figs. 2 and 3 and in Supplementary Movies 5-7.Similar setup can be used for hypocycloid rings, with inner and outer boundary curves each determined by the above parametric equations with two sets of a and b parameters (see some sample simulation snapshots given in Fig. 8 and Supplementary Movie 8).(v) Cycloid open channel : For a vertically aligned channel, the left boundary is a straight line located at x = −x 1 (so that r s = ±|x+x 1 | and n = (1, 0)), while the right boundary curve is of a cycloid or trochoid form described by It is a curtate cycloid if a > b, a prolate cycloid if a < b, and a cycloid when a = b which is used in our simulations.We choose −π ≤ θ ≤ π and set 2πa = l y to satisfy the periodic boundary condition along the y direction with open ends of the channel.The corresponding channel area is equal to π(2a 2 + b 2 ) + 2πa(x 0 + x 1 ).This channel configuration is implemented in our simulations through numerical calculations of φ, r s , and n, with examples given in Fig. 6.
In principle any other types of boundary geometries, as long as the corresponding analytic or numerical expressions of boundary curves are available, can be described via similar procedure and thus implemented in our modeling and simulations.This approach that we introduce here, based on Eqs. ( 4) and ( 5) and the above implicit representation of domain boundary, allows us to apply the pseudospectral method with periodic boundary conditions in the whole system to numerically solve the active PFC equations subjected to the confinement of various types of cavity or channel geometry.

Algorithm for defect detection
To identify the topological defects (dislocations, disclinations, and grain boundaries) in the simulated smectic pattern, we use an algorithm based on the combination of two methods given in Refs.[84,85], with some modifications and extension.The implementation steps are described below.
Given the local stripe orientation ns = ∇ψ/|∇ψ| = (cos θ s , sin θ s ) with θ s the local orientation angle of the smectic layer, we can calculate at each spatial location r = (x, y) To detect the locations of defect cores, at each grid point the local orientation gradient is calculated [84], i.e., A s = |∇θ s | 2 .If A s exceeds a threshold value A 0s (e.g., A 0s = 0.2/(∆x) 2 , with ∆x = π/4 the numerical grid spacing), the corresponding grid point is considered to be in a defect core region.To obtain the specific location of each individual defect core, first the individual cluster of sites for each defect core region is identified by using the Hoshen-Kopelman (Union-Find) algorithm with raster scan to connect neighboring grid points of each cluster tree with large enough local orientation variation (A s > A 0s ).The cluster's center of mass then gives the position r CM of the corresponding defect core, with r CM = j r j A s (r j ) where r j is the spatial coordinate of site j within the cluster.
To reduce the artifacts or ambiguities caused by the choice of threshold A 0s , if the size of a cluster is larger than a limit (e.g., 20 grid sites) this part is then re-clustered through the Union-Find algorithm to divide it into smaller sub-clusters by increasing its threshold value A 0s by a percentage (e.g., 1/8) of max(A s ) − min(A s ) of that cluster.In addition, if the distance between the centers of mass r CM of any two clusters is less than another threshold value (e.g., 5.5∆x), they will be merged if the merged/connected cluster size would not exceed an upper limit (e.g., 18 sites).This reclustering-merging process is conducted only once, and the corresponding defect core locations (i.e., cluster centers of mass) will be recalculated.
To identify the specific type of each individual defect, we follow the standard procedure of calculating the topological charge (winding number) of each defect core by performing a closed-path integral of θ s over a counterclockwise square loop around the position of defect core [84,85].The defect type (charge-0 dislocation vs ±1/2 disclination) is determined via the calculated value of topological charge.It is noted that all the above calculations are for the orientation of stripes (determined by the apolar density field ψ) and the corresponding topological charges, but not for the polar vector field P which would yield different topological charges via a similar procedure of calculation.A boundary defect is labeled if the location of its defect core is within a certain distance (e.g., 8 grid points) to the cavity or channel boundary.We can also identify the defect cores (clusters) belonging to a grain boundary (cluster chain), via the Union-Find algorithm again (but not merging them), if the distance between the centers of mass (r CM ) of any two clusters (defect cores) is less than or equal to a value (e.g., 25∆x) and if there are at least N GB (e.g., = 4) of such clusters (cores).
There would still be some ambiguities/uncertainties of defect identification, which are unavoidable for any detection algorithm particularly for the cases of close or crowded defect cores.We have checked the results by varying different parameters of the algorithm and comparing with some manual spot checks to identify the close-to-optimal or compromised choices of parameters, and to ensure the results are consistent statistically.

Fig. 1
Fig. 1 Topological defects in active smectics.Typical topological defects found in simulations of active smectic pattern as determined by the density field profile, including (a) +1/2 and (b) −1/2 disclinations, and (c) a single and a doublet of edge dislocations.Some local normal directions of smectic layers are indicated by arrows, and an integral of their orientational angle θ s over a counterclockwise closed-path loop around the defect core (white-circled) determines its topological charge via dθ s /2π = ±1/2, 0 for disclinations and dislocations respectively.(Note that angles θ s with a π difference are equivalent to each other.)

Fig. 2
Fig. 2 Three dynamical regimes of active smectics defect evolution in closed cavities.(a)-(c) Transitions between pinned, highly fluctuating (annihilation-nucleation), and self-rotating defect states with the increase of self-propulsion strength v 0 , in (a) circular, (b) epicycloid, and (c) hypocycloid cavities.Each regime is represented by sample snapshots of the spatial profiles of density field ψ obtained from simulations, with the ψ scale labelled by the color bars.In the mid panels the circled regions highlight the time-evolving process of boundary-induced defect generation.The bulk defects inside are labeled by white symbols, and boundary defects by black ones.Among them the square symbols indicate dislocations and the up or down triangles indicate ±1/2 disclinations respectively.(d) Sample time variation of defect number density in the fluctuation regime at v 0 = 0.31 (for system size of 512 × 512 grid points).(e) The corresponding normalized time correlation C n (τ ) of defect density, calculated over t = 10 5 − 10 6 and averaged over 80 simulations for each cavity.For a better illustration only the error-bar band for epicycloid of n = 6 cusps is shown, while those for other cases are of similar range.

Fig. 3
Fig. 3 Time correlation of defect density in closed cavities.Autocorrelation function C n (τ ) of defect number density for various epicycloid and hypocycloid closed cavities at v 0 = 0.31.Also shown are two sample snapshots of the ψ field spatial profile, with the scale of ψ indicated by the color bar and the same symbol labeling of defects as that in Fig. 2.

Fig. 4
Fig. 4 Power spectrum for time variation of defect density in closed cavities.Power spectrum S ω of the defect number density for circular, 6-cusp epicycloid, and 4-and 6-cusp hypocycloid cavities at v 0 = 0.31, each averaged over 80 simulations for t = 10 5 − 10 6 .

Fig. 5
Fig. 5 Spatial profiles of ψ and S fields and defect trajectories in hypocycloid cavity.(a) Sample simulation snapshot of the density field ψ profile for a 4-cusp hypocycloid (astroid) cavity at v 0 = 0.31, with the labeling of defects.(b) Trajectories of defects in the upper part of the cavity, starting from a single defect (circled on the right) nucleated at the boundary at time t 0 = 198650 up to the snapshot of (a) at t = 2 × 10 5 .Larger symbol size corresponds to later time, and samples for different time are also color-coded according to the color bar given on right.The hypocycloid boundary curve is also plotted.The corresponding time evolution of the pattern and defects can be found in Supplementary Movie 5. (c) Enlarged portion of the circled region in (a), with black arrows indicating the polarization vector field P and the color bar showing the scale of density field ψ.(d) The corresponding spatial profile of the divergence field S = ∇ • P, with P vector also shown as arrows and the scale of S field indicated by the color bar.Note a small phase shift of this S field profile with respect to the ψ profile in (c), which corresponds to a nonzero migration velocity.

Fig. 6
Fig. 6 Time evolution of defects in open channels.(a) Time evolution of defect number density in two types of open channels at v 0 = 0.31 (with periodic boundary condition along the vertical open ends).(b)-(e) Snapshots of ψ profiles at different times for (b) S-shaped and (c)-(e) cycloid channels.

Fig. 7
Fig. 7 Defect density and time correlation for channels and rings.(a) Sample time evolution of defect number density for annulus and single-cusp epicycloid rings of two different simulation system sizes at v 0 = 0.31.(b) Autocorrelation C n (τ ) corresponding to the sample simulation of cycloid channel given in Fig. 6 and for n = 1 epicycloid rings (averaged over 80 simulation runs at t = 10 5 -10 6 for each system size).(c) Sample snapshots of the annulus and one-cusp epicycloid ring simulated.

Fig. 8
Fig. 8 Defect density and time correlation for hypocycloid rings.Sample time evolution of defect number density and the results of autocorrelation C n (τ ) (each averaged over 80 runs for t = 10 5 -10 6 with system grid size 512 × 512) for various hypocycloid rings at v 0 = 0.31.Sample snapshots of the ψ field spatial profile are also shown, with the white-circled region corresponding to the local defect trapping near a cusp of n = 6 hypocycloid ring as seen in Supplementary Movie 8.