text:
# Classical technical analysis of Latin American market indices. Correlations in Latin American currencies (𝐴𝑅𝑆, 𝐶𝐿𝑃, 𝑀𝑋𝑃) exchange rates with respect to 𝐷𝐸𝑀, 𝐺𝐵𝑃, 𝐽𝑃𝑌 and 𝑈𝑆𝐷11footnote 1Happy Birthday, Dietrich; by now you should be rich !

## 1 Introduction

The buoyancy of the US dollar is a reproach to stagnant Japan, recessing Europe economy and troubled developing countries like Brazil or Argentina. Econophysics aims at introducing statistical physics techniques and physics models in order to improve the understanding of financial and economic matters. Thus when this understanding is established, econophysics might $`later`$ help in the well being of humanity. In so doing several techniques have been developed to analyze the correlations of the fluctuations of stocks or currency exchange rates. It is of interest to examine cases pertaining to rich or developing economies.

In the first sections of this report we recall the classical technical analysis methods of stock evolution. We recall the notion of moving averages and (classical) momentum. The case of IBM and Latin American market indices serve as illustrations.

In 1969 the International Monetary Fund created the special drawing rights $`SDR`$, an artificial currency defined as a basket of national currencies $`DEM`$, $`FRF`$, $`USD`$, $`GBP`$ and $`JPY`$. The $`SDR`$ is used as an international reserve asset, to supplement members existing reserve assets (official holdings of gold, foreign exchange, and reserve positions in the IMF). The $`SDR`$ is the IMF’s unit of account. Four countries maintain a currency peg against the $`SDR`$. Some private financial instruments are also denominated in $`SDR`$s. Because of the close connections between the developing countries and the IMF, we search for correlations between the fluctuations of $`ARS`$, $`CLP`$ and $`MXP`$ exchange rates with respect to $`SDR`$ and the currencies that form this artificial money. In the latest sections of this report, we compare the correlations of such fluctuations as we did in our previous results on $`EUR`$ exchange rates fluctuations with respect to $`USD`$, $`GBP`$ and $`JPY`$ .

## 2 Technical Analysis: IBM and Latin America Markets

Technical indicators as moving average and $`momentum`$ are part of the classical technical analysis and much used in efforts to predict market movements . One question is whether these techniques provide adequate ways to read the trends.

Consider a time series $`x(t)`$ given at N discrete times $`t`$. The series (or signal) moving average $`M_\tau (t)`$ over a time interval $`\tau `$ is defined as

$$M_\tau (t)=\frac{1}{\tau }\underset{i=t}{\overset{t+\tau 1}{}}x(i\tau )t=\tau +1,\mathrm{},N$$
(1)
i.e. the average of $`x`$ over the last $`\tau `$ data points. One can easily show that if the signal $`x(t)`$ increases (decreases) with time, $`M_\tau (t)<x(t)`$ ($`M_\tau (t)>x(t)`$). Thus, the moving average captures the trend of the signal given the period of time $`\tau `$. The IBM daily closing value price signal between Jan 01, 1990 and Dec 31, 2000 is shown in Fig. 1 (top figure) together with Yahoo moving average taken for $`\tau =50`$ days . The bottom figure shows the daily $`volume`$ in millions.

There can be as many trends as moving averages as $`\tau `$ intervals. The shorter the $`\tau `$ interval the more sensitive the moving average. However, a too short moving average may give false messages about the long time trend of the signal. In Fig. 2(a) two moving averages of the IBM signal for $`\tau `$=5 days (i.e. 1 week) and 21 days (i.e. 1 month) are compared for illustration.

The intersections of the price signal with a moving average can define so-called lines of resistance or support . A line of resistance is observed when the local maximum of the price signal $`x(t)`$ crosses a moving average $`M_\tau (t)`$. A support line is defined if the local minimum of $`x(t)`$ crosses $`M_\tau (t)`$. In Fig. 2(a) lines of resistance happen around May 1993 and lines of support around Sept 1993. Support levels indicate the price where the majority of investors believe that prices will move higher, and resistance levels indicate the price at which a majority of investors feel prices will move lower. Other features of the trends are the intersections between $`two`$ moving averages $`M_{\tau _1}`$ and $`M_{\tau _2}`$ which are usually due to drastic changes in the trend of $`x(t)`$ . Consider two moving averages of IBM price signal for $`\tau _1=5`$ days and $`\tau _2=21`$ days (Fig. 2(a)). If $`x(t)`$ increases for a long period of time before decreasing rapidly, $`M_{\tau _1}`$ will cross $`M_{\tau _2}`$ from above. This event is called a ”death cross” in empirical finance . In contrast, when $`M_{\tau _1}`$ crosses $`M_{\tau _2}`$ from below, the crossing point coincide with an upsurge of the signal $`x(t)`$. This event is called a ”gold cross”. Therefore, it is of interest to study the density of crossing points between two moving averages as a function of the size difference of the $`\tau `$’s defining the moving averages. Based on this idea, a new and efficient approach has been suggested in Ref. in order to estimate an exponent that characterizes the roughness of a signal.

The so called $`momentum`$ is another instrument of the technical analysis and we will refer to it here as the classical momentum, in contrast to the generalized momentum . The classical momentum of a stock is defined over a time interval $`\tau `$ as

$$R_\tau (t)=\frac{x(t)x(t\tau )}{\tau }=\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}t=\tau +1,\mathrm{},N$$
(2)
The momentum $`R_\tau `$ for three time intervals, $`\tau =5,21`$ and 250 days, i.e. one week, one month and one year, are shown in Fig. 2(b) for IBM. The longer the period the smoother the momentum signal. Much information on the price trend turns is usually considered to be found in a moving average of the momentum

$$R_\tau ^\mathrm{\Sigma }(t)=\underset{i=t}{\overset{t+\tau 1}{}}\frac{x(i)x(i\tau )}{\tau }t=\tau +1,\mathrm{},N$$
(3)
Moving averages of the classical momentum over 1 week, 1 month and 1 year for the IBM price difference over the same time intervals, $`R_\tau ^\mathrm{\Sigma }`$ are shown in Fig. 2(c). In Fig. 2(d) the IBM signal and its weekly (short-term), monthly (medium-term) and yearly (long-term) moving averages are compared to the weekly (short-term), monthly (medium-term) and yearly (long-term) momentum indicators in order to better observe the bullish and bearish trends in 1999.

The message that is coming out of reading the combination of the these six indicators states that one could start buying at the momentum bottom, as it is for both monthly and weekly momentum indicators around mid February 1999 and buy the rest of the position when the price confirms the momentum uptrend and rises above the monthly moving average which is around March 1999. The first selling signal is given during the second half of July 1999 by the death cross between short and medium term moving averages and by the maximum of the monthly momentum, which indicates the start of a selling. At the beginning of October 1999, occurs the maximum of the long-term momentum. It is recommended that one can sell the rest of the position since the price is falling down below the moving average. Hence, it is said that momentum indicators lead the price trend. They give signals before the price trend turns over.

Along the lines of the above for IBM, we analyze three Latin America financial indices, Argentinian MerVal, Brazilian Bovespa and Mexican IPC (Indice de Precios y Cotizaciones) applying the moving averages and the classical momentum concepts. In Fig. 3 the time evolution of the MerVal stock over the time interval mid 1996 \- mid 2001, is plotted with a simple moving average for $`\tau =`$50 days showing the medium range trend of the price. The moving average and classical momentum of the MerVal stock prices for time horizons equal to one week, one month and one year are shown Fig. 3.

The cases of the Brazilian Bovespa and Mexican IPC (Indice de Precios y Cotizaciones) are shown in Figs. 4 and 5 respectively.

## 3 Exchange Rates and Special Drawing Rights

In 1969 the International Monetary Fund (IMF) created the Special Drawing Rights (ticker symbol $`SDR`$; sometimes $`XDR`$ is also used), an artificial currency unit defined a basket of national currencies. The $`SDR`$ is used as an international reserve asset, to supplement members’ existing reserve assets (official holdings of gold, foreign exchange, and reserve positions in the IMF). The $`SDR`$ is the IMF’s unit of account: IMF voting shares and loans are all denominated in $`SDR`$s. The $`SDR`$ serves as the unit of account for a number of other international organizations, including the WB. Four countries maintain a currency peg against the $`SDR`$. Some private financial instruments are also denominated in $`SDR`$s.

The basket is reviewed every five years to ensure that the currencies included in the basket are representative of those used in international transactions and that the weights assigned to the currencies reflect their relative importance in the world’s trading and financial systems. Following the completion of the most recent regular review of $`SDR`$ valuation on October 11, 2000, the IMF’s Executive Board agreed on changes in the method of valuation of the $`SDR`$ and the determination of the $`SDR`$ interest rate, effective Jan. 01, 2001.

The $`SDR`$ artificial currency can be represented as an weighted sum of the five currencies $`C_i`$, $`i=1,5`$:

$$1SDR=\underset{i=1}{\overset{5}{}}\gamma _iC_i$$
(4)
where $`\gamma _i`$ are the currencies weights in percentage (Table 1) and $`C_i`$ denote the respective currencies, U.S. Dollar ($`USD`$), German Mark ($`DEM`$), French Franc ($`FRF`$), Japanese Yen ($`JPY`$), British Pound ($`GBP`$).

On January 1, 1999, the German Mark and French Franc in the $`SDR`$ basket were replaced by equivalent amounts of $`EUR`$. The relevant exchange rates (ExR) are shown in Figs. 6-8.

### 3.1 The DFA technique

The DFA technique  is often used to study the correlations in the fluctuations of stochastic time series like the currency exchange rates. Recall briefly that the DFA technique consists in dividing a time series $`y(t)`$ of length $`N`$ into $`N/\tau `$ nonoverlapping boxes (called also windows), each containing $`\tau `$ points . The local trend $`z(n)`$ in each box is defined to be the ordinate of a linear least-square fit of the data points in that box. The detrended fluctuation function $`F^2(\tau )`$ is then calculated following:

$$F^2(\tau )=\frac{1}{\tau }\underset{n=k\tau +1}{\overset{(k+1)\tau }{}}\left[y(n)z(n)\right]^2k=0,1,2,\mathrm{},\left(\frac{N}{\tau }1\right)$$
(5)
Averaging $`F^2(\tau )`$ over the $`N/\tau `$ intervals gives the mean-square fluctuations

$$\varphi (\tau )=<F^2(\tau )>^{1/2}\tau ^\alpha $$
(6)
The exponent $`\alpha `$ value implies the existence or not of long-range correlations, and is assumed to be identical to the Hurst exponent when the data is stationary. Moreover, $`\alpha `$ is an accurate measure of the most characteristic (maximum) dimension of a multifractal process . Since only the slopes and scaling ranges are of interest the various DFA-functions $`\varphi (\tau )`$ have been arbitrarily displaced for readability in Figs. 9-11. The $`\alpha `$ values are summarized in Table 2. It can be noted that the scaling ranges are usually from 5 days till 170 days for $`ARS`$ exchange rates, from 5 days to about 1 year for $`MXP`$ and $`CLP`$ exchange rates, with the exponent $`\alpha `$ close to 0.5 in that range. Crossover at 80 days from Brownian like to persistent correlations is obtained for $`CLP/JPY`$ and $`CLP/SDR`$.

### 3.2 Local scaling with DFA and Intercorrelations between fluctuations

The time derivative of $`\alpha `$ can usually be correlated to an entropy production  through market information exchanges. As done elsewhere, in order to probe the existence of locally correlated and decorrelated sequences, we have constructed an observation box, i.e. a 500 days wide window probe placed at the beginning of the data, calculated $`\alpha `$ for the data in that box. Moving this box by one day toward the right along the signal sequence and again calculating $`\alpha `$, a local $`\alpha `$ exponent is found (but not displayed here).

We eliminate the time between these data sets and construct a graphical correlation matrix of the time-dependent $`\alpha `$ exponents for the various exchange rates of interest (Fig.12-14). We show $`\alpha _{C_i/B_j}`$ vs. another $`\alpha _{C_i/B_j}`$, where a $`C_i`$ is a developing country currency while $`B_j`$ is $`USD`$, $`GBP`$, $`JPY`$, $`DEM`$, and $`SDR`$ for the available data. In so doing a so-called correlation matrix is displayed for the time interval of interest. Such bilateral correlations between different $`\alpha `$ exponents can be considered in order to estimate the strength and some nature of the correlations. As described elsewhere , such a correlation diagram can be divided into main sectors through a horizontal, a vertical and diagonal lines crossing at (0.5,0.5). If the correlation is strong the cloud of points should fall along the slope $`=+1`$ line. If there is no correlation the cloud should be rather symmetrical. The lack of symmetry of the plots and wide spreading of points outside expected clouds (see e.g. Figs. 13(c), 14(b,c) - mainly containing $`MXP`$) indicate highly speculative situations. Notice the marked imbalance of some plots, mainly involving $`ARS`$. It is fair to say that other techniques are also of great interest to observe correlations between financial markets .

## 4 Conclusions

The classical technical analysis methods of financial indices, stocks, futures, … are very puzzling. We have recalled them. Illustrations have used the IBM share price and Latin American financial indices. We have used the DFA method to search for scaling ranges and type of behavior of exchange rates between Latin American currencies ($`ARS`$, $`CLP`$, $`MXP`$) and other major currencies $`DEM`$, $`GBP`$, $`JPY`$ and $`USD`$, including $`SDR`$s. In all cases persistent to Brownian like behavior is obtained for scaling ranges from a week to about one year, with an exception of $`CLP/JPY`$ and $`CLP/SDR`$ for which there is a transition from Brownian like to persistent correlations with $`\alpha =0.70`$ and $`\alpha =0.74`$ for scaling ranges longer than 80 days. We have also sorted out to correlations and anticorrelations of such exchange rates with respect to currencies as $`DEM`$, $`GBP`$, $`JPY`$ and $`USD`$. They indicate a very complex or speculative behavior.

Acknowledgements

MA thanks to the organizers of the Stauffer 60th birthday symposium for their invitation and kind welcome.


format: text

text:
# Hardness of approximating the weight enumerator of a binary linear code

## 1 Matrix element problem and proof of Theorem 1

We choose the basis $`B=\{CNOT,H,T\}`$ for quantum circuits, where

$$CNOT:|a,b|a,ab;H=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right);T=\left(\begin{array}{cc}1& 0\\ 0& e^{i\pi /4}\end{array}\right).$$
(4)
It is well-known that this basis provides universal quantum computation. Solovay–Kitaev theorem guarantees the fast and effective approximation of any unitary operator by the operators in this basis.

Matrix element problem. Given a description of a quantum circuit in the basis $`B`$ and accuracy threshold $`\epsilon `$ find out the approximation $`\stackrel{~}{u}`$ of a matrix element $`0|U|0`$ with accuracy $`\epsilon `$: $`|\stackrel{~}{u}0|U|0|<\epsilon `$. Here $`U`$ stands for an operator realizing by the circuit.

Let’s denote an oracle answering the matrix element problem by $`ME`$.

###### Theorem 2.

$`\mathrm{P}^{ME}=\mathrm{P}^{\mathrm{GapP}}`$.

To prove Theorems 1 and 2 we construct three Turing reductions: a $`\mathrm{GapP}`$ oracle to the oracle $`ME`$ oracle; then the $`ME`$ oracle to the $`WE`$ oracle; then the $`WE`$ oracle to a $`\mathrm{GapP}`$ oracle.

### 1.1 Reduction an $`\mathrm{𝐆𝐚𝐩𝐏}`$-oracle to the $`𝑴𝑬`$ oracle

Let $`f`$ be a $`\mathrm{GapP}`$-function. We assume that $`f`$ is computed by a counting machine in the normal form and the machine makes only nondeterministic moves with branching degree 2. Lemma 4.3 from  guarantees that we do not loss the generality: the gaps computing by machines in normal form and general ones differ by factor 2.

Reduction of $`f`$-oracle to the $`ME`$ oracle reproduce the argument used in . Let $`q`$ be the number of nondeterministic moves along a computation path. The number does not depend on a choice of computation path due to the normal form condition. So, we may consider counting machine described above as a deterministic Turing machine $`M`$ supplied with an additional input: *a guess string* $`r\{0,1\}^q`$. Bits of the guess string determine the choices for the nondeterministic moves. This machine produces the same gap $`g_M(x)`$ as the initial one.

Now we transform the machine $`M`$ to a Boolean circuit $`C_M`$ in the basis $`\{,\}`$ computing 1 for pairs (input,guess) finishing at the rejecting state and computing 0 for pairs (input,guess) finishing at the accepting state. The size of the circuit $`C_M`$ is quadratic on running time of the $`M`$. For fixed input string $`x`$ the gap produced by $`M`$ is the difference between the number of guess strings $`u`$ that give the circuit output 0 and the number of guess strings $`u`$ that give the circuit output 1.

Hereinafter we will simplify notation and will use $`+`$ and $``$ instead of $``$ and $``$ respectively.

Let $`s`$ be the size of the circuit $`C_M`$. By $`z_k^{(j)}`$, $`1ks`$, we denote auxiliary variables of the $`C_M`$. We also assume that the circuit output is given by the value of the variable $`z_s^{(j)}`$. Each assignment in a circuit has the form $`z_k^{(j)}:=ab`$ where $`\{+,\}`$ and $`a,b`$ are either input or auxiliary variables. This assignment corresponds the equation $`Z_k^{(j)}=z_k^{(j)}+ab=0`$. Note that the values of input variables and the guess string ($`x,u`$) determine the values of all auxiliary variables. So, for each $`x`$ fixed the number of solutions of the system of equations $`Z_k^{(j)}=0`$, $`1k<s`$, $`z_s^{(j)}=\zeta `$ equals the number of rejecting computation paths of the counting machine $`M`$ for $`\zeta =1`$ and the number of accepting paths for $`\zeta =0`$.

For a polynomial $`p`$ over the two-element field $`𝔽_2`$ we denote by $`\mathrm{\Delta }p`$ the number

$$\mathrm{\Delta }p=\underset{x}{}(1)^{p(x)}.$$
(5)
Now let’s introduce polynomials over $`𝔽_2`$:

$$\begin{array}{cc}& F_x^\zeta (u,z,v)=\underset{k=1}{\overset{s1}{}}v_kZ_k^{(j)}+v_0(z_s^{(j)}+\zeta ),F_x=(1+w)F_x^{(0)}+w(1+F_x^{(1)}),\hfill \\ & \zeta \{0,1\},v\{0,1\}^s,w\{0,1\}.\hfill \end{array}$$
(6)
Let us check that

$$\mathrm{\Delta }F_x(u,z,v,w)=\mathrm{\Delta }F_x^0(u,z,v)\mathrm{\Delta }F_x^1(u,z,v)=2^sg_m(x).$$
(7)
The polynomials $`F_x^\zeta `$ are linear in variables $`v_j`$. So, if values $`x,u,z`$ does not satisfy a system $`Z_k^{(j)}=0`$, $`1k<s`$, $`z_s^{(j)}=\zeta `$ then all possible assignments to variables $`v`$ contribute the zero term to the sum (5) for the polynomial $`F_x^\zeta `$. If values $`x,u,z`$ satisfy the system then $`F_x^\zeta =0`$ for any assignment of $`v`$. Thus any accepting computation path of $`M`$ contributes $`2^s`$ to $`\mathrm{\Delta }F_x^0(u,z,v)`$ and any rejecting path contributes $`2^s`$ to $`\mathrm{\Delta }F_x^1(u,z,v)`$. The first equation in (7) follows immediately from the definition (5).

From the construction described above it is clear that function $`xF_x`$ is computable in polynomial time.

Now we relate the value $`\mathrm{\Delta }F_x`$ to a value $`0|U(F_x)|0`$ for some unitary operator acting in the space $`^N`$ of $`N`$ qubits where $`N`$ is the total number of variables in the polynomial $`F_x`$.

The operator $`S_J=\mathrm{\Lambda }^J(1)`$ (controlled phase shift) corresponds to the monomial $`x_J=_{jJ}x_j`$ of $`F_x`$. Controlled phase shift is defined as

$$\{\begin{array}{cc}\mathrm{\Lambda }^J(1)|x_1\mathrm{}x_n=|x_1\mathrm{}x_n,\hfill & x_J=1,\hfill \\ \mathrm{\Lambda }^J(1)|x_1\mathrm{}x_n=|x_1\mathrm{}x_n,\hfill & \text{otherwise}.\hfill \end{array}$$
(8)
Let

$$U(F_x)=\underset{j=1}{\overset{n}{}}H[j]\underset{JM(F_x)}{}S_J\underset{j=1}{\overset{n}{}}H[j]$$
(9)
where $`M(F_x)`$ is the set of the monomials of $`F_x`$.

We have

$$\begin{array}{c}0|U(F_x)|0=\hfill \\ \hfill =\frac{1}{2^N}\underset{u_1,\mathrm{},u_N}{}u_1,\mathrm{},u_N|\underset{JM(f)}{}S_J|u_1,\mathrm{},u_N=\frac{1}{2^N}\underset{u_1,\mathrm{},u_N}{}(1)^{_{JM(f)}u_J}=\\ \hfill =\frac{1}{2^N}\underset{u_1,\mathrm{},u_N}{}(1)^{F_x(u)}=2^N\mathrm{\Delta }F_x.\end{array}$$
(10)
Now, we approximate each factor in (9) $`U(f)`$ with accuracy $`\delta =2^{N2}/M`$ by a quantum circuit in the basis $`B`$. Here $`M=poly(n)`$ denotes the number of factors in (9). A matrix element $`0|\stackrel{~}{U}|0`$ of a product $`\stackrel{~}{U}`$ of these approximations is very close to $`0|U(F_x)|0`$:

$$|0|\stackrel{~}{U}U(F_x)|0|<2^{N2}.$$
(11)
So, we may compute the exact value of $`\mathrm{\Delta }F_x`$ if an approximate evaluation of $`0|\stackrel{~}{U}|0`$ with accuracy $`2^{N2}`$ is given.

To complete the proof of this reduction we note that each factor in (9) acts on 4 bits at most. So, by the Solovay – Kitaev theorem  the required efficient approximation $`\stackrel{~}{U}`$ can be constructed in $`poly(n)`$ time.

### 1.2 Reduction the $`𝑴𝑬`$ oracle to the $`𝑾𝑬`$ oracle

Consider a quantum circuit in the basis $`B`$ that realize an operator $`U`$. All computational paths of this circuit have the same scale factor $`1/2^{N/2}`$ where $`N`$ is the total number of Hadamard gates in the circuit. So, to compute the value $`0|U|0`$ we need to sum up phases along all computational paths that start at $`|0`$ and finish at $`|0`$. Note that an Hadamard gate that applies last to the particular qubit doesn’t affect the number of $`|0`$ to $`|0`$ paths. Let $`d`$ be the number of other Hadamard gates. Then there are $`2^d`$ summands in the expression of $`0|U|0`$. Each summand is indexed by a binary vector $`u𝔽_2^d`$ consisted of qubit values just after application of the corresponding Hadamard gate.

A phase shift along a computational path consists of two factors. The first describes the action of $`T`$ operators. It has a form

$$\underset{k=1}{\overset{t}{}}\omega ^{\mathrm{}_k(u)},\text{where}\omega =e^{i\pi /4},$$

and $`\mathrm{}_k(u)`$ is a $`𝔽_2`$-linear form. The second factor describes phase shifts induced by Hadamard gates. It has a form

$$(1)^{B(u)},\text{where}B(u)=\underset{j=1}{\overset{d}{}}b_j(u)u_j.$$

$`𝔽_2`$-linear form $`b_j`$ in this expression is the value of a qubit just before application of $`j`$th Hadamard gate. Using equations $`i^{ab}=i^ai^b(1)^{ab}`$ and $`i=\omega ^2`$ the second factor can be rewritten in the form similar to the first factor.

Finally, we get the following expression of the matrix element $`0|U|0`$:

$$0|U|0=\frac{1}{2^{N/2}}\underset{u𝔽_2^d}{}\underset{k=1}{\overset{n}{}}\omega ^{\beta _k(u)},$$
(12)
where $`\beta _k`$ are $`𝔽_2`$-linear forms.

It is clear that up to the factor $`1/2^{N/2}`$ this expression equals $`w_C(\omega )`$—the value of weight enumerator at the point $`\omega `$ for some binary linear code $`C`$. Coefficients of linear forms $`\beta _k(u)`$ are columns of the generator matrix for this code $`C`$.

So, we reduce the matrix element problem to the weight enumerator evaluation problem.

### 1.3 Reduction the $`𝑾𝑬`$ oracle to a $`\mathrm{𝐆𝐚𝐩𝐏}`$ oracle

Let $`C`$ be a linear binary code. It is a $`d`$-dimensional subspace of an $`n`$-dimensional space $`𝔽_2^n`$. By usual properties of $`\mathrm{GapP}`$ functions  the function $`E(D(C),q)=w_C(q)\mathrm{GapP}`$, where $`D(C)`$ is a description of the code $`C`$ and $`q`$ is integer. The required reduction follows from the fact that all coefficients of $`w_C(q)`$ can be restored efficiently from the value $`E(D(C),2^n)`$.

## 2 An approximate evaluation

Theorem 1 establishes the hardness of high-accuracy approximation of the weight enumerator.

Now we will prove two results concerning the approximation with moderate accuracy. Our choice of accuracy bounds is motivated by observation of Knill and Laflamme  that the power of quantum computation is just the ability to produce a Bernoulli sequences with a parameter $`p`$ given by a real or imaginary part of a matrix element $`0|U|0`$ of an operator realized by a quantum circuit. So, value of the matrix element can be estimated with an accuracy $`n^{O(1)}`$ and exponentially small error probability.

The formula (12) shows that the value of weight enumerator $`w_C(\omega )`$ can be evaluated in some cases with rather high accuracy $`n^{O(1)}2^{N/2}`$. Interesting quantum computation that does not reduced directly to the probabilistic one corresponds the case $`N<2d/\mathrm{log}d`$.

Corollaries 1 and 2 show that the general case differs from the case of these special codes. In this section we give the proofs of these corollaries.

### 2.1 Combining codes: direct sums and wreath sums

In this section we introduce two operations with linear codes and describe corresponding transformations of weight enumerators.

Consider an $`[n_1,d_1]`$ linear code $`A`$ and an $`[n_1,d_1]`$ linear code $`B`$. Let $`G_A`$ and $`G_B`$ are the generator matrices of codes $`A`$ and $`B`$ respectively.

A direct sum $`AB`$ is a $`[n_1+n_2,d_1+d_2]`$ code defined by generator matrix as follows:

$$\left(\begin{array}{cc}G_A& 0\\ 0& G_B\end{array}\right).$$
(13)
The weight enumerator of the direct sum is just the product of weight enumerators of summands:

$$w_{AB}(q)=w_A(q)w_B(q).$$
(14)
A wreath sum $`AB`$ is a $`[n_1n_2,d_1+d_2]`$ code. To define elements of the generator matrix $`G_{AB}`$ for $`AB`$ we assume that rows of $`G_{AB}`$ are indexed by pairs $`(k_1,0)`$, $`1k_1d_1`$, and $`(0,k_2)`$, $`1k_2d_2`$, and columns are indexed by pairs $`(l_1,l_2)`$, $`1l_1n_1`$, $`1l_2n_2`$. In this setting the elements of $`G_{AB}`$ are given by

$`(G_{AB})_{(k,0)(l_1,l_2)}=(G_A)_{k_1l_1}`$ (15)
$`(G_{AB})_{(0,k)(l_1,l_2)}=(G_B)_{k_2l_2}.`$ (16)
Each pair of words $`(a,b)`$, $`aA,bB`$ determines in natural way a word $`c_{ab}AB`$. Each bit $`\alpha `$ in the codeword $`a`$ is extended to a string $`\alpha ^{l_2}`$. This string is added by the codeword $`b`$ to produce the $`l_2`$-block of the codeword $`c_{a,b}`$. So, for weights of these codewords the following relation holds

$$|c_{ab}|=|a|(n_2|b|)+|b|(n_1|a|).$$

From this equation we get the following relation between weight enumerators:

$$w_{AB}(q)=\underset{aA}{}q^{n_2|a|}w_B(1+q^{n_12|a|}).$$
(17)
### 2.2 Proof of Corollary 1

We will construct a Turing reduction of the weight enumerator evaluation problem to the $`\mathrm{WE}(\alpha )`$ problem.

Namely, we will evaluate weight enumerators of a code $`C_k=CI_k`$ at several point $`q_j`$. Here $`I_k`$ means the trivial code given by identity generator matrix. The parameter $`k`$ will be chosen later.

Choose $`r`$ such that $`\beta =2^\alpha /(2r)<1`$. Note that $`(2^{\alpha (d+k)})`$-evaluation of $`w_{C_k}`$ implies much stronger estimation of $`w_C`$: if $`|q1|<r`$ and $`|w_{C_k}(q)\stackrel{~}{w}(q)|<2^{\alpha (d+k)}`$ then

$$\left|\frac{w_C(q)(1+q)^k\stackrel{~}{w}}{(1+q)^k}\right|=\left|w_C(q)\frac{\stackrel{~}{w}}{(1+q)^k}\right|<2^{\alpha d}\beta ^k=2^{\mathrm{\Omega }(k)}.$$
(18)
(In asymptotic notations we assume that $`k\mathrm{}`$ while other parameters are fixed.)

To compute the coefficients of weight enumerator $`w_C(q)`$ we use this tight bound for $`n+1`$ points inside the circle $`|z1|<r`$. The distances between these points can be chosen greater than $`r/(2(n+1))`$. So, the inverse of Vandermonde determinant for these points is $`2^{O(n^2\mathrm{log}n)}`$.

Thus, for some $`k=poly(n,d)`$ we can compute the coefficients of the weight enumerator with accuracy $`1/2`$.

### 2.3 Proof of Corollary 2

Note that reduction in the subsection 1.2 uses the evaluation of the weight enumerator at the point $`\omega `$ only.

To construct a Turing reduction of the weight enumerator evaluation problem to $`\mathrm{WE}_\omega (\alpha )`$ problem we will obtain an evaluation of the weight enumerator of a code $`C`$ at point $`\omega `$ using the ($`2^{\alpha (n+a+k)}`$)-estimations of weight enumerators for codes $`C(a,k)=C_aI_k`$ at the same point. Here $`I_k`$ stands for a trivial code of dimension $`k`$ and code $`C_a`$ is produced from the code $`C`$ by padding $`a`$ zeroes to each codeword. From (17) we get for these codes:

$$w_{C(a,k)}=\underset{wC}{}\omega ^{k|w|}(1+\omega ^{n+a2|w|})^k.$$
(19)
We assume that $`k1(mod8)`$.

Divide the weight enumerator sum for code $`C`$ in four parts:

$$w_C(\omega )=M_0+M_1+M_2+M_3;M_j=\underset{wC,|w|j(mod4)}{}\omega ^j.$$
(20)
It is clear that $`|M_j|<2^n`$.

Now rewrite (19) as

$$w_{C(a,k)}=M_0(1+\omega ^{n+a})^k+M_1(1i\omega ^{n+a})^k+M_2(1\omega ^{n+a})^k+M_3(1+i\omega ^{n+a})^k$$
(21)
Take $`a`$ such that $`n+a=1(mod8)`$. Let $`\stackrel{~}{w}_1`$ be an ($`2^{\alpha (n+a+k)}`$)-evaluation of the $`w_{C(a,k)}`$ and $`N=|1+\omega |`$. From (21) we have

$$|M_0+M_1\omega ^1+M_2\frac{(1\omega )^k}{N^k}+M_3\frac{(1\omega ^1)^k}{N^k}\frac{\stackrel{~}{w}_1}{N^k}|<2^{\alpha (n+a+k)}/N^k=2^{\alpha (n+a)}\left(\frac{2^\alpha }{N}\right)^k.$$
(22)
The last factor is $`2^{\mathrm{\Omega }(k)}`$. Note that

$$\frac{(1\omega )^k}{N^k}=2^{\mathrm{\Omega }(k)}\text{and}\frac{(1\omega ^1)^k}{N^k}=2^{\mathrm{\Omega }(k)}.$$

So, for some $`k=poly(n,\mathrm{log}\epsilon )`$ we obtain an $`\epsilon `$-evaluation $`\mu _0`$ of $`M_0+M_1\omega ^1`$.

Incrementing $`a`$ by 2 in the above argument we get $`\epsilon `$-evaluations $`\mu _j`$ of $`M_j+M_{j+1}\omega ^1`$. Solving the linear system

$$\begin{array}{cc}& x_0+x_1\omega ^1=\mu _0,\hfill \\ & x_1+x_2\omega ^1=\mu _1,\hfill \\ & x_2+x_3\omega ^1=\mu _2,\hfill \\ & x_3+x_0\omega ^1=\mu _3,\hfill \end{array}$$
(23)
we get the $`O(\epsilon )`$-approximations for $`M_j`$. It is enough to solve the weight enumerator evaluation problem for code $`C`$.


text:
# Precise Simulation of Near-critical Fluid Coexistence

\[

## Abstract

We present a novel method to derive liquid-gas coexisting densities, $`\rho ^\pm (T)`$, from grand canonical simulations (without knowledge of $`T_\text{c}`$ or criticality class). The minima of $`Q_Lm^2_L^2/m^4_L`$ in an $`L`$$`\times `$$`L\times `$$`L`$ box with $`m=\rho \rho _L`$ are used to generate recursively an unbiased universal finite-size scaling function. Monte Carlo data for a hard-core square-well fluid and for the restricted primitive model electrolyte yield $`\rho ^\pm `$ to $`\pm 1`$-$`2\%`$ of $`\rho _\text{c}`$ down to 1 part in $`10^4`$-$`10^3`$ of $`T_\text{c}`$ (and confirm well Ising character). Pressure mixing in the scaling fields is unequivocally revealed and indicates Yang-Yang ratios $`R_\mu =0.04_4`$ and $`0.2_6`$ for the two models, respectively.
PACS numbers: 64.60.Fr, 02.70.Rr, 05.70.Jk, 64.70.Fx

\]

Determining phase boundaries, critical points, and universality classes for various models that lack a clear symmetry has presented a serious difficulty in computer simulations . To tackle this problem, understanding scaling behavior in systems of finite size is crucial. However, as recently stressed , an important issue arises for asymmetric fluid criticality, even in the thermodynamic limit, namely, the potential presence of a Yang-Yang anomaly, in which the second derivative of the chemical potential, $`\mu _\sigma (T)`$, on the gas-liquid phase boundary diverges when the critical point, $`T_\text{c}`$, is approached from below. To describe a Yang-Yang anomaly requires pressure mixing in the scaling fields . This also generates a term varying as $`|t|^{2\beta }`$ \[with $`t`$$``$$`(TT_\text{c})/T_\text{c}`$\] in the gas-liquid coexistence diameter, that dominates the previously recognized $`|t|^{1\alpha }`$ term  and further distorts coexistence curves near criticality.

Our aim here is to show how coexistence curves may be estimated precisely and reliably near asymmetric critical points using grand canonical simulations, and to check our current understanding of scaling in such cases . It transpires that a finite-size scaling analysis at $`T_\text{c}`$ also elucidates pressure mixing and allows us to measure its strength using simulation data.

Figure 1 presents our estimates of $`\mathrm{\Delta }\rho _{\mathrm{}}(T)`$$``$$`\rho ^+`$$``$$`\rho ^{}`$,

the density discontinuity across the phase boundary, for a hard-core square-well (HCSW) fluid and for the restricted primitive model (RPM) electrolyte, where $`\rho ^+(T)`$ and $`\rho ^{}(T)`$ are the coexisting densities of liquid and vapor. The crosses represent new estimates obtained as explained below, while the open circles were derived previously directly from the observed double-peaked structure of the density distribution function in a finite grand canonical ensemble . Evidently the new approach yields estimates of $`\rho ^+(T)`$ and $`\rho ^{}(T)`$ of precision $`\pm 1`$-$`2\%`$ of $`\rho _\text{c}`$ or better, for temperatures $`1.5`$ to $`2.5`$ decades closer to the critical point. These results confirm convincingly that both models belong (as now expected ) to the same $`(d`$=$`3)`$-dimensional Ising universality class: see below and the dashed line in Fig. 1.

To outline the established situation , recall that for $`T`$$`<`$$`T_\text{c}`$ the grand canonical equilibrium distribution of the density, $`𝒫_L(T;\rho )`$, in a finite system of dimensions $`L^d`$ with periodic boundary conditions, has two Gaussian peaks near $`\rho ^\pm (T)`$ when $`L`$$``$$`a`$, where $`a`$ measures the particle size. For large $`L`$ the two peaks are clearly separated and thus provide reasonable estimates for the coexisting densities via the equal-weight prescription  — the open circles in Fig. 1 . However, when $`T_\text{c}`$ is approached, finite-size effects, arising from the divergence of the correlation length, soon blur the distinction between the vapor and liquid states thereby seriously hampering the reliable estimation of the coexistence curve. An alternative procedure applicable near $`T_\text{c}`$ is thus imperative.

Accordingly, we study the finite-system parameter $`Q_L`$ defined  by the dimensionless moment-ratio

$$Q_L(T;\rho _L)m^2_L^2/m^4_L,m=\rho \rho _L,$$
(1)
where $`_L`$ denotes a grand canonical expectation value at fixed $`T`$ and $`\mu `$. As well known, $`Q_L`$$``$$`\frac{1}{3}`$ when $`L`$$``$$`\mathrm{}`$ in any single-phase region of the $`(\rho ,T)`$ plane while $`Q_L`$$``$$`1`$ on the coexistence diameter, $`\rho _{\text{diam}}(T)`$$``$$`\frac{1}{2}(\rho ^++\rho ^{})`$. At criticality, $`Q_L`$ rapidly approaches a universal value $`Q_\text{c}`$ , e.g., $`Q_\text{c}`$$`=`$$`0.623_6`$ for $`d`$=$`3`$ Ising systems. The $`Q`$-loci, $`\rho _Q(T;L)`$, on which $`Q_L`$ attains isothermal maxima, have recently provided a route to estimating $`T_\text{c}`$ and $`\rho _\text{c}`$ with unprecedented precision .

In the two-phase region it has been known, but little appreciated, for some time \[10(a),13\], that $`Q_L(T;\rho )`$ displays a surprising singular behavior when $`L`$$``$$`\mathrm{}`$ . This is illustrated by the dashed-line plots in Fig. 2,

which follow directly from the double-peaked structure of $`𝒫_L(T;\rho )`$ below $`T_\text{c}`$  (together with our estimates for $`\rho ^+`$ and $`\rho ^{}`$). Specifically, $`Q_{\mathrm{}}(T;\rho )`$ exhibits a discontinuous drop from $`Q_{\mathrm{}}`$$`=`$$`\frac{1}{3}`$ to $`Q_{\mathrm{}}`$$`=`$$`0`$ on the two-phase boundaries, $`\rho ^{}`$ and $`\rho ^+`$, and a continuous (but nonconvex ) form for $`\rho ^{}`$$`<`$$`\rho `$$`<`$$`\rho ^+`$. For finite systems, however, the singularities are rounded and, as seen from the histogram-reweighted Monte Carlo simulations presented in Fig. 2, $`Q_L(T;\rho _L)`$ displays two smooth isothermal minima close to $`\rho ^+(T)`$ and $`\rho ^{}(T)`$. It is notable that while the HCSW data are fairly symmetrical about $`\rho _{\text{diam}}`$, the RPM displays a remarkably strong asymmetry.

Clearly, it is tempting to extrapolate these minima in order to estimate $`\rho ^+(T)`$ and $`\rho ^{}(T)`$ . However, when $`T_\text{c}`$ is approached, naive extrapolation fails badly owing to the finite-size effects: indeed the graph of $`Q_L(T;\rho _L)`$ still exhibits two distinct minima at and above $`T=T_\text{c}`$. Hence, some more powerful approach is necessary.

The behavior of $`Q_L(T;\rho _L)`$ near criticality can be understood via a recently developed ‘complete’ scaling theory that explicitly encompasses pressure mixing . Specifically, the full thermodynamics of a one-component fluid near the bulk critical point $`(p_\text{c},T_\text{c},\mu _\text{c})`$ can be described with three relevant scaling fields

$`\stackrel{~}{p}=\stackrel{ˇ}{p}k_0tl_0\stackrel{ˇ}{\mu }+\mathrm{},`$ (2)
$`\stackrel{~}{t}`$ $`=`$ $`tl_1\stackrel{ˇ}{\mu }j_1\stackrel{ˇ}{p}+\mathrm{},\stackrel{~}{h}=\stackrel{ˇ}{\mu }k_1tj_2\stackrel{ˇ}{p}+\mathrm{},`$ (3)
where the dimensionless deviations of the pressure and chemical potential from criticality are $`\stackrel{ˇ}{p}`$$``$$`(pp_\text{c})/\rho _\text{c}k_\text{B}T_\text{c},`$ and $`\stackrel{ˇ}{\mu }`$$``$$`(\mu \mu _\text{c})/k_\text{B}T_\text{c}`$: the coefficients $`j_1`$ and $`j_2`$ measure the degree of pressure mixing, the Yang-Yang ratio $`(`$$`T\mu _\sigma ^{\prime \prime }/C_V)`$ being fixed by $`R_\mu `$$`=`$$`j_2/(1j_2)`$ . For a finite box of dimensions $`L^d`$ with periodic boundary conditions, the finite-size scaling hypothesis now asserts

$$\rho _\text{c}\stackrel{~}{p}L^dY(x,z),x=D\stackrel{~}{t}L^{1/\nu },z=U\stackrel{~}{h}/|\stackrel{~}{t}|^\mathrm{\Delta },$$
(4)
where we have used the hyperscaling relation $`d\nu `$$`=`$$`2\alpha `$ (valid for $`d`$$`<`$$`4`$) and, for simplicity, neglected corrections to scaling. Note that $`D`$ and $`U`$ are nonuniversal amplitudes (of dimensions $`L^{1/\nu }`$ and $`L^0`$, respectively), while $`Y(x,z)`$ is a universal function that is even in $`z`$ and independent of microscopic details while depending on the geometry and the boundary conditions of the system.

It follows that the full scaling expression for $`Q_L`$ is

$$𝒬_Q(x,z)[1+A_jL^\kappa 𝒬_j(x,z)+A_lL^\lambda 𝒬_l(x,z)+\mathrm{}],$$
(5)
 with exponents and nonuniversal amplitudes

$`\kappa `$ $`=`$ $`\beta /\nu ,A_j=j_2D^\mathrm{\Delta }U/\rho _\text{c},`$ (6)
$`\lambda `$ $`=`$ $`(\mathrm{\Delta }1)/\nu ,A_l=(l_1+j_1)D^{1\mathrm{\Delta }}/(1j_2),`$ (7)
while the scaling functions $`𝒬_Q`$, $`𝒬_j`$, and $`𝒬_l`$ depend only on derivatives of $`Y(x,z)`$ thereby being universal. The symmetry of $`Y(x,z)`$ implies that $`𝒬_Q`$ is even in $`z`$ while $`𝒬_j`$ and $`𝒬_l`$ are odd. Notice that the pressure mixing coefficient $`j_2`$ provides the dominant asymmetric $`L`$-dependent correction (with Ising values $`\kappa `$$`=`$$`0.51_7`$$`<`$$`\lambda `$$`=`$$`0.89_6`$) which, indeed, describes the strong asymmetric behavior of $`Q_L(T;\rho _L)`$ for the RPM seen in Fig. 2(b).

Of course, the mean density $`\rho _L`$ also has a scaling form which we choose to write as

$`y(T;L)`$ $``$ $`2[\rho _L\rho _{\text{diam}}(T)]/\mathrm{\Delta }\rho _{\mathrm{}}(T)`$ (8)
$`=`$ $`𝒴\left[1+A_jL^\kappa 𝒴_j+A_lL^\lambda 𝒴_l+\mathrm{}\right],`$ (9)
where, again, the scaling functions $`𝒴(x,z)`$, $`𝒴_j(x,z)`$, and $`𝒴_l(x,z)`$ derive from $`Y(x,z)`$ and are universal, while $`𝒴`$ is odd in $`z`$, and $`𝒴_j`$ and $`𝒴_l`$ are even.

The crucial point here is that $`\rho _{\text{diam}}`$$``$$`(\rho ^++\rho ^{})`$ and $`\mathrm{\Delta }\rho _{\mathrm{}}`$$``$$`(\rho ^+\rho ^{})`$ embody the desired coexistence values $`\rho ^+(T)`$ and $`\rho ^{}(T)`$. Our strategy will be to determine values for $`\rho _{\text{diam}}`$ and $`\mathrm{\Delta }\rho _{\mathrm{}}`$ so that the minima of $`Q_L(T;\rho _L)`$, say, $`Q_\text{m}^+(T;L)`$ and $`Q_\text{m}^{}(T;L)`$, and their locations, $`\rho _\text{m}^+(T;L)`$$`>`$$`\rho _\text{m}^{}(T;L)`$, satisfy appropriate scaling relations. We focus first on $`\mathrm{\Delta }\rho _{\mathrm{}}`$ and, to minimize the effects of asymmetry (arising from the mixing coefficients $`j_2`$, $`j_1`$ and $`l_1`$), we examine the mean and difference

$$\overline{Q}_{\text{min}}\frac{1}{2}(Q_\text{m}^++Q_\text{m}^{}),\mathrm{\Delta }y_{\text{min}}\frac{1}{2}(y_\text{m}^+y_\text{m}^{}).$$
(10)
Now, on evaluating (4) and (6) at $`z_{\text{min}}^\pm `$ (which asymptotically fixes $`Q_\text{m}^\pm `$) and formally eliminating $`x`$$``$$`tL^{1/\nu }`$ between the resulting expressions, we see that $`\overline{Q}_{\text{min}}(T;L)`$ and $`\mathrm{\Delta }y_{\text{min}}(T;L)`$ should be related in a way that, to the orders displayed, is independent of $`T`$ and $`L`$ and (up to the neglected corrections to scaling) reflects only the universality class of the critical system under consideration. A priori this class is unknown — and, indeed, is to be determined. However, for any scalar order parameter the two-peaked, double-Gaussian structure of $`𝒫_L(T;\rho )`$ should be reproduced asymptotically when $`L`$$``$$`\mathrm{}`$ at fixed $`T`$$`<`$$`T_\text{c}`$. On this basis it is straightforward to calculate the universal relation for $`\overline{Q}_{\text{min}}`$$``$$`0`$: we find \[12(a)\]

$$\mathrm{\Delta }y_{\text{min}}(q)=1+\frac{1}{2}q+𝒪(q^2),q\overline{Q}_{\text{min}}\mathrm{ln}(4/e\overline{Q}_{\text{min}}),$$
(11)
which, to this order, is independent of any asymmetry.

Finally, we can employ our scaling analysis to generate the limiting coexistence curve recursively using finite-size simulation data for $`Q_L`$. Appropriate initial steps are: (i) Collect data sets $`\{Q_\text{m}^\pm (T;L_i),\rho _\text{m}^\pm (T;L_i)\}`$ for a range of values $`\{L_i\}_1^n`$ at fixed values of $`T`$$``$$`T_\text{c}`$. (ii) For a value $`T`$$`=`$$`T_0`$ sufficiently low that $`\overline{Q}_{\text{min}}`$$``$$`0.03`$ \[which corresponds to well separated peaks in $`𝒫_L(T_0;\rho )`$\], choose a density-jump value, say $`\mathrm{\Delta }\rho _{T_0}`$, independent of $`i`$, which leads to the best fit of $`\mathrm{\Delta }y_{\text{min}}^{(i)}`$$``$$`[\rho _\text{m}^+(T_0;L_i)\rho _\text{m}^{}(T_0;L_i)]/\mathrm{\Delta }\rho _{T_0}`$ vs. $`q_0^{(i)}`$$``$$`q(T_0;L_i)`$ to the relation (8) at small $`q`$: see the dashed lines in Fig. 3.

In light of the scaling relations (4) and (6), the parameter $`\mathrm{\Delta }\rho _{T_0}`$ can then be identified as an estimate for $`\mathrm{\Delta }\rho _{\mathrm{}}(T_0)`$. (iii) Increase $`T_0`$ to $`T_1`$$`=`$$`T_0+\mathrm{\Delta }T_0`$ by a small $`\mathrm{\Delta }T_0`$, chosen so that the new set $`\{q_1^{(i)}\}_1^n`$ overlaps the previous one. (iv) Determine a new value, $`\mathrm{\Delta }\rho _{T_1}`$, so that the plotted data display an optimal collapse that extends the previous numerical scaling function to larger values of $`q`$: see the gradual departure of the fits from the dashed lines in Fig. 3 as $`q`$ increases. In practice we have found that $`n`$$`=`$$`3`$ distinct box sizes with $`L_3`$$``$$`1.3L_1`$ may well suffice. (v) Repeat steps (iii) and (iv) generating successive estimates for $`\mathrm{\Delta }\rho _{\mathrm{}}(T_j)`$ for $`j`$$`=`$$`2`$, $`3`$, $`\mathrm{}`$. Smaller increments $`\mathrm{\Delta }T_j`$ are needed as $`T_j`$$``$$`T_\text{c}`$ and the $`q_j^{(i)}`$ increase to $`q_\text{c}`$$`=`$$`Q_{\text{min}}^\text{c}\mathrm{ln}(4/eQ_{\text{min}}^\text{c})`$ (see Fig. 3) so that histogram-reweighting procedures are crucial .

Figure 3 presents a scaling plot for the HCSW fluid constructed in this fashion: system sizes $`L^{}`$$``$$`L/a`$$`=`$$`9,10.5,`$ and $`12`$ were used and led to the estimates shown in Fig. 1 for $`\mathrm{\Delta }\rho _{\mathrm{}}(T)`$ from $`|t|`$$``$$`0.23`$ down to $`|t|`$$``$$`10^4`$. Purely for ease of presentation, Fig. 3 displays $`(\mathrm{\Delta }y_{\text{min}})^\psi `$ for selected values of $`\psi `$. In fact, the scaling analysis indicates that $`\mathrm{\Delta }y_{\text{min}}(q)`$ should diverge like $`(q_\text{c}q)^\beta `$ when $`q`$$``$$`q_\text{c}`$ as $`T`$$``$$`T_\text{c}`$, with $`q_\text{c}`$ a universal value (depending on geometry and boundary conditions) \[12(b)\]. For the HCSW fluid with periodic boundary conditions we find $`Q_{\text{min}}^\text{c}`$$`=`$$`0.110_2`$. To lower precision, the RPM data yield the same scaling plots and value of $`Q_{\text{min}}^\text{c}`$ \[12(b)\]. On the other hand, the approximate scaling form proposed by Tsypin and Blöte  for $`𝒫_L(T_\text{c};\rho )`$ for ($`d`$=$`3`$) Ising models gives $`Q_{\text{min}}^\text{c}`$$``$$`0.117`$, only $`6\%`$ higher than we observe. For ($`d`$=$`2`$) Ising models we estimate $`Q_{\text{min}}^\text{c}`$$``$$`0.28`$ using data in .

Evidently, the choice of $`\psi `$$`=`$$`1/\beta `$ should yield a plot that intersects the $`q`$ axis linearly; indeed, for the Ising value, $`\beta _{\text{Is}}`$$`=`$$`0.32_6`$, this is so. But, we emphasize that this observation plays no role in the calculation of Fig. 1.

Clearly, uncertainties in choosing $`\mathrm{\Delta }\rho _{T_j}`$, $`\mathrm{\Delta }\rho _{T_{j+1}}`$, $`\mathrm{}`$ in steps (ii) and (iv) will propagate. Well below $`T_\text{c}`$ (where care must be taken to ensure two-phase equilibrium) we can fit the limiting behavior (8) with a precision of $`\pm 1.0\%`$ or better in $`\mathrm{\Delta }\rho _T/\rho _\text{c}`$. The overall uncertainties then grow by factors of $`5`$-$`10`$ as $`|t|`$ decreases to $`10^4`$ for the HCSW fluid and $`10^3`$ for the RPM: see Fig. 1.

It is also remarkable that the $`\mathrm{\Delta }\rho _{\mathrm{}}(T)`$ estimates imply values for $`T_\text{c}`$. For the HCSW fluid we thus find $`T_\text{c}^{}`$$``$$`1.21821(2)`$ which lies close to the upper confidence limit of the previous estimate $`T_\text{c}^{}`$$``$$`1.2179(3)`$ : see also \[13(a)\] Eq. (5.6). For the RPM we obtain $`T_\text{c}^{}`$$``$$`0.05069(2)`$ which agrees precisely with Ref. . Explicit fits to $`\mathrm{\Delta }\rho _{\mathrm{}}(T)`$, that allow for the leading correction terms, yield $`\beta `$$`=`$$`0.324(10)`$ for the HCSW fluid and $`\beta `$$`=`$$`0.34(5)`$ for the RPM, so providing independent, albeit weaker confirmation of the Ising behavior established using data confined to $`TT_\text{c}`$ .

The scaling results (4) and (5) suggest that evidence for a pressure-mixing coefficient $`j_2`$ might be detected in finite-size data. Indeed, a detailed calculation \[12(b)\] of the asymmetry seen in the minima of $`Q_L`$ at $`T`$$`=`$$`T_\text{c}`$ yields

$$\left(𝒜_{\text{min}}\frac{Q_\text{m}^+Q_\text{m}^{}}{Q_\text{m}^++Q_\text{m}^{}}\right)_\text{c}=A_jc_jL^\kappa +A_lc_lL^\lambda +\mathrm{},$$
(12)
where $`c_j`$ and $`c_l`$ are universal numbers determined by expansion coefficients of $`Y(0,z)`$ about the minima at $`z_{\text{min}}^\pm `$. Recall from (5) that $`A_j`$ is proportional to $`j_2`$.

In Fig. 4 we present data for $`𝒜_{\text{min}}^\text{c}(L)`$ for the RPM and the HCSW fluid:

even by eye, the former strongly suggest a leading exponent closer to $`\kappa `$$`=`$$`0.51_7`$ than to $`\lambda `$$`=`$$`0.89_6`$. The fits in Fig. 4, using only the two leading terms in (9), support this but also indicate a weak $`j_2`$ contribution of opposite sign for the HCSW fluid. Further fairly elaborate analysis  yields $`j_2`$$`=`$$`0.35(7)`$, implying a strong, $`R_\mu `$$`=`$$`0.26(4)`$, Yang-Yang anomaly for the RPM, while $`j_2`$$`=`$$`0.042(3)`$ and $`R_\mu `$$`=`$$`0.044(3)`$ for the HCSW fluid. The latter result is consistent with the earlier, much less precise estimate $`R_\mu `$$``$$`0.08(12)`$ .

Finally, to determine the diameter $`\rho _{\text{diam}}(T)`$ we compare $`\overline{y}_{\text{min}}`$$``$$`\frac{1}{2}(y_\text{m}^++y_\text{m}^{})`$ and $`𝒜_{\text{min}}(T;L_i)`$. Analysis of the two-Gaussian limit \[12(b)\] yields $`\overline{y}_{\text{min}}/𝒜_{\text{min}}`$$`=`$$`\frac{1}{2}\overline{q}+𝒪(\overline{q}^2)`$ with $`\overline{q}`$$``$$`q\overline{Q}_{\text{min}}`$ which is again universal in leading order. Owing to the asymmetric terms in (4) and (6) the analogous scaling plots are now more sensitive to nonuniversal details and exhibit small, $`L`$-dependent corrections when $`q`$ approaches $`q_\text{c}`$. Nevertheless, the approach succeeds and the critical densities, $`\rho _\text{c}^{}`$, predicted from the diameters when $`T`$$``$$`T_\text{c}`$ are fully consistent with the previous, $`T`$$``$$`T_\text{c}`$ estimates . Details for both the RPM and the HCSW fluid will be presented elsewhere \[12(b)\].

In summary, we have shown how the finite-size scaling information hidden in precise simulation data can be systematically extracted via a novel $`Q`$-minima recursive algorithm to yield coexistence curves far closer to $`T_\text{c}`$ and with a much higher precision than previously appeared possible. As a byproduct, pressure mixing has been quantitatively resolved.

We are grateful to Dr. G. Orkoulas for help in generating $`Q_L`$ data for the hard-core square-well fluid, and to him, Prof. K. Binder, and Prof. A. Z. Panagiotopoulos for their interest. The support of the National Science Foundation (through Grant No. CHE 99-81772 to M.E.F.) and of the Department of Energy (for E.L. via the F. Seitz Materials Research Laboratory under contract No. DEFG02-91ER45439) is gratefully acknowledged.


text:
# Kink-Antikink Unbinding Transition in the Two Dimensional Fully Frustrated XY Model

## I Introduction

The two dimensional (2D) fully frustrated XY model (FFXY) Villain ; TJ  is one of the most intriguing of “simple” statistical mechanics models. The doubly degenerate checkerboard pattern of vortices in the ground state leads to an Ising-like discrete Z(2) symmetry in addition to the Kosterliz-Thouless-like continuous O(2) symmetry associated with the uniform rotation of all phase angles KT . It remains controversial whether there are two distinct phase transitions $`T_{KT}<T_I`$, with $`T_{KT}`$ marking the breaking of the O(2) symmetry and $`T_I`$ marking the breaking of the Z(2) symmetry double ; Olsson1 ; Olsson2 , or rather a single transition in which both symmetries are broken simultaneously single .

Recently, Korshunov Korshunov  presented an argument for a new interfacial transition $`T_w`$ in the 2D FFXY model, lying well below the bulk transition(s), arising from the unbinding of step excitations of unit height (“kink-antikink pairs”) on the domain walls associated with the Z(2) symmetry. Korshunov argued that the kink-antikink unbinding transition leads to a decoupling of phase coherence across domain boundaries, supporting the identification of the FFXY with the coupled XY-Ising model XYIsing . Korshunov further argued that this effect necessarily leads to the scenario of two distinct bulk transitions, $`T_{KT}<T_I`$.

Earlier, Lee and co-workers, first in simulations of the 2D FFXY model with Langevin dynamics LeeXY , and then in its dual Coulomb gas (CG) model with Monte Carlo (MC) dynamics LeeCG , found evidence for a transition in domain wall morphology in simulations of the ordering kinetics of domain growth following a sudden quench. They interpreted this as a finite temperature roughening transition of the Ising-like domain walls. Jeon et al. Jeon  made similar conclusions in simulations of the 2D FFXY with resistively shunted junction dynamics. Korshunov Korshunov , however, has argued that domain walls should be rough at all temperatures.

In this paper we present the first direct numerical evidence demonstrating the existence of the kink-antikink unbinding transition at a temperature $`T_w`$ below the bulk transition(s). In agreement with Korshunov’s predictions, we show that phase angles on opposite sides of the domain wall decouple above $`T_w`$. The numerical value we find for $`T_w`$ is comparable to that of the morphological transition found by Lee et al. LeeCG , however we explicitly demonstrate that domain walls are rough at temperatures well below $`T_w`$. This indicates that the transition seen by Lee and co-workers is really kink-antikink unbinding, rather than roughening.

## II The Model

### II.1 The fully frustrated XY model

We study the 2D FFXY model on a square lattice, given by the Hamiltonian Villain ; TJ ,

$$[\theta (𝐫_i)]=\underset{i,\mu }{}V\left(\theta (𝐫_i+\widehat{\mu })\theta (𝐫_i)A_\mu (𝐫_i)\right).$$
(1)
Here $`\theta (𝐫_i)`$ is the thermally fluctuating phase angle of the planar XY spin on site $`𝐫_i=x_i\widehat{x}+y_i\widehat{y}`$ ($`x_i`$, $`y_i`$ integers) of a $`L_x\times L_y`$ periodic square lattice, $`\widehat{\mu }=\widehat{x},\widehat{y}`$ label the bond directions of the lattice, and $`A_\mu (𝐫_i)`$ is the quenched gauge field on the bond leaving site $`𝐫_i`$ in direction $`\widehat{\mu }`$ (with $`A_\mu (𝐫_i+\widehat{\mu })A_\mu (𝐫_i)`$). For full frustration, the $`A_\mu (𝐫_i)`$ are constrained so that their directed sum going counterclockwise around any plaquette $`P`$ of the lattice is fixed (modulus $`2\pi `$) to,

$$\underset{P}{}A_\mu (𝐫_i)=\pi .$$
(2)
To implement the constraint of Eq. (2), we use the specific gauge choice,

$$A_x(𝐫_i)=0,A_y(𝐫_i)=(1)^{x_i}(\pi /2).$$
(3)
The interaction potential $`V(\varphi )`$ is periodic on $`[0,2\pi )`$, with a single quadratic minimum at $`\varphi =0`$. We will take for $`V(\varphi )`$ the commonly used Villain function Villain2 ,

$$V(\varphi )=T\mathrm{ln}\left[\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}e^{J(\varphi 2\pi m)^2/2T}\right].$$
(4)
The boundary conditions for the phase angles are, in the most general case, given by, twist1 ; twist2

$$\theta (𝐫_i+L_\mu \widehat{\mu })\theta (𝐫_i)=\mathrm{\Delta }_\mu ,$$
(5)
where $`\mathrm{\Delta }_\mu [0,2\pi )`$ is the total twist applied across the system in direction $`\widehat{\mu }`$. $`\mathrm{\Delta }_\mu =0`$ corresponds to periodic boundary conditions. Alternatively, if one makes the change of variables, $`\theta ^{}(𝐫_i)\theta (𝐫_i)𝐫_i𝐝`$, with $`d_\mu \mathrm{\Delta }_\mu /L_\mu `$, then the system has periodic boundary conditions in the $`\theta ^{}(𝐫_i)`$ and the applied twist appears as an additive constant to the gauge field, $`A_\mu (𝐫_i)A_\mu (𝐫_i)+\mathrm{\Delta }_\mu /L_\mu `$.

To study the behavior of the Ising-like domain walls we consider systems with sizes $`L_x=L`$, $`L_y=L+1`$, with $`L`$ even. The odd length $`L_y`$ forces into the ground state checkerboard pattern of vortices a single straight domain wall running the length of the system in the $`\widehat{x}`$ direction. This is illustrated in Fig. 1a, where a $`(+)`$ signifies a vortex in the phase angles $`\theta (𝐫_i)`$, and a $`()`$ signifies the absence of a vortex.

Phase coherence in the FFXY model is most conveniently studied by considering the dependence of the total free energy $`F`$ on the total twist $`\mathrm{\Delta }_\mu `$ applied across the system (see Eq. (5)). In a phase coherent ordered state, we expect that $`F(\mathrm{\Delta }_\mu )`$ varies with the twist $`\mathrm{\Delta }_\mu `$; in a phase incoherent disordered state, we expect $`F(\mathrm{\Delta }_\mu )`$ is independent of $`\mathrm{\Delta }_\mu `$ in the thermodynamic limit of $`L\mathrm{}`$. The dependence of the free energy on $`\mathrm{\Delta }_\mu `$ is readily obtained by using fluctuating twist boundary conditions twist2 , in which one treats the applied twist $`\mathrm{\Delta }_\mu `$ as a thermally fluctuating degree of freedom. If $`Z`$ is the partition function for this ensemble, then the probability $`P(\mathrm{\Delta }_\mu )`$ of finding a state with a particular twist $`\mathrm{\Delta }_\mu `$ is given by,

$$P(\mathrm{\Delta }_\mu )=\frac{\mathrm{e}^{F(\mathrm{\Delta }_\mu )/T}}{Z},$$
(6)
and so the free energy with respect to a reference twist $`\mathrm{\Delta }_{\mu 0}`$ is,

$$F(\mathrm{\Delta }_\mu )F(\mathrm{\Delta }_{\mu 0})=T\mathrm{ln}\left[P(\mathrm{\Delta }_\mu )/P(\mathrm{\Delta }_{\mu 0})\right].$$
(7)
The probability $`P(\mathrm{\Delta }_\mu )`$ is directly measured within our fluctuating twist Monte Carlo simulation. We choose the reference twist $`\mathrm{\Delta }_{\mu 0}`$ to be the value of the twist that minimizes the free energy $`F(\mathrm{\Delta }_\mu )`$. For the gauge choice of Eq. (3), it is straightforward to see to that the minimizing twist in the $`\widehat{x}`$ direction is at $`\mathrm{\Delta }_{x0}=0`$. In our simulations we keep a fixed twist $`\mathrm{\Delta }_x=0`$, and consider only the dependence of the free energy on the varying twist $`\mathrm{\Delta }_y`$, transverse to the Ising-like domain wall that is introduced in our $`L\times (L+1)`$ systems (see Fig. 1a). In Fig. 2 we show sample results from our simulations for $`F(\mathrm{\Delta }_y)F(0)`$ vs. $`\mathrm{\Delta }_y`$ at two different values of $`T<T_w`$, for a system of size $`L=128`$. We see that $`F(\mathrm{\Delta }_y)`$ has two equal minima at $`\mathrm{\Delta }_y=0`$ and $`\pi `$ (i.e. periodic and antiperiodic boundary conditions). One of these minima corresponds to states where the domain wall sits at even values of the height $`y`$, while the other corresponds to states where the domain wall sits at odd values of the height $`y`$. Below the kink-antikink unbinding transition $`T_w`$, the system is in a state of broken translational symmetry; due to the free energy barrier between the two minima, states in which the domain wall is at an even height cannot be reached from states in which the domain wall is at an odd height. As noted by Korshunov Korshunov , this broken symmetry is restored when phase coherence transverse to the wall is lost, i.e. when $`F(\mathrm{\Delta }_y)`$ becomes independent of $`\mathrm{\Delta }_y`$, and so the free energy barrier between $`\mathrm{\Delta }_y=0`$ and $`\mathrm{\Delta }_y=\pi `$ vanishes. Alternatively viewed, when the domain wall changes its height by an odd number, the system acquires an average twist of $`\pi `$ in the $`\widehat{y}`$ direction. Thus, restoring the symmetry of domain wall translations leads to phase angle fluctuations that destroy phase coherence transverse to the direction of the wall.

In our numerical work we will use two convenient measures of the variation of $`F(\mathrm{\Delta }_\mu )`$ with $`\mathrm{\Delta }_\mu `$. The first is the helicity modulusTJ ; twist1 , $`\mathrm{{\rm Y}}_\mu `$, which measures the curvature of $`F(\mathrm{\Delta }_\mu )`$ at its minimum,

$`\mathrm{{\rm Y}}_\mu (L_x,L_y)={\displaystyle \frac{L_\mu ^2}{L_xL_y}}{\displaystyle \frac{^2F}{\mathrm{\Delta }_\mu ^2}}|_{\mathrm{\Delta }_\mu =0}`$ (8)
$`=`$ $`{\displaystyle \frac{1}{L_xL_y}}\left\{{\displaystyle \underset{i}{}}V^{\prime \prime }(\varphi _{i\mu })_0{\displaystyle \frac{1}{T}}\left[{\displaystyle \underset{i}{}}V^{}(\varphi _{i\mu })\right]^2_0\right\},`$
where $`\varphi _{i\mu }\varphi _\mu (𝐫_i)\theta (𝐫_i+\widehat{\mu })\theta (𝐫_i)A_\mu (𝐫_i)`$, $`V^{}`$ and $`V^{\prime \prime }`$ are the first and second derivatives of the Villain function of Eq. (4), and $`\mathrm{}_0`$ indicates a thermodynamic average in the ensemble with fixed twist $`\mathrm{\Delta }_\mu =0`$. A second measure is,

$$\mathrm{\Delta }F=F_{\mathrm{max}}F_{\mathrm{min}}=F(\pi /2)F(0),$$
(9)
where $`F_{\mathrm{max}}`$ and $`F_{\mathrm{min}}`$ are the maximum and minimum values of $`F(\mathrm{\Delta }_y)`$, as $`\mathrm{\Delta }_y`$ is varied at fixed $`\mathrm{\Delta }_x=0`$. Since $`\mathrm{{\rm Y}}_\mu `$ is an intensive quantity, it should approach a value independent of system size as $`L\mathrm{}`$. The parameter $`\mathrm{\Delta }F`$ scales as $`\mathrm{{\rm Y}}(\mathrm{\Delta }/L)^2L^D`$ in $`D`$ dimensions, and so for $`D=2`$ it also becomes independent of system size as $`L\mathrm{}`$.

### II.2 The Coulomb gas

Although our simulations are carried out in the XY variables $`\theta (𝐫_i)`$, it is helpful to consider the situation from the viewpoint of the dual CG model of logarithmically interacting half integer charges Villain ; cg . For the case of a fixed total twist $`\mathrm{\Delta }_\mu `$, the XY Hamiltonian of Eq. (1) maps onto,

$$_{\mathrm{CG}}=_0+_1.$$
(10)
$`_0`$ is the logarithmic interaction of the charges,

$$_0=\frac{1}{2}(2\pi J)\underset{i,j}{}q_iG(𝐫_i𝐫_j)q_j,$$
(11)
where $`q_i=\pm 1/2`$ are the half integer charges, neutrality is imposed, $`_iq_i=0`$, and $`G(𝐫)`$ is the 2D periodic lattice Coulomb potentialVillain ; cg  with $`G(𝐫)\mathrm{ln}|𝐫|`$ for large $`1|𝐫|L/2`$. $`_1`$ arises from the fixed twist boundary condition and is given by Vallat ; Bokil ,

$$_1=V_x\left(\mathrm{\Delta }_xA_x^0+\frac{2\pi p_y}{L_y}\right)+V_y\left(\mathrm{\Delta }_yA_y^0\frac{2\pi p_x}{L_x}\right),$$
(12)
where $`V_x`$ and $`V_y`$ are Villain functions as in Eq. (4), but with couplings $`J_x=J(L_y/L_x)`$ and $`J_y=J(L_x/L_y)`$ respectively, $`𝐩`$ is the total dipole moment,

$$𝐩=\underset{i}{}q_i𝐫_i,$$
(13)
and $`A_x^0=_xA_x(x,y=0)`$, $`A_y^0=_yA_y(x=0,y)`$. For the gauge choice of Eq. (3) we have,

$$A_x^0=0,A_y^0=(L+1)\pi /2.$$
(14)
For the ground state as illustrated in Fig. 1a, one has $`p_y=0`$, and so again it is easy to see from Eq. (12) that the total ground state energy is minimized when $`\mathrm{\Delta }_x=0`$. However, for this ground state one has $`p_x=L/4`$, hence the ground state energy is minimized when $`\mathrm{\Delta }_y=(L+1)\pi /2\pi /2=L\pi /2`$. If the location of the domain wall was shifted by one unit in height, then the ground state would have $`p_x=L/4`$, and the energy would be minimized when $`\mathrm{\Delta }_y=(L+1)\pi /2+\pi /2=L\pi /2+\pi `$. For $`L`$ even, as we have required, these two values, modulus $`2\pi `$, are just equal to $`0`$ and $`\pi `$.

The helicity modulus $`\mathrm{{\rm Y}}_y/J`$ maps KT ; twist1 ; twist2  onto the inverse dielectric function $`ϵ_x^1`$ of the CG. As noted by Korshunov Korshunov , in order for the domain wall to move a unit lattice spacing in height, it must first form a unit step of finite width $`\mathrm{}`$; $`\mathrm{}`$ must be even to preserve charge neutrality (see Fig. 1b). We denote the left hand edge of the step as the kink, and the right hand edge as the antikink. As the kink and antikink separate out to infinity, the domain wall moves one unit in height. As shown by Halsey Halsey , a corner in a domain wall carries with it a net charge of $`\pm 1/4`$. The kink, consisting of two successive corners with equal $`+1/4`$ charge, carries a net charge of $`q_{\mathrm{kink}}=+1/2`$; the antikink carries a net charge of $`q_{\mathrm{kink}}=1/2`$. At low temperatures, the logarithmic attraction between the kink and antikink charges keeps them bound with a largest separation $`\mathrm{}_{\mathrm{max}}(T)`$. At higher temperatures, entropy wins out over energy, and there is a kink-antikink unbinding transition at $`T_w`$, where $`\mathrm{}_{\mathrm{max}}(T_w)\mathrm{}`$. Above $`T_w`$, the kink-antikink unbinding leads to diverging dipole fluctuations in the $`\widehat{x}`$ direction, driving $`ϵ_x^1`$ (and hence $`\mathrm{{\rm Y}}_y`$) to zero.

The problem of logarithmically interacting charges in one dimension (1D) has been treated by Bulgadaev Bulgadaev . KoshunovKorshunov  has applied these results to the unbinding transition of the kink-antikink pair along the one dimensional Ising-like domain wall. To include the screening effect of charge excitations in the bulk of the system on either side of the Ising-like domain wall, we take as the coupling between kink-antikink pairs separated at large distance to be the helicity modulus of the FFXY for an ordinary $`L\times L`$ system, $`\mathrm{{\rm Y}}(L,L)`$, in the limit of large enough $`L`$. Applying Bulgadaev’s exact result for the unbinding transition temperature, we conclude,

$$\frac{2\pi \mathrm{{\rm Y}}(L,L)q_{\mathrm{kink}}^2}{T_w}=2,\mathrm{or}T_w=\frac{\pi }{4}\mathrm{{\rm Y}}(L,L).$$
(15)
One can reproduce this result using a Kosterlitz-Thouless-like argument KT  as follows. In analogy with Lee et al. LeeCG , we consider the total free energy to have a single “free” (i.e. unbound) kink in the domain wall (see Fig. 1c). Fixing the kink at a given position on the domain wall, its free energy (averaging over all other fluctuations) is that of an isolated $`+1/2`$ vortex in a medium with phase stiffness $`\mathrm{{\rm Y}}(L,L)`$; here we use Bulgadaev’s result Bulgadaev  that kink-antikink pairs in 1D do not lead to a renormalization of the kink-antikink interaction, and so any screening of their interaction is due to charge excitations in the bulk on either side of the domain wall, and so accounted for by the large $`L`$ value of $`\mathrm{{\rm Y}}`$. As $`L\mathrm{}`$, the leading contribution to this energy is $`E=\pi q_{\mathrm{kink}}^2\mathrm{{\rm Y}}\mathrm{ln}L`$. The entropy of the kink is just that associated with its position along the domain wall, $`S=\mathrm{ln}L`$. Combining gives $`F_{\mathrm{kink}}=ETS=(\pi \mathrm{{\rm Y}}/4T)\mathrm{ln}L`$, which as $`L\mathrm{}`$ gives the instability temperature for the formation of free kinks as $`T_w=\pi \mathrm{{\rm Y}}/4`$, in agreement with Eq. (15).

## III Numerical Results

### III.1 Helicity modulus

We now present our numerical results. At each temperature, our simulations consist of typically $`10^810^9`$ ordinary MC passes through the entire lattice for the largest system sizes. In Figs. 3 and 4 we plot the helicity moduli $`\mathrm{{\rm Y}}_x`$ and $`\mathrm{{\rm Y}}_y`$ vs. $`T`$, as computed by Eq. (8) in an ensemble with fixed twists $`\mathrm{\Delta }_x=\mathrm{\Delta }_y=0`$. We show an “ordinary” case (no domain wall at $`T=0`$) of size $`64\times 64`$, which is large enough that any finite size effects are negligible for the temperatures shown. In comparison, we also show several “anomalous” cases (percolating domain wall at $`T=0`$) of sizes $`L\times (L+1)`$. For the ordinary case, $`\mathrm{{\rm Y}}_x=\mathrm{{\rm Y}}_y`$ and the bulk transition (where the $`\mathrm{{\rm Y}}_\mu `$ jump discontinuously to zero in the thermodynamic limit) is Olsson2  at $`T_{KT}0.81J`$. In comparison, as $`L`$ increases in the anomalous case, $`\mathrm{{\rm Y}}_x`$ (parallel to the domain wall) in Fig. 3 approaches the value of the ordinary case, and so presumably vanishes at the same $`T_{KT}`$. However the curves of $`\mathrm{{\rm Y}}_y`$ (transverse to the domain wall) in Fig. 4 clearly decrease below that of the ordinary case, and presumably vanish in the thermodynamic limit at a lower $`T_w`$.

The reduction seen in $`\mathrm{{\rm Y}}_y`$ for the $`L\times (L+1)`$ systems, as compared to the $`L\times L`$ system, shown in Fig. 4, is due to the kinks at the Ising-like domain wall. To explicitly see this, we can consider the helicity modulus at finite wavevector, $`\mathrm{{\rm Y}}_y(k_x)`$, defined as the response to a small sinusoidal perturbation in the vector potential $`A_y(𝐫_i)`$. If we take,

$$A_y(𝐫_i)A_y(𝐫_i)+\underset{k_x}{}\delta A_{k_x}\mathrm{e}^{ik_xx_i}$$
(16)
then $`\mathrm{{\rm Y}}_y(k_x)`$ is defined by Olsson2 ; Upsk ,

$$\mathrm{{\rm Y}}_y(k_x)=\frac{1}{L_xL_y}\frac{^2F}{\delta A_{k_x}\delta A_{k_x}}|_{\delta A_{k_x}=0}.$$
(17)
In view of the discussion following Eq. (5), equating the application of a uniform twist $`\mathrm{\Delta }_\mu `$ to the addition of a constant to the gauge field $`A_\mu `$, the helicity modulus $`\mathrm{{\rm Y}}_y`$ of Eq. (8) can also be viewed as the zero wavevector helicity $`\mathrm{{\rm Y}}_y(k_x=0)`$. In the CG representation, $`\mathrm{{\rm Y}}_y(k_x)`$ becomes the usual formula for the wavevector dependent inverse dielectric function twist1 ; twist2 ,

$$\mathrm{{\rm Y}}_y(k_x)/J=1\frac{4\pi ^2J}{T}\frac{q(k_x)q(k_x)}{L_xL_yk_x^2},$$
(18)
where $`q(𝐤)=_i\mathrm{e}^{i𝐤𝐫_i}q_i`$ is the Fourier transform of the charge distribution.

Unlike $`\mathrm{{\rm Y}}_y`$ of Eq. (8), which measures the response to a uniform twist applied at the boundaries, $`\mathrm{{\rm Y}}_y(k_x)`$ measures the response to a spatially varying twist applied throughout the bulk of the system. For a homogeneous system with periodic boundary conditions, one in general expects $`\mathrm{{\rm Y}}_y=lim_{k_x0}\mathrm{{\rm Y}}_y(k_x)`$, since the spatially varying twist becomes uniform as $`k_x0`$, and $`\delta A_{k_x}\mathrm{\Delta }_y/L_y`$. For free boundary conditions however, where the phase angle $`\theta (x,L_y)`$ is not coupled to the phase angle $`\theta (x,0)`$, this equality does not hold. For free boundary conditions, the absence of any constraint (such as in Eq. (5)) relating $`\theta (x,L_y)`$ to $`\theta (x,0)`$ means that the phase angles are free to untwist any additive constant to the gauge field, $`A_y(𝐫_i)A_y(𝐫_i)+\mathrm{\Delta }_y/L_y`$, by choosing $`\theta (x,y+1)\theta (x,y)=\mathrm{\Delta }_y/L_y`$; hence if one computes $`\mathrm{{\rm Y}}_y`$ by Eq. (8) in a free boundary ensemble, one necessarily has $`\mathrm{{\rm Y}}_y=0`$ at any temperature. For the spatially varying twist of Eq. (16) however, no such transformation is possible since the perturbing twist is a strictly transverse vector function, while the phase angle differences give a strictly longitudinal vector function. In this case one finds that $`lim_{k_x0}\mathrm{{\rm Y}}_y(k_x)`$ has the same value, as $`L\mathrm{}`$, that one has for the system with periodic boundary conditions.

We expect a similar effect to be true in our present case. The kink-antikink pairs confined to the one dimensional Ising-like domain wall can be viewed as a relaxation of the boundary condition. They can unwind, or soften the energy of a uniform twist $`\mathrm{\Delta }_y`$ applied at the boundary, but cannot unwind a spatially varying twist $`\delta A_{k_x}`$ applied throughout the bulk of the system. We therefore expect that, as $`L\mathrm{}`$, $`lim_{k_x0}\mathrm{{\rm Y}}_y(k_x)`$ will equal the value of $`\mathrm{{\rm Y}}_y`$ obtained for an ordinary $`L\times L`$ system, representing the stiffness of the bulk of the system on either side of the domain wall; $`\mathrm{{\rm Y}}_y`$, however, will be a lower value including effects due to the polarization of the kink-antikink pairs localized to the domain wall.

In Fig. 5 we plot, at several different temperatures around $`T_w`$, $`\mathrm{{\rm Y}}_y(k_x)`$ for finite $`k_x`$ for $`L\times (L+1)`$ systems with sizes $`L=64`$ and $`128`$. When we plot the results versus the scaled axis of $`Lk_x`$, we see that the data for the two different system sizes collapse to essentially a common curve, $`u(Lk_x)`$, at each temperature. As $`L\mathrm{}`$, such scaling implies that $`lim_{k_x0}\left[lim_L\mathrm{}\mathrm{{\rm Y}}_y(k_x)\right]=lim_\kappa \mathrm{}u(\kappa )`$, is different from $`\mathrm{{\rm Y}}_y`$, just as we have argued. We may estimate $`lim_\kappa \mathrm{}u(\kappa )`$ by taking the maximum value of $`\mathrm{{\rm Y}}_y(k_x)`$ for each system size $`L`$. In Fig. 6 we plot $`\mathrm{{\rm Y}}_y^{\mathrm{max}}\mathrm{max}_{k_x}[\mathrm{{\rm Y}}_y(k_x)]`$ versus $`1/L`$ for the same temperatures as shown in Fig. 5. We give results for $`L=32,64`$ and $`128`$. At $`1/L=0`$ we plot the value of $`\mathrm{{\rm Y}}_y(L,L)`$ for $`L=64`$, representing the large $`L`$ limit for the helicity modulus of the ordinary system without the Ising domain wall. We see that the values of $`\mathrm{{\rm Y}}_y^{\mathrm{max}}`$ extrapolate perfectly to $`\mathrm{{\rm Y}}_y(L,L)`$ as $`1/L0`$. This confirms the following conclusion. As $`L\mathrm{}`$ in an $`L\times (L+1)`$ system, the finite wavevector helicity $`\mathrm{{\rm Y}}_y(k_x)`$ is equal to the corresponding helicity modulus of an ordinary $`L\times L`$ system for all finite values of $`k_x`$; this measures the stiffness of the bulk of the system on either side of the Ising domain wall, and is unaffected by the kinks on the wall. However the zero wavevector response $`\mathrm{{\rm Y}}_y`$ to a uniform twist $`\mathrm{\Delta }_y`$, shown in Fig. 4, is softened, and above $`T_w`$ reduced to zero, by the polarization of the kinks on the domain wall.

### III.2 Finite size dependence

Returning to Fig. 4, we have indicated the kink-antikink unbinding transition temperature $`T_w`$, as predicted in Eq. (15), by the intersection of the line $`4T/\pi J`$ with the helicity modulus of the ordinary $`L\times L`$ system. This gives an estimate of $`T_w0.71J`$. This value occurs noticeably above the point where many of the curves $`\mathrm{{\rm Y}}_y(L,L+1)`$ appear to cross. Such a crossing point, if remaining constant as $`L`$ increases, is generally taken as an estimate for the phase transition in a 2D XY system. To explain the difference between this crossing point at $`0.60J`$ and the above estimate $`T_w0.71J`$, we need to examine the finite size dependence of $`\mathrm{{\rm Y}}_y(L,L+1)`$ more carefully.

To get the most accurate results, we have found it better to work in the fluctuating twist ensemble and compute the phase coherence parameter $`\mathrm{\Delta }F`$ of Eq. (9), rather than work with periodic boundary conditions and measure $`\mathrm{{\rm Y}}_y`$. In Fig. 7 we plot our results for $`\mathrm{\Delta }F`$ versus system size $`1/L`$, for various temperatures $`T`$ above and below $`T_w`$. We use $`L\times (L+1)`$ systems with $`L`$ ranging from $`32`$ up to $`512`$. We note that at low $`T`$, the behavior of $`\mathrm{\Delta }F`$ is non-monotonic; as $`L`$ increases, the system first softens with a decreasing $`\mathrm{\Delta }F`$, but then stiffens again as $`\mathrm{\Delta }F`$ reaches a minimum and then increases. We denote the system size at this minimum by $`\xi _k`$. Since the system becomes stiffer on length scales $`L>\xi _k`$, we assume that $`\xi _k`$ determines the size of the largest kink-antikink pairs on the domain wall. At higher $`T`$, the minimum in $`\mathrm{\Delta }F`$ disappears, and $`\mathrm{\Delta }F`$ continues to decrease as $`L`$ increases. The temperature that separates the two behaviors is somewhere between $`0.68J`$ and $`0.74J`$, in good agreement with the result $`T_w0.71J`$ from Fig. 4, based on the theoretical prediction of Eq. (15).

As a further check on our results, we have also directly simulated a 1D neutral system of logarithmically interacting charges $`q_{\mathrm{kink}}=\pm 1/2`$. The 1D interaction potential is taken as the 2D $`L\times L`$ lattice Coulomb potential, as in Eq.(11), but with height separation fixed at $`y=0`$. We use an interaction coupling constant of $`2\pi \mathrm{{\rm Y}}_y(L,L)`$, with $`\mathrm{{\rm Y}}_y(L,L)`$ obtained from our simulations of the ordinary $`L\times L`$ 2D FFXY for $`L=64`$, in order to model as closely as possible the interaction between kinks on the domain wall in the true $`L\times (L+1)`$ FFXY system. Within this 1D simulation we measure the normalized histogram of the total dipole moment, $`P(p_x)`$, and use it to construct Gupta  what would be the free energy $`F(\mathrm{\Delta }_y)`$ of the corresponding $`L\times (L+1)`$ 2D FFXY system,

$$F(\mathrm{\Delta }_y)=T\mathrm{ln}\left[\underset{p_x}{}P(p_x)\mathrm{e}^{V_y\left(\mathrm{\Delta }_yA_y^0\frac{2\pi p_x}{L}\right)/T}\right]+\mathrm{const}.,$$
(19)
where $`V_y`$ is the Villain function as in Eq. (12), and “const.” is a constant term independent of $`\mathrm{\Delta }_y`$.

In Fig. 8 we plot the resulting $`\mathrm{\Delta }F=F(\pi /2)F(0)`$ versus $`1/L`$ for the same temperatures and sizes $`L`$ as in Fig. 7. The agreement between Figs. 7 and 8 is not exact, since the coupling between kinks is only equal to $`2\pi \mathrm{{\rm Y}}_y(L,L)`$ on large length scales; the true screening of the kink interaction due to charge excitations in the bulk on either side of the domain wall is length scale dependent. Moreover, the domain wall in the 2D FFXY model is not a strictly straight one dimensional line; the roughness of the domain wall (see following section) means that height fluctuations can add to the distance of separation between kinks. Nevertheless, the agreement is qualitatively very good, indicating again that it is the polarization of kink-antikink pairs along the Ising domain wall that is responsible for the decrease in the phase stiffness transverse to the domain wall.

One evident feature of both Figs. 7 and 8 is the very large finite size effect. The asymptotic behavior of $`\mathrm{\Delta }F`$ only sets in at quite large length scales. Equivalently, the correlation length of the kink-antikink pairs, $`\xi _k`$, grows large well below $`T_w0.71J`$. Fitting the data of Fig. 7 to a quadratic in $`\mathrm{ln}L`$, and determining $`\xi _k`$ from the minima of these fitted curves, we plot the resulting $`\xi _k`$ versus $`T`$ in Fig. 9. The rapidly growing $`\xi _k`$ means that it is difficult to get a very precise estimate of $`T_w`$ directly from the data of Fig. 7, without going to prohibitively large system sizes $`L512`$. We believe that this is also the reason that Lee al. LeeCG  report a lower valuefoot3  of $`T/J=2\pi (0.09\pm 0.01)=0.57\pm 0.06`$ for the “roughening transition” temperature in their dual 2D Coulomb gas. Lee et al.’s simulations are on a lattice of size $`L=64`$. From Fig. 9 we see that $`\xi _k32`$ at $`T/J0.58`$, hence kinks already look unbound above this temperature for such a small system size. The large values of $`\xi _k`$ can also be compared to other length scales in the system. The correlation length of the Ising-like order parameter, $`\xi _I`$, gives the typical size of an Ising-like domain excitation in an ordinary $`L\times L`$ FFXY model. From Ref. Olsson2,  (see Fig. 15) we find that $`\xi _I`$ is quite small below $`T_w0.71J`$; in particular $`\xi _I<2.5`$ for $`T<0.77J`$. Thus we expect that kink-antikink pairs on the domain wall enclosing a typical thermally excited Ising-like domain are always effectively unbound. The energetics of kink-antikink unbinding will only effect the morphology of domains that are much larger than those due to typical thermal excitations.

### III.3 Roughening

Finally, we consider the roughness of the domain wall. As noted by Koshunov Korshunov , an isolated step on a domain wall of height two (see Fig. 1d) carries no net charge, each corner giving opposite $`1/4`$ charges that therefore cancel. Such isolated height of two steps therefore cost finite energy, and should roughen the domain wall at any finite temperature. To verify this we have explicitly measured the domain wall width squared,

$$W^2=\frac{1}{N_s}\underset{s}{}(y(s)\overline{y})^2,$$
(20)
where $`s`$ is an index that counts horizontal length traveled along the domain wall ($`s=x`$ if there are no overhangs), $`y(s)`$ is the height of the domain wall at index $`s`$, $`\overline{y}`$ is the average height of the domain wall in the particular configuration, and $`N_s`$ is the number of horizontal links in the domain wall ($`N_s=L`$ if there are no overhangs). To avoid problems with periodic boundary conditions, we always measure the domain wall height as relative to some initial starting position. In Fig. 10 we plot $`W^2`$ vs. $`L`$ for various temperatures. The linear growth in $`W^2`$ as $`L`$ increases indicates that the domain wall is rough for all temperatures shown, including $`T<T_w`$. The domain wall diffusion constant, $`dW^2/dL`$, shows no observable singularity at $`T_w0.71J`$ (see inset to Fig. 10).

## IV Discussion

Our numerical results demonstrate the existence of the unbinding transition for kink-antikink pairs along the domain wall of Ising-like excitations in the 2D FFXY model, and that the behavior of this transition is in good agreement with the theoretical predictions of KorshunovKorshunov . We show that domain walls are rough at all temperatures, and therefore argue that the “roughening transition” claimed by Lee and co-workersLeeXY ; LeeCG  is really the kink-antikink unbinding transition. We show that the effects of kink-antikink pairs are not readily apparent for Ising-like domains such as are typically present due to thermal excitation; because of the large length $`\xi _k`$, such effects are important only for much larger domains.

One of Korshunov’s main motivations for investigating the kink-antikink unbinding transition was to argue that such a transition necessarily implies the existence of two separate bulk transitions, $`T_{KT}<T_I`$. We now present our own thoughts on this issue. In the original paper by Teitel and Jayaprakash TJ , two possibilities were considered, $`T_{KT}<T_I`$ and $`T_{KT}=T_I`$. In discussing the first case, Teitel and Jayaprakash presented the following scenario. In terms of the dual CG model, the helicity modulus $`\mathrm{{\rm Y}}`$ gets reduced from its $`T=0`$ value by fluctuations that produce dipole moments. If Ising domains of typical size $`\xi _I`$ carried a total dipole moment proportional to their size, they would drive $`\mathrm{{\rm Y}}0`$ continuously due to the diverging $`\xi _I`$ as $`TT_I`$ from below. However, in addition to these domain excitations, there are also pair excitations. The original Kosterlitz-Thouless instability criteria would imply that pairs unbind once $`\mathrm{{\rm Y}}`$ falls below the critical value $`\mathrm{{\rm Y}}(T)=2T/\pi `$, which must happen at some $`T_{KT}`$ below $`T_I`$. However, if one computes domain energies at $`T=0`$ (and presumably the same holds for domain free energies at low $`T`$), one finds that it is only the domains with vanishing total dipole moment that have energies which scale with the perimeter; i.e. the only domains which are “Ising-like” are those which carry no dipole moment and so cannot give any reduction in $`\mathrm{{\rm Y}}`$! Fortunately, once $`T`$ increases above $`T_w`$, kink-antikink pairs on the boundary of the domain are free to unbind, and the Ising domains can now acquire large dipole moments at no cost in free energy. The scenario of Teitel and Jayaprakash is now restored. This conclusion is in complete agreement with the numerical work of Olsson Olsson1 ; Olsson2 , who finds two distinct transitions at $`T_{KT}<T_I`$, and argues that the non Ising-like critical behavior claimed by some single  at $`T_I`$ is in fact an artifact of finite size effects.

The kink-antikink unbinding transition may also have implications for the non-equilibrium steady state behavior of the system when the vortices are driven by a uniform force, such as is the case for a fully frustrated Josephson junction array in an applied uniform d.c. current $`I`$. Since $`I`$ couples linearly to the total dipole moment of a domain $`p`$, the force will can can lead to an instability MT  in domain growth provided the free energy of exciting the domain scales less than linearly with $`p`$. For $`T<T_w`$, the binding of kinks to antikinks prevents large dipole moments from building up on domains. The Ising-like domains, whose free energy scales like the domain length $`\mathrm{}`$, are only those domains whose total dipole moment vanishes, and hence these domains remain stable, and the Ising-like order should persist at small drives. For $`T>T_w`$, kink-antikink pairs can unbind, resulting in domains whose dipole moment may scale at least proportional to their length $`\mathrm{}`$. In this case, when the current exceeds an amount proportional to the Ising domain surface tension, domains will become unstable to growth and the Ising-like order will be destroyed.

## Acknowledgements

We wish to thank S. E. Korshunov for his very helpful correspondence and comments on an earlier version of this manuscript. This work was supported by the Engineering Research Program of the Office of Basic Energy Sciences at the Department of Energy grant DE-FG02-89ER14017 and the Swedish Research Council Contract No. 2002-3975. Travel between Rochester and Umeå was supported by grants NSF INT-9901379 and STINT 99/976(00).