many text and figure changes

5 files changed, 223 insertions, 127 deletions

diff --git a/figs/fig_clusters_ising2d.gplot b/figs/fig_clusters_ising2d.gplot
index 61172db..6f1a75d 100644
--- a/figs/fig_clusters_ising2d.gplot
+++ b/figs/fig_clusters_ising2d.gplot
@@ -1,140 +1,127 @@
 
-set terminal epslatex size 6.834, 2.81576
-set output "fig_clusters_ising2d.tex"
-
-cc1 = "#5e81b5"
-cc2 = "#e19c24"
-cc3 = "#8fb032"
-cc4 = "#eb6235"
-cc5 = "#8778b3"
-cc6 = "#c56e1a"
-
-set palette model RGB defined ( 8 cc1, 16 cc2, 32 cc3, 64 cc4, 128 cc5, 256 cc6  )
-
-unset colorbox
-set nokey
-set format x '$10^{%T}$'
-set logscale xy
-set ylabel '$\avg{s_{\text{\sc 1c}}}L^{-\gamma/\nu}$'
-set xlabel '$hL^{\beta\delta/\nu}$'
-
-set xrange [0.000005:12000000]
-
-
-set multiplot layout 2,3
+# data and constants, sorted by column
 
-set size 0.38, 0.5
-
-dataIsing2D = "data/cluster-size/cluster-size_ising2d.dat"
-
-set yrange [0.8:8]
-
-unset xlabel
-
-set label "2D Ising" at graph 0.1, 0.75
+isingRange = "[0.0000005:5000000]"
 
+dataIsing2D = "data/correlation-times/ising-2d.dat"
 betaIsing2D = 1.0 / 8.0
 nuIsing2D = 1.0
 deltaIsing2D = 15.0
 gammaIsing2D = 7.0 / 4.0
 
-plot dataIsing2D using ($3/$2 * exp(log($1) * betaIsing2D * deltaIsing2D / nuIsing2D)):($4 * exp(-log($1) * gammaIsing2D / nuIsing2D)):($5 * exp(-log($1) * gammaIsing2D/nuIsing2D)):1 with yerrorbars pt 0 lw 2 palette, dataIsing2D using ($3/$2 * exp(log($1) * betaIsing2D * deltaIsing2D / nuIsing2D)):($4 * exp(-log($1) * gammaIsing2D / nuIsing2D)):1 with lines palette, 1.5 * exp(log(x) * 2 / deltaIsing2D) with lines lc black
-
-unset label
-set origin 0.35, 0.5
-set size 0.35, 0.5
-
-data3Potts2D = "data/cluster-size/cluster-size_3potts2d.dat"
-
-set xrange [0.000005:9000000]
-set yrange [0.4:7]
-
-set label "2D 3-State Potts" at graph 0.1, 0.75
+dataIsing3D = "data/correlation-times/ising-3d.dat"
+betaIsing3D = 0.326419
+nuIsing3D = 0.629971
+deltaIsing3D = 4.78984
+gammaIsing3D = 1.237075
 
-unset xlabel
-unset ylabel
+pottsRange="[0.000001:5000000]"
 
+data3Potts2D = "data/correlation-times/potts-3.dat"
 beta3Potts2D = 1.0 / 9.0
 nu3Potts2D = 5.0 / 6.0
 delta3Potts2D = 14.0
 gamma3Potts2D = 13.0 / 9.0
 
-plot data3Potts2D using ($3/$2 * exp(log($1) * beta3Potts2D * delta3Potts2D / nu3Potts2D)):($4 * exp(-log($1) * gamma3Potts2D / nu3Potts2D)):($5 * exp(-log($1) * gamma3Potts2D/nu3Potts2D)):1 with yerrorbars pt 0 lw 2 palette, data3Potts2D using ($3/$2 * exp(log($1) * beta3Potts2D * delta3Potts2D / nu3Potts2D)):($4 * exp(-log($1) * gamma3Potts2D / nu3Potts2D)):1 with lines palette, 1.2 * exp(log(x) * 2 / delta3Potts2D) with lines lc black
-
-unset label
-set size 0.38, 0.5
-
-data2Vector3D = "data/cluster-size/cluster-size_2vector3d.dat"
-
-set xrange [0.000005:9000000]
-set yrange [0.6:200]
+data4Potts2D = "data/correlation-times/potts-4.dat"
+beta4Potts2D = 1.0 / 12.0
+nu4Potts2D = 2.0 / 3.0
+delta4Potts2D = 15.0
+gamma4Potts2D = 7.0 / 6.0
 
-set label "3D O(2)" at graph 0.1, 0.75
+vectorRange="[0.000001:10000000]"
 
+data2Vector3D = "data/correlation-times/planar-pert.dat"
 beta2Vector3D = 0.3470
 nu2Vector3D = 0.6703
 delta2Vector3D = 4.79539
 gamma2Vector3D = 1.3169
 
-plot data2Vector3D using ($3/$2 * exp(log($1) * beta2Vector3D * delta2Vector3D / nu2Vector3D)):($4 * exp(-log($1) * gamma2Vector3D / nu2Vector3D)):($5 * exp(-log($1) * gamma2Vector3D/nu2Vector3D)):1 with yerrorbars pt 0 lw 2 palette, data2Vector3D using ($3/$2 * exp(log($1) * beta2Vector3D * delta2Vector3D / nu2Vector3D)):($4 * exp(-log($1) * gamma2Vector3D / nu2Vector3D)):1 with lines palette, 3 * exp(log(x) * 2 / delta2Vector3D) with lines lc black
+data3Vector3D = "data/correlation-times/heisenberg-pert.dat"
+beta3Vector3D = 0.3662
+nu3Vector3D = 0.7073
+delta3Vector3D = 4.79465
+gamma3Vector3D = 1.3895
 
-unset label
+# define colors
 
-set origin -0.05,0.0
-set size 0.43, 0.55
+cc1 = "#5e81b5"
+cc2 = "#e19c24"
+cc3 = "#8fb032"
+cc4 = "#eb6235"
+cc5 = "#8778b3"
+cc6 = "#c56e1a"
+cc7 = "#5d9ec7"
 
-dataIsing3D = "data/cluster-size/cluster-size_ising3d.dat"
+set palette model RGB defined ( 8 cc1, 16 cc2, 32 cc3, 64 cc4, 128 cc5, 256 cc6 , 512 cc7 )
 
-set xrange [0.000005:9000000]
-set yrange [0.9:400]
+# global plotting settings
 
-set label "3D Ising" at graph 0.1, 0.75
+set terminal epslatex size 6.834, 2.81576
+set output "fig_clusters_ising2d.tex"
 
-set ylabel '$\avg{s_{\text{\sc 1c}}}L^{-\gamma/\nu}$'
-set xlabel '$hL^{\beta\delta/\nu}$'
+set multiplot layout 2,3 margins 0.05, 0.995, 0.15, 0.99 spacing 0.05, 0.015
+set size ratio 1 / 1.61803398875
+set key autotitle columnhead
+unset colorbox
+set nokey
+set logscale xy
 
-betaIsing3D = 0.326419
-nuIsing3D = 0.629971
-deltaIsing3D = 4.78984
-gammaIsing3D = 1.237075
+scaledylabel = '"$\\avg{s_{\\text{\\sc 1c}}}L^{-\\gamma/\\nu}$" offset 2,0'
 
-plot dataIsing3D using ($3/$2 * exp(log($1) * betaIsing3D * deltaIsing3D / nuIsing3D)):($4 * exp(-log($1) * gammaIsing3D / nuIsing3D)):($5 * exp(-log($1) * gammaIsing3D/nuIsing3D)):1 with yerrorbars pt 0 lw 2 palette, dataIsing3D using ($3/$2 * exp(log($1) * betaIsing3D * deltaIsing3D / nuIsing3D)):($4 * exp(-log($1) * gammaIsing3D / nuIsing3D)):1 with lines palette, 5 * exp(log(x) * 2 / deltaIsing3D) with lines lc black
+# plotting the first row
+set format x ""
+unset xlabel
+set ylabel @scaledylabel
 
+# ploting Ising 2D
+set xrange @isingRange
+#set yrange [0.8:10]
+
+plot dataIsing2D using (column("H")/column("T") * exp(log(column("L")) * betaIsing2D * deltaIsing2D / nuIsing2D)):(column("S") * exp(-log(column("L")) * (gammaIsing2D / nuIsing2D))):(column("dS") * exp(-log(column("L")) * (gammaIsing2D / nuIsing2D))):1 with yerrorbars pt 0 lw 2 palette,\
+  dataIsing2D using (column("H")/column("T") * exp(log(column("L")) * betaIsing2D * deltaIsing2D / nuIsing2D)):(column("S") * exp(-log(column("L")) * (gammaIsing2D / nuIsing2D))):1 with lines palette
 
-unset label
 unset ylabel
-set origin 0.35, 0
-set size 0.35, 0.55
 
-data4Potts2D = "data/cluster-size/cluster-size_4potts2d.dat"
+# plotting 3 Potts 2D
+set xrange @pottsRange
+set yrange [0.5:10]
 
-set xrange [0.000005:9000000]
-set yrange [0.6:4]
+plot data3Potts2D using (column("H")/column("T") * exp(log(column("L")) * beta3Potts2D * delta3Potts2D / nu3Potts2D)):(column("S") * exp(-log(column("L")) * (gamma3Potts2D / nu3Potts2D))):(column("dS") * exp(-log(column("L")) * (gamma3Potts2D / nu3Potts2D))):1 with yerrorbars pt 0 lw 2 palette,\
+  data3Potts2D using (column("H")/column("T") * exp(log(column("L")) * beta3Potts2D * delta3Potts2D / nu3Potts2D)):(column("S") * exp(-log(column("L")) * (gamma3Potts2D / nu3Potts2D))):1 with lines palette
 
-set label "2D 4-State Potts" at graph 0.1, 0.75
+# plotting 2 Vector 3D
+set xrange @vectorRange
+set yrange [0.4:250]
 
-beta4Potts2D = 1.0 / 12.0
-nu4Potts2D = 2.0 / 3.0
-delta4Potts2D = 15.0
-gamma4Potts2D = 7.0 / 6.0
+plot data2Vector3D using (column("H")/column("T") * exp(log(column("L")) * beta2Vector3D * delta2Vector3D / nu2Vector3D)):(column("S") * exp(-log(column("L")) * (gamma2Vector3D / nu2Vector3D))):(column("dS") * exp(-log(column("L")) * (gamma2Vector3D / nu2Vector3D))):1 with yerrorbars pt 0 lw 2 palette,\
+  data2Vector3D using (column("H")/column("T") * exp(log(column("L")) * beta2Vector3D * delta2Vector3D / nu2Vector3D)):(column("S") * exp(-log(column("L")) * (gamma2Vector3D / nu2Vector3D))):1 with lines palette
 
-plot data4Potts2D using ($3/$2 * exp(log($1) * beta4Potts2D * delta4Potts2D / nu4Potts2D)):($4 * exp(-log($1) * gamma4Potts2D / nu4Potts2D)):($5 * exp(-log($1) * gamma4Potts2D/nu4Potts2D)):1 with yerrorbars pt 0 lw 2 palette, data4Potts2D using ($3/$2 * exp(log($1) * beta4Potts2D * delta4Potts2D / nu4Potts2D)):($4 * exp(-log($1) * gamma4Potts2D / nu4Potts2D)):1 with lines palette, 1.2 * exp(log(x) * 2 / delta4Potts2D) with lines lc black
+# plotting the second row
+set format x '$10^{%T}$'
+set xlabel '$hL^{\beta\delta/\nu}$'
+set ylabel @scaledylabel
 
-unset label
-set size 0.38, 0.55
+# plotting Ising 3D
+set xrange @isingRange
+set yrange [0.5:500]
 
-data3Vector3D = "data/cluster-size/cluster-size_3vector3d.dat"
+plot dataIsing3D using (column("H")/column("T") * exp(log(column("L")) * betaIsing3D * deltaIsing3D / nuIsing3D)):(column("S") * exp(-log(column("L")) * (gammaIsing3D / nuIsing3D))):(column("dS") * exp(-log(column("L")) * (gammaIsing3D / nuIsing3D))):1 with yerrorbars pt 0 lw 2 palette,\
+  dataIsing3D using (column("H")/column("T") * exp(log(column("L")) * betaIsing3D * deltaIsing3D / nuIsing3D)):(column("S") * exp(-log(column("L")) * (gammaIsing3D / nuIsing3D))):1 with lines palette
 
-set xrange [0.000005:9000000]
-set yrange [0.6:200]
+unset ylabel
 
-set label "3D O(3)" at graph 0.1, 0.75
+# plotting 4 Potts 2D
+set xrange @pottsRange
+set yrange [0.8:4]
 
-beta3Vector3D = 0.3662
-nu3Vector3D = 0.7073
-delta3Vector3D = 4.79465
-gamma3Vector3D = 1.3895
+plot data4Potts2D using (column("H")/column("T") * exp(log(column("L")) * beta4Potts2D * delta4Potts2D / nu4Potts2D)):(column("S") * exp(-log(column("L")) * (gamma4Potts2D / nu4Potts2D))):(column("dS") * exp(-log(column("L")) * (gamma4Potts2D / nu4Potts2D))):1 with yerrorbars pt 0 lw 2 palette,\
+  data4Potts2D using (column("H")/column("T") * exp(log(column("L")) * beta4Potts2D * delta4Potts2D / nu4Potts2D)):(column("S") * exp(-log(column("L")) * (gamma4Potts2D / nu4Potts2D))):1 with lines palette
+
+# plotting 3 Vector 3D
+set xrange @vectorRange
+set yrange [0.1:100]
 
-plot data3Vector3D using ($3/$2 * exp(log($1) * beta3Vector3D * delta3Vector3D / nu3Vector3D)):($4 * exp(-log($1) * gamma3Vector3D / nu3Vector3D)):($5 * exp(-log($1) * gamma3Vector3D/nu3Vector3D)):1 with yerrorbars pt 0 lw 2 palette, data3Vector3D using ($3/$2 * exp(log($1) * beta3Vector3D * delta3Vector3D / nu3Vector3D)):($4 * exp(-log($1) * gamma3Vector3D / nu3Vector3D)):1 with lines palette, 2 * exp(log(x) * 2 / delta3Vector3D) with lines lc black
+plot data3Vector3D using (column("H")/column("T") * exp(log(column("L")) * beta3Vector3D * delta3Vector3D / nu3Vector3D)):(column("S") * exp(-log(column("L")) * (gamma3Vector3D / nu3Vector3D))):(column("dS") * exp(-log(column("L")) * (gamma3Vector3D / nu3Vector3D))):1 with yerrorbars pt 0 lw 2 palette,\
+  data3Vector3D using (column("H")/column("T") * exp(log(column("L")) * beta3Vector3D * delta3Vector3D / nu3Vector3D)):(column("S") * exp(-log(column("L")) * (gamma3Vector3D / nu3Vector3D))):1 with lines palette
 
-unset multiplot
diff --git a/figs/fig_correlation-times.gplot b/figs/fig_correlation-times.gplot
index 7dac658..db62a2d 100644
--- a/figs/fig_correlation-times.gplot
+++ b/figs/fig_correlation-times.gplot
@@ -73,7 +73,7 @@ unset colorbox
 set nokey
 set logscale xy
 
-scaledylabel = '"$\\avg{s_{\\text{\\sc 1c}}}L^{-\\gamma/\\nu}$" offset 2,0'
+scaledylabel = '"$\\tau\\avg{s_{\\text{\\sc 1c}}}L^{-z}$" offset 2,0'
 
 # plotting the first row
 set format x ""
diff --git a/figs/fig_harmonic-susceptibilities.gplot b/figs/fig_harmonic-susceptibilities.gplot
index 85d474f..6129723 100644
--- a/figs/fig_harmonic-susceptibilities.gplot
+++ b/figs/fig_harmonic-susceptibilities.gplot
@@ -1,5 +1,5 @@
 
-set terminal epslatex size 3.417, 5.3
+set terminal epslatex size 3.417, 4.8
 set output "fig_harmonic-susceptibilities.tex"
 
 data4 = "data/harmonic-susceptibilities/order-4.dat"
@@ -11,12 +11,12 @@ set cbrange[STATS_min:STATS_max]
 set logscale xycb
 set nokey
 
-set xrange [1:5000]
-set yrange [0.01:10000]
+set xrange [3:5000]
+set yrange [0.005:20000]
 
 unset colorbox
 
-set multiplot layout 2,1 margins 0.15, 0.99, 0.225, 0.99 spacing 0,0.015
+set multiplot layout 2,1 margins 0.15, 0.99, 0.1, 0.88 spacing 0,0
 set size ratio 1 / 1.61803398875
 
 unset xlabel
@@ -27,9 +27,9 @@ plot data4 using "L":"X":"dX":"H" with yerrorbars pt 0 lw 2 palette,\
      data4 using "L":"X":"H"      with lines palette
 
 set xlabel '$L$' 
-set colorbox horizontal user origin 0.15, 0.1 size 0.84, 0.03
-set cblabel '(Top) $h_4$ or (bottom) $h_6$'
-set cbtics format '$10^{%T}$'
+set colorbox horizontal user origin 0.15, 0.89 size 0.84, 0.03
+set cblabel '(Top) $h_4$ or (bottom) $h_6$' offset 0,6.5
+set cbtics format '$10^{%T}$' offset 0,2.8
 set format x
 
 plot data6 using "L":"X":"dX":"H" with yerrorbars pt 0 lw 2 palette,\
diff --git a/monte-carlo.bib b/monte-carlo.bib
index 70f3f49..bdca017 100644
--- a/monte-carlo.bib
+++ b/monte-carlo.bib
@@ -562,4 +562,68 @@ random field Ising model and finally of quantum spin glasses.},
 	keywords = {Autocorrelation time, Cluster algorithm, Dynamic critical exponent, Ising model, Monte Carlo, Potts model, Swendsen–Wang algorithm},
 	pages = {259--291},
 	file = {ScienceDirect Full Text PDF:/home/pants/.zotero/data/storage/MKA8WYZZ/Ossola and Sokal - 2004 - Dynamic critical behavior of the Swendsen–Wang alg.pdf:application/pdf;ScienceDirect Snapshot:/home/pants/.zotero/data/storage/YHGX7CDT/S0550321304003098.html:text/html}
+}
+
+@article{ala-nissila_numerical_1994,
+	title = {Numerical studies of the two-dimensional {XY} model with symmetry-breaking fields},
+	volume = {50},
+	url = {https://link.aps.org/doi/10.1103/PhysRevB.50.12692},
+	doi = {10.1103/PhysRevB.50.12692},
+	abstract = {We present results of numerical studies of the two-dimensional XY model with four- and eightfold symmetry-breaking fields. This model has recently been shown to describe hydrogen-induced reconstruction on the W(100) surface. Based on mean-field and renormalization-group arguments, we first show how the interplay between the anisotropy fields can give rise to different phase transitions in the model. When the fields are compatible with each other there is a continuous phase transition when the fourth-order field is varied from negative to positive values. This transition becomes discontinuous at low temperatures. These two regimes are separated by a multicritical point. In the case of competing four- and eightfold fields, the first-order transition at low temperatures opens up into two Ising transitions. We then use numerical methods to accurately locate the position of the multicritical point, and to verify the nature of the transitions. The different techniques used include Monte Carlo histogram methods combined with finite-size scaling analysis, the real-space Monte Carlo renormalization-group method, and the Monte Carlo transfer-matrix method. Our numerical results are in good agreement with the theoretical arguments.},
+	number = {17},
+	urldate = {2018-09-25},
+	journal = {Physical Review B},
+	author = {Ala-Nissila, T. and Granato, E. and Kankaala, K. and Kosterlitz, J. M. and Ying, S.-C.},
+	month = nov,
+	year = {1994},
+	pages = {12692--12701},
+	file = {Ala-Nissila et al. - 1994 - Numerical studies of the two-dimensional XY model .pdf:/home/pants/.zotero/data/storage/2C3RYCHG/Ala-Nissila et al. - 1994 - Numerical studies of the two-dimensional XY model .pdf:application/pdf;APS Snapshot:/home/pants/.zotero/data/storage/YU2G86C5/PhysRevB.50.html:text/html}
+}
+
+@article{kankaala_theory_1993,
+	title = {Theory of adsorbate-induced surface reconstruction on {W}(100)},
+	volume = {47},
+	url = {https://link.aps.org/doi/10.1103/PhysRevB.47.2333},
+	doi = {10.1103/PhysRevB.47.2333},
+	abstract = {We report results of a theoretical study on an adsorbate-induced surface reconstruction. Hydrogen adsorption on a W(100) surface causes a switching transition in the symmetry of the displacements of the W atoms within the ordered c(2×2) phase. This transition is modeled by an effective Hamiltonian, where the hydrogen degrees of freedom are integrated out. Based on extensive Monte Carlo renormalization-group calculations we show that the switching transition is of second order at high temperatures and of first order at low temperatures. This behavior is qualitatively explained in terms of an XY model where there is an interplay between fourfold and eightfold anisotropy fields. We also compare the calculated phase diagrams with a simple mean-field theory.},
+	number = {4},
+	urldate = {2018-09-25},
+	journal = {Physical Review B},
+	author = {Kankaala, Kari and Ala-Nissila, Tapio and Ying, See-Chen},
+	month = jan,
+	year = {1993},
+	pages = {2333--2343},
+	file = {APS Snapshot:/home/pants/.zotero/data/storage/EQCT6K8H/PhysRevB.47.html:text/html;Kankaala et al. - 1993 - Theory of adsorbate-induced surface reconstruction.pdf:/home/pants/.zotero/data/storage/8DGX9H9Y/Kankaala et al. - 1993 - Theory of adsorbate-induced surface reconstruction.pdf:application/pdf}
+}
+
+@article{selinger_theory_1988,
+	title = {Theory of {Hexatic}-to-{Hexatic} {Transitions}},
+	volume = {61},
+	url = {https://link.aps.org/doi/10.1103/PhysRevLett.61.416},
+	doi = {10.1103/PhysRevLett.61.416},
+	abstract = {A theory of transitions between tilted hexatic-I and -F phases in liquid-crystal films is developed. A renormalization-group analysis leads to four tilted hexatic phases: the hexatic-I and -F phases, an intermediate hexatic-L phase, and an unlocked phase. All these phases except the unlocked phase also exist if the films are crystalline rather than hexatic. These results are consistent with recent experiments on thermotropic and lyotropic liquid crystals.},
+	number = {4},
+	urldate = {2018-09-25},
+	journal = {Physical Review Letters},
+	author = {Selinger, Jonathan V. and Nelson, David R.},
+	month = jul,
+	year = {1988},
+	pages = {416--419},
+	file = {APS Snapshot:/home/pants/.zotero/data/storage/59CAGAAM/PhysRevLett.61.html:text/html;Selinger and Nelson - 1988 - Theory of Hexatic-to-Hexatic Transitions.pdf:/home/pants/.zotero/data/storage/EHI5X23X/Selinger and Nelson - 1988 - Theory of Hexatic-to-Hexatic Transitions.pdf:application/pdf}
+}
+
+@article{dierker_consequences_1986,
+	title = {Consequences of {Bond}-{Orientational} {Order} on the {Macroscopic} {Orientation} {Patterns} of {Thin} {Tilted} {Hexatic} {Liquid}-{Crystal} {Films}},
+	volume = {56},
+	url = {https://link.aps.org/doi/10.1103/PhysRevLett.56.1819},
+	doi = {10.1103/PhysRevLett.56.1819},
+	abstract = {By use of depolarized laser relection microscopy, unique textural defects have been observed in thin tilted hexatic liquid-crystal films which result from their bond-orientational order. The temperature dependence of the bond-orientational elasticity can be deduced from these textural patterns and is in qualitative agreement with the theory of two-dimensional melting. The effects of molecular chirality and film thickness on the defect textures is also reported along with observations on the polygonization of the two-dimensional crystalline phase.},
+	number = {17},
+	urldate = {2018-09-25},
+	journal = {Physical Review Letters},
+	author = {Dierker, S. B. and Pindak, R. and Meyer, R. B.},
+	month = apr,
+	year = {1986},
+	pages = {1819--1822},
+	file = {APS Snapshot:/home/pants/.zotero/data/storage/A8T7IR63/PhysRevLett.56.html:text/html;Dierker et al. - 1986 - Consequences of Bond-Orientational Order on the Ma.pdf:/home/pants/.zotero/data/storage/9XEH86XP/Dierker et al. - 1986 - Consequences of Bond-Orientational Order on the Ma.pdf:application/pdf}
 }
\ No newline at end of file
diff --git a/monte-carlo.tex b/monte-carlo.tex
index 370ac0f..b24fa9d 100644
--- a/monte-carlo.tex
+++ b/monte-carlo.tex
@@ -89,7 +89,7 @@ extra degree of freedom, allowing the method to be used in a subcategory of
 interesting fields \cite{alexandrowicz_swendsen-wang_1989, wang_clusters_1989,
 ray_metastability_1990}. Static fields have also been applied by including a
 separate metropolis or heat bath update step after cluster formation
-\cite{destri_swendsen-wang_1992, lauwers_critical_1989}, and other categories
+\cite{destri_swendsen-wang_1992, lauwers_critical_1989, ala-nissila_numerical_1994}, and other categories
 of fields have been applied using replica methods
 \cite{redner_graphical_1998,chayes_graphical_1998,machta_replica-exchange_2000}.
 We show that the scaling of correlation time near the critical point of
@@ -192,9 +192,12 @@ in the following way.
 \begin{enumerate}
   \item Pick a random site $m_0$ and add it to the stack.
 
-  \item Pick a random rotation $r\in R$ of order two from the set of such
-    rotations with a probability distribution of the form
-    $P_{m_0}(r)=f(\J(m_0,r\cdot m_0))$.
+  \item Select a rotation $r\in R$ of order two. For detailed balance it is
+    sufficient to ensure that probability distribution $f_{m_0}(r\mid \set s)$ the
+    rotation is sampled from depend only on $Z(s,r\cdot s)$ for every site,
+    since these numbers are invariant under a cluster flip, e.g., $Z(s',r\cdot
+    s')=Z(r\cdot s,r\cdot(r\cdot s))=Z(r\cdot s,(rr)\cdot s)=Z(r\cdot
+    s,s)=Z(s,r\cdot s)$.
   \item While the stack isn't empty,
     \begin{enumerate}
       \item pop site $m$ from the stack.
@@ -204,7 +207,7 @@ in the following way.
           \item For every $j$ such that $\{m,j\}\in E$, add site $j$ to the
             stack with probability
             \[
-              p_r(s_m,s_j)=\min\{0,1-e^{\beta(\J(r\cdot s_m,s_j)-\J(s_m,s_j))}\}.
+              p_r(s_m,s_j\mid q)=\min\{0,1-qe^{\beta(\J(r\cdot s_m,s_j)-\J(s_m,s_j))}\}.
             \]
           \item Take $s_m\mapsto r\cdot s_m$.
       \end{enumerate}
@@ -225,12 +228,12 @@ C\subset E$ is related to the probability of the reverse process $P(\set{s'}\to\
 \[
   \begin{aligned}
     \frac{P(\set s\to\set{s'})}{P(\set{s'}\to\set s)}
-    &=\frac{P_{m_0}(r)}{P_{m_0}(r^{-1})}\prod_{\{i,j\}\in
+    &=\frac{f_{\set s}(r\mid m_0)}{f_{\set s'}(r^{-1}\mid m_0)}\prod_{\{i,j\}\in
   C}\frac{p_r(s_i,s_j)}{p_{r^{-1}}(s_i',s_j')}\prod_{\{i,j\}\in\partial
   C}\frac{1-p_r(s_i,s_j)}{1-p_{r^{-1}}(s'_i,s'_j)}\\
   &=\prod_{\{i,j\}\in\partial
   C}e^{\beta(\J(r\cdot s_i,s_j)-\J(s_i,s_j))}
-  =\frac{p_r(s_i,s_j)}{p_{r}(s_i,s_j)}\frac{e^{-\beta\H(\set
+  =\frac{e^{-\beta\H(\set
   s)}}{e^{-\beta\H(\set{s'})}},
 \end{aligned}
 \]
@@ -340,6 +343,7 @@ representation, this form of the Hamiltonian might be considered the `ghost
 transformation' representation.
 
 \section{Examples}
+\label{sec:examples}
 
 Several specific examples from Table~\ref{table:models} are described in the
 following.
@@ -357,7 +361,10 @@ wang_clusters_1989, ray_metastability_1990}. The algorithm has been
 implemented by one of the authors in an existing interactive Ising
 simulator at \texttt{https://mattbierbaum.github.io/ising.js} \cite{bierbaum_ising.js_nodate}.
 
-\subsection{The $\mathrm O(n)$ model} In the $\mathrm O(n)$ model spins are described by vectors on the
+\subsection{The $\mathrm O(n)$ model}
+\label{sec:examples:on}
+
+In the $\mathrm O(n)$ model spins are described by vectors on the
 $(n-1)$-sphere $S^{n-1}$. Its symmetry group is $O(n)$, $n\times n$ orthogonal
 matrices, which act on the spins by matrix multiplication. The elements of
 $O(n)$ of order two are reflections about hyperplanes through the origin and
@@ -394,6 +401,11 @@ of large $L$ the expected square angle between neighbors being
 \[
   \avg{\theta^2}\simeq\frac{(n-1)T}{D+H/2}
 \]
+We shall in in the numeric experiments of Section \ref{sec:performance} that
+this choice ameliorates the problem but probably is not the best. A more
+subtle technique---for instance by matching statistics of the effective
+Gaussian model that results in these circumstances to the cluster
+statistics---may result in better performance.
 
 \subsection{The Potts model} In the $q$-state Potts model spins are described
 by elements of $\{1,\ldots,q\}$. Its symmetry group is the symmetric group
@@ -428,6 +440,7 @@ of the system.  A variant of the algorithm has been applied without a field
 \cite{evertz_stochastic_1991}.
 
 \section{Performance}
+\label{sec:performance}
 
 No algorithm is worthwhile if it doesn't run efficiently. This algorithm,
 being an extension of the Wolff algorithm into a new domain, should be
@@ -454,24 +467,19 @@ one should expect its natural extension in the presence of a field to scale
 roughly like $h^{-z\nu/\beta\delta}$ and collapse appropriately as a function
 of $hL^{\beta\delta/\nu}$.  We measured the autocorrelation time for the $D=2$
 square-lattice model at a variety of system sizes, temperatures, and fields
-$B(s)=hs/\beta$ using standard methods \cite{geyer_practical_1992}. The
+$B(s)=hs/\beta$ using standard methods \cite{ossola_dynamic_2004}. The
 resulting scaling behavior, plotted in
 Fig.~\ref{fig:correlation_time-collapse}, is indeed consistent with an
 extension to finite field of the behavior at zero field.
 
-\begin{figure}
-  \centering
-  \input{fig_correlation_collapse-hL}
-  \caption{Collapse of the correlation time $\tau$ of the 2D square lattice
-    Ising model along the critical isotherm at various systems sizes
-    $N=L\times L$ for $L=8$, $16$, $32$, $64$, $128$, and $256$ as a function
-    of the renormalization invariant $hL^{\beta\delta/\nu}$. The exponent
-    $z=0.30$ is taken from recent measurements at zero field
-  \cite{liu_dynamic_2014}. The solid black line shows a plot of
-  $(hL^{\beta\delta/\nu})^{-z\nu/\beta\delta}$.
+\begin{figure*}
+  \include{fig_correlation-times}
+  \caption{
+    Correlation times $\tau$ scaled by the average cluster size as a function
+    of external field for various models at various system sizes. 
   }
   \label{fig:correlation_time-collapse}
-\end{figure}
+\end{figure*}
 
 Since the formation and flipping of clusters is the hallmark of Wolff
 dynamics, another way to ensure that the dynamics with field scale like those
@@ -538,10 +546,50 @@ conjecture.
   \include{fig_generator-times}
 \end{figure}
 
+\section{Applying Nonlinear Fields to the xy Model}
+
+This far our numeric work has quantified the performance of existing
+techniques. Here, we apply our general framework in a new way:
+harmonic perturbations to the low-temperature {\sc xy}, or \twodee O(2),
+model. We consider fields of the form $B_n(s)=h_n\cos(n\theta(s))$, where
+$\theta$ is the angle made between $s$ and the $x$-axis, say. Corrections of
+these types are expected to appear in realistic models of systems na\"ively
+expected to exhibit Kosterlitz--Thouless critical behavior due to the presence
+of the lattice or substrate. Whether these fields are relevant or irrelevant
+in the renormalization group sense determines whether those systems spoil or
+admit that critical behaviour. Among many fascinating
+\cite{jose_renormalization_1977, kankaala_theory_1993,
+ala-nissila_numerical_1994, dierker_consequences_1986, selinger_theory_1988}
+results that emerge from systems with one or more of these fields applied,
+it is predicted that $h_4$ is relevant while $h_6$ is not at some
+sufficiently high temperatures below the Kosterlitz--Thouless point
+\cite{jose_renormalization_1977}.
+
 \begin{figure}
   \include{fig_harmonic-susceptibilities}
+  \caption{Susceptibilities as a function of system size for a \twodee O(2)
+    model at $T=0.7$ and with (top) fourfold symmetric and (bottom) eightfold
+    symmetric perturbing fields. Different field strengths are shown in
+    different colors.
+  }
+  \label{fig:harmonic-susceptibilities}
 \end{figure}
 
+We made a basic investigation of this result using our algorithm. Since we ran
+the algorithm at fairly high fields we did not choose reflections though the
+origin uniformly. Instead, we choose the planes of reflection first by
+rotating our starting spin by $2\pi m/n$ degrees for $m$ uniformly taken from
+$1,\ldots,n$ and generating a normal to the plane from that direction as
+described in Section \ref{sec:examples:on}. The resulting susceptibilities as
+a function of system size are shown for various field strengths in
+Fig.~\ref{fig:harmonic-susceptibilities}. In the fourfold case, for each field
+strength there is a system size at which the divergence in the susceptibility
+is cut off, while for the sixfold case we measured no such cutoff, even up to
+strong fields. This confirms the expected result, that even in a strong field
+the sixfold perturbations preserve the critical behavior.
+
+\section{Conclusions}
+
 We have taken several disparate extensions of cluster methods to spin models
 in an external field and generalized them to work for any model of a broad
 class.  The resulting representation involves the introduction of not a ghost
@@ -560,7 +608,7 @@ perturbations on spin models can be tested numerically
 \cite{jose_renormalization_1977, blankschtein_fluctuation-induced_1982,
 bruce_coupled_1975, manuel_carmona_$n$-component_2000}.
 
-<<<<<<< HEAD
+
 \appendix
 
 \section{$\mathrm O(n)$ model at high field}
@@ -629,9 +677,6 @@ whence
 \end{align}
 
 \section{Calculating autocorrelation time}
-\begin{figure*}
-  \include{fig_correlation-times}
-\end{figure*}
 
 \begin{acknowledgments}
   This work was supported by NSF grant NSF DMR-1719490.