Collective steady-state patterns of swarmalators with finite-cutoff interaction distance

For finite $r_{\rm c}>2R_{\rm s}$, the system with positive coupling strength $K$ is found to be in SS. It is interesting to note that the steady state with multiple discs in SS is discovered (see Fig. 3). Multiple discs of SS have not been found in the system with $r_{\rm c}=\infty$. Figure 3 is the numerical results of $N=800$ system for $r_{\rm c}=1.5$ when $K=0.1$ and $J=1$ (we below fix $J=1$), displaying the process how the multiple discs of SS are formed.

The radius of the discs is, in common, $\approx 1/\sqrt{2}$ that reminds us of $R_{\rm s}$ in Eq. (8) with $J=1$. A simple explanation of this numerical observation is following; the swarmalators are divided into three groups apart from others farther than $r_{\rm c}$ during transient period, and each group evolves into its own steady state. The snapshots are shown in Fig. 3, from initial $T=0$ via transient $T=50,~{}100$ to final $T=1000$.

Once a number of swarms is included in a circular region of radius $r_{\rm c}$ and this region is separated from the others farther than $r_{\rm c}$, the dynamics taking place thereafter in the region is identical to that by $r_{\rm c}=\infty$. If there are several regions separated that way (the number of regions depends on initial condition), the same number of patterns finally appears. We note, when $r_{\rm c}>2R_{\rm s}$, the self-consistent equation for density [Eq. (5)] does not change from that for $r_{\rm c}=\infty$. Then, each final pattern becomes no more than the SS of Eqs. (7) and (8). As expected, the discs in the last snapshot are, in common, of radius $1/\sqrt{2}$ with their own uniform densities. The minimal separation between the discs is $2R_{\rm s}$ as a consequence of the condition $r_{\rm c}>2R_{\rm s}$.

With no loss of generality, the condition $r_{\rm c}>2R_{\rm s}$ considered for SS case is generalized to

where $D_{\infty}$ is the diameter of a pattern appearing in the $r_{\rm c}\to\infty$ limit. It can be SA or SPW beside SS depending on $K$. That is, any of them is expected in the same reason when Eq. (20) holds and, when repeated, the minimal separation is given by the diameter.

In case of static async, $D_{\infty}=2R_{\rm a}$. Thus by Eq. (20), SA should appear for finite $r_{\rm c}>2R_{\rm a}$ when the negatively strong coupling $(K<K_{c}<0)$ is provided. As expected, the system is found to be in SA where the phases of all the swarmalators are random in the full range of $[-\pi,\pi)$. Also, similar to the previous sync case, the steady state composed of multiple discs of SA is discovered, which has not been found in the system with $r_{\rm{c}}=\infty$. Figure 4 shows multiple SA discs for $r_{\rm c}=2.1$ (note $D_{\infty}=2R_{\rm a}=2$ for SA).

Next, we find there appear anomalous steady-state patterns for $r_{\rm c}<2R_{\rm s}$ when $K>0$, as shown in Fig. 5(a) and (b).

Because Eq. (20) does not hold this time for $D_{\infty}=2R_{\rm s}$, what a different pattern may form was in question, and the shown anomalous disc is the result to be understood. The radius becomes larger to $R_{\rm s}+\Delta$. The density looks uniform inside the region of radius $R_{\rm s}^{\prime}(<R_{\rm s})$ while decreases thereafter. These observations can be explained with Eq. (5) as follows.

For finite $r_{\rm c}$, it is given from Eq. (6) that

where $H(x)$ is the Heaviside step-function assigned with $1$ for $x\geq 0$, or $0$ otherwise. Plugging Eq. (21) into Eq. (5) with the numerical observation $\rho({\bf r})=\rho(r)$, one may consider the geometric situations by $R_{\rm s}+\Delta$, $R_{\rm s}^{\prime}$, and $r_{\rm c}$. Figure 5(d) shows a situation when the interaction range of radius $r_{\rm c}$ marginally covers the whole distribution of radius $R_{\rm s}+\Delta$; see the blue arc and the dashed circle tangentially meet. Then, since $\nabla_{\bf r}^{2}U_{\rm att}({\bf r},{\bf r}^{\prime})$ is constant in the interaction range, Eq. (5) gives $\rho(r=R_{\rm s}^{\prime})=(2\pi)^{-1}2(1+J)\int d{\bf r}^{\prime}H(r_{\rm c}-|{\bf r}-{\bf r}^{\prime}|)\rho(r^{\prime})=((1+J)/\pi)\int d{\bf r}^{\prime}\rho(r^{\prime})=(1+J)/\pi$. This value does not change as long as the blue circle covers the dashed one, and this is the case when $r\leq R_{\rm s}^{\prime}$. Note $(1+J)/\pi$ is the very $(\pi R_{\rm s}^{2})^{-1}$ of Eq. (7).

When $r>R_{\rm s}^{\prime}$, the blue circle cannot include the dashed one, and the covered region decreases as $r$ increases. So that, in Eq. (5), the integration range becomes smaller for $r$, and this leads to decreasing $\rho(r)$ in the left hand side. Then, one may write

where $\rho_{0}\equiv(\pi R_{\rm s}^{2})^{-1}$ and $g(r)$ is a proper function that decreases for $r$ from $g(R_{\rm s}^{\prime})=\rho_{0}$. Figure 5(e) shows the schematic of Eq. (22). The density begins to decrease from $\rho_{0}$ following $O(\epsilon^{3/2})$ for small $\epsilon=r-R_{\rm s}^{\prime}>0$ as the area $\sim\epsilon^{3/2}$ is excluded from the integration in Eq. (5). A lower bound of $g(R_{\rm s}+\Delta)$ is simply arguable in the $R_{\rm s}^{\prime}=0^{+}$ case, which is accompanied with $r_{\rm c}=R_{\rm s}+\Delta$. In this case, when the blue circle is moved to locate its center at $r=R_{\rm s}+\Delta$, Eq. (5) gives $g(R_{\rm s}+\Delta)>\rho_{0}/4$ as the integration range covers, at least, a quarter sector of the disc in the dashed circle. Note the coverage is minimized when $R_{\rm s}^{\prime}=0^{+}$, so that $\rho_{0}/4$ can be a lower bound of $g(R_{\rm s}+\Delta)$. As the truncated $U_{\rm att}({\bf r},{\bf r}^{\prime})$ by $r_{\rm c}$ still shows $U_{\rm att}({\bf r},{\bf r}^{\prime})=U_{\rm att}(|{\bf r}-{\bf r}^{\prime}|)$ and the density in Eqs. (21) and (22) also does $\rho({\bf r})=\rho(r)$, $|V(0)|<\infty$ follows. Then, a consequent constant $V({\bf r})$ gives ${\bf v}({\bf r})=-\nabla_{\bf r}V({\bf r})=0$. Hence, the density in Eq. (22) can be a steady-state configuration. We call this anomalous static sync (aSS).

The aSS explained so far is conditioned in $r_{\rm c}<2R_{\rm s}$. Meanwhile, the numerical data manifest there should be a lower bound on $r_{\rm c}$. For a survey on lower bound, we revisit Fig. 5(d). It states

with $R_{\rm s}^{\prime}=R_{\rm s}^{\prime}(r_{\rm c})$ and $\Delta=\Delta(r_{\rm c})$. As $R_{\rm s}$ is independent of $r_{\rm c}$, the lower bound of $r_{\rm c}$ is reduced to that of $R_{\rm s}^{\prime}(r_{\rm c})+\Delta(r_{\rm c})$. Also in the independence, Eq. (23) shows $R_{\rm s}^{\prime}(r_{\rm c})+\Delta(r_{\rm c})$ is decreasing for $r_{\rm c}$. Since $R_{\rm s}^{\prime}$ approaches $0$ as $r_{\rm c}$ decreases, the lower bound $l$ of $r_{\rm c}$ is $l=R_{\rm s}+\Delta(l)$. Thus when $r_{\rm c}=l$, it reads $\rho(r)=g(r)$ in the absence of the central constant region. For a simple estimation of $l$, we regard $g(r)$ is linear and $g(l)\approx\rho_{0}/4$. From normalization $\int d{\bf r}\rho(r)=1$, it follows $l\approx\sqrt{2}R_{\rm s}$. Thus considering the above-mentioned upper bound, one finds an interval

for the formation of aSS, and this is, interestingly, consistent with our numerical data. In the other side in $r_{\rm c}\lesssim\sqrt{2}R_{\rm s}=1$ (note we set $J=1$ in numerics), we numerically observe the pattern becomes amorphous. $l\approx\sqrt{2}R_{\rm s}$ is by itself the upper bound of $R_{\rm s}+\Delta$, the radius of aSS. Depending on initial condition, there can appear multiple aSS. After transient period, each region separated farther than $r_{\rm c}$ from the others has its own disc (not shown here).

Our argument so far is also applicable for the change of the asynchronized-swam disc for $K<K_{\rm c}<0$ while $R_{\rm s}$ is replaced with $R_{\rm a}$. We call the pattern anomalous static async (aSA). An instance is shown in Fig. 5(c). We remark the async-counterpart of Eq. (24), $\sqrt{2}R_{\rm a}\lesssim r_{\rm c}<2R_{\rm a}$ is valid for the smaller criterion than $K_{\rm c}$. In the numerical data obtained for $N=400$, the lower bound shows deviation till $K\gtrsim-2$ (see the upper-left part of the aSA region in Fig. 8). There is more deviation for larger $K$ and smaller $r_{\rm c}$ in the corner. This observation is consistent with the reported stable-region of SA in Ref. [3] for $r_{\rm c}=\infty$, $K\lesssim K_{\rm c}\approx-1.2J$, in that the criterion $K_{\rm c}$ is a consequence of linear stability and a smaller $r_{\rm c}$ means a stronger perturbation. Different from aSS case, aSA does not become amorphous but a non-stationary pattern when it disappears by the decrease of $r_{\rm c}$.

The diameter of SPW that appears when $K=0$ for $r_{\rm c}=\infty$ is $D_{\infty}=2R_{2}$. Thus by Eq. (20), for $r_{\rm c}>2R_{2}$, a formation of the same SPW(s) is expected, and this is numerically observed. Next, for $r_{\rm c}<2R_{2}$, the results become different, also as expected. However, this time, the difference is not that apparent in shape; the inner and outer radii change a little though the cutoff value is decreased a lot to $r_{\rm c}\approx 1.6$ from $2.5\approx 2R_{2}$. Figure 6(a) is a pattern at $r_{\rm c}=1.6$, whose inner and outer radii are not discriminated so clearly from $R_{1}$ and $R_{2}$, respectively, known for $r_{\rm c}=\infty$ in Eq. (III).

Instead, we have found that a substantial change in density takes place while $r_{\rm c}$ decreases a lot. In Fig. 6(b), the red pluses come from the numerical data obtained at $r_{\rm c}=2.5$ that is slightly smaller than $2R_{2}=2.53..\,$. Each of them is the count of swarms in a narrow circular region between radii $r-\delta$ and $r+\delta$, where $r$ is discrete and small $\delta>0$ is a few times of the resolution of $r$. Their overall profile is guided by the bold red curve from Eqs. (13), (18), and (III). This guidance is plausible as $r_{\rm c}=2.5\approx 2R_{2}$. The blue crosses in Fig. 6(b) come from the numerical data at $r_{\rm c}=1.6$, and the blue bold curve is their overall shape as a guide to eyes. The substantial change in density is supposed to compensate for the small change in shape despite the considerable decrease of $r_{\rm c}$. We also observed the small change of shape mainly occurs below $r_{\rm c}\approx 1.8$. The inner and outer radii remain apparently same for $1.8\lesssim r_{\rm c}<2R_{2}$ while the density change takes place.

A heuristic argument for this observation is following. Let $r_{\rm c}$ be smaller that $2R_{2}$, and consider a circle of radius $r_{\rm c}$ whose center is in somewhere in SPW. If the center is near the outer boundary of SPW, the circle cannot cover the SPW. When this incomplete coverage is considered in Eq. (5), the reduced integral region would lead to density decrease at the center. Next, consider the center is near the inner boundary. In this case, the $r_{\rm c}$-radius circle covers the SPW when $r_{\rm c}-R_{2}>R_{1}$ is provided. For the density of such center, the integral region of Eq. (5) is conserved, and there is no explicit reason to decrease the density. Instead, an opposite reason was already due because the expected density decrease near the outer region requires mass conservation in the system. We consider the cross between the density profiles shown in Fig. 6(b) is the consequence.

If $r_{\rm c}-R_{2}>R_{1}$ does not hold, the integral region of Eq. (5) is always reduced whatever points of SPW is used for the center of $r_{\rm c}$-radius circle. This may indicate a termination of the shape-preserving density change. One then considers an interval

where such SPWs of negligible/considerable change in radius/density for $r_{\rm c}$ are respectively expected. The lower bound suggested this way is rather comparable with the numerical lower bound $r_{\rm c}\approx 1.8$ in that $R_{1}+R_{2}=1.74..\approx 1.8$ (see Eq. (III) with $J=1$ for the $R_{1}+R_{2}$ value). Thus when SPW is specified by the shape only, the condition for its formation is relaxed from $r_{\rm c}>2R_{2}$ to $r_{\rm c}\gtrsim R_{1}+R_{2}$.

We remark that some cases do not give SPW when $r_{\rm c}$ begins to be smaller than $1.8$, depending on initial conditions. No appearance of SPW below $r_{\rm c}\approx 1.8$ becomes more frequent for the smaller $r_{\rm c}$. Furthermore, SPW does not appear when $r_{\rm c}\lesssim 1.5$, at least, numerically. In brief, the numerical data suggests that the basin of attractor [34] for SPW begins to decrease from $r_{\rm c}\approx 1.8$, and then disappears below $r_{\rm c}\approx 1.5$.

It is worthy of noting that, when the pattern’s shape is not circular below $r_{\rm c}\approx 1.8$, it looks linear. According to the numerical data, linear shapes begin to appear for $r_{\rm c}\approx 1.8$, and disappear when $r_{\rm c}\lesssim 1$. Figure 7 shows a few linear bar-like patterns (Bar).

It is inappropriate we consider to view such patterns as the deformed APW by finite $r_{\rm c}$. As a first step to the understanding of such bar-like patterns, we below argue that each of them could be a deformed shape, by its own finite $r_{\rm c}$, of the bar given as the solution of Eq. (5) in the $r_{\rm c}\to\infty$ limit. An ironical situation here is no bar-like pattern has been observed yet in our numerical study neither reported in literature as far as we know. We thus conjecture that a solution bar exists for Eq. (5), of no stability, but it becomes stable for finite $1\lesssim r_{\rm c}\lesssim 1.8$ accompanied with proper deformation. If our scenario is valid, the pattern of Fig. 7(a) is most close to the solution bar among the figures therein because its $r_{\rm c}=1.8$ is larger than the others’.

We choose $x$-axis as the direction of bar with no loss of generality. In the numerical data, phase smoothly increases along the $x$-axis. The positive direction of the $x$-axis is selected in line with the phase increase. We thus introduce phase field $\theta({\bf r})=\theta(x)$ that is differentiable for $x$, and then write

in rectangular coordinate. No dependence on $r_{\rm c}$ shows Eq. (26) is that for the infinite cutoff.

Consider a bar of length $2l$ and of width $2w$, whose center is the origin. By the symmetry of $\rho(x,y)=\rho(|x|,|y|)$, we choose $\theta(0)=0$ for simplicity. We then give a relation between $\rho(x,y)$ and $\theta(x)$ as

in the consideration of i) $\theta_{i}$s are initialized with uniform random variables in our numerics and ii) any of them does not change in time in Eq. (2) with $K=0$. By Eq. (27), $\theta(x)$ increases from $\theta(0)=0$ to $\theta(l)=\pi$, and then remains same thereafter, i.e., $\theta(x)=\pi$ for $x>l$. By differentiabilty, $d\theta(l)/dx=0$ is given. By the symmetry mentioned above, it reads that $\theta(-x)=-\theta(x)$.

Plugging Eqs. (26) and (27) into Eq. (5), after a little algebra, one can obtain

for $c(l)=l+\int_{0}^{l}dx\cos\theta$ and $s(l)=\int_{0}^{l}dxx\sin\theta$. Obviously, $\rho(x,y)$ of Eq. (28) is considered for $|x|<l$ and $|y|<w$ while $\rho(x,y)=0$ elsewhere.

One integrates Eq. (28) to know the area $S$ of the solution bar is

Interestingly, the pattern in Fig. 7(a) covers a bar of area $\pi$ by length $\pi$ and width $1$, with a rather tolerable margin (see the guiding black rectangle). This numerical observation suggests $w\approx 1/2$ and $l\approx\pi/2$. We also observed the pattern becomes elongated and broadening as $r_{\rm c}$ decreases. Thus the patterns in Fig. 7(b)-(d) for the smaller $r_{\rm c}$ also cover the rectangle with the more margins (one can imagine the rectangle in each ones without difficulty). This is consistent with our scenario for the birth of bar-like patterns with finite $r_{\rm c}$, proposed above. We also observe in our numerical data that the bar-like pattern is so deformed for $r_{\rm c}\lesssim 1$ that it becomes amorphous.

The rectangle we have specified is a consequence of numerically motivated speculation encapsulated in Eq. (28). This equation is again converted with Eqs. (27) and (29) to

Its solution will lead to the self-consistent equations for $s(l)$ and $c(l)$, which can fix the value of $l$ and next that of $w$ in use of Eq. (29). This way, finding the size of the bar is dual to knowing its density profile. The shooting method looks one of the practical approaches for finding $l$. We leave related interesting task for future work as remarking $l\approx\pi/2$ and $w\approx 1/2$ is expected with the numerical data up to now.

We add that an 1-dimension result could be relevant to the size of the bar. In 1-dimension, $U_{\rm rep}(r,r^{\prime})=-|r-r^{\prime}|$ giving $\nabla_{r}^{2}U_{\rm rep}(r,r^{\prime})=2\delta(r-r^{\prime})$ is responsible for the repulsion part. Thus when $U_{\rm att}(r,r^{\prime})=(1/2)(r-r^{\prime})^{2}(1+J)$ is considered, uniform distribution of width $W(J)=2/(1+J)$ follows. We here consider the swarmalators on $y$-direction segment in the bar are effectively in 1-dimension in that their phases are same so that $U_{\rm att}({\bf r},{\bf r^{\prime}})=(1/2)(y-y^{\prime})^{2}(1+J)$, the distribution is uniform in the segment as written in Eq. (28), and its width in numerics with $J=1$ is $2w\approx 1=W(1)$. We leave expected explicit connection between the 1-dimension result and the bar for future work.

The steady-state condition Eq. (3) is a sufficient condition for Eq. (5). That is, the solution of Eq. (5) is not necessarily a steady state. We note $\rho(x)$ of Eq. (28) is not a steady state. For this, one can check the $y$-component of $v_{\rm att}({\bf r})$ given by $\int d{\bf r}^{\prime}(-\partial_{y}U_{\rm att}({\bf r},{\bf r}^{\prime}))\rho({\bf r}^{\prime})=-y$ cannot balance with its repulsion-counterpart for such $\rho(x)$ that vanishes at $x=l$ (note $\rho(l)\propto d\theta(l)/dx=0$). However, with $r_{\rm c}=1.8$, we observe a steady-state pattern that is rather close to the rectangle expected with Eq. (28), as depicted in Fig. 7(a). We interpret that the bar we speculated with $r_{\rm c}=\infty$ is properly deformed for $1\lesssim r_{\rm c}\lesssim 1.8$, and then it becomes stable while still keeping linear shape more or less. The patterns in Fig. 7 are the instances.