Phase transitions and pattern formation in ensembles of phase-amplitude solitons in quasi-one-dimensional electronic systems

TL;DR

This paper investigates phase transitions and pattern formation in quasi-one-dimensional electronic systems with degenerate phase and amplitude states, revealing complex topological configurations like vortices and soliton walls through analytical and Monte Carlo methods.

Contribution

It introduces an analytical and simulation-based study of topological defects and phase transitions in a two-fields XY-Ising model relevant to electronic systems.

Findings

Identification of temperature-driven confinement of topological defects

Discovery of soliton walls forming at lower temperatures

Potential realization of half-integer vortices with amplitude kinks

Abstract

Most common types of symmetry breaking in quasi-one-dimensional electronic systems possess a combined manifold of states degenerate with respect to both the phase $\theta$ and the amplitude $A$ sign of the order parameter $A\exp(i\theta)$. These degrees of freedom can be controlled or accessed independently via either the spin polarization or the charge densities. To understand statistical properties and the phase diagram in the course of cooling under the controlled parameters, we present here an analytical treatment supported by Monte Carlo simulations for a generic coarse-grained two-fields model of XY-Ising type. The degeneracies give rise to two coexisting types of topologically nontrivial configurations: phase vortices and amplitude kinks -- the solitons. In 2D, 3D states with long-range (or BKT type) orders, the topological confinement sets in at a temperature $T=T_1$ which binds together the kinks and unusual half-integer vortices. At a lower $T=T_2$, the solitons start to aggregate into walls formed as rods of amplitude kinks which are ultimately terminated by half-integer vortices. With lowering $T$, the walls multiply passing sequentially across the sample. The presented results indicate a possible physical realization of a peculiar system of half-integer vortices with rods of amplitude kinks connecting their cores. Its experimental realization becomes feasible in view of recent successes in real space observations and even manipulations of domain walls in correlated electronic systems.

Figures

Equations

$$H_{1D}=$$

$$-V\cos (2\theta )+(\hslash /4\pi \gamma )[v({\partial }_x\phi {)}^2+({\partial }_t\phi {)}^2/v]$$

$$\int dx\sum _{\alpha ,\beta }{J}_{\alpha ,\beta }{\mathcal{O}}_{\alpha }{\mathcal{O}}_{\beta }^{\ast }=$$

$$=\int dx\sum _{\alpha ,\beta }{J}_{\alpha ,\beta }\cos {\theta }_{\alpha }\cos {\theta }_{\beta }\exp (i({\phi }_{\alpha }-{\phi }_{\beta }))$$

$$H_c(\phi )$$

$$H_s(\theta )$$

$$\int dx\sum _{\alpha ,\beta }{J}_{\alpha ,\beta }\cos {\phi }_{\alpha }\cos {\phi }_{\beta }\exp (i({\theta }_{\alpha }-{\theta }_{\beta }))$$

$$H=$$

$$-A_{\perp }\sum _r{s}_r{s}_{r+\hat{y}}\cos ({\theta }_r-{\theta }_{r+\hat{y}})-J_{||}\sum _r{s}_r{s}_{r+\hat{x}}$$

$$\displaystyle\begin{cases}s_{x,y}\rightarrow-s_{x,y}\qquad&\mbox{for all $s,\theta$}\\ \theta\rightarrow\theta+\pi\qquad&\quad\mbox{at a given chain $y$}\end{cases}$$

$$H_0(A_{||},A_{\perp })=$$

$$-A_{\perp }\sum _r{s}_r{s}_{r+\hat{y}}\cos ({\theta }_r-{\theta }_{r+\hat{y}})$$

$$H_0=-A_{||}\sum _r{s}_r{s}_{r+\hat{x}}\cos ({\chi }_r-{\chi }_{r+\hat{x}})-A_{\perp }\sum _r\cos ({\chi }_r-{\chi }_{r+\hat{y}}),$$

$$\mathrm{Tr}\{s_r{s}_{r+\hat{x}}{e}^{-\beta H(A_{||},A_{\perp })}\}=\mathrm{Tr}\{s_r{s}_{r+\hat{y}}{e}^{-\beta H(A_{\perp },A_{||})}\}.$$

$${\tilde{H}}_0=$$

$$-\frac{1}{\beta }\sum _r\ln \cosh (\beta A_{\perp }\cos ({\theta }_r-{\theta }_{r+\hat{y}}))$$

$${\tilde{H}}_0=-A_{||}\sum _r\cos ({\theta }_r-{\theta }_{r+\hat{x}})-A_{\perp }\sum _r\left|\cos ({\theta }_r-{\theta }_{r+\hat{y}})\right|$$

$${\tilde{H}}_0=-A_{||}\sum _r\cos ({\theta }_r-{\theta }_{r+\hat{x}})-\frac{A_{\perp }}{4}\sum _r\cos 2({\theta }_r-{\theta }_{r+\hat{y}})$$

$$Z$$

$$\exp (\pi \beta A\sum _{R\ne R^{\prime }}(m_R+{\nu }_R)\ln \frac{|\mathbf{R}-{\mathbf{R}}^{\prime }|}{a}(m_{R^{\prime }}+{\nu }_{R^{\prime }})$$

$$+(\ln y_0)\sum _R(m_R+{\nu }_R{)}^2+\beta J_{||}\sum _r{s}_r{s}_{r+\hat{x}}).$$

$$F=E-TS\simeq J_{||}H-T\ln (HL/a^2)$$

$$\frac{T}{J_{||}}=\frac{H}{\ln (HL/a^2)}$$

$$\nu \approx \exp \left(-\frac{2J_{||}+2{\tilde{E}}_{c}A}{T}\right)$$

$$J_{||}(T,\nu )\approx \frac{T}{2}\ln \frac{1}{\nu }-{\tilde{E}}_{c}A$$

$$T_2\approx \frac{2{\tilde{E}}_{c}A}{\ln (1/\nu )}$$

$$E(l)=E_{s}l+2\widetilde{E_c}A+\frac{1}{2}\pi A\ln l.$$

$$n(l)$$

$$=\exp \left(\frac{\mu l-E_{s}l-2{\tilde{E}}_{c}A-\frac{1}{2}\pi A\ln l}{T}\right).$$

$$\nu =\sum _{l=1}^{\infty }ln(l)$$

$$=e^{-2E_c/T}{\mathrm{Li}}_{\alpha -1}(e^{-m})$$

$$\nu \le e^{-2{\tilde{E}}_{c}A/T}\zeta (\alpha -1)\approx e^{-2{\tilde{E}}_{c}A/T}$$

$$E(R)=\frac{1}{2}{\pi }^{2}AR\ln R+2\pi R{\tilde{E}}_{c}A+\pi R^2{E}_s$$

$$n(R)$$

$$=\exp \left(\frac{(\mu -E_s)\pi R^2-2\pi E_{c}R-\frac{1}{2}{\pi }^{2}AR\ln R}{T}\right)$$

$$\nu$$

Full text

The presented results indicate a possible physical realization of a peculiar system of half-integer vortices with rods of amplitude kinks connecting their cores. Its experimental realization becomes feasible in view of recent successes in real space observations and even manipulations of domain walls in correlated electronic systems.

topological defect, soliton, XY model, kink, stripe, confinement, aggregation

I Introduction

Symmetry broken electronic states give rise to topological defects: from extended objects like vortices or domain walls as solitonic lattices to microscopic solitons as anomalous quasi-particles and instantons in their dynamics. Role of solitons in electronic properties was appreciated in theories since mid 70’s (see reviews BK:84 ; YuLu , and also braz-moriond ; JSSS08 ; SB-Landau100:2009 for a more recent development). The solitons were firstly accessed in experiments on conducting polymers of early 80’s Heeger-RMP . New motivations came in early 2000’s from discoveries of the ferroelectric charge ordering (see reviews brazov-springer ; Monceau-springer ; Monceau-advances ) in organic conductors lebed-springer , from new accesses via nano-scale tunneling experiments latyshev-prl:05 ; latyshev-prl:06 in materials with Charge Density Waves (CDW), from optics of new conducting polymers vardeni .

Most clearly the existence of solitons was confirmed by recent observations in CDWs where the individual solitons have been visually captured by STM experiments Brun-soliton ; Kim:2012 . There were amplitude solitons (through which the order parameter changes the sign) corresponding to the long-sought special quasi-particles – the spinons or the holons. New opportunities are coming from the latest facilities of manipulating cooperative electronic states by short optical pulses and by a local STM injection. On this way, the networks of domain walls – the line aggregates of solitons – have been created in a layered CDW compound $\mathrm{TaS_{2}}$ Stojchevska:2014 ; Yoshida:2014 ; Vaskivskyi:2016 ; Ma:2016 ; Cho:2016 ; Svetin:2016 ; Gerasimenko:2017 ; Karpov:2018 .

A major puzzle, as well as the inspiration, coming from the experimental observations, is that the solitons were observed Brun-soliton within the low temperature ($T<T_{c}$) phase with the long-range 3D order. The hidden obstacle is the effect of the confinement appearing in higher dimensions $D>1$ JSSS08 ; braz-moriond . Commuting between degenerate minima on only one chain would lead to a loss of the interchain ordering energy proportional to $|x|$ – the length along the chain till the next soliton. In case of discrete symmetries (when solitons are the amplitude kinks as in the $Z_{2}$ case of the dimerization) the solitons are bound in topologically trivial pairs with an option for a subsequent phase transition to form cross-sample domain walls Bohr:1983 ; Teber:2001 ; Teber:2002 ; Karpov:2016 . But for a continuous symmetry (the phase degeneracy in superconductors (SC) and incommensurate charge density waves (ICDW) or the directional degeneracy in spin-isotropic antiferromagnets (AFM)), both phase and spin degeneracy in spin density waves (SDW)) the gapless mode can cure the interruption from the amplitude kink.

As we shall explain below, continuously broken symmetries allow for individual solitons entering the low-temperature phases with long- range ordered states. The solitons adapt by forming topologically bound complexes with half-integer vortices of gapless modes: $\pi$-rotons braz-moriond ; JSSS08 ; SB-Landau100:2009 . For cases of repulsion and attraction correspondingly, that results in the spin- or charge- roton configurations with charge- or spin- amplitude kinks localized in the core. Beyond the quasi-1D electronic systems, the interference of amplitude and phase (or angle) topological defects has been noticed already for superfluid phases of surface layers of $\mathrm{He_{3}}$ [Korshunov:1985, ] which provoked the studies of different generalizations of the XY model Lee:1985 ; Korshunov:1986 ; Shi:2011 ; Serna:2017 .

In this paper, we construct and study the coarse-grained model of XY-Ising type, which describes local and the long-range ordering under the monitored mean concentration of solitons. The paper is organized as follows. In Sec. II we set the playground: we introduce the general concept of solitons and, in particular, of solitons in quasi-1D electronic cooperative. In Sec. III we present the discretized version of the model and perform its qualitative analysis in various regimes. Section IV describes the results of numerical simulations for 2D and 3D systems. Finally, in Sections V and VI we present the Discussion and draw the Conclusions.

II Solitons in quasi-1D systems.

At the 1D level, for systems with a complex order parameter $\mathcal{O}=A\exp(i\varphi)$, the amplitude soliton $A(x=-\infty)\Rightarrow-A(x=+\infty)$ performs the sign change of the order parameter $\mathcal{O}$ at an arbitrary $\varphi=const$. At higher space and/or space-time dimensions, the requirement appears that the order parameter is mapped onto itself. Even being favorable in energy in comparison with the band electron, the soliton cannot be created dynamically already in $D=1$ space dimension, and is prohibited to exist even stationary at $D>1$. The resolution is to invoke the combined symmetry: the amplitude kink $A\Rightarrow-A$ coupled with the half-integer $\varphi\Rightarrow\varphi\pm\pi$ vortex of the phase rotation which compensates for the amplitude sign change. The resulting kink-roton complex, illustrated in Figs. 1, 2 allows for several interpretations in applications to ICDW or SC: 2D view is a pair of $\pi$-vortices sharing the common core bearing one unpaired spin which stabilizes the state; 3D view is a ring of a half-flux vortex line, its center confines the spin; at any $D>1$ this is a nucleus of the melted Fulde-Ferrel-Larkin-Ovchinnikov phase in the spin-polarized superconductors Fulde:1964 ; Larkin:1964 , or its CDW analogue Buzdin ; SB-Landau100:2009 .

II.1 Solitons in incommensurate CDWs

II.1.1 Amplitude solitons in 1D.

Here we quote some basic facts on microscopic origination of solitons in ICDW, see BK:84 ; braz-moriond ; JSSS08 ; SB-Landau100:2009 , based upon the Peierls-Froelich model with generalization to BCS case for filamentary superconductors Buzdin ; Machida . The ICDW is a crystal $A\cos(2K_{f}x+\varphi)$ of electron-hole pairs which order parameter is $\mathcal{O}_{ICDW}=A\exp(i\varphi)$. The singlet pair can be broken into spin $1/2$ components, but that will not be an expectedly liberated electron-hole pair at the gap rims $\pm\Delta_{0}$. Rather, there will be spin-carrying “amplitude solitons” (AS) – nodes of the order parameter distributed over the length $\xi_{0}=\hbar v_{F}/\Delta_{0}$, see BK:84 ; braz-moriond ; JSSS08 ; SB-Landau100:2009 . The total energy of the soliton is $\approx 2\Delta_{0}/3<\Delta_{0}$ making it energetically favorable with respect to the energy $\Delta_{0}$ of the undressed unpaired electron which now is trapped at the midgap state associated with the amplitude soliton. This electron brings the spin $1/2$, but the charge of the AS is compensated to zero by the local dilatation in the filled manifold of singlet vacuum states BK:84 . That makes the AS be a CDW realization of the “spinon”, with a direct generalization to a superconductor kwon:02 .

II.1.2 An unpaired spin in the ICDW or SC environment at $D>1$.

As a nontrivial topological object (the order parameter $O_{cdw}(x)=A(x)\exp(i\varphi(x))$ does not map onto itself), the pure AS (with $\varphi(x)=const$) is prohibited in $D>1$ environment. Nevertheless, the AS becomes allowed if it acquires phase tails with the total increment $\delta\varphi=\pm\pi$. The length of these tails $\xi_{\varphi}\gg\xi_{0}=\hbar v_{F}/\Delta_{0}$ is determined by the weak interchain coupling. As in 1D, the sign of $A(x)$ changes within the scale $\xi_{0}$ but further on, at the scale $\xi_{\varphi}$, the factor $\exp(i\varphi)$ also changes the sign, thus leaving the product in $O_{cdw}$ to be invariant. As a result, the 3D allowed particle is formed with the AS core $\xi_{0}$ carrying the spin, and the two $\pi/2$-twisting wings stretched over $\xi_{\varphi}$, each carrying the charge $e/2$. This picture can be directly reformulated for quasi-1D superconductors by redefining the meaning of the phase. The phase tails form now the elementary $\pi$-junction braz-moriond ; kwon:02 .

II.2 Solitons within the bosonisation language

II.2.1 Solitons in spin-singlet systems.

Recall a universal microscopic insight to excitations in 1D electronic systems, first in application to the spin-gap cases – SC or CDW. The starting single chain level is well described by the bosonization language, see Giamarchi:book ; Tsvelik:book . The Hamiltonian

[TABLE]

is written in terms of phases $\theta$ and $\varphi$ for spin and charge degrees of freedom correspondingly. The energy $\sim V$ comes from the backward exchange scattering $V\sim g_{1}$ of electrons. The pair-breaking excitation - the $s=1/2$ spinon, is the soliton connecting the degenerate minima of $H_{1D}$: $\theta\Rightarrow\theta\pm\pi$. The singlet order parameter, for either SC or CDW (with the inversion of the parameter $\gamma$ depending on a definition of the charge phase $\varphi$) is $\mathcal{O}_{SC,CDW}\sim\cos\theta\exp(i\varphi)$. Its amplitude $A=\cos\theta$ changes the sign across the allowed $\pi$ soliton in ${\theta}$, hence the spinon is an alternative description of the same amplitude soliton which appears in BCS-Peierls type models.

For the ensemble of chains, labeled below by indices $\alpha,\beta$, the relevant interchain coupling reads

[TABLE]

II.2.2 Solitons in the Mott insulator: SDW route.

Consider the quasi-1D system with repulsion at a nearly half-filled band, which is the SDW route to a general doped antiferromagnetic Mott-Hubbard insulator SB-Landau100:2009 . The single-chain state has no magnetic long-range order, gapeless spin degrees of freedom are described by the same phase $\theta$ as in (1), (2), charge degrees of freedom are described by the chiral phase $\varphi(x,t)$ corresponding to the CDW version of (1), (2). The 1D bosonized Hamiltonian Giamarchi:book ; Tsvelik:book becomes $H_{1D}=H_{c}+H_{s}$

[TABLE]

Here $v_{c}$ and $v_{s}$ are the charge and spin phase velocities, $U$ is the Umklapp scattering amplitude $U\sim g_{3}$ which is responsible for the charge gap formation DL ; Luther+Emery . The constant $\gamma$ (commonly called the Luttinger liquid parameter – $K_{\rho}$, dependent on interactions) determines several regimes: $\gamma<1$ at any repulsion, then $U$ is preserved if it is already build-in - this is the generic Mott state case for a half-filled band; $\gamma<1/2$ at a stronger repulsion, then $U$ is spontaneously generated even away from the bare half-filling - this is the Charge Ordering case brazov-springer .

In 1D, the excitations are the holon/doublon pairs seen as the $\mp\pi$ solitons in $\varphi$, and the spin sound in $\theta$, which are decoupled. In a quasi-1D system, multiple forms of the interchain coupling appear among which the relevant one is the coupling via the SDW order parameter: the staggered magnetization $\mathcal{O}_{sdw}\propto S_{x}+iS_{y}\propto\cos\varphi\exp(i\theta)$.

For the ensemble of chains, analogously to (2), the relevant interchain coupling reads

[TABLE]

We see that the $\pi$-soliton in $\varphi$ becomes the amplitude kink ($\cos\varphi\rightarrow-\cos\varphi$), and to survive in $D>1$ it should enforce the $\pi$ rotation in $\theta$, then the sign changes in both of the two factors composing $\mathcal{O}_{sdw}$ compensate and the configuration becomes allowed.

This construction can be generalized beyond quasi-1D systems by considering a vortex configuration bound to an unpaired electron. Extending to AFMs like the pristine $\mathrm{CuO_{2}}$ planes, the SDW becomes a staggered magnetization; the soliton becomes a hole, which motion leaves the string of reversed AFM sublattices Bulaevskii:1968 ; Brinkman:1970 ; the $\pi$-wings become the magnetic semi-vortices. The resulting configuration is a half-integer vortex ring of the staggered magnetization (a semi-roton) with the holon confined in its center (see Fig. 1b – in 2D when the vortex ring is reduced to the pair of vortices.) Such a combined semi-vortex – Fig. 1 might be significant also in sliding incommensurate SDWs which possesses an even richer order parameter NK+SB:00 .

II.3 Inverse rout: from stripes to solitons and fractional vortices.

At low temperature and at large concentration, the amplitude solitons tend to organize themselves macroscopically, as a periodic array of stripes crossing the whole sample in the transverse direction. The arrays of spin solitons have been identified experimentally e.g. in spin-Peierls systems Horvatic:1999 . No confinement energy is lost in this perfect state since the order parameter alternates coherently across all chains. But when the solitonic lattice melts, the defects appear such as termination points of lines of the amplitude nodes. To be allowed, each such a point should be complemented by the pair of $\pi$-vortices, Fig. 2.

III Formulation of the effective model and its qualitative analysis

A discrete version of Hamiltonians (1), (2) or (3), (4) can be written as Bohr:1983

[TABLE]

Here the Ising variable $s=\pm 1$ corresponds to the normalized real amplitude of the order parameter and the angle $\theta$ corresponds to the phase of the complex order parameter; $r$ runs over sites of a square lattice. The last term fixes the chemical potential $\mu_{s}$ of amplitude solitons (and therefore their average number in the system, see Karpov:2016 for the details). For the last constraint to be held, the constant $J_{||}=J_{||}(\nu,T)=(E_{s}-\mu_{s}(\nu,T))/2$ has to be found self-consistently by maintaining the average concentration (per site) of solitons $\nu=-\partial\Omega/\partial\mu_{s}=T\partial\ln Z/\partial\mu_{s}$. The discretization length $a$ for integrals (2) or (4) is chosen to be equal to the minimal in-chain distance between solitons, which is determined by the short-range repulsion of solitons at a given chain.

The Hamiltonian (5) is invariant under a symmetry (gauge) transformation of simultaneous flipping of all Ising and all XY spins within any given chain:

[TABLE]

which implies vanishing of macroscopical averages of Ising and XY fields separately: $\langle s_{r}\rangle=0$, $\langle e^{i\theta_{r}}\rangle=0$. Still a combined variable $\Delta_{r}=s_{r}e^{i\theta_{r}}$ can have nontrivial average and correlation functions.

Now we give a brief “dictionary” for different terms we use throughout the paper from the point of view of the discretized Hamiltonian (5) the Ising variables $s_{r}$.

An amplitude soliton (of just a soliton for brevity) is an interruption of the Ising order along the chains’ ($x$) direction, i.e. the soliton is present when $s_{r}s_{r+\hat{x}}=-1$.

A bisoliton consists of two solitons at the minimal in-chain distance $a$ between them; in the language of the Ising variables $s_{r}$ it is a reversed spin in the otherwise ordered part of the chain.

A solitonic rod (or a solitonic line) is transverse to the chains aggregation of solitons, when solitons at several neighboring chains have the same $x$-coordinate. Generalized to the 3D case such a structure would be called a solitonic disk. For the transverse aggregate of bisolitons instead of solitons, we accordingly use terms bisolitonic rod (in 2D) or bisolitonic disk in 3D.

A solitonic wall can appear when a solitonic line grows and crosses the whole transverse section of the sample. Likewise, a bisolitonic wall appears, when a bisolitonic rod grows across the sample.

III.1 Zero chemical potential of amplitude solitons

We start studying the model (5) from the simple, but important case $J_{||}=0$. That corresponds to the zero chemical potential of solitons (with respect to their activation energy), maintained by a solitonic reservoir, with which they can exchange the solitons with no energy cost. For an isolated system, the role of such a reservoir is played by the solitonic domain walls Bohr:1983 ; Teber:2001 ; Karpov:2016 .

For $J_{||}=0$, the model becomes an anisotropic version of the XY-Ising model Yosefin:1985 :

[TABLE]

with the XY interaction only along the chains’ $x$-direction and XY-Ising interaction only along the inter-chain $y$-direction.

Using the identity $s_{r}\cos(\theta)=\cos(\theta+(1-s_{r})\pi/2)$ we can show that the model is self-dual with respect to transformation of the angle variable $\chi_{r}=\theta_{r}+(1-s_{r})\pi/2$. Indeed, the Hamiltonian (7) transforms to

[TABLE]

which is equivalent to (7) under the substitutions $\theta\leftrightarrow\chi$, $A_{||}\leftrightarrow A_{\perp}$, $x\leftrightarrow y$. This means that all thermodynamic functions have to be symmetric with respect to the interchange $A_{||}\leftrightarrow A_{\perp}$. Also, the correlation functions have to be symmetric with respect to the simultaneous change $A_{||}\leftrightarrow A_{\perp}$, $x\leftrightarrow y$, for example

[TABLE]

Particularly, in the limit $A_{\perp}=A_{||}$ we have $\langle s_{r}s_{r+\hat{x}}\rangle=\langle s_{r}s_{r+\hat{y}}\rangle$ .

The summation over the Ising variables in the partition function generated by the Hamiltonian (7) can be performed exactly. For a fixed configuration of the XY field $\theta_{r}$ in (7), the Hamiltonian describes a set of non-interacting transverse Ising chains with random interactions between the neighboring spins $s_{r}$ across the chain.

Let $s_{r}s_{r+\hat{y}}=\tilde{s}_{r}=\pm 1$, then the Hamiltonian (7) becomes $H_{0}=-A_{||}\sum_{r}\cos(\theta_{r}-\theta_{r+\hat{x}})-A_{\perp}\sum_{r}\tilde{s}_{r}\cos(\theta_{r}-\theta_{r+\hat{y}})$, where $\tilde{s}_{r}$ are independent Ising variables. Summing out $\tilde{s}_{r}$ we arrive at the effective Hamiltonian for the modified XY-field

[TABLE]

The interaction is fundamentally anisotropic: it is $2\pi$-periodic in the chains’ ($x$) direction, and only $\pi$-periodic in the transverse ($y$) direction. That differs from the isotropic Korshunov model Korshunov:1985 ; Lee:1985 ; Korshunov:1986 , where both $\pi$ and $2\pi$ interactions coexist in all directions. In the low-temperature limit $\beta A_{\perp}\gg 1$, the Hamiltonian (8) becomes

[TABLE]

For slowly varying configurations (apart from the transverse $\pi$-jumps in $\theta_{r}$), keeping only the lowest order periodicities in increments of $\theta$, the Hamiltonians (8), (9) can be mapped upon the form

[TABLE]

The isotropic version of this model (i.e. when both $\theta$- and $2\theta$-terms are present in both $x$ and $y$ directions simultaneously) has been studied in Korshunov:1985 ; Lee:1985 ; Korshunov:1986 . For $A_{||}\ll A_{\perp}$ it gives rise to half-integer $\pi$-vortices (allowed by the $\sim A_{\perp}$ term) which are bound in pairs by “strings” generated by the $2\theta$ term (here, the strings are borders separating $\pi$-mismatches of $\theta$ Korshunov:1985 ; Lee:1985 ; Korshunov:1986 ). But for the model we study here, where each type of terms is associated with a particular direction, the situation is different. Even when $A_{||}\simeq A_{\perp}$, the model allows for $\pi$-vortices with strings attached in $x$-direction (so that $\pi$-mismatches in phase $\theta$ take place only in $y$-direction, for an example see Fig. 3). These $\pi$-vortices take the role of the amplitude solitons which are hiddenly present in the excluded Ising variables. The $A_{||}$ term still provides the confinement of $\pi$-vortices into pairs, but only in the $y$-direction.

III.2 Fixed non-zero chemical potential of solitons

Now consider the complete Hamiltonian (5) with $J_{||}\neq 0$, putting for simplicity $A_{||}=A_{\perp}\equiv A$. Using the Villain approximation we can transform (5) to the Coulomb gas representation as described in Appendix A:

[TABLE]

Here $a$ is a lattice spacing cutoff, $y_{0}=\exp(-\pi^{2}\beta A/2)$, summation goes over all integer values of $m_{R}$ and all configurations of spins $s_{r}=\pm 1$, the index $r$ corresponds to the original lattice and $R$ corresponds to the dual one; half-integers $\nu_{R}=\frac{1}{4}(s_{r+\hat{x}}s_{r+\hat{x}+\hat{y}}-s_{r}s_{r+\hat{y}})$. The vortex “charge” $m_{R}+\nu_{R}$ becomes dependent on the local configuration of spins. Also apart from the long-range “Coulomb” energy associated with half-vortices, there is an additional local Ising-term $\sim J_{||}$. This term is important, because as long as $J_{||}\neq 0$ it leads to confinement of the transverse half-vortex pairs, situated at the ends of solitonic rods (see also Sec. III.3.3).

Since the fugacity $y_{0}$ is small, we can take into account only the leading terms with minimal vorticities which are $m_{R}=0$ if $\nu_{R}=0$, $m_{R}=0,-1$ if $\nu_{R}=+1/2$, $m_{R}=0,+1$ if $\nu_{R}=-1/2$.

III.3 Qualitative phase diagram

In this section we reconstruct the phase diagram of the system, based on the qualitative interpretation of partition function (11) corresponding to the Hamiltonian (5). Again, for simplicity, we put $A_{||}=A_{\perp}=A$.

We start construction of the phase diagram with the limiting cases: the weak coupling limit $\beta J_{||}\rightarrow 0$ (corresponding to the zero chemical potential of solitons, when the Hamiltonian is given by (8)) and the strong coupling limit $\beta J_{||}\rightarrow\infty$ (corresponding to a strongly negative chemical potential of solitons $\mu_{s}-E_{s}\rightarrow-\infty$, hence to their low concentration).

III.3.1 Weak coupling limit for Ising spins $\beta J_{||}\rightarrow 0$

Consider, first, the weak coupling limit for the Ising spins, $\beta J_{||}\rightarrow 0$, which corresponds to the case when the chemical potential of amplitude solitons is trapped at zero (this happens e.g. below the temperature of solitons’ walls formation in 3D or in the limit $T\rightarrow 0$ in 2D, see Secs. III.5, G below). This limit was already discussed in Sec. III.1. Integer vortices disintegrate into pairs of half-integer ones, which lowers the system energy. A string of XY phase mismatch is attached to every half-integer vortex (Fig. 3). We note that unlike the model of Korshunov Korshunov:1985 ; Lee:1985 ; Korshunov:1986 , where the strings did cost some energy and thus were uniquely physically defined, in our case, only the solitonic rod is uniquely defined, while for the string only its $y$-coordinate is fixed. Still there is an ambiguity from which side, from the left or from the right, the string is attached (for example, applying the gauge transformation (6) to all the chains in the upper half-plane in configuration in Fig. 3 we will get the string attached to the half-vortex from the left side). This differs from conventional XY-vortices, where the direction of the ray of $2\pi$-jump of the phase can be chosen completely arbitrarily. Since the string energy density is zero, then the string does not affect the statistics of the associated vortices. Then, analogously to BKT-case, we expect that at low temperatures all half-vortices are bound into neutral pairs and the transition proceeds via their dissociation. In comparison to the usual BKT-transition, this one has the temperature approximately 4 times lower $T^{half}_{BKT}=T_{BKT}/4$ and the 4 times larger jump of helicity modulus Korshunov:1986-FFXY .

III.3.2 Strong coupling limit for Ising spins $\beta J_{||}\rightarrow\infty$

Now consider the opposite, strong coupling limit for Ising spins $\beta J_{||}\rightarrow\infty$; in terms of the original solitons this means that their chemical potential goes to minus infinity i.e. their concentration vanishes. If all Ising variables were equal (showing the perfect Ising long-range order, i.e. the total absence of amplitude solitons), even then XY spins would show only the algebraic order $\langle e^{i(\theta_{r}-\theta_{r^{\prime}})}\rangle\sim|{\mathbf{r}}-{\mathbf{r}^{\prime}}|^{-1/2\pi K}$ at low temperatures.

But taking into account the gauge symmetry (6) it follows that $\langle S_{r}S_{r^{\prime}}\rangle=0$ and $\langle e^{i(\theta_{r}-\theta_{r^{\prime}})}\rangle=0$ for any $r$ and $r^{\prime}$ lying at different chains, even at $T=0$, i.e. for $\beta J_{||}\rightarrow\infty$. The natural question then arises: is there Ising order at one given chain at $T\neq 0$? The answer is no: the combined topological defect of a unit height and an arbitrary length (Fig. 4) has a finite energy $E$. So applying the known argument for 1D systems Landau-argument ; Landau5 we see that such defects break the in-chain Ising long-range order.

Nevertheless, the long-range correlations can be present for the physically significant combined order parameter $\Delta_{r}=S_{r}e^{i\theta_{r}}$. That can already be seen from the fact that for $J_{||}\rightarrow\infty$ all the half-vortices are suppressed, and we have the usual Berezinskii-Kosterlitz-Thouless (BKT) transition at the temperature $T_{BKT}\approx 0.9A$ [Shugard:1978, , Olsson:1995, ], with $2/\pi$-jump of the helicity modulus; for big values $J_{||}\gg A$ the transition should also be of the BKT type. Then we expect that in the low-temperature (ordered) phase the correlation function $\langle\Delta_{r}\Delta_{r^{\prime}}\rangle$ decays algebraically, and in the high-temperature disordered phase it decays exponentially.

III.3.3 Intermediate $\beta J_{||}$ and the phase diagram

For intermediate $J_{||}$ we expect a crossover between the two analyzed limiting cases: strong and weak coupling limits for the Ising subsystem, when integer or half-integer vortices dominate correspondingly. Fig. 5 shows the expected phase diagram of the system.

In the thermodynamic limit, an arbitrarily small $J_{||}$ prevents formation of single half-vortices (Fig. 3). In the course of heating, the order-disorder transition is of the BKT type with the dissociation of integer vortex-antivortex (VA) pairs. Half-vortex pairs may also be present (Fig. 6a), but they are bound by a constant (“confinement”) force, so they always have a finite size and are irrelevant under rescaling.

But if $J_{||}$ is exactly zero, then the confinement force between half-antivortices vanishes. Therefore, they can dissociate, and this indeed takes place at the temperature $T_{BKT}^{half}\approx T_{BKT}/4$.

The crossover between these two regimes ($J_{||}>0$ and $J_{||}\equiv 0$) occurs at some small $J_{||}$, which is related to the finite size of the system. In order to find this crossover temperature we can recast Kosterlitz-Thouless argument Kosterlitz:1973 . The free energy of a single half-vortex is

[TABLE]

where $L$ and $H$ are the system dimensions in the chains’ and in the transverse directions respectively.

If half-integer VA pairs can escape from the linear confinement at $T<T_{BKT}/4$, then the BKT transition will occur via dissociation of these pairs. If half-VA pairs escape from the linear confinement at $T_{BKT}/4<T<T_{BKT}$, then the transition from the algebraically ordered to the disordered state occurs via Ising mechanism (such a possibility was recently explored by Shi et.al.Shi:2011 ).

From (12) we find the critical value of $T/J_{||}$, when half-VA pairs become important:

[TABLE]

This crossover temperature (the dashed line in Fig. 5) goes to infinity as $H\rightarrow\infty$.

III.4 Pattern formation

In this subsection we describe the pattern formation in relation to the phase diagram. Consider the cooling trajectory across the phase diagram when the number of amplitude solitons is fixed (red line in Fig. 5). In the high-temperature phase, the vortices exist in a plasma state and half-vortices are strongly confined via the linear potential, typically only the smallest half-VA pairs with the transverse distance $1$ are present. At the temperature $T_{1}$ the system enters the region of algebraic order and the vortices become confined in VA pairs. The earlier free amplitude solitons cannot survive anymore and they start to be dressed by half-VA pairs. Upon further cooling at some crossover temperature $T_{2}$ (an estimate is given below in (16)) the solitons start to aggregate in the transverse direction, thus elongating the half-VA pairs. At some size-dependent temperature $T_{F}\sim T_{2}$, $J_{||}(T_{2})$ becomes so small (such that (13) is fulfilled) that the solitonic walls grow across the whole sample (Fig. 6b,c), then the confinement force does not play any further role, half-VA pairs are weakly confined only by 2D Coulomb log-interaction. Since at such low temperatures half-VA pairs are already very rare, then the regular “stripe phase” sets in.

The picture of the two phase transitions resembles the previously studied case of Ising model Karpov:2016 . For that case, with lowering the temperature, firstly the long-range order (versus quasi-long-range algebraic order described here) is established, and then at the lower temperature $T_{2}$ the solitonic walls grow across the whole system.

III.5 Estimates for $T_{1},T_{2}$.

In this subsection we find estimates for the transition temperatures $T_{1},T_{2}$.

$T_{1}$ is the BKT transition temperature, $T_{1}\approx T_{BKT}\approx 0.9A$ [Shugard:1978, , Olsson:1995, ], and corrections in soliton concentration $\nu$ can be found perturbatively.

$T_{2}$ can be estimated by extrapolating of $J_{||}(T)$ from high to low temperatures in the approximation of non-interacting solitons (assuming their small concentration: $\nu\ll 1$). Each soliton bears the energy $2J_{||}$ via Ising subsystem; XY subsystem tries to “cure” the defect, then the half-VA pair forms with the minimal transverse size of $1$ site and the energy of the half-vortex cores $2E_{core}\equiv 2\tilde{E}_{c}A$ ($\tilde{E}_{c}$ is the dimensionless core energy of a half-vortex). The total energy of such combined defect is $2J_{||}+2\tilde{E}_{c}A$. Hence their concentration is

[TABLE]

Inverting this relation for $J_{||}$, we get

[TABLE]

From the condition $J_{||}(T,\nu)\approx 0$, when the approximation of non-interacting solitons breaks and they start to aggregate in transverse “rods”, we find an estimation for $T_{2}$:

[TABLE]

At lower $T\lesssim T_{2}$ we cannot further neglect the interaction of solitons. Two scenarios may be envisaged then: (i) $T_{2}$ can be a crossover temperature, below which solitons start to aggregate in transverse lines that are gradually growing as $T\rightarrow 0$; (ii) $T_{2}$ can be a phase transition temperature, at which infinite lines appear, crossing the whole sample. In the next section we test these two scenarios numerically with results in favor of the scenario (i).

The estimate for $T_{2}$ can also be done in another way – using the argument of solitonic transverse “rods” Bohr:1983 . When small transverse aggregates become favorable, they can further stick together, growing into infinite walls. This can be supported by the following argument. Let transverse lines of lengths $l>2$ became favorable at $T\rightarrow 0$. The energy of two such lines is $\sim 2\ln l$. Merging them into a single line of length $2l$, the energy becomes $\sim\ln(2l)=\ln l+\ln 2$. Therefore for $l>2$ such merging always lowers the energy, which provides the possibility for the rods to grow to infinity.

Here we consider a thermodynamically equilibrium system of such transverse segments of solitonic linesBohr:1983 (the solitonic “rods”). We shall work with the grand canonical ensemble of solitons with the chemical potential $\mu$ to be determined a posteriori by the fixed average concentration $\nu$.

Let the amplitude solitons be aggregated in rods of variable lengths $l$. Because of exchange among the rods, their chemical potentials are $\mu(l)=\mu l$. Their energies $E(l)$ can be estimated as:

[TABLE]

Here $E_{s}$ is the core energy of the amplitude soliton; $2E_{c}\equiv cA$ is the energy of two half-vortex cores and $\tfrac{1}{2}\pi A\ln l$ is the interaction between them.

The distribution of rods over their lengths Bohr:1983

[TABLE]

Denote $m=(E_{s}-\mu)/T>0$, $\alpha=\pi A/2T$, so that $n(l)=e^{-2\tilde{E}_{c}A/T}l^{-\alpha}e^{-ml}$. Then the total concentration of solitons is

[TABLE]

where $\operatorname{Li}_{s}(z)=\sum_{n=1}^{\infty}z^{n}/n^{s}$ is the polylogarithm function; its maximum value is reached at $z=1$, when the series is just Riemann zeta function: $\operatorname{Li}_{s}(1)=\zeta(s)$.

As $T\rightarrow 0$, $e^{-2\tilde{E}_{c}A/T}\rightarrow 0$, so to keep $\nu$ fixed, the factor $\operatorname{Li}_{\alpha-1}(e^{-m})$ need to grow to infinity. However for $\alpha>2$ ($T<\pi A/4\simeq T_{BKT}$) it is bounded from above; therefore the ensemble of finite lines cannot further accommodate all the solitons in the system. This means that infinite (crossing the sample width $H$) domain walls have to grow.

Assuming that the corresponding critical temperature $T_{2}\ll T_{BKT}\approx 0.9A$ (which we then check for self-consistency), so that $\alpha\gg 1$ and $\zeta(\alpha-1)\approx 1$, we get

[TABLE]

Therefore, the characteristic temperature for transverse walls nucleation is $T_{2}\approx 2\tilde{E}_{c}A/\ln(1/\nu)$; this estimate is consistent with the one based on the $J_{||}=0$ argument (see (16)). The assumption $T_{2}\ll T_{BKT}$ is satisfied when $\nu\ll 1$.

When the temperature decreases towards $T_{2}$, the solitonic rods begin to elongate, then above some critical length $l^{*}$ it may become energetically favorable for them to group into pairs and form the bisolitonic rods. In the next section we estimate this characteristic length $l^{*}$ and discuss how it is affected by the system’s anisotropy.

III.6 Role of the anisotropy

In this subsection we discuss the effect of the anisotropy $A_{\perp}/A_{||}\neq 1$ on the growth of the rods and the $T_{2}$ transition, and also find the characteristic length $l^{*}$, above which solitonic rods combine to bisolitonic ones.

Simple scaling arguments show that BKT-transition temperature $T_{1}$ depends only on $A=\sqrt{A_{||}A_{\perp}}$ [Pereira:1995, ], so the expression $T_{1}\approx 0.9A$ holds. Likewise, the vortex-antivortex interaction energy also depends only on $A$, therefore the estimate (16) for $T_{2}$ also holds. However, as we argue below, anisotropy can affect the character of the $T_{2}$ transition: it can facilitate the growth of either solitonic or bisolitonic rods.

Consider a low-temperature regime of a system with just two solitons (or one bisoliton) confined to a given chain. Depending on the ratio of $A_{\perp}/A_{||}$ it is energetically favorable for the system to either have separate solitons or a tight bisoliton (Fig. 7). Configuration (a) with two kinks in $\theta$-phase has the energy $4A_{||}l$, configuration (b) with two half-VA pairs has the energy $(4\tilde{E}_{c}+\pi\ln l)\sqrt{A_{||}A_{\perp}}$, and configuration (c) with a tight bisoliton has the energy $4A_{\perp}$ (for generality we also took into account the transverse length of the rods $l$; in case of Fig. 7 we have $l=1$). Depending on the ratio $A_{\perp}/A_{||}$ one of the three configurations is favored for a given $l$.

Energy of the configuration of Fig. 7 (a) strongly depends on $l$ and it becomes unstable towards configuration (b) already for rather small values $l\gtrsim\tilde{E}_{c}\sqrt{A_{\perp}/A_{||}}$. Since the energy of the configuration (b) depends on $l$ only logarithmically, then such type of configurations can be stabilized for a broad range of rods lengths $l<l^{*}\sim\exp(\frac{4}{\pi}\sqrt{A_{\perp}/A_{||}})$, which exponentially depends on the anisotropy of the system.

This means that for $A_{\perp}/A_{||}\gtrsim\tilde{E}_{c}^{2}$, when stable solitonic rods exist, their growth occurs in two stages: first, solitonic rods grow to the characteristic size $l^{*}$, after which it becomes favorable to glue two solitonic rods (with the energy $\sim A\ln l$) to one bisolitonic rod with compensating pairs of half-vortices at their terminals; then only the energy $\sim A_{\perp}$ is left. After this gluing, the system enters the regime of slow further growth of bisolitonic rods, similarly to the previously studied Ising case Bohr:1983 ; Teber:2001 ; Karpov:2016 . This means that here also in the 2D case $T_{2}$ is only a crossover and not a phase transition temperature.

We note, however, that for non-equilibrium processes, like supercooling the system with a sudden quench (as it is implemented in our Monte Carlo modeling, see Sec. IV.4), this first stage can be drastically promoted, since a growing solitonic wall creates uncompensated half-vortices at the ends which attract mobile smaller rods, while the slow process of gluing of low-mobile big rods can be kinetically suppressed.

The above analysis does not take into account the screening effects from multiple interactions of vortices Kosterlitz:1973 , which limits the literal applicability to the case of distant and non-interacting solitonic or bisolitonic rods. That leads to a restriction upon their characteristic length $l^{*}\ll 1/\nu$ ($l=1/\nu$ is such a length when every square $l\times l$ contains one rod of length $l$, which means that $\nu=l/l^{2}=1/l$). When the rods become too long i.e. for $l^{*}\gtrsim 1/\nu$, the interaction among them cannot be further neglected and the system enters the strongly interacting regime of plasma of half-vortices and antivortices.

We conclude that in 2D the walls’ formation temperature $T_{2}$ is a crossover, becoming for finite samples effectively a phase transition with the transverse size dependent temperature $T_{F}(H)$. Depending on the ratio of the interchain couplings $A_{\perp}/A_{||}$, either solitonic or bisolitonic walls grow across the sample at $T<T_{F}(H)$. Contrarily, only topologically trivial bisolitonic walls could be formed for the case of the real order parameter described by the Ising model Bohr:1983 ; Teber:2001 ; Karpov:2016 .

III.7 Extension to a 3D system

For a 3D system we expect the phase diagram similar to 2D one (Fig. 5) except for a character of phase transitions. In the low-temperature ordered phase, there is still no Ising order. This can be seen using the same argument as in Sec. III.3.2): the combined XY-Ising defect has a finite energy, so the entropy gain in the thermodynamic limit will always favor such defects at all $T>0$, thus breaking the Ising order. This leads to only a “nematic” $2\theta$-order in XY subsystem (i.e. nonzero expectation value $\langle\cos(2(\theta_{r}-\theta_{r^{\prime}}))\rangle$ opposed to a pure 3D XY model that possesses true long-range $\theta$-order $\langle\cos(\theta_{r}-\theta_{r^{\prime}})\rangle\neq 0$ for $|\mathbf{r}-\mathbf{r^{\prime}}|\rightarrow\infty$). In terms of the original picture of solitons this means that an isolated soliton with a finite energy is allowed to exist, since the interchain order which it breaks can be cured by the adjusted phase degree of freedom; such isolated solitons break the Ising order at any non-zero temperature. This “nematic” order can be also described by the combined XY-Ising variable $\Delta_{r}=s_{r}e^{i\theta_{r}}$, which acquires a non-zero expectation value below the ordering-transition temperature.

Repeating the arguments of the Sec. III.5, we obtain the estimate (16) for the walls’ formation temperature $T_{2}$, with a difference that now it is the true transition rather than a crossover.

Instead of the vortex pairs, the lowest energy topological excitations in 3D XY model are the closed vortex loops. Proliferation of finite vortex loops into infinite ones triggers the order-disorder phase transition Kleinert:book-multivalued .

The energy of the loop of radius $R$ is:

[TABLE]

where the first term describes the loop vorticity energy Williams:1987 , the second term is the core energy of the vortex loop (distributed along the length of the loop; $\tilde{E}_{c}$ is a numerical constant, which is generally different from the analogous one in 2D), and the third term is the energy of the enclosed solitonic disk. The chemical potential of the disk is $\mu(R)=\pi R^{2}\mu$, therefore the concentration of disks with a radius $R$ is:

[TABLE]

Hence the total concentration of solitons is:

[TABLE]

Clearly, even if $\mu-E_{s}=0$, as $T\rightarrow 0$ the sum in the RHS goes to [math]. This means that the disk-phase cannot accommodate all the solitons, and at some $T=T_{2}$ domain walls need to be produced via a phase transition.

Expressions (21)-(23) hold also for the anisotropic case with $A=A_{\perp}^{2/3}A_{||}^{1/3}$ [Pereira:1995, ]. Analogously to the 2D case the effect of anisotropy can play an important role for the mechanism of growth of the solitonic disks across the sample. If $A_{\perp}/A_{||}\gtrsim\tilde{E}_{c}^{3}$, then solitonic disks are stable and can grow up to the characteristic radius $R^{*}\sim\exp(\frac{4}{\pi^{2}}\sqrt[3]{A_{\perp}/A_{||}})$, after which they glue to bisolitonic disks.

IV Numerical simulations

In this section we study numerically the phase transitions, described in the previous section qualitatively. We model the cooling behavior of a system described by the Hamiltonian (5) using Metropolis Monte Carlo (MC) method. The numerical algorithm preserves the number of amplitude solitons (so the term $J_{||}\sum_{r}s_{r}s_{r+\hat{x}}=const$, where $r$ runs through all sites of the lattice) which substantially simplifies the analysis. There is no restrictions upon the phase degree of freedom $\theta$. An ambiguity of the visual representation related to the gauge transformation (6) is resolved by maximizing the amount of $s=+1$ Ising spins when drawing the figures.

IV.1 Observed topological defects

We start with the description of topological defects observed in the modeling. Figure 8 shows the half-vortex enforced by the termination of the line of amplitude solitons (the unabridged version of the image is shown in Fig. 11d). Traveling counterclockwise around the core over the complete $2\pi$-revolution, the arrow (corresponding to the XY degree of freedom $\theta$) rotates by $\pi$ in the clockwise direction – in some similarity to the field of directors in nematic liquid crystals.

Figure 9 shows a configuration with 4 solitons (i.e. $s\rightarrow-s$ kinks breaking the interchain order) which are “dressed” by the phase adjustments $\theta\rightarrow\theta+\pi$, restoring the interchain order. Here we observe two neutral half-VA pairs (outer ones) and also two charged half vortex-vortex and antivortex-antivortex pairs (inner ones). The presence of such configurations with combined phase-amplitude topological defects grows with the increase of $A_{\perp}$. The sharper versions of such combined defects can be also found in a great number in Fig. 11.

Figure 10 shows several snapshots of the MC evolution of the system. From (a) to (b) we see the attraction of vortex-vortex to antivortex-antivortex pair and their annihilation with the formation of a segment of solitonic wall. From (b) to (c,d) we observe how the long bisolitonic pair shrinks, presumably because of the weak dipole-dipole interaction of the half-VA pairs at its ends. The final configuration (d) constitutes a segment of a bisolitonic wall, which no longer possesses a topological defect in $\theta$, with all the energy loss $\sim A_{\perp}$ being concentrated in two transverse $s$-misfits (coming from the term $-A_{\perp}s_{r}s_{r+\hat{y}}\cos(\theta_{r}-\theta_{r+\hat{y}})$ in (5)). The tendency to formation of such bisoliton walls becomes more pronounced with decreasing of $A_{\perp}$ (see also Fig. 13).

IV.2 Crossover and finite-size phase transitions in 2D

Here we study in detail the limiting cases of strong ($A_{\perp}=4A_{||}$) and weak ($A_{\perp}=A_{||}/4$) interchain couplings, and also the intermediate symmetric case ($A_{\perp}=A_{||}$). In the thermodynamic limit, the properties of these three cases qualitatively match each other, but for finite systems there is a difference among them as described below. For all three cases we take the same value of $A=\sqrt{A_{\perp}A_{||}}$. Concentration of solitons is taken to be $\nu=0.05$ (per site) throughout this subsection. The size $L\times H$ of the studied system is $25\times 100$; periodic boundary conditions are imposed. The system is cooled down from $T=1.5A>T_{BKT}$ to $T=0.01A$.

IV.2.1 Strong interchain coupling: $A_{\perp}=2A$, $A_{||}=A/2$.

We start from the strong interchain coupling: $A_{\perp}=4A_{||}$. Cooling the system down starting from the high temperatures, first we observe only short transverse solitonic rods (Fig. 11a). At temperature $T_{2}\approx 0.3A$ we observe a sudden change: one of the solitonic rods starts to grow until it reaches the borders of the sample (Fig. 11b). With further cooling more solitonic walls appear (Fig. 11c) until the material for their building – free solitons – is exhausted. All the walls are installed through the grows of solitonic rods.

The transition can also be detected using integrated characteristics. Fig. 12 shows the dependence of the average number of solitons’ transverse neighbors versus temperature. With lowering the temperature, the condensation of the first wall occurs at $T\simeq 0.3A$. With further cooling, more walls condense sequentially giving rise to the step-like features on the plot.

IV.2.2 Weak interchain coupling: $A_{\perp}=A/2$, $A_{||}=2A$.

Here we consider the case of weak interchain coupling $A_{\perp}=0.25A_{||}$. The walls’ formation here is similar to the previously studied case of the pure Ising model Bohr:1983 ; Teber:2001 ; Karpov:2016 . At $T\gg A_{\perp}$ different chains are only weakly correlated with each other. For $T\sim A_{\perp}\ll A_{||}$ the in-chain variation of the phase is strongly suppressed, so the phase degree of freedom cannot any more “cure” the intrachain misfit created by the soliton, thus a (pseudo-) confinement force between solitons is established, which binds them into pairs. With further cooling, bisolitons aggregate into growing rods which finally cross the whole sample; after that the paired bisolitonic walls can freely diverge to separated solitonic ones (Fig. 13).

Fig. 14 shows the dependence of the average number of solitons’ transverse neighbors versus temperature: after the first bisolitonic rod crossing the sample, the fast growth of the number of neighbors with lowering the temperature (starting at $T\sim 0.4A$) finally saturates to almost flat behaviour at $T\approx 0.25A$, since the system runs out of the new building blocks (free bisolitons) for the further growth of the remaining rods.

IV.2.3 Intermediate interchain coupling: $A_{\perp}=A$, $A_{||}=A$.

The symmetric case $A_{\perp}=A_{||}$ is somewhat intermediate between the cases of strong and weak interchain couplings and combines some traits of both.

Unlike the case of strong interchain coupling, there are stable solitonic rods with transverse lengths down to the smallest size $l=1$ dressed by half-VA pairs, thus solitonic rods elongate gradually starting from the smallest ones. When $T$ approaches $T_{2}\approx 0.41A$ from above, there is already a considerable amount of solitonic rods with $l>1$ (Fig. 15a). First walls grow through solitonic rods (Fig. 15a), but then the system starts to run out of the solitons, hence the screening of vortices becomes inefficient and solitonic rods recombine to bisolitonic ones, so the rest of the walls grows through the bisolitonic rod mechanism (Fig. 15b). The temperature dependence of the number of the transverse neighbors of solitons (Fig. 16) is qualitatively similar to the case of the weak interchain coupling (Fig. 12), which reflects the gradual rods’ elongation.

We conclude this subsection by noting that all the cases of weak, intermediate, and strong couplings probably match each other in the thermodynamic limit. Supplementary Figure 2 (Appendix B) shows the results of the simulations for the system $100\times 100$ for the case of strong coupling. We observe that solitonic rods can be favorable up to some maximal length, however the domain walls grow through the bisolitonic mechanism. This suggests that even in the case of strong interchain coupling in 2D the walls’ growth proceeds through a crossover and not a phase transition.

IV.3 Phase transitions in 3D

In this subsection we describe results of the numerical modeling for the three-dimensional case. Here we consider the system with the size $10\times 10\times 100$ with concentration of solitons $\nu=0.06$ and the symmetric coupling $A_{\perp}=A_{||}\equiv A$. System is cooled from $T=3A$ down to $T=0.02A$.

In the three-dimensional case, unlike 2D, there is a true long-range ordering phase transition happening at $T_{1}\approx 2A$ (see Fig. 17; the small concentration of solitons only slightly shifts down the 3D XY model transition temperature $T_{XY}^{3D}\approx 2.2A$ [Pereira:1995, ]). With lowering the temperature, the first solitonic wall grows across the sample at $T\approx 0.24A$. For the periodic boundary conditions used in the simulation, actually, a pair of solitonic walls has to grow – in Fig. 18b we can see the fully grown wall at the left side of the figure and an emerging wall at the right side. The system decomposes to domains with alternating values of the order parameter $\Delta=\sum s_{i}e^{i\theta_{i}}$. In a macroscopic sample we expect an equal amount of $+\Delta$ and $-\Delta$ domains, so the combined XY-Ising order parameter $\Delta$ should vanish; for a finite system it just drops to a smaller value. Still, the system possesses a non-zero value of the nematic order parameter $\langle\cos(2\theta)\rangle$. With lowering the temperature more walls are installed through the bisolitonic mechanism – Fig. 18c.

Temperature dependence of the average number of transverse solitons’ neighbors (Fig. 19) indicates that there is a pronounced start of the transverse aggregation of solitons at $T\approx 0.26A$, which eventually leads to solitonic walls’ formation.

IV.4 Quench

In this subsection we analyze the scenario of a sudden quench from high to low temperatures $T<T_{2}$. Such a scenario is relevant to photo- and electrical pulse induced phase transitions PIPT-book:2004 ; impact:12 and in particular to formation of domain walls in CDW systems Stojchevska:2014 ; Yoshida:2014 ; Vaskivskyi:2016 ; Ma:2016 ; Cho:2016 ; Svetin:2016 ; Gerasimenko:2017 ; Karpov:2018 . We show that even for the case of weak interchain coupling the solitonic mechanism of domain walls growth can be important.

Figure 20 shows configuration of the system after a sudden quench from a high temperature $T>T_{2}$ to low temperature $T=0.25A<T_{2}$ for a system with weak interchain coupling $A_{\perp}=A_{||}/4$. We observe that the domain wall was created by the solitonic mechanism. When the solitonic rod initially starts to grow, half-vortices with opposite charges (vorticities) are created at its ends. These charges attract smaller solitonic rods via logarithmic interaction mediated by the XY subsystem, promoting the further growth of the solitonic rods. On the contrary, the growth of bisolitonic rods goes via attaching of new bisolitons at the ends, which is essentially a much slower random walk process. At the same time, longer solitonic rods become low-mobile, therefore their recombination to bisolitonic rods is kinetically suppressed.

V Discussion

In the course of the cooling the solitonic lines (“rods”) are growing in two stages. The first stage lasts until the rods reach some maximal characteristic size $l=l^{*}\sim\exp(\frac{4}{\pi}\sqrt{A_{\perp}/A_{||}})$ when the bisolitonic rods become energetically favorable with respect to the solitonic ones (or until the rods’ lengths become comparable to the distance between them, which happens at $l\sim 1/\nu$). At this first stage of the growth the energy gain from merging of two rods is $\sim A\ln l$, so the growth process is self-accelerating and for finite samples the solitonic rods can even grow across the whole sample. For big enough samples with the transverse size $H>l^{*}$, the single-soliton rods become eventually attracted to each other and glue into bisolitonic rods. Then the system enters the second stage of the growth, when bisolitonic rods continue to gradually elongate, but the energy gain from merging of two bisolitonic rods is only $\sim A$. So at the second stage the rods grow much more “slowly” with lowering the temperature. The bisolitonic rods grow to infinity as $T\rightarrow 0$, so the second stage is similar to the previously studied case of the pure Ising model Bohr:1983 ; Teber:2001 ; Karpov:2016 . For finite systems there should be some transverse-size-dependent temperature $T_{F}(H)$, at which either solitonic or bisolitonic walls grow across the whole sample.

An important role is played here by the combined amplitude-phase topological defects. The interchain order broken by an amplitude soliton $s\rightarrow-s$ (with a small core size) can be cured over the long distances by the phase degree of freedom $\theta\rightarrow\theta+\pi$. This allows the existence of single solitons at any non-zero temperature (without confinement to bisolitonic pairs), hence there is no Ising order: $\langle s\rangle=0$.

VI Conclusions

The motivation of presented studies was to notice that most common types of symmetry breaking (superconductivity, spin density waves, antiferromagnetic Mott state, incommensurate charge density waves) in quasi-one-dimensional electronic systems possess a combined manifold of degenerate states. Beyond the standard continuous XY-type degeneracy with respect to the phase $θ$ degree of freedom (the Goldstone mode), there is also an Ising type discrete degeneracy with respect to the sign of the amplitude $A$ of the order parameter $A\exp(i\theta)$. These degrees of freedom can be controlled or accessed independently via either the spin polarization or via the charge or current densities. The degeneracies give rise to two coexisting types of topologically nontrivial configurations: phase vortices or amplitude kinks – the solitons. Being decoupled at the 1D level of isolated chains, the two degrees of freedom experience confinement which may be enforced or screened depending on the concentration of solitons (given by charge or spin decompensation) and on the temperature. In 2D, 3D states with long-range (or quasi-long-range BKT type) orders, the topological confinement takes place which binds together kinks and half-integer vortices. These combined amplitude-phase solitons are the lowest-energy excitations of either spin or charge taking that role from conventional electrons.

With lowering the temperature at a finite concentration of solitons (controlled by the spin polarization or the charge doping) the combined solitonic complexes experience further aggregation resulting in a pattern formation in the course of passing through temperatures of specific phase transitions or crossovers. To understand these transformations, we have given an analytical treatment of two-fields system supported by Monte Carlo simulations for a generic coarse-grained model. In some similarity Bohr:1983 ; Teber:2001 ; Karpov:2016 with systems possessing only the discrete degree of freedom (CDW of a bond or site dimerization type, charge ordering) here there are also two phase transitions $T_{1}$, $T_{2}$. At the higher $T_{1}$ the (quasi) long-range order of the XY type sets in leading to the confinement which dresses the kinks by half-integer vortices. At the lower $T_{2}$ the liquid of combined solitons starts to structure by aggregating the solitons into walls growing across the sample. In a 2D system, $T_{2}$ is only a crossover temperature below which the walls’ lengths grow; but for finite samples the modeling yields sequential events of walls crossing the whole sample separating it into order parameter alternating domains. Depending on the sample width and the degree of the interchain coupling, that proceeds via either solitonic or bi-solitonic domain walls. The growing solitonic walls are formed as rods of transversely correlated amplitude kinks, with rods’ termination points being ultimately accomplished by half-integer vortices. Attractions of vortices from different walls may force the walls to glue together forming a topologically trivial bi-solitonic wall. Estimations and modeling results for the 3D system indicate that $T_{2}$ is a true phase transition.

For non-equilibrium processes, like supercooling the system with a sudden quench, the solitonic-wall scenario can be drastically promoted, because a growing solitonic rod has uncompensated half-vortex “charges” at the ends, which attract mobile smaller rods, while the slow process of gluing of low-mobile long solitonic rods into bisolitonic ones can be kinetically suppressed.

In view of recent successes in real space observations and even manipulations Stojchevska:2014 ; Yoshida:2014 ; Vaskivskyi:2016 ; Ma:2016 ; Cho:2016 ; Svetin:2016 ; Gerasimenko:2017 ; Karpov:2018 of domain walls in correlated electronic systems, we are calling for the attention to described here intriguing opportunities.

Acknowledgements

We acknowledge the financial support of the Ministry of Education and Science of the Russian Federation in the framework of Increase Competitiveness Program of NUST MISiS (N K2-2017-085). SB acknowledges funding from the ERC AdG 320602 “Trajectory”.

Appendix A Transformation from XY-picture to Coulomb gas

Consider the simplified version of the Hamiltonian (5) with $J_{||}=0$ and $A_{||}=A_{\perp}\equiv A$ when the Hamiltonian acquires the form

[TABLE]

The partition function of the system is ($K=\beta A$):

[TABLE]

We introduce “solitonic” bond variables $\sigma$ showing whether there is a flip of spins along a given bond in $y$-direction ($\sigma_{r,r+\hat{y}}=1$ if there are different spins at sites $r$ and $r+\hat{y}$ and $\sigma_{r,r+\hat{y}}=0$ if they are identical), and which are identically zero for $x$-direction

[TABLE]

Denote also $\prod_{r}\int_{-\pi}^{\pi}d\theta_{r}/2\pi=\int d[\theta]$, and the index $\hat{\epsilon}=\hat{x},\hat{y}$, then

[TABLE]

We shall employ the Villain approximation Villain:1975 : $\exp(K\cos\theta)\approx\sum_{n}\exp(-K(\theta-2\pi n)^{2}/2)$ and spread the integration over $\theta$ from ${-\infty}$ to ${\infty}$. In this approximation spin-wave and vortical degrees of freedom are completely decoupled.

[TABLE]

Now we make dual transformations, adapting the general scheme outlined in Jose:1977 ; Chaikin:book ; Lee:1985 :

1. Make the Fourier transform for variables $\theta$ which gives us discrete gaussian variables $h$.

Using the identity

[TABLE]

we arrive at

[TABLE]

2. Next we integrate over $\theta$ which yield the factor $\prod_{r}\delta(h_{rx}-h_{r-\hat{x},x}+h_{ry}-h_{r-\hat{y},y})$

3. Introduce integer variables $h_{R}$ on the dual lattice in order to automatically satisfy the above constraint (Fig. A):

[TABLE]

Denoting $\sigma_{rx}\equiv\sigma_{Ry}\equiv 0$, $\sigma_{ry}\equiv\sigma_{Rx}$, this yields to

[TABLE]

If $\sigma\equiv 0$ we get the so-called “discrete gaussian model”, which is a particular case of SOS model (with a quadratic potential). Here the terms with $\sigma\neq 0$ contribute to the partition function with “negative weights” – a similar approach was explicitly used in Korshunov:1986 .

4. Using the Poisson summation formulaKadanoff:book

[TABLE]

we trade the $h$-summation to $m$-summation and $\phi$-integration; then rescaling the field $\phi_{new}=\phi/K$ we obtain

[TABLE]

5. Finally we can integrate out Gaussian variables $\phi$ and arrive at the Coulomb gas representation.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $a$ is the cutoff of the order of the lattice spacing. Note that the Ising variables $\sigma$ do not appear in the matrix $A$; they enter only the “current” term

[TABLE]

Here we denoted half-integer contributions to vorticity as $\nu_{R}=(\sigma_{R,x}-\sigma_{R+\hat{x},x}+\sigma_{R,y}-\sigma_{R+\hat{y},y})/2$. We see that the standard XY-vorticity is just shifted by $\nu_{R}=0,\pm 1/2$, which defines the Coulomb gas with a mixture of integer and half-integer charges. Completing the Gaussian integration over $\phi$ and substituting back $K=\beta A$ we obtain

[TABLE]

Here $y_{0}=\exp(-\pi^{2}\beta A/2)$, the sum $\sum_{\{m_{R}\}}$ contains only neutral configurations $\sum_{R}(m_{R}+\nu_{R})=0$, the sum $\sum_{\{\nu(s)\}}$ counts all possible configurations of spins $\{s\}$ with corresponding values of $\nu(\{s\})$ (recall the definition of $\sigma_{ry}$ (26))

[TABLE]

So, we have constructed a dual model with two types of vorticity variables $m_{R}$ and $\nu_{R}$, making the vorticities effectively half-integer. We use the partition function (33) in Sec. III for constructing the strong coupling part of the phase diagram corresponding to $J_{||}=0$.

Appendix B Modeling for 2D system $100\times 100$ with strong coupling $A_{\perp}=4A_{||}$