Defect interactions through periodic boundaries in two-dimensional p-atics

Abstract

Periodic boundary conditions are a common theoretical and computational tool used to emulate effectively infinite domains. However, two-dimensional periodic domains are topologically distinct from the infinite plane, eliciting the question: How do periodic boundaries affect systems with topological properties themselves? In this work, I derive an analytical expression for the orientation fields of two-dimensional p-atic liquid crystals, systems with p-fold rotational symmetry, with topological defects in a flat domain subject to periodic boundary conditions. I show that this orientation field leads to an anomalous interaction between defects that deviates from the usual Coulomb interaction, which is confirmed through continuum simulations of nematic liquid crystals (p = 2). The interaction is understood as being mediated by non-singular topological solitons in the director field which are stabilized by the periodic boundary conditions. The results show the importance of considering domain topology, not only geometry, when analyzing interactions between topological defects.

---

I. Introduction

Topological defects are topologically protected excitations that arise in systems with spontaneous symmetry breaking.^1,2 Examples of topological defects manifest in biological, material, high-energy, and cosmological physics and include monopoles, disclinations, dislocations, vortices, Skyrmions, Hopfions, and cosmic strings.^3–13 In biological physics, defects may execute biological functions such as cell apoptosis and morphogenesis;^12,14–16 in material physics, defects may be manipulated to alter a material's optical and mechanical properties, facilitate self-assembly, actuate surfaces, or perform computations;^17–23 while in cosmological physics, defects are predicted to interact with light and matter, but have yet to be directly observed.^7,24 Thus, understanding the dynamics of defects is of fundamental importance to many areas of physics.

One of the simplest systems that admits topological defects are the p-atic liquid crystals, which describe systems with p-fold rotational symmetry. In two dimensions, the most common p-atics are those with p = 1, 2, 6, corresponding to the continuous XY-model, nematic liquid crystals, and hexatic liquid crystals. While these three models have been heavily studied,² there have been recent experimental realizations of other p-atics,^25–27 as well as recent theoretical work involving the hydrodynamics, anyonic braiding, and interaction with curvature of arbitrary p-atic defects.^28–31

The dynamics of p-atics, both in equilibrium and non-equilibrium settings (and particularly when topological defects are involved), are often studied using computational models. A common choice for boundary conditions in computations is periodic boundary conditions (PBC), which are often employed to emulate an effectively infinite system. In the infinite plane, it is known that topological defects interact via Coulomb forces, akin to electrostatic charges, for systems with a single elastic constant.² The Coulomb interaction between defects has also been observed experimentally in passive nematic liquid crystals.³² However, since domains with PBC are topologically distinct from infinite systems, it is not obvious that the Coulomb interactions in the infinite plane should be the same as interactions in PBC. For example, in the context of fluid dynamics, it has been shown that nontrivial domain topology can lead to incomplete solutions of the stream function-formulated incompressible Navier–Stokes equations.^33,34 There have been several prior analytical and computational studies of topological defects on the torus, a domain topologically identical to PBC, which have found anomalous interactions between defects,^35–39 though this is typically attributed to curvature, and to my knowledge the role of topology (as opposed to curvature) in these interactions has not yet been clearly ascertained.

Here, I study p-atic topological defect interactions in flat systems with PBC. In Sections II I lay the theoretical foundation regarding defect configurations in two dimensional p-atics and review the relevant results regarding defects in the infinite plane. Then, in Sections III, I attempt to derive the p-atic configuration with defects subject to PBC. I show that if one assumes that defects behave like electrostatic particles (as they do in the infinite plane) the resulting configuration will not be periodic. Instead, to conform to the boundary conditions, the configuration must be modified, which alters the interaction between defects to be highly asymmetric. In Sections IV, I show how this asymmetric interaction may be interpreted as being mediated by non-singular topological solitons that are stabilized by the PBC and how these results impact the analysis of computational studies. Finally, in Sections V I conclude by positing how these asymmetric effects may play a role in physically realizable domains.

II. p-Atic topological defects in the plane

p-Atic liquid crystals are characterized by local orientational order with p-fold rotational symmetry. In two dimensions, the orientational configuration of a p-atic may be summarized with an angle field θ_p(x,y), which represents the angle of the orientation at each point in space and it is understood that θ_p is defined modulo 2π/p. Variations in θ_p are penalized through an elastic free energy:where K is the elastic constant. For p = 1 or p = 2 there may be two elastic constants, but I will take the one-constant approximation here for analytical tractability. Thus, energy minimizing configurations must satisfy ∇²θ_p = 0.

Topological defects in two-dimensional p-atics are points where the orientation is singular. They carry a winding number, or topological charge, defined as

where n is an integer and C is a closed loop. Configurations that contain defects will thus be energy stabilizing if they satisfy ∇²θ_p = 0 with eqn (2) as a constraint. For the rest of the paper I will assume, for simplicity, that only ±1/p charge defects populate the system. Higher charge defects are energetically unstable, unless stabilized by confinement or external fields.

It is often useful to describe the configuration using a complex-valued, analytic function Ω(z) = ϕ(z,z̄) − ipθ_p(z,z̄) where z = x + iy and ϕ is the harmonic conjugate of θ_p.^40–44 Then ∇²θ_p = 0 is automatically satisfied and topological defects are represented as poles of Ω. It will also be useful to identify ϕ(z,z̄) = Re[Ω(z)] as the electrostatic potential of an equivalent configuration in which the poles of Ω are electrostatic charges.

A single defect in the plane has the configuration

θ_p(z,z̄) = k Im[log(z − z₀)] + θ₀ (3)

where k is the defect charge, z₀ = x₀ + iy₀ is the location of the defect, and θ₀ is a global phase. For a system with many defects, the configuration is given by the sum of single defect configurations, assuming all relative defect orientations minimize the free energy.^42,45,46 To analyze the interaction between defects I will follow the work of Halperin and Mazenko,^47–50 later adapted for liquid crystals,^51,52 which gives the defect velocity as a function of derivatives of the order parameter. If defect motion is induced only by elastic relaxation (i.e. there are no flows or external fields) the velocity, v = v_x + iv_y, of the jth defect is⁵¹where ∂_z ≡ (1/2)(∂_x − i∂_y), [small theta, Greek, tilde]_p is the non-singular part of the configuration at z_j, k_j is the defect charge, K is the elastic constant, and γ is the rotational viscosity (for the remainder of the paper I will work in units so that K = γ = 1). If the multi-defect orientationis inserted into the expression, the result is a Coulomb interaction between defects:where the complex conjugate, v̄, is used to directly compare with the electrostatic interaction between point charges (the complex electric field is often defined by E = E_x − iE_y).⁵³ This well-established result can also be derived using other methods.^2,8,54 The result also follows conceptually from applying the Cauchy–Riemann equations to Ω(z), namely ip∂_zθ_p = −∂_zϕ. As discussed above, identifying ϕ as the electrostatic potential of the equivalent charge configuration reinterprets this as v̄ ∼ E. A natural question to ask is: can we treat p-atic defects as Coulomb charges in PBC?

III. Periodic boundary conditions

To handle PBC, the goal is to find a configuration θ_p that has the property θ_p(z + 2[small script l]ω₁ + 2mω₂) = θ_p(z)mod 2π/p where [small script l] and m are integers and ω₁ and ω₂ are complex numbers such that Im[ω₂/ω₁] > 0 which describe the shape and size of the periodic domain. Note that, in general, the domain is a parallelogram [see Fig. 1], but if ω₂/ω₁ is purely imaginary the domain is rectangular. There are two immediate consequences of the periodicity: first, given a configuration, we may modify it to have any integer number of 2π/p rotations, n₁, n₂, along the ω₁, ω₂ cycles. That is, given a periodic θ_p,where A = 4 Im[[small omega, Greek, macron]₁ω₂] is the area of the domain, will also be a periodic configuration. The second consequence of the PBC is that the sum of all defect charges must be zero by the Poincaré–Hopf theorem. Then, the defects populate the domain in pairs and a configuration may be constructed by taking a periodic tiling of the center cell, shown schematically in Fig. 1.

Fig. 1 Schematic of the periodic system. Positive defects are located at z⁺, negative defects at z⁻, and a test defect at z. ω₁ and ω₂ give the principal directions of the periodic lattice.

A. Failure of the electrostatic analogy

To find a periodic configuration, I will first attempt to make an analogy with electrostatics, since this analogy holds in the case of the infinite plane. To determine the configuration, I place a “test defect” at location z, which is meant to measure the force induced by the defect configuration already present in the domain, θ_p(z,z̄). Following the expectation that p-atic defects interact like Coulomb charges, I assume that the force on the test defect is equivalent to the electric field of the corresponding point charge configuration, that is, I will assume v̄(z) = E(z).⁵³ Additionally, a useful identity is

∂_zΩ = −2ip∂_zθ_p (8)

which follows from the Cauchy–Riemann equations. Using eqn (4) and (8), Ω(z) satisfieswhere Δz_j = z_j⁺ − z_j⁻ is the coordinate difference between the jth positive and negative charge pair, ζ(z) is the Weierstrass zeta function with quasi-periods 2ω₁ and 2ω₂,⁵⁵ and f(Δz_j) is a function that only depends on the charges' positions and the center cell geometry, and, hence, may be regarded as a constant for a given configuration. For a square domain with system length L it reduces to f(Δz_j) = (π/L²)Δz̄_j. For completeness, a derivation of the right hand side of eqn (9) may be found in Appendix B.

A solution to eqn (9) is given by

where σ(z) is the Weierstrass sigma function [defined by d log σ/dz = ζ(z)]⁵⁵ and C is an integration constant. Thus, the imaginary part of eqn (10) represents a candidate analytical solution for the angle of the p-atic orientation field. However, even though its derivative is periodic, eqn (10) is only quasi-periodic:where η_k = ζ(ω_k).⁵⁵

For electrostatic charges, quasi-periodicity of the potential is not an issue, as long as its derivative, the electric field, is periodic. For p-atics, however, the physically relevant object is the orientation field θ_p, therefore, it must be periodic. I show examples of θ_p = −(1/p)Im[Ω(z)] computed from eqn (10) for p = 1 (XY-model) in Fig. 2(a) and p = 2 (nematic liquid crystal) in Fig. 2(b), where it is clear that the periodic boundary conditions are not satisfied.

Fig. 2 Orientation fields compatible with defects interacting as Coulomb charges, eqn (10), for a system with ω₁ = 0.5, ω₂ = 0.5i, z⁺ = 0.2 + 0.1i, z⁻ = −0.2 − 0.1i and (a) p = 1 or (b) p = 2. In this case, periodic boundary conditions are not satisfied.

B. Periodic configuration and interaction between defects

It follows from the above analysis that, in general, topological defects in p-atics do not share the same interaction as Coulomb charges in periodic systems. What, then, is the interaction between defects? To answer this, we first need a periodic configuration; fortunately, the imaginary part of eqn (10) can be forced to be periodic by adding terms proportional to z. The angle of the orientation field, θ_p(z,z̄) = −(1/p)Im[Ω(z)], is thenwhere A is the domain area, the second to last term comes from the invariance under rotations by 2π/p (n₁ and n₂ are integers), and θ₀ is a global phase (all examples shown here will have θ₀ = 0). Note that the corresponding analytic potential Ω(z) will still only be quasi-periodic in general. Indeed, it is not possible for a function with arbitrary poles to be both doubly periodic and analytic.⁵⁵

The orientation field corresponding to eqn (12) is plotted in Fig. 3(a) and (b) for the cases shown previously in Fig. 2(a) and (b) (with n₁ = n₂ = 0) and it is clear that the boundary conditions are now satisfied. Additionally, example configurations for p = 3, 4, 5, 6 are shown in Fig. S1.⁵⁶ To predict the defects' velocities we may use eqn (4), which gives the velocity of the kth positive defect:

where the notation j ≠ k⁺ indicates a sum that excludes only the kth positive defect. The velocity of the negative defects may be found by inverting the sign and excluding the kth negative defect from the first sum. In deriving these equations I have relied on complex-valued functions as the notation is more compact; however, it is possible to derive the same results using only real-valued functions, which I show in Appendix C.

Fig. 3 Orientation fields corresponding to eqn (12) for the same system of defects as in Fig. 2 for (a) p = 1 or (b) p = 2.

Let us now analyze the predicted defect velocity, eqn (13), for the simple case of two defects in a square domain of length L = 1 and n₁ = n₂ = 0. In this case, the velocity simplifies to

v_defect⁺ = −2[ζ(Δz̄) + πΔz]. (14)

Compare this to the force (f = f_x + if_y) experienced by the positive Coulomb charge in the same setup:

f_Coulomb⁺ = −2[ζ(Δz̄) − πΔz]. (15)

The two expressions differ only by a sign, yet the sign difference drastically alters the interaction. In Fig. 4(a)–(d) I have plotted eqn (14) and (15) as complex color plots and as vector plots. The color plots, Fig. 4(a) and (b), clearly show the zeros and poles of the interaction. Poles occur when defects are on top of one another (they are attracted with infinite force) while zeros occur due to force balance through the periodic boundaries. For electrostatic charges, the zeros are located on a square half the size of the domain since the interaction is symmetric. If the separation between defects is larger than half the system size, the electrostatic charges will attract through the periodic boundary instead of towards the interior of the domain. On the other hand, the zeros of the p-atic defect interaction are no longer located on a square and occur very close to the periodic boundary. Additionally, p-atic defects will only attract through the periodic boundaries in a small locus of separations near the zeros of the interaction. The vector plot in Fig. 4(c) shows that p-atic defects will almost always attract toward the center of the domain.

Fig. 4 (a) and (b) Complex color plots of the velocity of the positive p-atic defect, eqn (14) and force on the positive Coulomb charge, eqn (15) for a configuration with two defects or charges in a square domain with ω₁ = 0.5 and ω₂ = 0.5i as a function of defect separation Δz = Δx + iΔy. The color represents the argument of the function while the brightness indicates the magnitude. The black line indicates the subset of the plot that is shown in (e). (c) and (d) Corresponding vector plots of eqn (14) and (15) where the arrows indicate the argument of the function while the color indicates the magnitude. (e) Velocity of the positive defect for the case Δy = 0 (green solid line). The purple dashed line shows the force on the positive Coulomb charge in the same configuration. The black dot indicates the initial separation for the simulation shown in (f). (f) Time snapshots of a simulation of two annihilating defects in a nematic liquid crystal (p = 2) with initial separation Δx = 0.81. The color is the scalar order parameter S and the white lines are the nematic director n̂.

In Fig. 4(e), I have plotted eqn (14) and (15) for the case where Δy = 0. In this one-dimensional cut of Fig. 4(a) and (b) it is easier to compare the two directly. If the defects are close to each other, either at the center of the domain (Δx = 0) or at the edge of the domain (Δx = 1), the defect velocity is similar to the Coulomb interaction. For two electrostatic charges, the force balance point is Δx = 1/2 which is required by the symmetry of the interaction. However, for p-atic defects, eqn (14) predicts that the force balance point is much closer to the boundary, Δx ≈ 0.82. I have tested this prediction by simulating the annihilation of two defects in a nematic liquid crystal (p = 2) with the initial configuration given by eqn (12) with Δx = 0.81. The simulation minimizes the Landau–de Gennes free energy⁵⁴ of the nematic order parameter and is subject to periodic boundary conditions (further simulation details are provided in Appendix D). Fig. 4(f) shows several time snapshots of the nematic scalar order parameter S and director n̂, indicating that the defects indeed annihilate in the center even though they are initially closer to each other through the periodic boundary. In the SI, I have provided additional examples of predicted velocities and corresponding simulations for two or more defects that highlight the qualitative difference between defect interactions and Coulomb interactions in PBC.⁵⁶

For systems with more than two defects, a systematic analysis of defect behaviors quickly becomes intractable. However, it can be seen from eqn (13) that if and n₁ = n₂ = 0 then the defect velocities will coincide with the forces on electrostatic charges in the same configuration. Thus, there is a special case in which defects are expected to interact the same as electrostatic charges: when the topological dipole moment vanishes. This cannot occur for two defects, but can lead to interesting states in which the force on defects is zero. An example of this is shown in the SI.⁵⁶

IV. Topological solitons

The comparison between eqn (14) and (15) suggests that one could understand the interaction asymmetry between defects as Coulomb charges with an effective uniform electric field proportional to the distance between defects. How can we understand the emergence of this effective field physically? Consider a periodic domain with just two p-atic defects and consider a curve C that runs from the bottom to the top of the domain such that the end points are at the same x coordinate. C is a closed curve in PBC, so we may measure the topological charge, eqn (2), by this curve. If we shift the curve past a defect the value of the measured topological charge must change by ±1/p by definition of a topological defect. Fig. 5(a) shows an example of this for a defect pair in the XY-model (p = 1). The curve between the defects measures a charge k = −1, while the curve on the other side of the defects measures a charge k = 0. Such a distortion, which I will call a topological soliton, is required for each defect pair in the system and are stabilized by the periodic boundary condition. Because topological solitons disrupt the p-atic configuration, they cost elastic energy and hence induce additional forces on the defects, thus producing an asymmetric interaction between defect pairs akin to a uniform electric field.

Fig. 5 (a) Defects in the XY-model subject to periodic boundary conditions. The curve between the defects measures a charge [eqn (2)] k = −1, indicating a topological excitation. If the curve is moved past either defect, it measures k = 0, as shown by the curve to the right of the positive defect. (b) Time snapshots of a nematic liquid crystal simulation with defects initially separated by Δx = 0.2 and n₂ = 1. Color is the scalar order parameter S while the white lines are the orientation field. (c) Time snapshots of a simulation similar to (b) but with initial separation Δx = 0.1.

I emphasize that these solitons are “topological” because they carry a distinct, nonzero topological charge along a closed loop. Of course, it is not topologically forbidden for a distortion that preserves the topological charge to occur. However, such a distortion would no longer be a free energy minimum, and hence the configuration would no longer be harmonic, which is the focus of this work. While non-topological distortions may be important for applications in non-equilibrium p-atics, the topological solitons described here are required when defects are present in PBC for both equilibrium and non-equilibrium conditions.

For every pair of ±1/p defects in the system, there must be a topological soliton between them. The sign of the topological soliton [as measured from the bottom (left) of the domain to the top (right)] depends on the relative position of the defect pair to the measuring curve. For example, the defect pair shown in Fig. 5(a) has the positive defect to the right of the purple curve and the negative defect to the left which results in the curve measuring k = −1. If the defect positions were switched so that the positive defect was on the left of the measuring curve, one would measure k = +1 instead. The green curve would still measure k = 0 since it does not straddle the defects on the interior of the domain. If there are multiple pairs of defects, the solitons induced by the pairs may either add or subtract to one another, either increasing or reducing the interaction strength between all other defects (a few examples of this are shown in the SI⁵⁶). Indeed, the cancellation of topological solitons in this manner leads to the coincidence of the force on p-atic defects and electrostatic force when the topological dipole moment goes to zero.

One might worry that this analysis depends on the location of the periodic boundaries since, technically, a periodic domain is not bounded (one should think of a torus or donut). If I change the location of the “boundaries” in Fig. 5(a) so that the green curve straddles the defects, I should still measure k = 0 since I have only changed my field of view, and not anything within the system. The resolution to this is that, once one has set the “boundaries” of the system, they should be consistently held, and then one has to then take into account an additional, key feature that I have not yet analyzed: the effect of “global” topological solitons.

A. Global topological solitons

Topological solitons may exist without defects. Indeed, this is the physical interpretation of the integers n₁ and n₂ in eqn (7) and (12): in the absence of topological defects, n₁ (n₂) is the number of windings the orientation field makes along the ω₁ (ω₂) direction, and hence counts the number of “global” topological solitons in that direction. The term “global” here refers to the fact that n₁ windings in the ω₁ direction occur for all coordinates along ω₂, and vice versa. When defects are present, these global topological solitons may add or subtract to the local topological solitons which appear between defect pairs, potentially changing the direction of the interaction. Fig. 5(b) shows a simulation of two defects initially separated by distance Δx = 0.2 with a global topological soliton added along the y-direction, n₂ = 1, which has the effect of removing the topological soliton on the inner domain, and adding a topological soliton on the outer domain. For example, compare the first panel of the simulation in Fig. 4(f) and the first panel in Fig. 5(b). The former has a local topological soliton in the middle of the domain between the defect, while the latter has moved the topological soliton to the outer boundary. As the simulation snapshots show, the defects now move to annihilate on the outer boundary, even though they are initially closer to each other inside the domain.

It is not possible for a uniform configuration to spontaneously shift to one with a global topological soliton. That is, it is not possible to spontaneously change the global (in the sense discussed above) topological charge. It is possible, however, to go from a uniform configuration to a topological solitonic configuration by nucleation and then subsequent annihilation of a defect pair. Fig. 5(c) shows an example of a pair of simulated defects annihilating to stabilize a global topological soliton. Here the simulation is initialized in the same way as in Fig. 5(b), except the initial separation is Δx = 0.1 so the defects are close enough to attract towards the center of the domain. In general, defects that are interacting via elastic relaxation will move towards the local topological solitons, unless they are within the locus of points near the zeros of the interaction shown in Fig. 4(a).

Transitioning from a uniform configuration to a global soliton costs energy, and so would not occur spontaneously in gradient flow simulations, though could occur in simulated quenches and defect coarsening or non-equilibrium contexts such as applied external fields or active liquid crystals.^57–59 In such systems topological defect unbinding and annihilation events may change the number of global topological solitons. To analyze the elastic interaction between defects in these contexts, one needs to keep track of the topological indices n₁ and n₂ as well as defect positions. To measure n₁ (n₂) for an arbitrary configuration one could measure the winding number on a loop that goes along ω₂ (ω₁) and add the winding of the local topological solitons that are required by the topological defects, as discussed above. For temporal data, if the indices n₁ and n₂ are known at a given time, the following rule may be used to keep track of n₁ and n₂ for future times: if a positive defect moves through the periodic boundary in the positive ω₂ direction, then n₁ → n₁ + 1 and if a positive defect moves through the periodic boundary in the positive ω₁ direction, then n₂ → n₂ − 1. The rule for the negative direction is the opposite, and rules for negative defects apply opposite to those for positive defects. Additionally, if one also wants to keep track of the configuration θ_p, then every time a defect passes a periodic boundary one should add π/p to θ₀. This will ensure that the orientation of defects remain continuous as they move across the periodic boundary. Considering the effects of global solitons, it is interesting to note that the defect behaviors not only depend on their current positions, but also on the history of the defect dynamics in the system.

V. Physically realizable domains

Since PBC are reserved to theoretical or computational studies, a natural question is whether these effects would be observable in experimentally relevant domains. Nematic textures with defects separating topological solitons have been observed in channel geometries with fixed boundary conditions when a strong flow is applied to the channel. In this setup, the topologically nontrivial “Dowser” state is effectively a two-dimensional topological soliton.^60,61 Topological defects then mark boundaries between topological solitons and topologically trivial configurations. The analysis of defect dynamics presented above may be useful in understanding defect dynamics in this system. While one should note that hydrodynamics play an important role here as well, it has been shown that aspects of the “Dowser” textures can be recast as coming from effective equilibrium field theories.^60–62

A physically realizable geometry that has explicit periodicity is the annulus. If the inner and outer boundary conditions are the same, the argument for the existence of solitons between defects still applies, and configurations from a square geometry with periodic boundaries on one edge may be mapped to the annulus. In this case, even equilibrium systems (i.e. no applied flow) should, at least qualitatively, permit the above analysis of their defect dynamics. In Fig. 6 I show simulation snapshots of a nematic liquid crystal in an annular geometry (radii r₁ = 1 and r₂ = 2) with two defects initialized with angular separation Δϕ = 2π/5 but such that the topological soliton between them is on the longer side of the domain. In this case the defects annihilate along the longer path even though they are initially close to each other in the annulus. Note that here the defect trajectories are not symmetric due to the influence of the boundary curvature which interacts differently with positive and negative defects. Even without a quantitative prediction for the velocity we can estimate the direction the defects will move based on the location of the topological solitons between them.

Fig. 6 Time snapshots of a nematic liquid crystal simulation in an annulus with radii r₁ = 1, r₂ = 2 and initial angular defect separation Δϕ = 2π/5 with the local topological soliton on the longer side of the annulus. Color is the scalar order parameter S while the white lines are the orientation field.

Another geometry with these properties that may be physically realized is a torus, since this can be mapped directly to a periodic domain. The concepts of local and global topological solitons would apply to the torus, so it may be possible to observe interactions similar to those shown in Fig. 6 for the annulus where defects annihilate over the longer part of the domain. Simulations on the torus are beyond the scope of this paper, and I note that curvature plays an important role in this system as well, since it is known that Gaussian curvature attracts defects of like sign, while extrinsic curvature induces a preferred orientation,^36,37,39 so it is unclear what the relative strength of the topological forces may be.

VI. Conclusion

The effect of boundary conditions on phases with broken rotational symmetry has been of great interest for decades. While fixed (Dirichlet) boundary conditions are often the subject of such studies, I have shown analytically that PBC bestow a different interaction between defects than in infinite domains. The interaction is due to the domain topology, requiring non-singular elastic excitations that must accompany defect pairs. I stress that the defect interactions are qualitatively different than those of Coulomb charges, and that this should be recognized when comparing computational data to experimental data. Further, additional topological indices (n₁ and n₂) which represent overall windings along the periodic lattice directions, must be taken into account when considering the elastic forces on defects, and these windings depend on the history of the defect dynamics.

A future challenge will be to adapt previous work regarding arbitrary defect orientations^42,45,46,63 to the case of periodic bounded systems. Here I have assumed that the relative orientations between defects have the lowest energy. It would also be interesting to test the predictions on toroidal geometries where curvature will also play a role in defect dynamics and to posit topological interactions for more complex domains such as the n-genus surface. Additionally, generalizing these results to three dimensions will be important for analyzing the interactions between defects in three-dimensional ferromagnets or nematics, which no longer share the same topology and types of defects. There has been recent interest in topological soliton-mediated interactions between topological defects in 3D nematic liquid crystals,^64,65 and it would be interesting to investigate how such interactions manifest in topologically nontrivial domains since the 2D topological solitons discussed here are the conceptual equivalent of Skyrmions or merons in 3D.⁶⁶ The methods and results of this work will prove useful in these future endeavors.

Conflicts of interest

There are no conflicts to declare.

Data availability

MATLAB scripts for computing the periodic director structure given defect positions are available at https://doi.org/10.5281/zenodo.15690961. Simulation data and scripts are available from the author upon request.

Supplementary information is available. The supplementary information includes figures of defect configurations for p = 3--6, an analysis of a few examples of two or four defect interactions in periodic boundary conditions, and movies of the simulations. See DOI: https://doi.org/10.1039/d5sm00773a.

Appendices

A Halperin–Mazenko formalism for p-atic defect velocities

For completeness, I derive here the equation for the velocity of defects in a p-atic liquid crystal using the formalism of Halperin and Mazenko.^47–50 This will mirror the derivation for active nematic liquid crystals in ref. 51, except that in this case the order parameter has p-fold rotational symmetry and the only time dependence I will consider comes from the elastic interaction.

The p-atic order parameter may be written as a complex-valued function

ψ_p(z,z̄) = S_p(z,z̄)e^{ipθ_p(z,z̄)} (A1)

where S_p is the local degree of orientational order and θ_p is the local orientational angle. As a side note, the reader should not confuse the order parameter with the local director n_p(z,z̄) = e^iθ_p which gives the local orientation of the p-atic phase and is implicitly understood to be p-fold degenerate under rotations by 2π/p. ψ_p is defined so that it is explicitly invariant under θ_p → θ_p + 2π/p.

The magnitude of the order parameter S_p is constant and S_p ≈ 1 when changes in θ_p are small, but S_p → 0 if θ_p becomes discontinuous to preserve the continuity of the order parameter. This is precisely the case for the core of topological defects, and the Halperin–Mazenko formalism takes advantage of this fact. The topological defect density may be written in two ways:

where j indexes individual defects and k_j is the charge of the jth defect. The first equality is the usual definition for defect charge density while the second equality comes from a coordinate change from the domain to order parameter space. The Jacobian determinant D_p is given by

D_p(z,z̄) = ∂_zψ_p∂_z̄[small psi, Greek, macron]_p − ∂_z[small psi, Greek, macron]_p∂_z̄ψ_p (A3)

where ∂_z = (1/2)(∂/∂x − i∂/∂y) and . To get the defect velocity, one takes advantage of two conservation laws:

∂_tρ_p + ∇·J^(ρ) = 0 (A4)

∂_tD_p + ∇·J^(ψ)_p = 0 (A5)

where the J are respective currents. The first comes from topological charge conservation, the second can be derived from taking a time derivative of eqn (A3). It is given by

J^(ψ)_p = −∂_tψ_p∂_z̄[small psi, Greek, macron]_p + ∂_t[small psi, Greek, macron]_p∂_z̄ψ_p. (A6)

Using the conservation laws and eqn (A2), one may write

where the first equality can be derived by taking a time derivative of eqn (A2) and defining v_j = ∂z_j/∂t and the second equality combines the conservation of D_p with the definition of the defect density. Comparing the equations, and noting that the δ-functions require this only hold at z = z_j, yieldswhere ∂_t = ∂/∂t.

Eqn (A8) gives a kinematic equation for the velocity of a defect as a function of the local order parameter. To get a prediction for the defect velocity one must approximate both the space and time derivatives of ψ_p. To handle the time derivative, I will assume a Landau–Ginzburg model for the free energy of the p-atic phase:

where A and C are phenomenological parameters (set so that the ordered phase is stable) and K is the elastic constant. The time dependence of the order parameter is further assumed to be given by free energy relaxationwhere γ is a rotational viscosity.

Eqn (A10) can be substituted into eqn (A8) to yield an expression that depends only on spatial derivatives of ψ_p (terms proportional to ψ_p go to zero at the defect core). To handle the spatial derivatives I use the fact that near the defect core S_p ∼ |z − z_j| so that

where the + (−) is used for a positive (negative) defect and [small theta, Greek, tilde]_p is the non-singular orientation field (i.e. the orientation field due to all other defects). Substituting this into eqn (A8) yields eqn (4) above:

It is interesting to note that since the charges are ±1/p the dependence of p drops out of the velocity, predicting that defects elastically interact at the same strength independent of p. This is initially surprising, since the energy of a defected configuration is smaller for larger p (the field rotates by less, hence the overall gradient is smaller), so the rotation rate of the relaxing p-atic phase should be smaller. However, for increasing p, the required rotation to move a defect a given distance is smaller, exactly accounting for the slower rotation rate.

B The Coulomb interaction and Weierstrass functions

Eqn (9) of the main text asserts that the velocity of defects that interact as Coulomb charges in a periodic domain will be given by Weierstrass ζ functions and an additional constant that depends on the domain shape and size and the locations of defects. I show this here for completeness.

The electric field of a point charge in two-dimensions may be written in complex coordinates as

where z₀ is the location of the point charge and I am using units so that e/2πε₀ = 1. In the periodic configuration described in the main text, with half periods ω₁ and ω₂, the electric field will be the superposition of all fields created by the periodic copies:

Compare this with the definition of the Weierstrass zeta function:⁵⁵

where Ω_{[small script l],m} = 2ω₁[small script l] + 2ω₂m. Thus, evidently, eqn (B2) may be written as

There is, however, an issue with eqn (B4): it is only quasi-periodic. The quasi-periodicity is a result of the distribution of fields for a periodic array of charges.⁵³ A periodic function can be constructed by adding terms proportional to z − z₀ and z̄ − z̄₀.⁵³ The result is

where A is the area of the domain. One may worry that the field in eqn (B5) is no longer analytic, and hence does not obey Maxwell's equations for a single point charge. The resolution is that the non-analytic part (proportional to z̄) actually adds a constant uniform density electric charge that balances the point charge. This is consistent with the requirement of zero net electric flux through the periodic boundaries, i.e. the total charge on a torus must be zero.

Finally, to yield eqn (9) in the main text, I note that the electric field produced by two charges of opposite sign is

E = ζ(z − z⁺) − ζ(z − z⁻) + f(Δz) (B6)

which can be shown from eqn (B5). Since the charge is balanced in this configuration, the resulting electric field is analytic.

C Real-valued orientation configuration and defect velocity

The arguments in the main text using complex-valued functions can be made using only real-valued functions. I derive the resulting expressions here. Eqn (4) can be equivalently derived for real-valued order parameters (see e.g. ref. 52). For p-atics it takes the formwhere (x_j,y_j) is the location of the jth defect.

As in the main text, I will first assume that the interaction is that of Coulomb charges, so that the orientation field should satisfy

where i indexes the defect pair and [Doublestruck Z] denotes the set of integers. Here I have restricted the domain to be rectangular for simplicity; the following can be easily generalized to arbitrary parallelogram domains.⁶⁷Eqn (C2) and (C3) can be simplified by summing on one of the indices first using a method from the computational electrostatics literature known as Lekner summation.^67–69 Summing on [small script l] and then integrating the result giveswhere Δy_i = y_i⁺ − y_i⁻.

Eqn (C4) is equivalent to eqn (10) in the main text and still has the problem that it is only quasi-periodic. Indeed, it satisfies

As with the complex-valued version, we only have to slightly modify eqn (C4) to force periodicity: . The full expression for the orientation angle is

where θ₀ is a global phase and the second to last term comes from the invariance under rotations by 2π/p (n₁ and n₂ are integers). Eqn (C6) is the real-valued equivalent of eqn (12) in the main text (when ω₂/ω₁ = iL_y/L_x). I note that there is another equivalent expression if m is summed over in eqn (C2) and (C3):

The velocity for a positive defect can be obtained by substituting either eqn (C6) or eqn (C7) into eqn (C1). Using eqn (C6),

where v_x,j⁺ denotes the x component of the velocity of the ith positive defect and the notation i ≠ j⁺ indicates summation on all defect positions that are not the ith positive defect. Eqn (C8) and (C9) are the real-valued equivalents of eqn (8) in the main text.

Finally, even though eqn (C6)–(C9) include infinite sums, they are actually quickly converging sums due to the hyperbolic trigonometric functions. The infinite sums appearing in the equations will converge to machine precision [error ∼ O(10⁻¹⁶)] after summing up to |m|, |[small script l]| = 5.⁶⁹ I have provided a simple MATLAB script which computes the orientation angle θ_p(x,y) and the predicted velocities of defects given the positions of defects to facilitate future studies.⁷⁰

D Computational details

Here I give details regarding the simulations described in the main text and SI. The simulations solve the time evolution equations for the nematic tensor order parameter Q = S[n̂ ⊗ n̂ − (1/2)I] where S is the local degree of order and n̂ is the director. Q may be mapped to the complex order parameter by , hence the procedure described here may also be regarded as a simulation of the complex-valued order parameter.

The free energy of a nematic configuration is given by the Landau–de Gennes free energy:⁵⁴

and for all simulations I set K = 10⁻³, C = 4, and A = −1 which sets the value of S in the ordered phase to . To minimize the free energy for a given starting configuration, the time dependence of Q iswhere γ is a rotational viscosity that is set to γ = 1 for all simulations.

Eqn (D2) is discretized in space on a square mesh of length L = 1 with 20 004–29 041 vertices or an annular mesh with radii r₁ = 1 and r₂ = 2 with 20 004 vertices using the MATLAB/C++ finite element package FELICITY.⁷¹ The time dependence is solved using a backwards Euler method with time step Δt = 0.5 and a Newton–Raphson sub-routine to solve the resulting nonlinear equation. In the square domain, periodic boundary conditions are used while in the annular domain Dirichlet boundary conditions such that S = S* and θ_p = ϕ (where ϕ is the azimuthal angle) are used on the inner and outer boundary. Defect configurations are initialized using the periodic configuration eqn (C6) with p = 2.

Acknowledgements

I am indebted to Oleg Tchernyshyov for pointing out the connection between the orientation field and the electrostatic potential, as well as for the use of the Weierstrass functions to describe the electric field of periodic point charges. This work also benefited greatly from conversations with Brandon Klein and Daniel Beller. I gratefully acknowledge support from the Department of Physics and Astronomy at Johns Hopkins University through the William H. Miller III postdoctoral fellowship.

References

1. N. D. Mermin, The topological theory of defects in ordered media, Rev. Mod. Phys., 1979, 51, 591.
2. P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, 1995.
3. A. Abrikosov, The magnetic properties of superconducting alloys, J. Phys. Chem. Solids, 1957, 2, 199.
4. B. I. Halperin and D. R. Nelson, Theory of twodimensional melting, Phys. Rev. Lett., 1978, 41, 121.
5. Formation and Interactions of Topological Defects, ed. A.-C. Davis and R. Brandenberger, Springer, New York, NY, 1994.
6. G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin and V. M. Vinokur, Vortices in high-temperature superconductors, Rev. Mod. Phys., 1994, 66, 1125.
7. T. W. B. Kibble, Phase transitions and topological defects in the early universe, Aust. J. Phys., 1997, 50, 697.
8. L. M. Pismen, Vortices in Nonlinear Fields, Oxford University Press, 1999.
9. M. Kleman and J. Friedel, Disclinations, dislocations, and continuous defects: A reappraisal, Rev. Mod. Phys., 2008, 80, 61.
10. P. J. Ackerman and I. I. Smalyukh, Diversity of knot solitons in liquid crystals manifested by linking of preimages in torons and hopfions, Phys. Rev. X, 2017, 7, 011006.
11. G. Duclos, R. Adkins, D. Banerjee, M. S. E. Peterson, M. Varghese, I. Kolvin, A. Baskaran, R. A. Pelcovits, T. R. Powers, A. Baskaran, F. Toschi, M. F. Hagan, S. J. Streichan, V. Vitelli, D. A. Beller and Z. Dogic, Topological structure and dynamics of three-dimensional active nematics, Science, 2020, 367, 112.
12. Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun and K. Keren, Topological defects in the nematic order of actin fibres as organization centers of Hydra morphogenesis, Nat. Phys., 2021, 17, 251.
13. J.-S. B. Tai, A. J. Hess, J.-S. Wu and I. I. Smalyukh, Field-controlled dynamics of skyrmions and monopoles, Sci. Adv., 2024, 10, eadj9373.
14. A. Livshits, L. Shani-Zerbib, Y. Maroudas-Sacks, E. Braun and K. Keren, Structural inheritance of the actin cytoskeletal organization determines the body axis in regenerating hydra, Cell Rep., 2017, 18, 1410.
15. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Topological defects in epithelia govern cell death and extrusion, Nature, 2017, 544, 212.
16. L. A. Hoffmann, L. N. Carenza, J. Eckert and L. Giomi, Theory of defect-mediated morphogenesis, Sci. Adv., 2022, 8, eabk2712.
17. M. Peach and J. S. Koehler, The forces exerted on dislocations and the stress fields produced by them, Phys. Rev., 1950, 80, 436.
18. B. Escaig, Dissociation and mechanical properties. dislocation splitting and the plastic glide process in crystals, J. Phys. Colloq., 1974, 35, C7.
19. M. Ravnik, M. Škarabot, S. Žumer, U. Tkalec, I. Poberaj, D. Babič, N. Osterman and I. Muševič, Entangled nematic colloidal dimers and wires, Phys. Rev. Lett., 2007, 99, 247801.
20. X. Wang, D. S. Miller, E. Bukusoglu, J. J. de Pablo and N. L. Abbott, Topological defects in liquid crystals as templates for molecular self-assembly, Nat. Mater., 2016, 15, 106.
21. G. Babakhanova, T. Turiv, Y. Guo, M. Hendrikx, Q.-H. Wei, A. P. Schenning, D. J. Broer and O. D. Lavrentovich, Liquid crystal elastomer coatings with programmed response of surface profile, Nat. Commun., 2018, 9, 456.
22. Y. Guo, M. Jiang, S. Afghah, C. Peng, R. L. B. Selinger, O. D. Lavrentovich and Q.-H. Wei, Photopatterned designer disclination networks in nematic liquid crystals, Adv. Opt. Mater., 2021, 9, 2100181.
23. R. Zhang, A. Mozaffari and J. J. de Pablo, Logic operations with active topological defects, Sci. Adv., 2022, 8, eabg9060.
24. G. Servant and P. Simakachorn, Ultrahigh frequency primordial gravitational waves beyond the khz: The case of cosmic strings, Phys. Rev. D, 2024, 109, 103538.
25. K. Zhao, R. Bruinsma and T. G. Mason, Local chiral symmetry breaking in triatic liquid crystals, Nat. Commun., 2012, 3, 801.
26. P.-Y. Wang and T. G. Mason, A brownian quasi-crystal of pre-assembled colloidal penrose tiles, Nature, 2018, 561, 94.
27. J.-M. Meijer, V. Meester, F. Hagemans, H. N. W. Lekkerkerker, A. P. Philipse and A. V. Petukhov, Convectively assembled monolayers of colloidal cubes: Evidence of optimal packings, Langmuir, 2019, 35, 4946.
28. L. Giomi, J. Toner and N. Sarkar, Long-ranged order and flow alignment in sheared p-atic liquid crystals, Phys. Rev. Lett., 2022, 129, 067801.
29. L. Giomi, J. Toner and N. Sarkar, Hydrodynamic theory of p-atic liquid crystals, Phys. Rev. E, 2022, 106, 024701.
30. A. Mietke and J. Dunkel, Anyonic defect braiding and spontaneous chiral symmetry breaking in dihedral liquid crystals, Phys. Rev. X, 2022, 12, 011027.
31. F. Vafa, G. H. Zhang and D. R. Nelson, Defect absorption and emission for p-atic liquid crystals on cones, Phys. Rev. E, 2022, 106, 024704.
32. Y. Zushi and K. A. Takeuchi, Scaling and spontaneous symmetry restoring of topological defect dynamics in liquid crystal, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2207349119.
33. I. Nitschke, S. Reuther and A. Voigt, Vorticity-stream function approaches are inappropriate to solve the surface Navier–Stokes equation on a torus, PAMM, 2021, 20, e202000006.
34. H. Yin, M. S. Nabizadeh, B. Wu, S. Wang and A. Chern, Fluid cohomology, ACM Trans. Graphics, 2023, 42, 1–25.
35. M. Bowick, D. R. Nelson and A. Travesset, Curvatureinduced defect unbinding in toroidal geometries, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 041102.
36. R. L. B. Selinger, A. Konya, A. Travesset and J. V. Selinger, Monte carlo studies of the XY model on Two-Dimensional curved surfaces, J. Phys. Chem. B, 2011, 115, 13989.
37. P. W. Ellis, D. J. G. Pearce, Y.-W. Chang, G. Goldsztein, L. Giomi and A. Fernandez-Nieves, Curvature-induced defect unbinding and dynamics in active nematic toroids, Nat. Phys., 2018, 14, 85.
38. D. J. G. Pearce, P. W. Ellis, A. Fernandez-Nieves and L. Giomi, Geometrical control of active turbulence in curved topographies, Phys. Rev. Lett., 2019, 122, 168002.
39. J. Rojo-González, L. N. Carenza, A. de la Cotte, L. A. Hoffmann, L. Giomi and A. Fernandez-Nieves, Defectpopulated configurations in nematic solid tori and cylinders, Phys. Rev. Res., 2024, 6, L012065.
40. T. G. Chandler and S. E. Spagnolie, A nematic liquid crystal with an immersed body: equilibrium, stress and paradox, J. Fluid Mech., 2023, 967, A19.
41. T. G. J. Chandler and S. E. Spagnolie, Exact and approximate solutions for elastic interactions in a nematic liquid crystal, SIAM J. Appl. Math., 2024, 84, 2476.
42. S. Čopar and Ž. Kos, Many-defect solutions in planar nematics: interactions, spiral textures and boundary conditions, Soft Matter, 2024, 20, 6894.
43. H. Miyoshi, H. Miyazako and T. Nara, Free energy formulae for confined nematic liquid crystals based on analogies with Kirchhoff–routh theory in vortex dynamics, Proc. R. Soc. A, 2024, 480, 20240405.
44. H. Miyoshi, D. G. Crowdy, H. Miyazako and T. Nara, Equivalent local force conditions minimizing the frank free energy for topological defect equilibria in nematic liquid crystals, Proc. R. Soc. A, 2025, 481, 20240569.
45. A. J. Vromans and L. Giomi, Orientational properties of nematic disclinations, Soft Matter, 2016, 12, 6490.
46. X. Tang and J. V. Selinger, Orientation of topological defects in 2d nematic liquid crystals, Soft Matter, 2017, 13, 5481.
47. B. I. Halperin, in Physics of Defects, ed. R. Balian, North-Holland Pub. Co., 1981.
48. F. Liu and G. F. Mazenko, Defect-defect correlation in the dynamics of first-order phase transitions, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 5963.
49. G. F. Mazenko and R. A. Wickham, Ordering kinetics of defect structures, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 57, 2539.
50. G. F. Mazenko, Velocity distribution for strings in phaseordering kinetics, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 1574.
51. L. Angheluta, Z. Chen, M. C. Marchetti and M. J. Bowick, The role of fluid flow in the dynamics of active nematic defects, New J. Phys., 2021, 23, 033009.
52. C. D. Schimming and J. Viñals, Kinematics and dynamics of disclination lines in three-dimensional nematics, Proc. R. Soc. A, 2023, 479, 20230042.
53. D. Hill, V. Slastikov and O. Tchernyshyov, Chiral magnetism: a geometric perspective, SciPost Phys., 2021, 10, 078.
54. P. G. de Gennes, The Physics of Liquid Crystals, Oxford University Press, 1975.
55. E. T. Whittaker and G. N. Watson, A Course of Modern Analysis, Cambridge Mathematical Library, Cambridge University Press, 4th edn, 1996.
56. See SI for further examples of the interaction between defects and movies of the simulations.
57. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Defect annihilation and proliferation in active nematics, Phys. Rev. Lett., 2013, 110, 228101.
58. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Active nematics, Nat. Commun., 2018, 9, 3246.
59. J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas and L. Giomi, Epithelia are multiscale active liquid crystals, Nat. Phys., 2023, 19, 1773.
60. T. Emeršič, R. Zhang, Ž. Kos, S. Čopar, N. Osterman, J. J. de Pablo and U. Tkalec, Sculpting stable structures in pure liquids, Sci. Adv., 2019, 5, eaav4283.
61. S. Čopar, Ž. Kos, T. Emeršič and U. Tkalec, Microfluidic control over topological states in channel-confined nematic flows, Nat. Commun., 2020, 11, 59.
62. X. Tang and J. V. Selinger, Minimization principle for shear alignment of liquid crystals, Phys. Rev. E, 2020, 101, 032701.
63. D. J. G. Pearce and K. Kruse, Properties of twisted topological defects in 2d nematic liquid crystals, Soft Matter, 2021, 17, 7408.
64. J. Pollard and R. G. Morris, Morse theory and meronmediated interactions between disclination lines in nematic materials, Phys. Rev. X, 2025, 15, 021099.
65. J. Pollard and R. G. Morris, Defect dynamics in cholesterics: beyond the peach–koehler force, Soft Matter, 2025, 21, 2300.
66. T. Machon and G. P. Alexander, Global defect topology in nematic liquid crystals, Proc. R. Soc. A, 2016, 472, 20160265.
67. N. Grønbech-Jensen, Summation of logarithmic interactions in nonrectangular periodic media, Comput. Phys. Commun., 1999, 119, 115.
68. J. Lekner, Summation of coulomb fields in computersimulated disordered systems, Phys. A, 1991, 176, 485.
69. N. Grønbech-Jensen, Summation of logarithmic interactions in periodic media, Int. J. Mod. Phys. C, 1996, 07, 873.
70. C. Schimming, schim111/periodicpaticconfig: Periodicpaticconfig, 2025.
71. S. W. Walker, Felicity: A matlab/c++ toolbox for developing finite element methods and simulation modeling, SIAM J. Sci. Comput., 2018, 40, C234.

---