Simulated Annealing for Topological Solitons
Mark Hale
and
Oliver Schwindt and Tom Weidig
{Dept.~of Physics, UMIST, Manchester M60 1QD, UK}
The search for solutions of field theories allowing for topological solitons requires that we find the field configuration with the lowest energy in a given sector of topological charge. The standard approach is based on the numerical solution of the static Euler-Lagrange differential equation following from the field energy. As an alternative, we propose to use a simulated annealing algorithm to minimize the energy functional directly. We have applied simulated annealing to several nonlinear classical field theories: the sine-Gordon model in one dimension, the baby Skyrme model in two dimensions and the nuclear Skyrme model in three dimensions. We describe in detail the implementation of the simulated annealing algorithm, present our results and get independent confirmation of the studies which have used standard minimization techniques.
Cooling1D.tex
% GNUPLOT: LaTeX picture
\setlength{\unitlength}{0.240900pt}
\ifx\plotpoint\undefined\newsavebox{\plotpoint}\fi
\sbox{\plotpoint}{\rule[-0.200pt]{0.400pt}{0.400pt}}%
\begin{picture}(1500,900)(0,0)
\font\gnuplot=cmr10 at 10pt
\gnuplot
\sbox{\plotpoint}{\rule[-0.200pt]{0.400pt}{0.400pt}}%
\put(100.0,82.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,82){\makebox(0,0)[r]{-1}}
\put(1419.0,82.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,160.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,160){\makebox(0,0)[r]{0}}
\put(1419.0,160.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,238.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,238){\makebox(0,0)[r]{1}}
\put(1419.0,238.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,315.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,315){\makebox(0,0)[r]{2}}
\put(1419.0,315.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,393.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,393){\makebox(0,0)[r]{3}}
\put(1419.0,393.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,471.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,471){\makebox(0,0)[r]{4}}
\put(1419.0,471.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,549.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,549){\makebox(0,0)[r]{5}}
\put(1419.0,549.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,627.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,627){\makebox(0,0)[r]{6}}
\put(1419.0,627.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,704.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,704){\makebox(0,0)[r]{7}}
\put(1419.0,704.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,782.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,782){\makebox(0,0)[r]{8}}
\put(1419.0,782.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,860.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,860){\makebox(0,0)[r]{9}}
\put(1419.0,860.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(100,41){\makebox(0,0){0}}
\put(100.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(291.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(291,41){\makebox(0,0){2}}
\put(291.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(483.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(483,41){\makebox(0,0){4}}
\put(483.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(674.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(674,41){\makebox(0,0){6}}
\put(674.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(865.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(865,41){\makebox(0,0){8}}
\put(865.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1056.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1056,41){\makebox(0,0){10}}
\put(1056.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1248.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1248,41){\makebox(0,0){12}}
\put(1248.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1439.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1439,41){\makebox(0,0){14}}
\put(1439.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(100.0,82.0){\rule[-0.200pt]{322.565pt}{0.400pt}}
\put(1439.0,82.0){\rule[-0.200pt]{0.400pt}{187.420pt}}
\put(100.0,860.0){\rule[-0.200pt]{322.565pt}{0.400pt}}
\put(100.0,82.0){\rule[-0.200pt]{0.400pt}{187.420pt}}
\put(1279,820){\makebox(0,0)[r]{Minimized Configuration}}
\put(800,0){\makebox(0,0)[r]{x}}
\put(120,900){\makebox(0,0)[r]{$\phi(x)$}}
\put(1299.0,820.0){\rule[-0.200pt]{24.090pt}{0.400pt}}
\put(100,160){\usebox{\plotpoint}}
\put(194,159.67){\rule{3.132pt}{0.400pt}}
\multiput(194.00,159.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(100.0,160.0){\rule[-0.200pt]{22.645pt}{0.400pt}}
\put(261,160.67){\rule{3.132pt}{0.400pt}}
\multiput(261.00,160.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(207.0,161.0){\rule[-0.200pt]{13.009pt}{0.400pt}}
\put(301,161.67){\rule{3.132pt}{0.400pt}}
\multiput(301.00,161.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(274.0,162.0){\rule[-0.200pt]{6.504pt}{0.400pt}}
\put(328,162.67){\rule{3.132pt}{0.400pt}}
\multiput(328.00,162.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(314.0,163.0){\rule[-0.200pt]{3.373pt}{0.400pt}}
\put(354,163.67){\rule{3.373pt}{0.400pt}}
\multiput(354.00,163.17)(7.000,1.000){2}{\rule{1.686pt}{0.400pt}}
\put(368,164.67){\rule{3.132pt}{0.400pt}}
\multiput(368.00,164.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(341.0,164.0){\rule[-0.200pt]{3.132pt}{0.400pt}}
\put(395,165.67){\rule{3.132pt}{0.400pt}}
\multiput(395.00,165.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(408,167.17){\rule{2.700pt}{0.400pt}}
\multiput(408.00,166.17)(7.396,2.000){2}{\rule{1.350pt}{0.400pt}}
\put(421,168.67){\rule{3.373pt}{0.400pt}}
\multiput(421.00,168.17)(7.000,1.000){2}{\rule{1.686pt}{0.400pt}}
\put(435,169.67){\rule{3.132pt}{0.400pt}}
\multiput(435.00,169.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(448,171.17){\rule{2.900pt}{0.400pt}}
\multiput(448.00,170.17)(7.981,2.000){2}{\rule{1.450pt}{0.400pt}}
\put(462,173.17){\rule{2.700pt}{0.400pt}}
\multiput(462.00,172.17)(7.396,2.000){2}{\rule{1.350pt}{0.400pt}}
\put(475,175.17){\rule{2.700pt}{0.400pt}}
\multiput(475.00,174.17)(7.396,2.000){2}{\rule{1.350pt}{0.400pt}}
\multiput(488.00,177.61)(2.918,0.447){3}{\rule{1.967pt}{0.108pt}}
\multiput(488.00,176.17)(9.918,3.000){2}{\rule{0.983pt}{0.400pt}}
\multiput(502.00,180.61)(2.695,0.447){3}{\rule{1.833pt}{0.108pt}}
\multiput(502.00,179.17)(9.195,3.000){2}{\rule{0.917pt}{0.400pt}}
\multiput(515.00,183.60)(1.797,0.468){5}{\rule{1.400pt}{0.113pt}}
\multiput(515.00,182.17)(10.094,4.000){2}{\rule{0.700pt}{0.400pt}}
\multiput(528.00,187.60)(1.943,0.468){5}{\rule{1.500pt}{0.113pt}}
\multiput(528.00,186.17)(10.887,4.000){2}{\rule{0.750pt}{0.400pt}}
\multiput(542.00,191.60)(1.797,0.468){5}{\rule{1.400pt}{0.113pt}}
\multiput(542.00,190.17)(10.094,4.000){2}{\rule{0.700pt}{0.400pt}}
\multiput(555.00,195.59)(1.489,0.477){7}{\rule{1.220pt}{0.115pt}}
\multiput(555.00,194.17)(11.468,5.000){2}{\rule{0.610pt}{0.400pt}}
\multiput(569.00,200.59)(1.123,0.482){9}{\rule{0.967pt}{0.116pt}}
\multiput(569.00,199.17)(10.994,6.000){2}{\rule{0.483pt}{0.400pt}}
\multiput(582.00,206.59)(0.950,0.485){11}{\rule{0.843pt}{0.117pt}}
\multiput(582.00,205.17)(11.251,7.000){2}{\rule{0.421pt}{0.400pt}}
\multiput(595.00,213.59)(0.890,0.488){13}{\rule{0.800pt}{0.117pt}}
\multiput(595.00,212.17)(12.340,8.000){2}{\rule{0.400pt}{0.400pt}}
\multiput(609.00,221.59)(0.728,0.489){15}{\rule{0.678pt}{0.118pt}}
\multiput(609.00,220.17)(11.593,9.000){2}{\rule{0.339pt}{0.400pt}}
\multiput(622.00,230.58)(0.704,0.491){17}{\rule{0.660pt}{0.118pt}}
\multiput(622.00,229.17)(12.630,10.000){2}{\rule{0.330pt}{0.400pt}}
\multiput(636.00,240.58)(0.539,0.492){21}{\rule{0.533pt}{0.119pt}}
\multiput(636.00,239.17)(11.893,12.000){2}{\rule{0.267pt}{0.400pt}}
\multiput(649.00,252.58)(0.497,0.493){23}{\rule{0.500pt}{0.119pt}}
\multiput(649.00,251.17)(11.962,13.000){2}{\rule{0.250pt}{0.400pt}}
\multiput(662.00,265.58)(0.497,0.494){25}{\rule{0.500pt}{0.119pt}}
\multiput(662.00,264.17)(12.962,14.000){2}{\rule{0.250pt}{0.400pt}}
\multiput(676.58,279.00)(0.493,0.616){23}{\rule{0.119pt}{0.592pt}}
\multiput(675.17,279.00)(13.000,14.771){2}{\rule{0.400pt}{0.296pt}}
\multiput(689.58,295.00)(0.494,0.607){25}{\rule{0.119pt}{0.586pt}}
\multiput(688.17,295.00)(14.000,15.784){2}{\rule{0.400pt}{0.293pt}}
\multiput(703.58,312.00)(0.493,0.734){23}{\rule{0.119pt}{0.685pt}}
\multiput(702.17,312.00)(13.000,17.579){2}{\rule{0.400pt}{0.342pt}}
\multiput(716.58,331.00)(0.493,0.774){23}{\rule{0.119pt}{0.715pt}}
\multiput(715.17,331.00)(13.000,18.515){2}{\rule{0.400pt}{0.358pt}}
\multiput(729.58,351.00)(0.494,0.754){25}{\rule{0.119pt}{0.700pt}}
\multiput(728.17,351.00)(14.000,19.547){2}{\rule{0.400pt}{0.350pt}}
\multiput(743.58,372.00)(0.493,0.853){23}{\rule{0.119pt}{0.777pt}}
\multiput(742.17,372.00)(13.000,20.387){2}{\rule{0.400pt}{0.388pt}}
\multiput(756.58,394.00)(0.494,0.754){25}{\rule{0.119pt}{0.700pt}}
\multiput(755.17,394.00)(14.000,19.547){2}{\rule{0.400pt}{0.350pt}}
\multiput(770.58,415.00)(0.493,0.853){23}{\rule{0.119pt}{0.777pt}}
\multiput(769.17,415.00)(13.000,20.387){2}{\rule{0.400pt}{0.388pt}}
\multiput(783.58,437.00)(0.493,0.814){23}{\rule{0.119pt}{0.746pt}}
\multiput(782.17,437.00)(13.000,19.451){2}{\rule{0.400pt}{0.373pt}}
\multiput(796.58,458.00)(0.494,0.717){25}{\rule{0.119pt}{0.671pt}}
\multiput(795.17,458.00)(14.000,18.606){2}{\rule{0.400pt}{0.336pt}}
\multiput(810.58,478.00)(0.493,0.734){23}{\rule{0.119pt}{0.685pt}}
\multiput(809.17,478.00)(13.000,17.579){2}{\rule{0.400pt}{0.342pt}}
\multiput(823.58,497.00)(0.493,0.655){23}{\rule{0.119pt}{0.623pt}}
\multiput(822.17,497.00)(13.000,15.707){2}{\rule{0.400pt}{0.312pt}}
\multiput(836.58,514.00)(0.494,0.570){25}{\rule{0.119pt}{0.557pt}}
\multiput(835.17,514.00)(14.000,14.844){2}{\rule{0.400pt}{0.279pt}}
\multiput(850.58,530.00)(0.493,0.536){23}{\rule{0.119pt}{0.531pt}}
\multiput(849.17,530.00)(13.000,12.898){2}{\rule{0.400pt}{0.265pt}}
\multiput(863.00,544.58)(0.536,0.493){23}{\rule{0.531pt}{0.119pt}}
\multiput(863.00,543.17)(12.898,13.000){2}{\rule{0.265pt}{0.400pt}}
\multiput(877.00,557.58)(0.539,0.492){21}{\rule{0.533pt}{0.119pt}}
\multiput(877.00,556.17)(11.893,12.000){2}{\rule{0.267pt}{0.400pt}}
\multiput(890.00,569.58)(0.652,0.491){17}{\rule{0.620pt}{0.118pt}}
\multiput(890.00,568.17)(11.713,10.000){2}{\rule{0.310pt}{0.400pt}}
\multiput(903.00,579.59)(0.786,0.489){15}{\rule{0.722pt}{0.118pt}}
\multiput(903.00,578.17)(12.501,9.000){2}{\rule{0.361pt}{0.400pt}}
\multiput(917.00,588.59)(0.950,0.485){11}{\rule{0.843pt}{0.117pt}}
\multiput(917.00,587.17)(11.251,7.000){2}{\rule{0.421pt}{0.400pt}}
\multiput(930.00,595.59)(1.026,0.485){11}{\rule{0.900pt}{0.117pt}}
\multiput(930.00,594.17)(12.132,7.000){2}{\rule{0.450pt}{0.400pt}}
\multiput(944.00,602.59)(1.123,0.482){9}{\rule{0.967pt}{0.116pt}}
\multiput(944.00,601.17)(10.994,6.000){2}{\rule{0.483pt}{0.400pt}}
\multiput(957.00,608.59)(1.378,0.477){7}{\rule{1.140pt}{0.115pt}}
\multiput(957.00,607.17)(10.634,5.000){2}{\rule{0.570pt}{0.400pt}}
\multiput(970.00,613.59)(1.489,0.477){7}{\rule{1.220pt}{0.115pt}}
\multiput(970.00,612.17)(11.468,5.000){2}{\rule{0.610pt}{0.400pt}}
\multiput(984.00,618.60)(1.797,0.468){5}{\rule{1.400pt}{0.113pt}}
\multiput(984.00,617.17)(10.094,4.000){2}{\rule{0.700pt}{0.400pt}}
\multiput(997.00,622.61)(2.918,0.447){3}{\rule{1.967pt}{0.108pt}}
\multiput(997.00,621.17)(9.918,3.000){2}{\rule{0.983pt}{0.400pt}}
\multiput(1011.00,625.61)(2.695,0.447){3}{\rule{1.833pt}{0.108pt}}
\multiput(1011.00,624.17)(9.195,3.000){2}{\rule{0.917pt}{0.400pt}}
\multiput(1024.00,628.61)(2.695,0.447){3}{\rule{1.833pt}{0.108pt}}
\multiput(1024.00,627.17)(9.195,3.000){2}{\rule{0.917pt}{0.400pt}}
\put(1037,631.17){\rule{2.900pt}{0.400pt}}
\multiput(1037.00,630.17)(7.981,2.000){2}{\rule{1.450pt}{0.400pt}}
\put(1051,633.17){\rule{2.700pt}{0.400pt}}
\multiput(1051.00,632.17)(7.396,2.000){2}{\rule{1.350pt}{0.400pt}}
\put(1064,635.17){\rule{2.700pt}{0.400pt}}
\multiput(1064.00,634.17)(7.396,2.000){2}{\rule{1.350pt}{0.400pt}}
\put(1077,637.17){\rule{2.900pt}{0.400pt}}
\multiput(1077.00,636.17)(7.981,2.000){2}{\rule{1.450pt}{0.400pt}}
\put(1091,638.67){\rule{3.132pt}{0.400pt}}
\multiput(1091.00,638.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(1104,639.67){\rule{3.373pt}{0.400pt}}
\multiput(1104.00,639.17)(7.000,1.000){2}{\rule{1.686pt}{0.400pt}}
\put(1118,640.67){\rule{3.132pt}{0.400pt}}
\multiput(1118.00,640.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(1131,641.67){\rule{3.132pt}{0.400pt}}
\multiput(1131.00,641.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(1144,642.67){\rule{3.373pt}{0.400pt}}
\multiput(1144.00,642.17)(7.000,1.000){2}{\rule{1.686pt}{0.400pt}}
\put(381.0,166.0){\rule[-0.200pt]{3.373pt}{0.400pt}}
\put(1171,643.67){\rule{3.373pt}{0.400pt}}
\multiput(1171.00,643.17)(7.000,1.000){2}{\rule{1.686pt}{0.400pt}}
\put(1158.0,644.0){\rule[-0.200pt]{3.132pt}{0.400pt}}
\put(1198,644.67){\rule{3.132pt}{0.400pt}}
\multiput(1198.00,644.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(1185.0,645.0){\rule[-0.200pt]{3.132pt}{0.400pt}}
\put(1225,645.67){\rule{3.132pt}{0.400pt}}
\multiput(1225.00,645.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(1211.0,646.0){\rule[-0.200pt]{3.373pt}{0.400pt}}
\put(1292,646.67){\rule{3.132pt}{0.400pt}}
\multiput(1292.00,646.17)(6.500,1.000){2}{\rule{1.566pt}{0.400pt}}
\put(1238.0,647.0){\rule[-0.200pt]{13.009pt}{0.400pt}}
\put(1412,647.67){\rule{3.373pt}{0.400pt}}
\multiput(1412.00,647.17)(7.000,1.000){2}{\rule{1.686pt}{0.400pt}}
\put(1305.0,648.0){\rule[-0.200pt]{25.776pt}{0.400pt}}
\put(1426.0,649.0){\rule[-0.200pt]{3.132pt}{0.400pt}}
\put(1279,779){\makebox(0,0)[r]{Initial Configuration}}
\multiput(1299,779)(20.756,0.000){5}{\usebox{\plotpoint}}
\put(1399,779){\usebox{\plotpoint}}
\put(100,160){\usebox{\plotpoint}}
\put(100.00,160.00){\usebox{\plotpoint}}
\put(119.43,167.30){\usebox{\plotpoint}}
\put(139.15,173.74){\usebox{\plotpoint}}
\put(158.67,180.80){\usebox{\plotpoint}}
\put(178.04,188.25){\usebox{\plotpoint}}
\put(197.54,195.36){\usebox{\plotpoint}}
\put(217.00,202.57){\usebox{\plotpoint}}
\put(236.41,209.93){\usebox{\plotpoint}}
\put(256.05,216.58){\usebox{\plotpoint}}
\put(275.56,223.60){\usebox{\plotpoint}}
\put(295.01,230.86){\usebox{\plotpoint}}
\put(314.44,238.16){\usebox{\plotpoint}}
\put(333.93,245.28){\usebox{\plotpoint}}
\put(353.30,252.73){\usebox{\plotpoint}}
\put(372.91,259.51){\usebox{\plotpoint}}
\put(392.58,266.14){\usebox{\plotpoint}}
\put(411.97,273.53){\usebox{\plotpoint}}
\put(431.44,280.73){\usebox{\plotpoint}}
\put(450.87,288.02){\usebox{\plotpoint}}
\put(470.34,295.21){\usebox{\plotpoint}}
\put(489.76,302.50){\usebox{\plotpoint}}
\put(509.49,308.88){\usebox{\plotpoint}}
\put(528.87,316.31){\usebox{\plotpoint}}
\put(548.36,323.45){\usebox{\plotpoint}}
\put(567.85,330.59){\usebox{\plotpoint}}
\put(587.23,338.01){\usebox{\plotpoint}}
\put(606.71,345.18){\usebox{\plotpoint}}
\put(626.44,351.59){\usebox{\plotpoint}}
\put(645.90,358.81){\usebox{\plotpoint}}
\put(665.30,366.18){\usebox{\plotpoint}}
\put(684.77,373.37){\usebox{\plotpoint}}
\put(704.27,380.49){\usebox{\plotpoint}}
\put(723.64,387.94){\usebox{\plotpoint}}
\put(743.42,394.16){\usebox{\plotpoint}}
\put(762.85,401.45){\usebox{\plotpoint}}
\put(782.29,408.73){\usebox{\plotpoint}}
\put(801.71,416.04){\usebox{\plotpoint}}
\put(821.16,423.29){\usebox{\plotpoint}}
\put(840.57,430.63){\usebox{\plotpoint}}
\put(860.27,437.16){\usebox{\plotpoint}}
\put(879.83,444.09){\usebox{\plotpoint}}
\put(899.20,451.54){\usebox{\plotpoint}}
\put(918.70,458.65){\usebox{\plotpoint}}
\put(938.14,465.91){\usebox{\plotpoint}}
\put(957.57,473.22){\usebox{\plotpoint}}
\put(977.15,480.04){\usebox{\plotpoint}}
\put(996.72,486.89){\usebox{\plotpoint}}
\put(1016.22,494.01){\usebox{\plotpoint}}
\put(1035.59,501.46){\usebox{\plotpoint}}
\put(1055.09,508.57){\usebox{\plotpoint}}
\put(1074.46,516.02){\usebox{\plotpoint}}
\put(1094.03,522.93){\usebox{\plotpoint}}
\put(1113.72,529.47){\usebox{\plotpoint}}
\put(1133.13,536.82){\usebox{\plotpoint}}
\put(1152.58,544.06){\usebox{\plotpoint}}
\put(1172.01,551.36){\usebox{\plotpoint}}
\put(1191.50,558.50){\usebox{\plotpoint}}
\put(1210.87,565.95){\usebox{\plotpoint}}
\put(1230.65,572.17){\usebox{\plotpoint}}
\put(1250.13,579.33){\usebox{\plotpoint}}
\put(1269.52,586.74){\usebox{\plotpoint}}
\put(1288.99,593.93){\usebox{\plotpoint}}
\put(1308.39,601.30){\usebox{\plotpoint}}
\put(1327.85,608.52){\usebox{\plotpoint}}
\put(1347.59,614.92){\usebox{\plotpoint}}
\put(1367.06,622.10){\usebox{\plotpoint}}
\put(1386.44,629.52){\usebox{\plotpoint}}
\put(1405.93,636.66){\usebox{\plotpoint}}
\put(1425.42,643.79){\usebox{\plotpoint}}
\put(1439,649){\usebox{\plotpoint}}
\sbox{\plotpoint}{\rule[-0.400pt]{0.800pt}{0.800pt}}%
\put(1279,738){\makebox(0,0)[r]{Heated Configuration}}
\put(1299.0,738.0){\rule[-0.400pt]{24.090pt}{0.800pt}}
\put(100,160){\usebox{\plotpoint}}
\multiput(101.41,150.23)(0.509,-1.402){19}{\rule{0.123pt}{2.354pt}}
\multiput(98.34,155.11)(13.000,-30.114){2}{\rule{0.800pt}{1.177pt}}
\multiput(113.00,123.08)(0.533,-0.509){19}{\rule{1.062pt}{0.123pt}}
\multiput(113.00,123.34)(11.797,-13.000){2}{\rule{0.531pt}{0.800pt}}
\multiput(128.41,107.08)(0.509,-0.616){19}{\rule{0.123pt}{1.185pt}}
\multiput(125.34,109.54)(13.000,-13.541){2}{\rule{0.800pt}{0.592pt}}
\multiput(140.00,97.41)(0.533,0.509){19}{\rule{1.062pt}{0.123pt}}
\multiput(140.00,94.34)(11.797,13.000){2}{\rule{0.531pt}{0.800pt}}
\multiput(155.41,109.00)(0.509,0.698){19}{\rule{0.123pt}{1.308pt}}
\multiput(152.34,109.00)(13.000,15.286){2}{\rule{0.800pt}{0.654pt}}
\multiput(168.41,127.00)(0.509,2.643){19}{\rule{0.123pt}{4.200pt}}
\multiput(165.34,127.00)(13.000,56.283){2}{\rule{0.800pt}{2.100pt}}
\multiput(180.00,190.09)(0.492,-0.509){21}{\rule{1.000pt}{0.123pt}}
\multiput(180.00,190.34)(11.924,-14.000){2}{\rule{0.500pt}{0.800pt}}
\multiput(195.41,178.00)(0.509,1.195){19}{\rule{0.123pt}{2.046pt}}
\multiput(192.34,178.00)(13.000,25.753){2}{\rule{0.800pt}{1.023pt}}
\multiput(208.41,201.24)(0.509,-0.912){21}{\rule{0.123pt}{1.629pt}}
\multiput(205.34,204.62)(14.000,-21.620){2}{\rule{0.800pt}{0.814pt}}
\multiput(222.41,172.21)(0.509,-1.567){19}{\rule{0.123pt}{2.600pt}}
\multiput(219.34,177.60)(13.000,-33.604){2}{\rule{0.800pt}{1.300pt}}
\multiput(235.41,131.42)(0.509,-1.857){19}{\rule{0.123pt}{3.031pt}}
\multiput(232.34,137.71)(13.000,-39.709){2}{\rule{0.800pt}{1.515pt}}
\multiput(248.41,98.00)(0.509,3.239){21}{\rule{0.123pt}{5.114pt}}
\multiput(245.34,98.00)(14.000,75.385){2}{\rule{0.800pt}{2.557pt}}
\put(261,180.34){\rule{2.800pt}{0.800pt}}
\multiput(261.00,182.34)(7.188,-4.000){2}{\rule{1.400pt}{0.800pt}}
\multiput(275.41,180.00)(0.509,1.360){19}{\rule{0.123pt}{2.292pt}}
\multiput(272.34,180.00)(13.000,29.242){2}{\rule{0.800pt}{1.146pt}}
\multiput(288.41,214.00)(0.509,0.531){21}{\rule{0.123pt}{1.057pt}}
\multiput(285.34,214.00)(14.000,12.806){2}{\rule{0.800pt}{0.529pt}}
\multiput(302.41,209.01)(0.509,-3.056){19}{\rule{0.123pt}{4.815pt}}
\multiput(299.34,219.01)(13.000,-65.005){2}{\rule{0.800pt}{2.408pt}}
\multiput(315.41,144.87)(0.509,-1.293){21}{\rule{0.123pt}{2.200pt}}
\multiput(312.34,149.43)(14.000,-30.434){2}{\rule{0.800pt}{1.100pt}}
\multiput(341.00,117.06)(1.768,-0.560){3}{\rule{2.280pt}{0.135pt}}
\multiput(341.00,117.34)(8.268,-5.000){2}{\rule{1.140pt}{0.800pt}}
\multiput(355.41,114.00)(0.509,0.645){21}{\rule{0.123pt}{1.229pt}}
\multiput(352.34,114.00)(14.000,15.450){2}{\rule{0.800pt}{0.614pt}}
\multiput(369.41,132.00)(0.509,2.601){19}{\rule{0.123pt}{4.138pt}}
\multiput(366.34,132.00)(13.000,55.410){2}{\rule{0.800pt}{2.069pt}}
\multiput(381.00,197.40)(0.710,0.514){13}{\rule{1.320pt}{0.124pt}}
\multiput(381.00,194.34)(11.260,10.000){2}{\rule{0.660pt}{0.800pt}}
\multiput(396.41,200.32)(0.509,-0.740){19}{\rule{0.123pt}{1.369pt}}
\multiput(393.34,203.16)(13.000,-16.158){2}{\rule{0.800pt}{0.685pt}}
\put(408,183.84){\rule{3.132pt}{0.800pt}}
\multiput(408.00,185.34)(6.500,-3.000){2}{\rule{1.566pt}{0.800pt}}
\multiput(422.41,184.00)(0.509,1.255){21}{\rule{0.123pt}{2.143pt}}
\multiput(419.34,184.00)(14.000,29.552){2}{\rule{0.800pt}{1.071pt}}
\multiput(436.41,212.57)(0.509,-0.698){19}{\rule{0.123pt}{1.308pt}}
\multiput(433.34,215.29)(13.000,-15.286){2}{\rule{0.800pt}{0.654pt}}
\multiput(448.00,201.41)(0.533,0.509){19}{\rule{1.062pt}{0.123pt}}
\multiput(448.00,198.34)(11.797,13.000){2}{\rule{0.531pt}{0.800pt}}
\multiput(462.00,211.08)(0.847,-0.520){9}{\rule{1.500pt}{0.125pt}}
\multiput(462.00,211.34)(9.887,-8.000){2}{\rule{0.750pt}{0.800pt}}
\multiput(475.00,203.08)(0.654,-0.514){13}{\rule{1.240pt}{0.124pt}}
\multiput(475.00,203.34)(10.426,-10.000){2}{\rule{0.620pt}{0.800pt}}
\multiput(489.41,195.00)(0.509,2.361){21}{\rule{0.123pt}{3.800pt}}
\multiput(486.34,195.00)(14.000,55.113){2}{\rule{0.800pt}{1.900pt}}
\multiput(502.00,256.08)(0.847,-0.520){9}{\rule{1.500pt}{0.125pt}}
\multiput(502.00,256.34)(9.887,-8.000){2}{\rule{0.750pt}{0.800pt}}
\multiput(516.41,250.00)(0.509,0.905){19}{\rule{0.123pt}{1.615pt}}
\multiput(513.34,250.00)(13.000,19.647){2}{\rule{0.800pt}{0.808pt}}
\multiput(529.41,273.00)(0.509,0.988){21}{\rule{0.123pt}{1.743pt}}
\multiput(526.34,273.00)(14.000,23.383){2}{\rule{0.800pt}{0.871pt}}
\multiput(542.00,298.08)(0.536,-0.511){17}{\rule{1.067pt}{0.123pt}}
\multiput(542.00,298.34)(10.786,-12.000){2}{\rule{0.533pt}{0.800pt}}
\multiput(556.41,276.73)(0.509,-1.637){21}{\rule{0.123pt}{2.714pt}}
\multiput(553.34,282.37)(14.000,-38.366){2}{\rule{0.800pt}{1.357pt}}
\multiput(569.00,242.08)(0.737,-0.516){11}{\rule{1.356pt}{0.124pt}}
\multiput(569.00,242.34)(10.186,-9.000){2}{\rule{0.678pt}{0.800pt}}
\multiput(583.41,235.00)(0.509,1.608){19}{\rule{0.123pt}{2.662pt}}
\multiput(580.34,235.00)(13.000,34.476){2}{\rule{0.800pt}{1.331pt}}
\put(328.0,119.0){\rule[-0.400pt]{3.132pt}{0.800pt}}
\multiput(610.41,275.00)(0.509,2.270){19}{\rule{0.123pt}{3.646pt}}
\multiput(607.34,275.00)(13.000,48.432){2}{\rule{0.800pt}{1.823pt}}
\multiput(623.41,318.78)(0.509,-1.789){21}{\rule{0.123pt}{2.943pt}}
\multiput(620.34,324.89)(14.000,-41.892){2}{\rule{0.800pt}{1.471pt}}
\multiput(636.00,284.41)(0.536,0.511){17}{\rule{1.067pt}{0.123pt}}
\multiput(636.00,281.34)(10.786,12.000){2}{\rule{0.533pt}{0.800pt}}
\multiput(650.41,295.00)(0.509,1.939){19}{\rule{0.123pt}{3.154pt}}
\multiput(647.34,295.00)(13.000,41.454){2}{\rule{0.800pt}{1.577pt}}
\multiput(663.41,343.00)(0.509,1.179){21}{\rule{0.123pt}{2.029pt}}
\multiput(660.34,343.00)(14.000,27.790){2}{\rule{0.800pt}{1.014pt}}
\multiput(677.41,365.23)(0.509,-1.402){19}{\rule{0.123pt}{2.354pt}}
\multiput(674.34,370.11)(13.000,-30.114){2}{\rule{0.800pt}{1.177pt}}
\multiput(690.41,340.00)(0.509,1.293){21}{\rule{0.123pt}{2.200pt}}
\multiput(687.34,340.00)(14.000,30.434){2}{\rule{0.800pt}{1.100pt}}
\multiput(703.00,376.41)(0.492,0.509){19}{\rule{1.000pt}{0.123pt}}
\multiput(703.00,373.34)(10.924,13.000){2}{\rule{0.500pt}{0.800pt}}
\multiput(717.41,378.74)(0.509,-1.319){19}{\rule{0.123pt}{2.231pt}}
\multiput(714.34,383.37)(13.000,-28.370){2}{\rule{0.800pt}{1.115pt}}
\multiput(730.41,355.00)(0.509,1.827){21}{\rule{0.123pt}{3.000pt}}
\multiput(727.34,355.00)(14.000,42.773){2}{\rule{0.800pt}{1.500pt}}
\multiput(743.00,402.08)(0.847,-0.520){9}{\rule{1.500pt}{0.125pt}}
\multiput(743.00,402.34)(9.887,-8.000){2}{\rule{0.750pt}{0.800pt}}
\multiput(757.41,396.00)(0.509,2.895){21}{\rule{0.123pt}{4.600pt}}
\multiput(754.34,396.00)(14.000,67.452){2}{\rule{0.800pt}{2.300pt}}
\multiput(770.00,474.40)(0.737,0.516){11}{\rule{1.356pt}{0.124pt}}
\multiput(770.00,471.34)(10.186,9.000){2}{\rule{0.678pt}{0.800pt}}
\multiput(784.41,482.00)(0.509,2.105){19}{\rule{0.123pt}{3.400pt}}
\multiput(781.34,482.00)(13.000,44.943){2}{\rule{0.800pt}{1.700pt}}
\multiput(797.41,534.00)(0.509,0.645){21}{\rule{0.123pt}{1.229pt}}
\multiput(794.34,534.00)(14.000,15.450){2}{\rule{0.800pt}{0.614pt}}
\multiput(811.41,552.00)(0.509,1.608){19}{\rule{0.123pt}{2.662pt}}
\multiput(808.34,552.00)(13.000,34.476){2}{\rule{0.800pt}{1.331pt}}
\multiput(823.00,593.41)(0.492,0.509){19}{\rule{1.000pt}{0.123pt}}
\multiput(823.00,590.34)(10.924,13.000){2}{\rule{0.500pt}{0.800pt}}
\multiput(836.00,603.07)(1.355,-0.536){5}{\rule{2.067pt}{0.129pt}}
\multiput(836.00,603.34)(9.711,-6.000){2}{\rule{1.033pt}{0.800pt}}
\multiput(851.41,599.00)(0.509,1.112){19}{\rule{0.123pt}{1.923pt}}
\multiput(848.34,599.00)(13.000,24.009){2}{\rule{0.800pt}{0.962pt}}
\multiput(864.41,622.37)(0.509,-0.569){21}{\rule{0.123pt}{1.114pt}}
\multiput(861.34,624.69)(14.000,-13.687){2}{\rule{0.800pt}{0.557pt}}
\multiput(878.41,611.00)(0.509,1.443){19}{\rule{0.123pt}{2.415pt}}
\multiput(875.34,611.00)(13.000,30.987){2}{\rule{0.800pt}{1.208pt}}
\multiput(891.41,638.00)(0.509,-1.278){19}{\rule{0.123pt}{2.169pt}}
\multiput(888.34,642.50)(13.000,-27.498){2}{\rule{0.800pt}{1.085pt}}
\multiput(904.41,615.00)(0.509,0.950){21}{\rule{0.123pt}{1.686pt}}
\multiput(901.34,615.00)(14.000,22.501){2}{\rule{0.800pt}{0.843pt}}
\multiput(918.41,627.91)(0.509,-1.939){19}{\rule{0.123pt}{3.154pt}}
\multiput(915.34,634.45)(13.000,-41.454){2}{\rule{0.800pt}{1.577pt}}
\multiput(931.41,583.87)(0.509,-1.293){21}{\rule{0.123pt}{2.200pt}}
\multiput(928.34,588.43)(14.000,-30.434){2}{\rule{0.800pt}{1.100pt}}
\multiput(945.41,558.00)(0.509,2.725){19}{\rule{0.123pt}{4.323pt}}
\multiput(942.34,558.00)(13.000,58.027){2}{\rule{0.800pt}{2.162pt}}
\multiput(958.41,614.46)(0.509,-1.526){19}{\rule{0.123pt}{2.538pt}}
\multiput(955.34,619.73)(13.000,-32.731){2}{\rule{0.800pt}{1.269pt}}
\multiput(970.00,585.08)(1.088,-0.526){7}{\rule{1.800pt}{0.127pt}}
\multiput(970.00,585.34)(10.264,-7.000){2}{\rule{0.900pt}{0.800pt}}
\multiput(985.41,580.00)(0.509,2.518){19}{\rule{0.123pt}{4.015pt}}
\multiput(982.34,580.00)(13.000,53.666){2}{\rule{0.800pt}{2.008pt}}
\multiput(998.41,627.41)(0.509,-2.171){21}{\rule{0.123pt}{3.514pt}}
\multiput(995.34,634.71)(14.000,-50.706){2}{\rule{0.800pt}{1.757pt}}
\multiput(1012.41,567.59)(0.509,-2.477){19}{\rule{0.123pt}{3.954pt}}
\multiput(1009.34,575.79)(13.000,-52.794){2}{\rule{0.800pt}{1.977pt}}
\multiput(1025.41,516.29)(0.509,-0.905){19}{\rule{0.123pt}{1.615pt}}
\multiput(1022.34,519.65)(13.000,-19.647){2}{\rule{0.800pt}{0.808pt}}
\multiput(1038.41,492.29)(0.509,-1.065){21}{\rule{0.123pt}{1.857pt}}
\multiput(1035.34,496.15)(14.000,-25.145){2}{\rule{0.800pt}{0.929pt}}
\multiput(1052.41,463.53)(0.509,-1.029){19}{\rule{0.123pt}{1.800pt}}
\multiput(1049.34,467.26)(13.000,-22.264){2}{\rule{0.800pt}{0.900pt}}
\multiput(1065.41,440.59)(0.509,-0.533){19}{\rule{0.123pt}{1.062pt}}
\multiput(1062.34,442.80)(13.000,-11.797){2}{\rule{0.800pt}{0.531pt}}
\multiput(1077.00,432.40)(0.639,0.512){15}{\rule{1.218pt}{0.123pt}}
\multiput(1077.00,429.34)(11.472,11.000){2}{\rule{0.609pt}{0.800pt}}
\multiput(1092.41,442.00)(0.509,1.815){19}{\rule{0.123pt}{2.969pt}}
\multiput(1089.34,442.00)(13.000,38.837){2}{\rule{0.800pt}{1.485pt}}
\multiput(1105.41,487.00)(0.509,1.904){21}{\rule{0.123pt}{3.114pt}}
\multiput(1102.34,487.00)(14.000,44.536){2}{\rule{0.800pt}{1.557pt}}
\multiput(1119.41,538.00)(0.509,1.360){19}{\rule{0.123pt}{2.292pt}}
\multiput(1116.34,538.00)(13.000,29.242){2}{\rule{0.800pt}{1.146pt}}
\put(1131,569.84){\rule{3.132pt}{0.800pt}}
\multiput(1131.00,570.34)(6.500,-1.000){2}{\rule{1.566pt}{0.800pt}}
\multiput(1145.41,571.00)(0.509,1.865){21}{\rule{0.123pt}{3.057pt}}
\multiput(1142.34,571.00)(14.000,43.655){2}{\rule{0.800pt}{1.529pt}}
\multiput(1159.41,621.00)(0.509,1.857){19}{\rule{0.123pt}{3.031pt}}
\multiput(1156.34,621.00)(13.000,39.709){2}{\rule{0.800pt}{1.515pt}}
\multiput(1172.41,656.44)(0.509,-1.522){21}{\rule{0.123pt}{2.543pt}}
\multiput(1169.34,661.72)(14.000,-35.722){2}{\rule{0.800pt}{1.271pt}}
\multiput(1186.41,626.00)(0.509,1.278){19}{\rule{0.123pt}{2.169pt}}
\multiput(1183.34,626.00)(13.000,27.498){2}{\rule{0.800pt}{1.085pt}}
\multiput(1199.41,647.72)(0.509,-1.484){19}{\rule{0.123pt}{2.477pt}}
\multiput(1196.34,652.86)(13.000,-31.859){2}{\rule{0.800pt}{1.238pt}}
\multiput(1212.41,613.05)(0.509,-1.103){21}{\rule{0.123pt}{1.914pt}}
\multiput(1209.34,617.03)(14.000,-26.027){2}{\rule{0.800pt}{0.957pt}}
\multiput(1226.41,579.19)(0.509,-1.733){19}{\rule{0.123pt}{2.846pt}}
\multiput(1223.34,585.09)(13.000,-37.093){2}{\rule{0.800pt}{1.423pt}}
\multiput(1239.41,548.00)(0.509,1.179){21}{\rule{0.123pt}{2.029pt}}
\multiput(1236.34,548.00)(14.000,27.790){2}{\rule{0.800pt}{1.014pt}}
\multiput(1253.41,580.00)(0.509,1.236){19}{\rule{0.123pt}{2.108pt}}
\multiput(1250.34,580.00)(13.000,26.625){2}{\rule{0.800pt}{1.054pt}}
\multiput(1265.00,609.08)(1.000,-0.526){7}{\rule{1.686pt}{0.127pt}}
\multiput(1265.00,609.34)(9.501,-7.000){2}{\rule{0.843pt}{0.800pt}}
\multiput(1279.41,604.00)(0.509,1.904){21}{\rule{0.123pt}{3.114pt}}
\multiput(1276.34,604.00)(14.000,44.536){2}{\rule{0.800pt}{1.557pt}}
\multiput(1293.41,655.00)(0.509,1.236){19}{\rule{0.123pt}{2.108pt}}
\multiput(1290.34,655.00)(13.000,26.625){2}{\rule{0.800pt}{1.054pt}}
\put(1305,685.34){\rule{3.132pt}{0.800pt}}
\multiput(1305.00,684.34)(6.500,2.000){2}{\rule{1.566pt}{0.800pt}}
\multiput(1319.41,680.77)(0.509,-0.988){21}{\rule{0.123pt}{1.743pt}}
\multiput(1316.34,684.38)(14.000,-23.383){2}{\rule{0.800pt}{0.871pt}}
\multiput(1333.41,651.74)(0.509,-1.319){19}{\rule{0.123pt}{2.231pt}}
\multiput(1330.34,656.37)(13.000,-28.370){2}{\rule{0.800pt}{1.115pt}}
\multiput(1346.41,615.55)(0.509,-1.827){21}{\rule{0.123pt}{3.000pt}}
\multiput(1343.34,621.77)(14.000,-42.773){2}{\rule{0.800pt}{1.500pt}}
\multiput(1359.00,580.40)(0.654,0.514){13}{\rule{1.240pt}{0.124pt}}
\multiput(1359.00,577.34)(10.426,10.000){2}{\rule{0.620pt}{0.800pt}}
\multiput(1373.41,589.00)(0.509,0.533){19}{\rule{0.123pt}{1.062pt}}
\multiput(1370.34,589.00)(13.000,11.797){2}{\rule{0.800pt}{0.531pt}}
\multiput(1385.00,604.41)(0.581,0.511){17}{\rule{1.133pt}{0.123pt}}
\multiput(1385.00,601.34)(11.648,12.000){2}{\rule{0.567pt}{0.800pt}}
\put(595.0,275.0){\rule[-0.400pt]{3.373pt}{0.800pt}}
\multiput(1413.41,615.00)(0.509,0.836){21}{\rule{0.123pt}{1.514pt}}
\multiput(1410.34,615.00)(14.000,19.857){2}{\rule{0.800pt}{0.757pt}}
\multiput(1426.00,639.40)(0.589,0.512){15}{\rule{1.145pt}{0.123pt}}
\multiput(1426.00,636.34)(10.623,11.000){2}{\rule{0.573pt}{0.800pt}}
\put(1399.0,615.0){\rule[-0.400pt]{3.132pt}{0.800pt}}
\end{picture}

CoolingCurve.tex
% GNUPLOT: LaTeX picture
\setlength{\unitlength}{0.240900pt}
\ifx\plotpoint\undefined\newsavebox{\plotpoint}\fi
\sbox{\plotpoint}{\rule[-0.200pt]{0.400pt}{0.400pt}}%
\begin{picture}(1500,900)(0,0)
\font\gnuplot=cmr10 at 10pt
\gnuplot
\sbox{\plotpoint}{\rule[-0.200pt]{0.400pt}{0.400pt}}%
\put(100.0,82.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,82){\makebox(0,0)[r]{0}}
\put(1419.0,82.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,179.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,179){\makebox(0,0)[r]{10}}
\put(1419.0,179.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,277.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,277){\makebox(0,0)[r]{20}}
\put(1419.0,277.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,374.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,374){\makebox(0,0)[r]{30}}
\put(1419.0,374.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,471.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,471){\makebox(0,0)[r]{40}}
\put(1419.0,471.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,568.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,568){\makebox(0,0)[r]{50}}
\put(1419.0,568.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,666.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,666){\makebox(0,0)[r]{60}}
\put(1419.0,666.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,763.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,763){\makebox(0,0)[r]{70}}
\put(1419.0,763.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,860.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(80,860){\makebox(0,0)[r]{80}}
\put(1419.0,860.0){\rule[-0.200pt]{4.818pt}{0.400pt}}
\put(100.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(100,41){\makebox(0,0){0}}
\put(100.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(234.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(234,41){\makebox(0,0){5}}
\put(234.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(368.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(368,41){\makebox(0,0){10}}
\put(368.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(502.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(502,41){\makebox(0,0){15}}
\put(502.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(636.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(636,41){\makebox(0,0){20}}
\put(636.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(770.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(770,41){\makebox(0,0){25}}
\put(770.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(903.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(903,41){\makebox(0,0){30}}
\put(903.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1037.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1037,41){\makebox(0,0){35}}
\put(1037.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1171.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1171,41){\makebox(0,0){40}}
\put(1171.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1305.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1305,41){\makebox(0,0){45}}
\put(1305.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1439.0,82.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(1439,41){\makebox(0,0){50}}
\put(1439.0,840.0){\rule[-0.200pt]{0.400pt}{4.818pt}}
\put(100.0,82.0){\rule[-0.200pt]{322.565pt}{0.400pt}}
\put(1439.0,82.0){\rule[-0.200pt]{0.400pt}{187.420pt}}
\put(100.0,860.0){\rule[-0.200pt]{322.565pt}{0.400pt}}
\put(100.0,82.0){\rule[-0.200pt]{0.400pt}{187.420pt}}
\put(1279,820){\makebox(0,0)[r]{Cooling Curve}}
\put(1299.0,820.0){\rule[-0.200pt]{24.090pt}{0.400pt}}
\put(100,692){\usebox{\plotpoint}}
\multiput(100.00,692.59)(2.389,0.482){9}{\rule{1.900pt}{0.116pt}}
\multiput(100.00,691.17)(23.056,6.000){2}{\rule{0.950pt}{0.400pt}}
\multiput(127.58,692.54)(0.497,-1.531){51}{\rule{0.120pt}{1.315pt}}
\multiput(126.17,695.27)(27.000,-79.271){2}{\rule{0.400pt}{0.657pt}}
\multiput(154.58,616.00)(0.497,2.410){49}{\rule{0.120pt}{2.008pt}}
\multiput(153.17,616.00)(26.000,119.833){2}{\rule{0.400pt}{1.004pt}}
\multiput(180.58,740.00)(0.497,0.949){51}{\rule{0.120pt}{0.856pt}}
\multiput(179.17,740.00)(27.000,49.224){2}{\rule{0.400pt}{0.428pt}}
\multiput(207.58,776.81)(0.497,-4.196){51}{\rule{0.120pt}{3.419pt}}
\multiput(206.17,783.90)(27.000,-216.905){2}{\rule{0.400pt}{1.709pt}}
\multiput(234.00,565.92)(0.753,-0.495){33}{\rule{0.700pt}{0.119pt}}
\multiput(234.00,566.17)(25.547,-18.000){2}{\rule{0.350pt}{0.400pt}}
\multiput(261.00,549.58)(0.686,0.495){35}{\rule{0.647pt}{0.119pt}}
\multiput(261.00,548.17)(24.656,19.000){2}{\rule{0.324pt}{0.400pt}}
\multiput(287.58,562.48)(0.497,-1.550){51}{\rule{0.120pt}{1.330pt}}
\multiput(286.17,565.24)(27.000,-80.240){2}{\rule{0.400pt}{0.665pt}}
\multiput(314.00,485.61)(5.820,0.447){3}{\rule{3.700pt}{0.108pt}}
\multiput(314.00,484.17)(19.320,3.000){2}{\rule{1.850pt}{0.400pt}}
\multiput(341.00,486.94)(3.844,-0.468){5}{\rule{2.800pt}{0.113pt}}
\multiput(341.00,487.17)(21.188,-4.000){2}{\rule{1.400pt}{0.400pt}}
\multiput(368.58,484.00)(0.497,0.630){51}{\rule{0.120pt}{0.604pt}}
\multiput(367.17,484.00)(27.000,32.747){2}{\rule{0.400pt}{0.302pt}}
\multiput(395.58,510.11)(0.497,-2.273){49}{\rule{0.120pt}{1.900pt}}
\multiput(394.17,514.06)(26.000,-113.056){2}{\rule{0.400pt}{0.950pt}}
\multiput(421.58,401.00)(0.497,0.611){51}{\rule{0.120pt}{0.589pt}}
\multiput(420.17,401.00)(27.000,31.778){2}{\rule{0.400pt}{0.294pt}}
\multiput(448.58,429.71)(0.497,-1.174){51}{\rule{0.120pt}{1.033pt}}
\multiput(447.17,431.86)(27.000,-60.855){2}{\rule{0.400pt}{0.517pt}}
\multiput(475.58,365.67)(0.497,-1.493){51}{\rule{0.120pt}{1.285pt}}
\multiput(474.17,368.33)(27.000,-77.333){2}{\rule{0.400pt}{0.643pt}}
\multiput(502.00,289.92)(0.541,-0.496){45}{\rule{0.533pt}{0.120pt}}
\multiput(502.00,290.17)(24.893,-24.000){2}{\rule{0.267pt}{0.400pt}}
\multiput(528.00,267.58)(0.498,0.497){51}{\rule{0.500pt}{0.120pt}}
\multiput(528.00,266.17)(25.962,27.000){2}{\rule{0.250pt}{0.400pt}}
\multiput(555.58,291.19)(0.497,-0.724){51}{\rule{0.120pt}{0.678pt}}
\multiput(554.17,292.59)(27.000,-37.593){2}{\rule{0.400pt}{0.339pt}}
\put(582,253.17){\rule{5.500pt}{0.400pt}}
\multiput(582.00,254.17)(15.584,-2.000){2}{\rule{2.750pt}{0.400pt}}
\multiput(609.00,251.92)(0.798,-0.495){31}{\rule{0.735pt}{0.119pt}}
\multiput(609.00,252.17)(25.474,-17.000){2}{\rule{0.368pt}{0.400pt}}
\multiput(636.00,234.92)(0.619,-0.496){39}{\rule{0.595pt}{0.119pt}}
\multiput(636.00,235.17)(24.765,-21.000){2}{\rule{0.298pt}{0.400pt}}
\multiput(662.00,213.92)(0.753,-0.495){33}{\rule{0.700pt}{0.119pt}}
\multiput(662.00,214.17)(25.547,-18.000){2}{\rule{0.350pt}{0.400pt}}
\multiput(689.00,195.93)(1.748,-0.488){13}{\rule{1.450pt}{0.117pt}}
\multiput(689.00,196.17)(23.990,-8.000){2}{\rule{0.725pt}{0.400pt}}
\multiput(716.00,187.93)(1.543,-0.489){15}{\rule{1.300pt}{0.118pt}}
\multiput(716.00,188.17)(24.302,-9.000){2}{\rule{0.650pt}{0.400pt}}
\put(743,180.17){\rule{5.500pt}{0.400pt}}
\multiput(743.00,179.17)(15.584,2.000){2}{\rule{2.750pt}{0.400pt}}
\multiput(770.00,180.93)(1.682,-0.488){13}{\rule{1.400pt}{0.117pt}}
\multiput(770.00,181.17)(23.094,-8.000){2}{\rule{0.700pt}{0.400pt}}
\multiput(796.00,172.93)(2.936,-0.477){7}{\rule{2.260pt}{0.115pt}}
\multiput(796.00,173.17)(22.309,-5.000){2}{\rule{1.130pt}{0.400pt}}
\put(823,167.67){\rule{6.504pt}{0.400pt}}
\multiput(823.00,168.17)(13.500,-1.000){2}{\rule{3.252pt}{0.400pt}}
\put(850,166.17){\rule{5.500pt}{0.400pt}}
\multiput(850.00,167.17)(15.584,-2.000){2}{\rule{2.750pt}{0.400pt}}
\multiput(877.00,164.95)(5.597,-0.447){3}{\rule{3.567pt}{0.108pt}}
\multiput(877.00,165.17)(18.597,-3.000){2}{\rule{1.783pt}{0.400pt}}
\put(903,161.67){\rule{6.504pt}{0.400pt}}
\multiput(903.00,162.17)(13.500,-1.000){2}{\rule{3.252pt}{0.400pt}}
\put(957,160.67){\rule{6.504pt}{0.400pt}}
\multiput(957.00,161.17)(13.500,-1.000){2}{\rule{3.252pt}{0.400pt}}
\put(930.0,162.0){\rule[-0.200pt]{6.504pt}{0.400pt}}
\put(1011,159.67){\rule{6.263pt}{0.400pt}}
\multiput(1011.00,160.17)(13.000,-1.000){2}{\rule{3.132pt}{0.400pt}}
\put(984.0,161.0){\rule[-0.200pt]{6.504pt}{0.400pt}}
\put(1037.0,160.0){\rule[-0.200pt]{96.842pt}{0.400pt}}
\put(1279,779){\makebox(0,0)[r]{Minimal Energy is 8}}
\put(1000,0){\makebox(0,0)[r]{Temperature Decrease Step}}
\put(100,900){\makebox(0,0)[r]{Energy}}
\multiput(1299,779)(20.756,0.000){5}{\usebox{\plotpoint}}
\put(1399,779){\usebox{\plotpoint}}
\put(100,160){\usebox{\plotpoint}}
\put(100.00,160.00){\usebox{\plotpoint}}
\put(120.76,160.00){\usebox{\plotpoint}}
\put(141.51,160.00){\usebox{\plotpoint}}
\put(162.27,160.00){\usebox{\plotpoint}}
\put(183.02,160.00){\usebox{\plotpoint}}
\put(203.78,160.00){\usebox{\plotpoint}}
\put(224.53,160.00){\usebox{\plotpoint}}
\put(245.29,160.00){\usebox{\plotpoint}}
\put(266.04,160.00){\usebox{\plotpoint}}
\put(286.80,160.00){\usebox{\plotpoint}}
\put(307.55,160.00){\usebox{\plotpoint}}
\put(328.31,160.00){\usebox{\plotpoint}}
\put(349.07,160.00){\usebox{\plotpoint}}
\put(369.82,160.00){\usebox{\plotpoint}}
\put(390.58,160.00){\usebox{\plotpoint}}
\put(411.33,160.00){\usebox{\plotpoint}}
\put(432.09,160.00){\usebox{\plotpoint}}
\put(452.84,160.00){\usebox{\plotpoint}}
\put(473.60,160.00){\usebox{\plotpoint}}
\put(494.35,160.00){\usebox{\plotpoint}}
\put(515.11,160.00){\usebox{\plotpoint}}
\put(535.87,160.00){\usebox{\plotpoint}}
\put(556.62,160.00){\usebox{\plotpoint}}
\put(577.38,160.00){\usebox{\plotpoint}}
\put(598.13,160.00){\usebox{\plotpoint}}
\put(618.89,160.00){\usebox{\plotpoint}}
\put(639.64,160.00){\usebox{\plotpoint}}
\put(660.40,160.00){\usebox{\plotpoint}}
\put(681.15,160.00){\usebox{\plotpoint}}
\put(701.91,160.00){\usebox{\plotpoint}}
\put(722.66,160.00){\usebox{\plotpoint}}
\put(743.42,160.00){\usebox{\plotpoint}}
\put(764.18,160.00){\usebox{\plotpoint}}
\put(784.93,160.00){\usebox{\plotpoint}}
\put(805.69,160.00){\usebox{\plotpoint}}
\put(826.44,160.00){\usebox{\plotpoint}}
\put(847.20,160.00){\usebox{\plotpoint}}
\put(867.95,160.00){\usebox{\plotpoint}}
\put(888.71,160.00){\usebox{\plotpoint}}
\put(909.46,160.00){\usebox{\plotpoint}}
\put(930.22,160.00){\usebox{\plotpoint}}
\put(950.98,160.00){\usebox{\plotpoint}}
\put(971.73,160.00){\usebox{\plotpoint}}
\put(992.49,160.00){\usebox{\plotpoint}}
\put(1013.24,160.00){\usebox{\plotpoint}}
\put(1034.00,160.00){\usebox{\plotpoint}}
\put(1054.75,160.00){\usebox{\plotpoint}}
\put(1075.51,160.00){\usebox{\plotpoint}}
\put(1096.26,160.00){\usebox{\plotpoint}}
\put(1117.02,160.00){\usebox{\plotpoint}}
\put(1137.77,160.00){\usebox{\plotpoint}}
\put(1158.53,160.00){\usebox{\plotpoint}}
\put(1179.29,160.00){\usebox{\plotpoint}}
\put(1200.04,160.00){\usebox{\plotpoint}}
\put(1220.80,160.00){\usebox{\plotpoint}}
\put(1241.55,160.00){\usebox{\plotpoint}}
\put(1262.31,160.00){\usebox{\plotpoint}}
\put(1283.06,160.00){\usebox{\plotpoint}}
\put(1303.82,160.00){\usebox{\plotpoint}}
\put(1324.57,160.00){\usebox{\plotpoint}}
\put(1345.33,160.00){\usebox{\plotpoint}}
\put(1366.09,160.00){\usebox{\plotpoint}}
\put(1386.84,160.00){\usebox{\plotpoint}}
\put(1407.60,160.00){\usebox{\plotpoint}}
\put(1428.35,160.00){\usebox{\plotpoint}}
\put(1439,160){\usebox{\plotpoint}}
\end{picture}

sa.tex
% This LaTeX file was written by tom.weidig@physics.org (UMIST, UK)
% (in collaboration with Mark Hale (Durham, UK) and Oliver Schwindt (UMIST, UK))
% this is the revised version   February 1999
% needs 9 postscript(eps) files
% needs bibtex file called paper_ref.bib
% should be 31 pages long (single-spaced)
% 5/7/00 - Updated with corrections for publication. Mark Hale

\documentclass[11pt,a4paper]{article}
\usepackage{epsfig} %include package for eps-file
\usepackage{graphicx}
\usepackage{graphpap,geometry}

\setlength{\parindent}{0mm}
\setlength{\hoffset}{0cm}
\setlength{\textheight}{230mm}
\setlength{\textwidth}{170mm}
\renewcommand{\baselinestretch}{1.4} % Line Spacing
\newcommand*{\e}{\mathop{\mathrm{e}}\nolimits}
\newcommand*{\ph}{\vec{\phi}} \newcommand*{\pha}{\phi_{a}}
\renewcommand*{\d}{\mathrm{d}}
\newcommand*{\dmu}{\partial_{\mu}} \newcommand*{\umu}{\partial^{\mu}}
\newcommand*{\dnu}{\partial_{\nu}} \newcommand*{\unu}{\partial^{\nu}}
\newcommand*{\di}{\partial_{i}} \newcommand*{\ui}{\partial^{i}}
\renewcommand*{\dj}{\partial_{j}} \newcommand*{\uj}{\partial^{j}}

\title{Simulated Annealing for Topological Solitons}

\author{Mark Hale \\
{\it Centre for Particle Theory, University of Durham, Durham DH1 3LE, UK} \\
and \\
Oliver Schwindt and Tom Weidig \\
{\it Dept.~of Physics, UMIST, Manchester M60 1QD, UK}}

\begin{document}
\maketitle

\begin{abstract}
The search for solutions of field theories allowing for
topological solitons requires that we find the field configuration
with the lowest energy in a given sector of topological charge. The
standard approach is based on the numerical solution of the static
Euler-Lagrange differential equation following from the field
energy. As an alternative, we propose to use a simulated annealing
algorithm to minimize the energy functional directly. We have applied
simulated annealing to several nonlinear classical field theories: the
sine-Gordon model in one dimension, the baby Skyrme model in two
dimensions and the nuclear Skyrme model in three dimensions. We
describe in detail the implementation of the simulated annealing
algorithm, present our results and get independent confirmation of the
studies which have used standard minimization techniques.
\end{abstract}


\section{Introduction}

Solitons have played an ever increasing role in the description of
physical phenomena since their discovery by Russell \cite{russell} in
1834.  Generally speaking, a soliton is a stable localized solution of
a nonlinear partial differential equation which propagates at a
constant speed and stays localized even after an interaction with
another soliton, see Ref.~\cite{soliton:intro} for an introduction to
solitons. They are particle-like extended objects.  Well-known soliton
models in 1D are the Korteweg-de Vries (KdV) and the sine-Gordon
model.  The stability of the KdV soliton is due to the dynamical
balance between the nonlinear and the dispersive term in the KdV
equation. This differential equation belongs to the class of
integrable models which can be solved exactly. The sine-Gordon model
is also integrable and the stability of its soliton is also based on
the balance between nonlinearity and dispersion.  However, it also
belongs to the wider class of models whose solitons are stable by
conservation of a topological charge or winding number, as discussed
see in Sec.~\ref{sine}. In topology, the field is interpreted as a
mapping from physical space to field space and a field configuration
with given topological charge cannot dynamically change into a field
configuration with a different charge. In this paper, when we discuss
topological solitons, we shall consider only the field configuration
with lowest energy in a non-zero topological charge sector.

Topological solitons arise in many areas of physics: field theory
(e.g.~vortices, monopoles and instantons; as discussed in
Ref.~\cite{ryder}), condensed matter (e.g.~baby skyrmions \cite{qhe}),
nuclear physics (skyrmions, see Ref.~\cite{holzwarth:review}),
cosmology (e.g.~cosmic strings \cite{cosmic:string}) and string
theory/M-theory (e.g.~Olive-Montonen duality \cite{duff}).  They
possess many interesting properties.  The conservation of topological
charge can be used to model particle conservation and annihilation:
the number of solitons is conserved and two solitons with opposite
charge can annihilate. Since the stability of solitons is assured by
topology, there exits considerable freedom in constructing appropriate
lagrangians for physical systems.  The common constraint of Lorentz
invariance is easily imposed by using covariant terms. However, the
use of topological models has a major disadvantage. The models are
generally not integrable and few analytical techniques are available
(only Ans\"atze, topological bounds, etc.). It is therefore crucial to
study the models using reliable and efficient numerical methods.

We are looking for the (static) lowest energy field configuration in a
given topological sector. Thus, we need to minimize the energy
functional $E$ of the field theory, the integral of an energy density
$\mathcal{E}$ over a manifold $M$:
\begin{equation}
 E[C]=\int_{M}\!\d^nx~\mathcal{E}\left(x,f(x),f'(x)\right)
\label{min}
\end{equation}
where $f(x)$ is a field configuration $C$. The topology typically
imposes  some boundary conditions $f(\partial M)$. We are searching
for the function $f_{\rm{min}}(x)$ which gives the lowest value for
$E$. There are two possible  approaches: to solve the  Euler-Lagrange
equation of the functional $E$  with respect to the function $f(x)$ or
to minimize $E$ through some other means.

Hitherto only the first approach, via the Euler-Lagrange equation, has
been used with topological systems.  We shall review the standard
numerical techniques which apply shooting or relaxation methods and
discuss their reliability and ease of use.  In this paper, we show how
to minimize the energy functional directly by using the simulated
annealing (SA) algorithm, as proposed in Refs.~\cite{sc:second,tom}.
SA is based on the fact that a solid which is slowly cooled down,
assuring thermal equilibrium at each temperature, reaches its ground
state. The SA algorithm describes the cooling process and a Metropolis
subalgorithm brings a system into thermal equilibrium.  SA has been
applied to minimization problems in such diverse areas as
combinatorial optimization (such as the traveling salesman problem),
circuit design, finance, physics and military warfare: see
Ref.~\cite{ingber} and Ref.~\cite[Sec.~10.9]{num}. We give a detailed
introduction to SA and describe our implementations. We find the
topological solitons of the sine-Gordon model in 1D, the baby Skyrme
model in 2D and the nuclear Skyrme model in 3D and compare our results
to those obtained using standard minimization techniques.


\section{Minimization via Euler-Lagrange Equation}

The standard procedure uses the Euler-Lagrange equation resulting from
the variation of the functional $E$ with respect to the function
$f(x)$ to find the minimal energy solution $f_{\rm{min}}(x)$.  In the
1D case, for example, the problem is a two point boundary value
problem satisfying the differential equation
\begin{equation}
 \frac{\d}{\d x}\left(\frac{\d\mathcal{E}}{\d f'}\right)
 -\frac{\d\mathcal{E}}{\d f}=0.
\label{DE}
\end{equation}
It is a second order ODE, and a PDE in two or more dimensions, which
is equivalent to a set of first order ODEs. Let $f(a)$ and $f(b)$
represent the boundary conditions a field configuration has to satisfy
over the interval $[a,b]$. There are two standard approaches: the
shooting and the relaxation methods (as discussed in
Ref.~\cite[Chap.~17]{num}).

\subsection{The Shooting Method}

The shooting method is usually based on integration from one boundary
to the other.  The value of the function at the point $x=a$ is taken
to be $f(a)$ and an initial guess $\alpha$ for its derivative is
made. A numerical integration, for example with a Runge-Kutta method,
up to the other boundary point $x=b$ then gives an estimate
$f_{\alpha}(b)$ for $f$ at $b$.  This value is compared to the known
boundary value $f(b)$ and $\alpha$ is adjusted to match
$f_{\alpha}(b)$ closer to $f(b)$. This procedure is repeated until the
desired accuracy is achieved. The shooting method is unrivaled in
speed and accuracy, but only applicable in one dimension.

\subsection{Relaxation Methods}

Gauss-Seidel over-relaxation (SOR) is commonly used to solve the
boundary problem directly.  A time-dependent differential equation, a
diffusion equation, is constructed out of the 1D ODE (\ref{DE}),
\begin{equation}
 \frac{\d f(x,t)}{\d t}=\omega\, \d x^2 \left[ \frac{\d}{\d x}\left(\frac{\d\mathcal{E}} {\d f'}\right)-\frac{\d\mathcal{E}}{\d f} \right].
\label{omegaDE}
\end{equation}
If the system reaches equilibrium, i.e.~$\frac{\d f}{\d t}=0$, this
configuration is a solution to (\ref{DE}). One starts out with a
configuration satisfying the boundary conditions.  The coefficient
$\omega$ of the leading term, which has the form $\frac{\d^2f}{\d x^2}$,
is dimensionless and determines the speed of convergence.  The choice
of integration method is not very sensitive, we can use Euler
integration, the Runge-Kutta method or the Crank-Nicholson method.
The standard SOR uses the Euler method with updated information from
already computed field values at lattice points and ensures better
convergence.

In many cases one also studies the time-evolution of the
models. Computer codes developed for this purpose can be adapted to
find minimal energy solutions by adding a damping term to the equation
of motion.  For example in \cite{baby_pots}, two baby skyrmions are
put in an attractive configuration and form an oscillating bound
state. As the system has been made dissipative, the energy of the
system decreases with time until a minimum is reached. Thus, finding a
minimal energy solution is translated into a damped time-evolution.
The same effect can often be reached by working with a finite box and
absorbing any outward propagating radiation on the boundaries.

Relaxation techniques are well documented and applicable in any
dimension, see Ref.~\cite{num}.  The SOR method has theoretically the
best rate of convergence but might be less than optimal since the best
choice of $\omega$ can rarely be determined for a nonlinear system and
must be made by trial and error. Using damping in a time-evolution
problem is convenient, but one first has to set up the time-evolution
code. It is impossible to estimate the error on an integration step
and one needs to monitor conserved quantities. This is especially
important if, as is often the case, the field has to satisfy a
constraint. Furthermore, if the initial configuration is far from the
global minimum, we might end up in a local minimum.  Moreover, the
derivation of the corresponding Euler-Lagrange equation becomes
increasingly difficult when higher order terms are added to the
lagrangian or when complicated constraints on the field space are
present.

We have come to the conclusion that the weaker points of iterative
minimization techniques via the Euler-Lagrange equation are:
\begin{itemize}
\item Uncertainty about the global nature of the minimum obtained.
\item Lack of direct control over the integration errors (important
for constrained fields).
\item Tedious derivation for complicated lagrangians.
\end{itemize}


\section{Minimization via Simulated Annealing}

Minimizing the energy functional directly is a more straightforward
approach than solving the equations of motion and we propose to use
the flexible and easy-to-implement simulated annealing technique.

\subsection{Metropolis Principle}

In 1953 Metropolis et.~al.~\cite{me:equation} proposed an algorithm,
now called the {\it Metropolis} or {\it $M(RT)^{2}$} algorithm, that
can be used to bring a statistical system into thermal
equilibrium. The {\it $M(RT)^{2}$} is most commonly used to evaluate
thermal averages $\langle\mathcal{F}\rangle$ of a quantity $\mathcal{F}(C)$,
\begin{equation}
 \langle\mathcal{F}\rangle=\frac{\int\mathcal{F}(C)P(C)\,\d C}{\int P(C)\,\d C}.
\label{average}
\end{equation}
Here $P(C)$ is a probability distribution for configurations $C$. It
must satisfy $P(C)\!\ge\!0$ and $\int P(C)\,\d C\!<\!\infty$ in order to be
normalisable. For example, $\langle\mathcal{F}\rangle$ can stand for the thermal
average of our energy functional $E[C]$ in Eq.~(\ref{min}).  In fact,
the Metropolis algorithm is only one of the possible sampling methods
for the Monte Carlo evaluation of the integral, see
Ref.~\cite[Sec.~3.7]{ka:monte}.

For the system to reach thermal equilibrium it needs to satisfy the
condition of detailed balance,
\begin{equation}
 K(C_2\mid C_1)P(C_1)=K(C_1\mid C_2)P(C_2).
\end{equation}
Here $P(C)$ is the probability to find the system in the
configuration, or state, $C$, and $K(C_2\mid C_1)$ is the conditional
probability to move from $C_1$ to $C_2$.  The conditional probability
$K$ is usually decomposed as
\begin{equation}
 K(C_2\mid C_1)=A(C_2\mid C_1)T(C_2\mid C_1),
\end{equation}
where the transition probability $T(C_2\mid C_1)$ can be chosen to be
any normalized distribution. It is used to select a random trial move
from $C_1$ to $C_2$.  The complications are kept in $A(C_2\mid C_1)$,
which gives the probability of accepting this move and is the
correction to the arbitrarily chosen $T(C_2\mid C_1)$. The key element
of the algorithm is the evaluation of the function $A(C_2\mid C_1)$ by
a rejection technique. Thus the function $T(C_2\mid C_1)$ is sampled,
and the resulting configuration is accepted or rejected depending on
the value of $A(C_2\mid C_1)$. One usually defines
\begin{equation}
 q(C_2\mid C_1)=\frac{T(C_1\mid C_2)P(C_2)}{T(C_2\mid C_1)P(C_1)}\ge 0
\end{equation}
and
\begin{equation}
 A(C_2\mid C_1)={\rm min}(1,q(C_2\mid C_1)).
\label{accept}
\end{equation}
If the configuration $C_2$ has a lower energy than $C_1$, it is
accepted.  Otherwise, it is accepted with the probability $q(C_2\mid
C_1)$. This procedure is repeated a large number of times and
eventually the system reaches an equilibrium. Here we define an
equilibrium to be the ensemble of states where the average of the
energy does not show systematic changes.

After $L$ steps, equilibrium is established and the system fluctuates
around $\langle\mathcal{F}\rangle$. The thermal average is approximated
by the sum
\begin{equation}
 \langle\mathcal{F}\rangle=\frac{1}{N}\sum_{i=L+1}^{L+N}\mathcal{F}(C_{i}),
\end{equation}
where $C_{i}$ is a state at thermal equilibrium and $N$ is the number
of iterations over which we compute the average. We can interpret
every trial move to represent a unit of quasi-time having
passed. This can not be converted to real units of time, but it is
possible to average thermodynamic properties over quasi-time when a
system is in equilibrium. It is possible to compute the mean of the
quantity because the unit of measurement, which is the number of trial
moves, cancels.  If a trial move is rejected, the old configuration
has to be counted in any averages.

For all the examples to be discussed, $P(C)=\e^{\beta E[C]}$, where the
temperature is defined by $\beta=(k_B T)^{-1}$. If the transition
probability $T(C_{2}\mid C_{1})$ is chosen to be uniform, $q(C_{2}\mid
C_{1})=\e^{\beta\left(E[C_{2}]-E[C_{1}]\right)}$, or
\begin{equation}
 q(C\mid C_\mathrm{new})=\e^{\beta(E-E_\mathrm{new})}.
\label{prob:trans}
\end{equation}
Here, $E$ is the energy of the system in the present configuration $C$
and $E_\mathrm{new}$ is the energy of the new configuration $C_\mathrm{new}$,
that was obtained through a random change in the state of the
system (sampled from the distribution $T$). If the energy of the new
configuration $C_\mathrm{new}$ is lower, the change is accepted. If the
energy is higher, the system accepts this upward step with a
transition probability $q$.  Thus the system can escape local minima
and achieve thermal equilibrium.  In Fig.~\ref{metropolis} we show a
flow diagram representing the Metropolis algorithm. In this diagram,
$\mathcal{F}$ denotes the functional being minimized.


\begin{figure}
\begin{picture}(460,275)(-100,110)
\put(100,350){\framebox(200,20){Randomly PERTURB $C \longrightarrow
C_{\rm new}$}} \put(200,350){\vector(0,-1){30}}
\put(100,300){\framebox(200,20){Change $\Delta=\mathcal{F}(C)-\mathcal{F}(C_\mathrm{new})$}}
\put(200,300){\vector(0,-1){80}}
\put(205,240){negative} \put(-100,250){\framebox(150,20){ACCEPT
Change: $C=C_\mathrm{new}$}} \put(200,260){\vector(-1,0){150}}
\put(125,265){positive} \put(-25,270){\line(0,1){90}}
\put(-25,360){\vector(1,0){125}}
\put(100,200){\framebox(200,20){Probability of transition
$q=\e^{\beta\Delta}$ }} \put(200,200){\vector(0,-1){30}}
\put(100,150){\framebox(200,20){(Random Number $\in$ $[0,1])>q$?}}
\put(100,160){\line(-1,0){125}} \put(40,165){yes}
\put(-25,160){\vector(0,1){90}} \put(200,150){\line(0,-1){30}}
\put(200,120){\line(1,0){150}} \put(300,125){no}
\put(350,120){\line(0,1){240}} \put(350,360){\vector(-1,0){50}}
\end{picture}
\caption{The Metropolis algorithm: Scheme for thermal
equilibrium.}\label{metropolis}
\end{figure}

\subsection{Simulated Annealing}

In 1982, Kirkpatrick and others \cite{sim:science} observed a deep
analogy between annealing of solids and optimization or minimization
problems.  A solid that is sufficiently slowly cooled down, i.e.~at
thermal equilibrium at each temperature, will reach its ground state.
If the energy is the functional to be optimized or minimized, the SA
scheme should find the minimum energy function of this functional
according to a statistical proof by Geman and Geman \cite{geman}. One
starts out with a configuration at a high temperature and runs the
Metropolis algorithm. The size of change of the functional is
proportional to the temperature. Once we have reached thermal
equilibrium, the temperature is decreased according to a cooling
schedule and the procedure repeated as often as necessary.  SA is a
conceptually easy-to-understand minimization technique.

There are several varieties of SA algorithms, each designed to speed
up the minimization of a particular problem.  The application to a
minimization of a continuous problem deserves some reflection on the
discretization of the derivatives. The most important question is the
cooling schedule; a sufficiently slow cooling is crucial for the
``statistical proof of convergence''. Quite often, a slow cooling is
not needed to reach  the global minimum and a faster cooling schedule
can be used: this SA version  is called simulated quenching. For more
information on SA, we recommend reading Ref.~\cite{ingber} and
Ref.~\cite[Sec.~10.9]{num}.  There is no unique way of implementing
the SA scheme and there exits ample opportunity to improve the
code. However, every SA implementation faces the same issues.

{\bf Which initial guess?} Unlike the case of the relaxation method,
the initial configuration is not important as the system should be
able to jump out of local minima. However, an initial guess close to
the global minimum solution can lead to a reduction of the running
time.

{\bf What sampling method to use?} The changes to $C$ should be made
such that the configuration space is well sampled. As the temperature
decreases, so should the size of the changes. Usually the choices made
are random in the configuration space with a Gaussian  or Lorentzian
distribution. However, since we are dealing with functionals rather
than functions such steps are expensive to compute. Therefore, we
restrict  ourselves to changes at individual gridpoints.

{\bf At what initial temperature to start?} A high temperature puts
the system in thermal equilibrium quickly, but at too high temperature
the soliton unwinds. Too low a choice for the initial temperature can
leave the system in a local minimum. The best choice is found by trial
and error.

{\bf When is equilibrium reached?} The determination of the
equilibrium position is crucial and a statistical study on the changes
of the system is essential. Equilibrium is reached when the energy of
the system fluctuates, but does not show any systematic trend.
However, meta-stable states have been observed (like glasses), where
the state of the system changes so slowly that one runs the risk to
interpreting it as being constant.

{\bf The choice of cooling schedule.} The temperature should decrease
according to a log rule to assure convergence to the global minimum
(see Ref.~\cite{geman}).  However, it takes a long time to reach
thermal equilibrium with such a cooling schedule.  Often, the cooling
is speeded up by an exponential cooling schedule using big temperature
decreases or a weaker equilibrium condition.

{\bf Discretization of continuous functional $\mathcal{F}$?} It is
important not to use the central difference, because it does not
depend on the function at the center point. However, this problem can
be  overcome by computing the derivatives midway between two
gridpoints. We discuss this later.

{\bf Use of constraints?} Constraints are no problem, because we use
random changes that satisfy the constraints.



\section{Simulated Annealing in 1D}

We have used the sine-Gordon model for our 1D SA implementation,
because it is one of the simplest field theories exhibiting extended
structures and it is exactly solvable, see Ref.~\cite{raj}.  The
sine-Gordon model is also a very good toy model for solitonic quantum
field theories, for the quantum mass correction and the S-matrix are
exactly known. Further, Coleman \cite{coleman:sine} has shown that the
quantum sine-Gordon model and the massive Thirring model are dual to each
other:  the bosonic soliton in sine-Gordon is a fermion in the massive
Thirring model. Finally and most importantly for this paper, the
soliton solution is known exactly and we can compare it to our SA
results.

\subsection{The Sine-Gordon Model}\label{sine}

The sine-Gordon model is described by the lagrangian density
\begin{equation}
\mathcal{L}=\frac{1}{2}\dmu\phi\,\umu\phi-(1-\cos\phi).
\end{equation}
For simplicity, we have set the mass and the coupling to one.  The
lagrangian is invariant under $\phi\longrightarrow\phi+2\pi n$,
$n \in \mathbf{Z}$. Here $\phi(x,t)$ is an angle in field space, the circle
$\mathrm{S}^1$. The field has to go to the vacuum sufficiently fast for the
soliton to be localized and of finite energy. Therefore, we can
identify the spatial infinity in each direction with one single point
and compactify the one-dimensional space $\mathbf{R}^1$ to $\mathrm{S}^1$.
The field theory of the sine-Gordon model can be described by the map
\begin{equation}
 \phi(t): \mathrm{S}^1 \longrightarrow \mathrm{S}^1
\end{equation}
at a given time $t$. This non-trivial mapping gives us the possibility
to partition the space of all possible field configurations into
equivalence classes having the same topological charge or winding
number. We can visualize this concept with a belt. We can trivially
close it or we can twist one side by 180 degrees and close it or we
can anti-twist it by 180 degrees, i.e.~twist it by $-$180 degrees, and
close it. The twist in the belt cannot be undone unless one opens the
belt. Topological solitons can be thought of as twisted field
configurations. The homotopy group $\Pi_1(\mathrm{S}^1)=\mathbf{Z}$ describes the
twists in the map.  For example, if we twist the belt twice and then
anti-twist it twice, we get back to an untwisted belt: very much like
an annihilation process in particle physics.  The `twist', i.e.~the
topological charge, is fixed by boundary conditions and conserved.

The corresponding Euler-Lagrange equation for the sine-Gordon model is
\begin{equation}
\ddot{\phi}-\phi''+\sin\phi=0.
\end{equation}
One can find the minimal energy solution by solving the static version
using theoretical or numerical methods. The 1-soliton, i.e.~minimal
energy solution of topological charge one, can be derived from the
Bogomolnyi equation (see Ref.~\cite[Sec.~2.5]{raj}):
\begin{equation}
\phi'=\pm \sqrt{2(1-\cos\phi)}.
\end{equation}
Rewriting this in terms of $\sin(\phi/2)$, integrating and inverting
the resulting relation, one finds
\begin{equation}
\phi_{\rm st}(x)=4\arctan [\exp\left(x+x_0\right)].
\label{sine:soliton}
\end{equation}
We can derive the solitons with higher charge via a Backlund
transformation \cite{raj}.  The static minimal energy solution
satisfies the boundary condition $\phi(-\infty)=0$ and
$\phi(\infty)=2\pi$, the field winds around the field sphere
$\mathrm{S}^1$ once. The expression of the energy density is
\begin{equation}
\mathcal{E}(x)=4\sin^2\left[\phi_{\rm{st}}(x)\right].
\end{equation}
The energy goes to zero at spatial infinity and the integral is
finite.  The total energy is $\int\!\mathcal{E}(x)~\d x=8$. In the next
section, we discuss the 1-soliton, Eq.~(\ref{sine:soliton}), and
calculate its total energy using the SA scheme.

\subsection{Implementation of Simulated Annealing}

Three aspects of the SA implementation  are crucial for successful
minimization: the derivatives, the sampling method and the cooling
schedule with thermal equilibrium.

\subsubsection*{Derivatives}

The most accurate discretized derivative is the centered difference.
However, this causes problems with derivative terms as it does not
depend on the function at the point where the energy is being
evaluated.  This results in a decoupling between neighboring points
which gives rise to two independent sub-lattices.  The configuration
becomes spiky since the values jump between the two sub-lattices.  To
avoid this problem, the energy is computed between the gridpoints
rather than at the gridpoints.  The value of the function between the
gridpoints is taken to be the average of the values at the surrounding
points.
\begin{eqnarray}
\phi(x_{i+\frac{1}{2}}) & = & \frac{\phi(x_i)+\phi(x_{i+1})}{2}, \\
\label{central}
\frac{\partial\phi(x_{i+\frac{1}{2}})}{\partial x}
& = & \frac{\phi(x_{i+1})-\phi(x_{i})}{\d x}.
\end{eqnarray}

\subsubsection*{Sampling}

Typically, SA is used to minimize a function $f(x)$ with $x$ being a
vector.  The general form of a change to a configuration is
\[ x_i\rightarrow x_i+M_{ij}U_j \]
where $M_{ij}$ is a matrix, and $U_i$ is a vector of random numbers
satisfying an appropriate probability distribution, see
Ref.~\cite{sa:cont}.  The matrix $M_{ij}$ needs to be chosen such that
the configuration space is well-sampled.  Information from the cooling
process can be used to dynamically adjust $M_{ij}$.  We are interested
in minimizing energy functionals on a lattice of $N$ gridpoints so our
$x$ vector will have $N$ components.  This makes calculating a new
configuration quite an intensive process.  To simplify matters we
sweep across the grid changing individual points at a time.  The
random numbers $U_i$ are taken from a Lorentzian distribution, rather
than a Gaussian distribution.  This is a quite common modification to
the original SA algorithm as the Lorentzian has a longer tail.  The
mean and width of the distribution need to be chosen so that we get a
good sampling of configuration space.  A narrow distribution will only
sample the local neighborhood while a wide distribution will spend too
much time probing irrelevant configurations.  To achieve a good
balance the width is adjusted so that 50\% of all the proposed new
configurations are accepted (this is called the acceptance rate).  If
the mean is taken to be linearly dependent on the temperature then the
acceptance rate will remain roughly constant throughout the cooling.
This leaves the constant of proportionality to be determined at the
start of the cooling process.

\subsubsection*{Cooling Schedule and Thermal Equilibrium}

We use an exponential cooling schedule; the temperature is decreased
by a fixed ratio at each cooling step.  This violates Geman and
Geman's statistical guarantee of reaching the minimum solution.  Since
we do not expect many local minima in 1D this should not be a problem
here.

There are several approaches to determining whether the configuration
is in equilibrium.  A popular one is based on a sliding average, also
known as binning, where the mean energy calculated over a number of
iterations is monitored to see whether it has converged to a fixed
value.  We employ a simpler, related condition which monitors the
lowest energy obtained during a set sequence or chain of iterations
until no new low from one chain to the other is found.  Since
equilibrium is by its very nature statistical it is important to know
how many iterations need to be sampled.  A ballpark figure seems to be
10 to 15 samples per point on average.  We chose the total number of
points in the chain to be 10$N$, where $N$ is the number of gridpoints
and checked empirically that this was enough.


\subsection{Results of Simulated Annealing}

\begin{figure}
\begin{center}
\input{Cooling1D.tex}
\end{center}
\caption{1D SA cooling for the sine-Gordon model.}
\label{cooling1D}
\end{figure}

\begin{figure}
\begin{center}
\input{CoolingCurve.tex}
\end{center}
\caption{Typical cooling curve for 1D SA in the sine-Gordon model.}
\label{coolingcurve}
\end{figure}

Specifically, we look here at the sine-Gordon model where we need to
minimize the following energy functional
\begin{equation}
E[\phi]=\int_{x=-a}^{x=a}\!\d x\;
\left(\frac{1}{2}\partial_i\phi\,\partial_i\phi+1-\cos\phi\right).
\end{equation}
We impose a winding or topological charge of one by setting
$\phi(-a)=0$ and $\phi(a)=2\pi$. We could use a 2D constrained field
to represent $\mathrm{S}^1$, i.e.~$\ph=(\phi_1,\phi_2)$ with
$\ph\cdot\ph=1$. However, we opted for an angle representation,
because it allows us to use fixed boundaries and the soliton cannot
unwind. In the 2D constrained coordinates, the winding is a twist in
the configuration over the whole grid and a very big fluctuation
induced by a high temperature undoes the twist.

We use different grid sizes. In Figs.~\ref{cooling1D} and
\ref{coolingcurve}, we represent different aspects of the cooling of a
sine-Gordon soliton. We start out with an initial field configuration,
here a straight line satisfying the boundary conditions. We then heat
up the configuration until thermal equilibrium is reached (thick solid
line in Fig.~\ref{cooling1D}). We cool it down by slowing decreasing
the temperature after reaching thermal equilibrium. Finally, we obtain
a minimal energy solution close to Eq.~(\ref{sine:soliton}).

We set the acceptance rate to $50\%$ and take $10N$ for the length of
the monitored chain to test thermal equilibrium. A change in these
values effects the speed and quality of minimization. We have re-run
the minimization under the same setting and find the same energy
value. This is a good indicator that the monitored chain is long
enough to obtain thermal equilibrium.  Our box size is 20 units. The
more points we use the closer the result becomes to the exact soliton
energy, namely eight, see Table~\ref{1dtable}.

To conclude, we find our SA code to be a very convenient minimization
technique in 1D. We have successfully tested our 1D SA code on many
different models.  The implementation of the SA search for solitons is
faster, for we did not have to derive the  Euler-Lagrange equation.

\begin{table}
\begin{center}
\begin{tabular}{|c||c|c|} \hline
Points & Energy \\ \hline
51 &  8.035 \\
101 & 8.009 \\
201 & 8.002 \\
301 & 8.001 \\ \hline
\end{tabular}
\end{center}
\caption{Energy versus number of points used for the sine-Gordon model.}
\label{1dtable}
\end{table}


\section{Simulated Annealing in 2D}

We have used the baby Skyrme model \cite{baby_pots} for our 2D SA
implementation, because there are exact and numerical solutions
available which we can compare to our SA results. The baby Skyrme
model is used to study some aspects of the quantum Hall effect (see
Ref.~\cite{qhe}) and is a convenient (2+1)D toy model for the (3+1)D
nuclear Skyrme model (see Ref.~\cite{holzwarth:review}) which requires
much more computational resources.

\subsection{The Baby Skyrme Model}

The nonlinear $\sigma$-model is described by the lagrangian
\begin{equation}
\mathcal{L}=\frac{1}{2}\dmu\ph\cdot\umu\ph
\label{sigmamodel}
\end{equation}
where $\ph$ is a three-dimensional field vector on the sphere $\mathrm{S}^2$
i.e.~$\ph\cdot\ph=1$. The field at a time $t$ is a map
\begin{equation}
 \phi(t): \mathrm{S}^2 \longrightarrow \mathrm{S}^2
\end{equation}
and the associated homotopy group is $\Pi_2(\mathrm{S}^2)=\mathbf{Z}$. The
existence of the topological charge, i.e.~the twisted field
configuration representing a soliton, is ensured by topology. It is
given by
\begin{equation}
B=\frac{1}{8\pi}\int\!\d^2x~{\epsilon}^{\mu\nu}\ph\cdot(\dmu\ph\times\dnu\ph).
\end{equation}

However,
we also need to make sure that the soliton has a stable scale.  From
Derrick's theorem \cite{derrick}, the energy functional corresponding
to Eq.~(\ref{sigmamodel}) is scale invariant.  A change of scale does
not change the energy and therefore numerical errors can significantly
change the scale of the soliton. It is therefore necessary to add
extra terms to stabilize the soliton, fixing the scale.  If we were to
extend the model to (3+1)D, the $\sigma$-model term would lead to an
expanding soliton and a balancing term needs to be added to ensure
stability.

The baby Skyrme model is a modified version of the $\mathrm{S}^2$
$\sigma$-model and the lagrangian is
\begin{equation}
\mathcal{L}=\frac{1}{2}\partial_{\mu}\ph\cdot\partial^{\mu}\ph
-\theta_{S}\left[(\partial_{\mu}\vec{\phi} \cdot
\partial^{\mu}\vec{\phi})^{2} -(\partial_{\mu}\vec{\phi} \cdot
\partial_{\nu}\vec{\phi}) (\partial^{\mu}\vec{\phi} \cdot
\partial^{\nu}\vec{\phi})\right] - \theta_V V(\ph).
\end{equation}
The addition of a potential and a Skyrme term to the lagrangian
ensures stable solitonic solutions. The Skyrme term has its origin in
the nuclear Skyrme model and the baby Skyrme model can therefore be
viewed as its (2+1)D analogue.  Further, in (2+1)D, a potential term
is necessary in the baby Skyrme models to ensure stability of
skyrmions; this term is optional in the (3+1)D nuclear Skyrme model.
One drawback of the model is that the potential term is free for us to
choose. The most common choices are: $V=(1+\phi_3)^4$ (the holomorphic
model has an exact 1-skyrmion solution, see
Ref.~\cite{pz:93}), $V=(1+\phi_3)$ (a 1-vacuum potential
studied in Ref.~\cite{psz:95a}) and $V=(1-\phi_3)(1+\phi_3)$ (a
2-vacua potential studied in Ref.~\cite{baby_pots}). Except for the
first choice, no closed form minimal energy solutions are known. The
baby Skyrme model is a non-integrable system and explicit solutions to
its resulting differential equations are nearly impossible to
find. Numerical methods are the only way forward.


\subsection{Implementation}\label{ThePartitionFunction}

Our 2D and 3D SA implementations originate from a more general
framework, the study of phase transitions in topological systems at
finite temperature and densities \cite{os:finite}.  The thermodynamic
partition function describes a system at a given temperature and can
only be evaluated numerically in the Skyrme models. The evaluation of
thermal averages, as discussed before, can be done with Monte Carlo
techniques and the Metropolis principle is one of the possible
sampling techniques. Conveniently, the thermal average of the energy at
zero temperature is equivalent to the minimal energy of the energy
functional to be minimized.

We start out with the grand canonical partition function
\cite{ma:statistical}
\begin{equation}
\mathcal{Z}(\beta,V,\mu) =\int_{\rm all~x_1,\dots,x_N}
\int_{\rm all~p_1,\dots,p_N}
\d^nx_1 \, \d^np_1 \cdots \d^nx_N \, \d^np_N \sum_{i=0}^{N}\exp[\beta(\mu_i-E_i)],
\label{part}
\end{equation}
where $E_i$ is the energy of the $i$th particle system at temperature
$\beta=(k_B T)^{-1}$ and $V$ is the integration range.  The integral
ranges over all phase space and is $2nN$ dimensional, where $n$ is the
number of space dimensions.  The thermodynamic partition function for
the baby Skyrme model has the form
\begin{equation}
\mathcal{Z}(\beta,V,\mu)=\int\prod_{k}\delta(\ph_{k}\cdot\ph_{k}-1)\,\d^{3}\phi_{k}
\left(\sqrt{\frac{1+\partial_{i}\ph_{k}\cdot\partial_{i}\ph_{k}}{\det(\mathcal{M}_{k})}}\right)\exp[-\beta(\mathcal{V}_{k}-\mu\mathcal{B}_{k})].
\label{BabyThermoPart}
\end{equation}
Here $\mathcal{M}$ is the mass density matrix, $\mathcal{V}$ the potential
energy density, and $\mathcal{B}$ is the topological charge density. The input
parameters of the thermodynamic partition function are the temperature
$\beta$, the volume of the system $V$ and the chemical potential $\mu$.
The $\delta$-function is required due to the constraint on the $\ph$ field.

At zero temperature, the factor in front of the exponent becomes
irrelevant and if we further set $\mu=0$,  the thermodynamic partition
function reduces to  the integral
\begin{equation}
 Z=\int\prod_{k}\delta(\ph_{k}\cdot\ph_{k}-1) \, \d^{3}\phi_{k}\exp(-\beta\mathcal{V}_{k}),
\label{Zfunction}
\end{equation}
where $\mathcal{V}_k$ is the potential energy density.  The
implementation of $\mathcal{Z}$ is similar to the implementing of $Z$,
but $\mathcal{Z}$ contains information that is not necessary for finding
minimal energy solutions The value of $Z$ is not of interest to us,
but the $\ph$ distribution as $\beta\rightarrow\infty$ is.  The
probability density function which is sampled at every point $k$ when
applying Monte Carlo is the sum over all neighbors of $k$,
\begin{equation}
P_{k}=\exp\left(-\beta\sum^{\stackrel{\rm
{number~of}}{\rm{neighbors}}}_{i=1} \mathcal{V}_{i}\right).
\label{pfunction}
\end{equation}
Examples of Monte Carlo calculations using a grand canonical ensemble
are given in Refs.~\cite{bi:monte,fr:understanding}. A good discussion
of possible errors and how to deal with them is given in
Ref.~\cite{fr:understanding}.


\subsubsection*{Monte Carlo for the Baby Skyrme Model}
\label{MCBaby}\label{MCBabyImportance}

We apply the Metropolis principle in the simplest possible way and
select a new vector $\ph_{\rm new}(x_k)$ at a gridpoint $k$ from a uniform
probability distribution function over the unit sphere, as the
integration measure with the $\delta$-function implies. We choose each
of the components $\phi_{a}$ uniformly between $-1$ and $1$. If
the sum of the squares of the components
$\phi_{1}{}^{2}+\phi_{2}{}^{2}+\phi_{3}{}^{2}$ is larger than one the
sample is rejected. All accepted vectors are scaled to obtain unit
length \cite{ka:monte}. The transition probability of this simple
method is
\begin{equation}
 T(C_2\mid C_1)=\left\{\begin{array}{cl}
\frac{1}{\rm surface~area~of~sphere} & {\rm on~the~unit~sphere} \\
0 & {\rm otherwise}
\end{array}\right.,
\end{equation}
and therefore
\begin{equation}
 q(C_2\mid C_1)=\frac{P(C_2)}{P(C_1)},
\end{equation}
where the present vector $\ph_{\rm present}(x_k)\in C_1$ and the
newly selected $\ph_{\rm new}(x_k)\in C_2$. The quantity $q$ is the
new integrand of the partition function (\ref{pfunction}) divided by
the present one. It is easiest to test a trial move for one gridpoint
at a time, although other methods will be discussed. The acceptance
probability defined by (\ref{accept}) is calculated by
\begin{equation}
 A(C_2\mid C_1)=\min\left(1,\exp\left[-\beta\left(\mathcal{V}_{\rm new}-\mathcal{V}_{\rm present}\right)\right]\right).
\label{Zaccept}
\end{equation}
A change to the vector $\ph$ at lattice point $k$ modifies the
potential energy on the gridpoint $k$ and its neighbors only (using a
linear approximation for the derivatives). This is all the information
needed to apply the Metropolis method.  The quantities of interest to
measure are the potential energy $E$ and the topological charge $B$
which should be conserved and is a check on the numerics.

A uniform sampling of the distribution has an extremely low acceptance rate;
too many vectors are rejected. We therefore use a biased sampling
technique where a new vector $\ph_{\rm new}$ is sampled near $\ph_{\rm
present}$, the vector that it is supposed to replace. We sample
vectors in an intrinsic frame where the $z$-axis corresponds to the
present vector.  The vector
\begin{equation}
 \vec{n}^{\rm int} =(n_{1}^{\rm int},n_{2}^{\rm int},n_{3}^{\rm int})
=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),
\label{vecn}
\end{equation}
gives the components of the new vector in cartesian coordinates in the
intrinsic frame. The Euler angles that define the rotation from the
intrinsic frame to the laboratory frame are given by
\begin{equation}
 \vec{k}^{\rm lab}=(\cos\beta\sin\alpha,\sin\beta\sin\alpha,\cos\alpha).
\end{equation}

The angles $\alpha$ and $\beta$ are used to rotate the new vector
$\vec{n}^{\rm int}$ to the laboratory frame,
\begin{equation}
 \vec{n}^{\rm lab}=\left(\begin{array}{c}
n_{1}^{\rm lab} \\
n_{2}^{\rm lab} \\
n_{3}^{\rm lab}
\end{array}\right)=
\left(\begin{array}{ccc}
\cos\alpha\cos\beta & -\sin\beta & \sin\alpha\cos\beta \\
\cos\alpha\sin\beta & \cos\beta & \sin\alpha\sin\beta \\
-\sin\alpha & 0 & \cos\alpha
\end{array}\right)
\left(\begin{array}{c}
n_{1}^{\rm int} \\
n_{2}^{\rm int} \\
n_{3}^{\rm int}
\end{array}\right).
\label{BabyRotationMatrix}
\end{equation}
In the previous terminology, the new unit field vector at the
gridpoint $k$ is $\ph_{\rm new}(x_k)=\vec{n}^{\rm lab}$, and it is a
vector selected from a particular probability distribution function
which is rotationally symmetric about the present vector $\ph_{\rm
present}(x_k)~(=\vec{k}^{\rm lab})$. The vectors $\ph_{\rm new}(x_k)$
and $\ph_{\rm present}(x_k)$ are inserted into the acceptance
probability. The angle $\theta$ is sampled uniformly on $[0,A)$, where
$A\le\pi$, and $\phi$ is sampled uniformly on $[0,2\pi)$.  The
corresponding transition probability is
\begin{equation}
T(C_2\mid C_1)=\left\{\begin{array}{cl}
\frac{1}{2\pi A} & {\rm if}~1-\cos\theta<A \\
0 & {\rm otherwise}
\end{array}\right..
\label{trans}
\end{equation}
No importance sampling is imposed and the probability distribution
function is uniform, so that the acceptance rate (\ref{Zaccept}) can
be applied directly. The optimal value of $A$, which we are free to
choose, depends on the particular configuration, especially on the
topological charge and the temperature $\beta$. No $A$-dependency was
discovered in any ensemble averages other than the acceptance
rate. Therefore, we believe that this method is reliable and
efficient. We allow $A$ to vary automatically to achieve an acceptance
rate near 40\%. The choice of $A$ influences the rate at which
equilibrium is reached, which is defined as the absence of change in
the average energy over a large number of steps.  At each temperature,
the system was required to reach equilibrium before being cooled
further. This ensures that cooling does not occur too quickly.

The choice (\ref{trans}) only samples a portion, rather than the whole
of the unit sphere. This is a valid method because we are modeling a
continuous system, and therefore the vectors can reach any region in a
number of steps. As $A$ is varied automatically, the whole unit sphere
is sampled for high temperatures.  If the region which is unsampled
for low temperatures were sampled, the vectors selected there would
have virtually zero probability of being accepted. For generality, we
discuss a more rigorous method using importance sampling in the
appendix. There, new vectors are chosen from a gaussian (or other)
distribution centered around the present vector. Importance sampling
allows vectors from all over the unit sphere to be selected at any
temperature, and this might be necessary for some systems, especially
when calculating thermal averages.  The disadvantage of importance
sampling over restricting the transition probability is the increased
amount of computing time.

\subsubsection*{Calculation of Field Derivatives for a Field on $\mathbf{S^2}$}

We calculate derivatives in a similar way to the 1D implementation
where measurable quantities are calculated in the center of the
plaquettes.  An illustration of a plaquette is given in
Fig.~\ref{Plaqpic}, where the field vectors are evaluated at the
intersections of the lines and all measurable quantities are
calculated at the midpoints.
\begin{figure}
\begin{center}
\includegraphics[width=3cm]{Plaqpic.eps}
\caption{A picture of the plaquette, where the fields are evaluated at
the intersections of the lines and the measured quantities calculated
at the midpoints $\times$.}
\label{Plaqpic}
\end{center}
\end{figure}
If a field vector is altered, using the Monte Carlo method as
described above, the effect of a change has to be calculated on the
four surrounding midpoints.

All field vectors at each of the four gridpoints lie on a unit sphere
and the average of the four field vectors must also be of unit length,
because the topology requires unitarity everywhere.  A simple average
fails this criterion unless all vectors point in identical directions.
We still use the average, corrected by scaling it to unit length.

This also impacts the calculation of derivatives. For example, the
error in taking the $x$-derivatives is minimized if
$\frac{1}{2}\sum_{b=\{b_1,b_2\}}\phi_b (x_{n+1})$ and
$\frac{1}{2}\sum_{a=\{a_1,a_2\}}\phi_a (x_{n})$ are scaled to unit
length before the latter is subtracted from the former to obtain the
$x$-derivative. Here $\{a_1,a_2\}$ are the 2 gridpoints with the
coordinates $(x_{n},y_{m})$ and $(x_{n},y_{m+1})$, and $\{b_1,b_2\}$
are the 2 gridpoints with the coordinates $(x_{n+1},y_{m})$ and
$(x_{n+1},y_{m+1})$. These vectors are on the corners of the
plaquette, see Fig.~\ref{DerivPos}. Thus, the derivative is calculated
by
\begin{figure}
\begin{center}
\includegraphics[width=10cm]{derivpos.eps}
\caption{Illustration of the scaling of the center derivative.}
\label{DerivPos}
\end{center}
\end{figure}
\begin{equation}
 \left.\frac{\partial\phi}{\partial x}\right|_{x_{n+\frac{1}{2}}}={\rm
 Scaled}\left[\frac{1}{2}\sum_{b=\{b_1,b_2\}}\phi_b
 (x_{n+1})\right]-{\rm
 Scaled}\left[\frac{1}{2}\sum_{a=\{a_1,a_2\}}\phi_a
 (x_{n})\right].
\label{deriv}
\end{equation}
This derivative works very well in practice. At very high temperatures
the numerics may break down, because the derivative (\ref{deriv}) is
by definition an underestimate. The 3D analog of this formula is
obtained by replacing `$\frac{1}{2}$' by `$\frac{1}{4}$' and summing
over the four components of $a$ and the four components of $b$.

\subsubsection*{Updating Mechanisms}

In our 1D simulations, we have randomly selected which gridpoint
should be sampled.  In our 2D and 3D implementation, we sweep over the
gridpoints sequentially. The new vectors are stored and the changes to
the field are only updated after a complete sweep over the entire grid
to avoid unwanted sequential correlations. We have split the grid into
four independent subgrids, each labelled by a different symbol in
Fig.~\ref{SubgridN}.
\begin{figure}
\begin{center}
\includegraphics[width=2.5cm]{SubgridNorm.eps}
\caption{The four subgrids for single point changes, each labelled by
a symbol.}
\label{SubgridN}
\end{center}
\end{figure}
The subgrids are chosen at random and, at each sweep, only one of the
possible four composing vectors of the derivatives and field averages
at the midpoints are changed. This avoids the creation of fluctuations
between neighboring vectors for high acceptance rates, which produce
an unphysical increase in energy.  Unfortunately, the change of single
gridpoints at a time does not favor collective motion, where a
localized energy distribution moves in one direction. Changing regions
of gridpoints at a time has proven to be more efficient.

We have successfully used a plaquette updating mechanism. For a given
plaquette, we sample a vector in the intrinsic frame, see
Eq.~(\ref{vecn}). Then this frame is rotated as in Sec.~\ref{MCBaby}
for each of the four vectors on the plaquette, separately.
Fig.~\ref{Plaqpic} illustrates a plaquette surrounded by the affected
nine midpoints.  The total change in the energy density of all these
midpoints is now calculated and all four vectors are accepted or they
are all rejected. Again, we have split the grid into four subgrids,
each shaded differently in Fig.~\ref{SubgridP}. However, some
midpoints are affected by two or four plaquettes from the same
subgrid.  Therefore, we choose the subgrids randomly rather than
sequentially to avoid unwanted correlations.

\begin{figure}
\begin{center}
\includegraphics[width=2.55cm]{SubgridPlaq.eps}
\caption{The four subgrids for plaquette changes, each shaded
differently.}
\label{SubgridP}
\end{center}
\end{figure}

\subsection{Results and Comparison}

We have investigated the three different baby Skyrme models mentioned
before. To that end, we need to minimize the energy functional:
\begin{equation}
E\left[\ph\right]=\int\!\d^2x
\left[\frac{1}{2}\partial_{i}\ph\cdot\partial_{i}\ph
+\theta_{S}\left[(\partial_{i}\vec{\phi} \cdot
\partial_{i}\vec{\phi})^{2} -(\partial_{i}\vec{\phi} \cdot
\partial_{j}\vec{\phi}) (\partial_{i}\vec{\phi} \cdot
\partial_{j}\vec{\phi})\right] + \theta_V V(\ph)\right].
\end{equation}

First, we have looked at the simplest holomorphic model with the
potential $V=(1+\phi_3)^4$. There exists an explicit 1-soliton solution,
\begin{equation}
W= \sqrt[4]{\frac{\theta_V}{2\theta_S}}~(x+\mathrm{i}y),
\end{equation}
where the $W$ field is the stereographic projection of $\ph$ on the
complex plane, given by $W=2(\phi_1+\mathrm{i}\phi_2){(1-\phi_3)}^{-1}$.
We choose $\theta_S=\theta_V=\frac{1}{2}$ where its total energy equals
$4\pi(1+\frac{8}{3\sqrt{2}})\approx 36.2618$. Since the soliton profile
has a polynomial decay, we need a large lattice. With a $350\times350$
grid and lattice spacing $h=0.05$, we obtain $E=36.4890$ and
$B=0.9999$. Here, the energy is slightly higher than the exact
solution because of the finite lattice effects. The holomorphic baby
skyrmion has the slowest decay of any of the models discussed, and
therefore can be seen as the worst case scenario.

\begin{table}
\begin{center}
\begin{tabular}{|c||c|c|c||c|c|c|}
\hline & \multicolumn{3}{|c||}{1-vacuum model} & \multicolumn{3}{|c|}{2-vacua model} \\ \cline{2-7}
Topological & \multicolumn{2}{|c|}{SA} & Ref.~\cite{baby_pots} & \multicolumn{2}{|c|}{SA} & Ref.~\cite{baby_pots} \\
Charge & $B$ & $E/B$ & $E/B$ & $B$ & $E/B$ & $E/B$ \\ \hline
1 & 0.99978 & 19.6505 & 19.47 & 0.99989 & 19.6572 & 19.65 \\
2 & 1.99973 & 18.4452 & 18.27 & 1.99984 & 17.6530 & 17.65 \\
3 & 2.99962 & 18.5257 & 18.34 & 2.99983 & 17.2259 & 17.22 \\
4 & 3.99952 & 18.4014 & 18.22 & 3.99989 & 17.0677 & 17.07 \\ \hline
\end{tabular}
\caption{Baby Skyrme models: Comparison of our SA results with
Euler-Lagrange results.}
\label{table2d}
\end{center}
\end{table}

We have also looked at the baby Skyrme models with one vacuum, where
$V=1+\phi_3$, and two vacua, where $V=1-\phi_3{}^2$. The parameters for
the 1-vacuum model have been fixed to $\theta_S=\frac{1}{4}$ and
$\theta_V=0.1$ and for the 2-vacua to $\theta_S=0.44365$ and
$\theta_V=0.05$ in agreement with existing literature, see
Ref.~\cite{baby_pots}.  We use an $80\times 80$ grid with periodic
boundary conditions and lattice spacing $h=0.4$.  The minimal energy
solution in the first four topological sectors is shown in
Fig.~\ref{Baby1-4} and Fig.~\ref{NewBaby1-4} respectively.  We compare
the energies per charge with the calculations from
Ref.~\cite{baby_pots} in Table~\ref{table2d}.

The results for the 2-vacua model are the same for both studies within
an accuracy of a few parts in $10^{4}$. The results from
Ref.~\cite{baby_pots} have been obtained via the shooting method. We
can apply this accurate method, because the skyrmions are radially
symmetric and the minimization reduces to a 1D problem.  There is a
slight discrepancy in energy when comparing the 1-vacuum model. The
energies in Ref.~\cite{baby_pots} were obtained on a 2D lattice using
a damped time-evolution.  The energy of the 1-skyrmion calculated
using the shooting method is $E=19.65$. Our SA result agrees well with
this value.  The inaccuracies in Ref.~\cite{baby_pots} arise due to a
different derivative approximation, which tends to reduce the
energy. The errors due to the fixed boundaries are also greater in
Ref.~\cite{baby_pots}, although this tends to increase the energy of
the 1-skyrmion.  We believe that periodic boundary conditions, which
have been used for our SA result, have the advantage that the tails of
the skyrmions can spread out further. The finite lattice effects are
still present, as skyrmions could then interact with themselves over
the boundaries. For our SA model, a study \cite{sc:second} on the
1-skyrmion case shows that the gridsize used induces an error of the
order 0.01\%. The error due to not having relaxed the system
properly is 0.01\%. The largest error is due to finite difference
effects and has a possible error of 0.1\%. Because of the successful
agreement of our SA result for the 1-skyrmion with the result of the
shooting method, it is believed that the SA solutions with higher
topological charge are also more accurate than the results quoted in
Ref.~\cite{baby_pots}. Finally, we show shows an example of the
cooling schedule for its 3-skyrmion in Fig.~\ref{BCoolSeq}.
\begin{figure}
\begin{center}
\includegraphics[width=11cm]{Baby1-4Q.eps}
\caption{Energy density plot with identical vertical energy scale. We
show the 1-vacuum skyrmions of charge one (a), two (b), three (c) and
four (d). The range plotted is $2.4\times 1.2$ units in each figure.}
\label{Baby1-4}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=11cm]{NewBaby1-4Q.eps}
\caption{Energy density plot with identical vertical energy scale. We
show the 2-vacua skyrmions of charge one (a), two (b), three (c) and
four (d). The range plotted is $2.0\times 1.0$ units in each figure.}
\label{NewBaby1-4}
\end{center}
\end{figure}

We studied three different baby Skyrme models. Changing from one
potential to the other could not have been easier. In the case of the
iterative techniques, changing the differential equation is in itself
conceptually easy, but, in practice, a lot of time is spent on getting
the coefficients right and checking the derivation.  For future
research, we intend to use SA to do a systematic check on the
multi-skyrmion structure. Specifically, we are interested in the class
of potentials that leads to radially symmetric multi-skyrmions.

\begin{figure}
\begin{center}
\includegraphics[width=16cm]{BCoolSeqQ.eps}
\caption{Energy density plot of identical vertical energy scale. We
show a SA cooling for the 3-skyrmion in the 1-vacuum model: (a) the
starting configuration; (b) system is heated to $\beta=500$ and
skyrmions repel each other because their isospins are initially in the
same direction; (c) isospins rotate relative to each other and
skyrmions attract each other; (d) equilibrium has been reached for
$\beta=500$; (e) system is cooled to $\beta=5000$; (f) minimal energy
solution at $\beta=\infty$.  The range plotted is $1.8\times 1.6$
units in each figure.}
\label{BCoolSeq}
\end{center}
\end{figure}


\section{Simulated Annealing in 3D}

We have chosen the nuclear Skyrme model \cite{holzwarth:review} for
our 3D SA implementation, because we believe that SA is a flexible
tool for exploring the multi-skyrmion structure.

In the sixties, Skyrme constructed an effective  field theory of
mesons where the baryons are the topological solitons of the
theory. Research by 't Hooft and Witten  has established that the
nuclear Skyrme model shows important similarities to the low-energy
effective lagrangian of QCD \cite{thooft:2,witten}. The 1-skyrmion can
be interpreted as a nucleon with reasonable success \cite{anw}.  The
numerical work by Braaten, Townend, and Carson \cite{braaten}, and Battye and
Sutcliffe \cite{multi:num} on the structure of classical
multi-skyrmions supports the idea that an  appropriate quantization of
these minimal energy solutions for a given topological sector could
possibly lead to an effective description of atomic nuclei.  However,
the calculation of quantum properties of multi-skyrmions is very
difficult. Part of this is due to the fact that these minimal energy
solutions are not radially symmetric and the theory is
non-renormalisable.  This is rather frustrating, for the claim that
the Skyrme model, descending from a large $N$-QCD approximation models
mesons, baryons and higher nuclei is a very attractive one.  Numerical
methods are probably the only way forward and the SA scheme might be
useful in exploring further the multi-skyrmion structure for different
versions of the Skyrme model.


\subsection{The Nuclear Skyrme Model}

The nuclear Skyrme lagrangian
\begin{equation}
\mathcal{L}=\frac{1}{2}\partial_{\mu}\ph\cdot\partial^{\mu}\ph
-\frac{1}{4}\left[(\partial_{\mu}\vec{\phi} \cdot
\partial^{\mu}\vec{\phi})^{2} -(\partial_{\mu}\vec{\phi} \cdot
\partial_{\nu}\vec{\phi}) (\partial^{\mu}\vec{\phi} \cdot
\partial^{\nu}\vec{\phi})\right]
\label{FullSkLag}
\end{equation}
is a straightforward extension of the nonlinear $\sigma$ model
containing an additional fourth-order term called the Skyrme term. We
need to include this extra term to ensure stability of the
soliton. The mapping becomes
\begin{equation}
 \phi(t): \mathrm{S}^3 \longrightarrow \mathrm{S}^3.
\end{equation}
More realistic lagrangians should probably include higher order
correction terms. The SA scheme is especially useful, because extra
term can be included trivially.


\subsection{Implementation}

The 3D implementation is very similar to the 2D code and we will only
mention new issues relevant to the 3D case.  On a 4D unit sphere, the
integration measure is given by
\begin{eqnarray}
\int_{\rm unit~sphere}\!\d x~\d y~\d z~\d w &=&
\int_{\theta=0}^{2\pi}\int_{\phi=0}^{\pi}\int_{\chi=0}^{\pi}\!
\sin\theta\sin^{2}\chi ~\d\phi~\d\theta~\d\chi\nonumber\\
&=& \int_{\theta=0}^{2\pi}\int_{\cos\theta=-1}^{1}
\int_{\frac{\chi}{2}-\frac{1}{4}\sin(2\chi)=0}^{\frac{\pi}{2}}\!\d\phi
~\d\cos\theta~\d\left(\frac{\chi}{2}-\frac{1}{4}\sin(2\chi)\right)\nonumber\\
&=& \int_{\theta=0}^{2\pi}\int_{u=-1}^{1}\int_{v=0}^{\frac{\pi}{2}}\!\d\theta~\d u~\d v.
\end{eqnarray}
In order to rotate the new vector from the intrinsic frame to the
laboratory frame, the angle $\chi$ must be evaluated from $v$. The
equation
\begin{equation}
 v=\frac{\chi}{2}-\frac{1}{4}\sin(2\chi)
\end{equation}
can not be rewritten in terms of $\chi$, and therefore it must be
solved numerically. The observation that most solutions will be in the
region of small $v$ due to the importance sampling, implies that a
small $\chi$ approximation might be useful. For small $\chi$,
$v\approx\frac{\chi^{3}}{3}$. Therefore, the inversion is performed
numerically by tabulating $\chi$ uniformly on $[0,\pi]$ against
$v^{1/3}=\left(\frac{\chi}{2}-\frac{1}{4}\sin(2\chi)\right)^{\frac{1}{3}}$.
The $\chi$ value corresponding to the selected $v^{\frac{1}{3}}$ is
found in this table, and a linear interpolation is applied between the
two nearest values of $v^{\frac{1}{3}}$ to give a better approximation
to $\chi$. Using 1000 pre-calculated entries the error on the whole
region $\chi\in[0,\pi]$ is then less than $10^{-7}$. On the region
$\chi\in[0,0.1]$, where almost the entire selection of $\chi$ lies,
the error is less than $10^{-14}$, which corresponds to the precision
of double real numbers. This method does not create any significant
errors and is time-consuming. No faster method seems to be possible.

As the cosines and sines of $\phi$, $\theta$, and $\chi$ are known,
the rotation can be performed using a rotation matrix similar to
(\ref{BabyRotationMatrix}). The new choice of vector in the intrinsic
frame is given by
\begin{equation}
 \vec{n}^{\rm int} =(n_{1}^{\rm int},n_{2}^{\rm int},n_{3}^{\rm int},n_{4}^{\rm int})
 =(\sin\chi\sin\theta\cos\phi,\sin\chi\sin\theta\sin\phi,\sin\chi\cos\theta,
 \cos\chi).
\end{equation}
The $z$-axis in the intrinsic frame coincides with $\vec{k}^{\rm
lab}=\ph_{\rm present}(x_k)$ in the laboratory frame, similarly as in
the baby Skyrme model. The rotation angles $\gamma$, $\alpha$, and
$\beta$ (in that order) are defined by
\begin{equation}
 \vec{k}^{\rm lab}=(\cos\beta\sin\alpha\sin\gamma,\sin\beta\sin\alpha\sin\gamma,
\cos\alpha\sin\gamma,\cos\gamma).
\end{equation}

The transformation $\vec{n}^{\rm int}$ to $\vec{n}^{\rm lab}$ is
performed using the matrix
\begin{equation}
\left(\begin{array}{c}
n_{1}^{\rm lab} \\
n_{2}^{\rm lab} \\
n_{3}^{\rm lab} \\
n_{4}^{\rm lab}
\end{array}\right)=
\left(\begin{array}{cccc}
\cos\alpha\cos\beta & -\sin\beta & \sin\alpha\cos\beta\cos\gamma & \sin\alpha\cos\beta\sin\gamma \\
\cos\alpha\sin\beta & \cos\beta & \sin\alpha\sin\beta\cos\gamma & \sin\alpha\sin\beta\sin\gamma \\
-\sin\alpha & 0 & \cos\alpha\cos\gamma & \cos\alpha\sin\gamma \\
0 & 0 & -\sin\gamma & \cos\gamma
\end{array}\right)
\left(\begin{array}{c}
n_{1}^{\rm int} \\
n_{2}^{\rm int} \\
n_{3}^{\rm int} \\
n_{4}^{\rm int}
\end{array}\right).
\label{SkRotationMatrix}
\end{equation}
For similar reasons as already discussed previously, we use a
restricted transition probability to be uniform over $[0,A)$, where
$A\le\pi$. The method of importance sampling has also been
investigated for the nuclear Skyrme model and is given in the
appendix. The transition probability is therefore
\begin{equation}
T(Y\mid X)=\left\{\begin{array}{cl}
\frac{1}{4\pi A} & {\rm if}~\frac{\chi}{2}-\frac{1}{4}\sin(2\chi)<A \\
0 & {\rm otherwise.}
\end{array}\right.
\end{equation}
The angles are sampled by
\begin{eqnarray}
 v&=&\frac{\chi}{2}-\frac{1}{2}\sin(2\chi)=A\xi_{1},\nonumber\\
 u&=&\cos\theta=2\xi_{2}-1,\nonumber\\
 \phi&=&2\pi\xi_{3},
\label{3dsampling}
\end{eqnarray}
where $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ are three random variables
uniformly sampled from 0 to 1.  The new vector
$\ph_n(x_k)=\vec{n}^{\rm lab}$ is a vector which has been selected
from a uniform probability distribution centered around the previous
vector $\ph_p(x_k)~(=\vec{k}^{\rm lab})$ where only the angle $\chi$
between these vectors has an upper limit. Finally, $\ph_n(x_k)$ and
$\ph_p(x_k)$ are inserted into (\ref{Zaccept}) to find the acceptance
rate.  Similarly to the baby Skyrme models, the value of $A$ is
automatically chosen to have an acceptance rate near 40\%. The cooling
is also controlled in the same manner as described in 2D.

\subsection{Results and Comparison}

\begin{figure}
\begin{center}
\includegraphics[width=16cm]{FullSk1-4Q.ps}
\caption{Plot of the same constant energy density surface. We show the
first four multi-skyrmion solutions (from left to right). All plots
are to the same scale. The mesh spacing of the plotted objects is
0.12 units.}
\label{FullSk1-4}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=15cm]{FullBeq4CoolCutQ.ps}
\caption{Plot of constant energy density surface. We show a SA search
of the $B=4$ skyrmion: (a) the starting configuration of four
1-skyrmions; (b) the system heated to $\beta=500$ where the skyrmions
fuse into one; (c) the system in equilibrium at $\beta=500$ where the
structure emerges; (d) the minimal energy solution at $\beta=\infty$.
All plots are to the same scale. The mesh spacing of the plotted
objects is 0.1 units.}
\label{3Dcooling}
\end{center}
\end{figure}
The 3D implementation of SA is computationally much more intensive
than the 2D case.  The accuracy of our numerical simulations is
therefore reduced due to limited resources in computation time and
memory. The intergrid spacing used is 0.12 Skyrme units and is close
to the upper limit where the numerics break down (at a reasonably high
temperature required to do SA sufficiently fast). The maximum gridsize
that can be used to obtain results in a reasonable time is $80\times
80\times 80$. The finite volume causes an increase of energy for the
1-skyrmion because it repels itself over the periodic boundary, and
therefore induces an error of 1\% \cite{sc:second}. The error due to
not having relaxed the system properly is 0.1\%. The error due to
finite difference effects has a maximum of 0.3\%.
The skyrmions of topological charge one to four are shown in
Fig.~\ref{FullSk1-4}. We show an an example of the cooling schedule
for the 4-skyrmion in Fig.~\ref{3Dcooling}.  These energies per charge
are contrasted in table~\ref{table3d} with those results obtained by
Battye and Sutcliffe \cite{multi:num}. It is very difficult to compare
the results.  The 1-skyrmion solution gives more information for
comparison. It is spherically symmetric and the shooting method in the
hedgehog ansatz can be used.  The energy of a 1-skyrmion minimized in
the 1D SA code is $E=73.12$ (in the continuous limit). We are not sure
if the result for the 1-skyrmion in Ref.~\cite{multi:num} is truly
more accurate or just a coincidence. Their topological charge, an
indicator for the discretization error, is certainly less accurate.
The finite lattice effect increases the energy of the 1-skyrmion and
therefore the energy obtained should be larger than 73.12.
Unfortunately, we do not know which lattice parameters they have used,
making a good comparison impossible. However, we get the same minimal
energy structure.
\begin{table}
\begin{center}
\begin{tabular}{|c||c|c|c|c|} \hline
Topological & \multicolumn{2}{|c|}{SA} & \multicolumn{2}{|c|}{Ref.~\cite{multi:num}} \\
Charge & $B$ & $E/B$ & B & $E/B$ \\ \hline
1 & 1.0015 & 73.75 & 0.984 & 72.96 \\
2 & 2.0030 & 70.31 & 1.972 & 69.34 \\
3 & 3.0042 & 68.52 & 2.960 & 67.69 \\
4 & 4.0048 & 66.30 & 3.948 & 66.09 \\ \hline
\end{tabular}
\caption{The nuclear Skyrme model: SA results versus published results.}
\label{table3d}
\end{center}
\end{table}


\section{Conclusion}

We have shown that SA is an alternative way of finding the minimal
energy solution in a given topological charge sector. We
independently confirmed the validity of the studies using the standard
minimization techniques. It is very hard to objectively compare the
different approaches. However, we have found SA to be a more
convenient and flexible minimization technique. The implementation
and fine-tuning of our SA codes took a fair amount of time due to a
lack of prior research in this area. In comparison to other methods,
we are confident that future implementations will take us
considerably less time. We did not find any significant differences
in speed of minimization. The SA codes can be made faster by
fine-tuning the cooling parameters. We prefer SA minimization because
of its ease of use. Speed considerations are irrelevant in 1D and 2D
and we can use parallel computing in the 3D case.
There are several areas we want to look at next. First of all, we will
optimize SA by using more sophisticated update and cooling mechanisms
and by parallelization. We are also currently investigating the
possibility of doing time-evolution via SA minimization of the action.
At the same time, we intend to look at a wide range of models.  We
shall investigate the multi-skyrmion structure of several baby Skyrme
models. Research is also underway in the use of symmetry
breaking terms for the nuclear Skyrme model. Moreover, the 2D and 3D
code will be used to study phase transitions in the baby and nuclear
Skyrme model at finite temperatures and finite densities
\cite{os:finite}. To conclude, SA is a flexible tool; all we really
need is an energy functional to minimize!


\section*{Acknowledgements}

First of all, the authors would like to thank Niels Walet for his
support and advice. We would also like to acknowledge discussions with
Klaus Gernoth, Laur J\"arv and Bernard Piette. The authors acknowledge
PPARC (M.~H.) and EPSRC (O.~S.~\& T.~W.) for research grants.

\appendix
\subsection*{Appendix: Importance Sampling}

For importance sampling, the integral (\ref{average}) is rewritten as
\begin{equation}
 \langle\mathcal{F}\rangle=\int\left[\frac{\mathcal{F}(C)P(C)}{\tilde{P}(C)}\right]\tilde{P}(C)\,\d C.
\end{equation}
Here $P(C)$ can be normalized without loss of generalization, see
Ref.~\cite{ka:monte}, and $\tilde{P}(C)$ is a different probability
density function which satisfies
\begin{equation}
 \tilde{P}(C)\ge 0,~~~~~\int\tilde{P}(C)\,\d C=1,
\label{constrftilde}
\end{equation}
and
\begin{equation}
 \frac{\mathcal{F}(C)P(C)}{\tilde{P}(C)}<\infty
\end{equation}
except on a countable set of points. In this method  one chooses a
$\tilde{P}(C)$ that minimizes the variance which is
\begin{equation}
{\rm var}\{\langle\mathcal{F}\rangle\}=\int\frac{\mathcal{F}^{2}(C)P^{2}(C)}
{\tilde{P}^{2}(C)}\d C-{\langle\mathcal{F}\rangle}^{2}.
\end{equation}
The measurement of statistical accuracy is given by
${\rm var}\{\langle\mathcal{F}\rangle\}$.
More samples reduce the variance. Alternatively, the same
variance using fewer samples can be achieved with importance
sampling. In practice, the closer $\tilde{P}(C)$ is  to $\mathcal{F}(C)P(C)$,
the lesser the variance becomes.  It is known that if
\begin{equation}
 \tilde{P}(C)=\frac{\mathcal{F}(C)P(C)}{\langle\mathcal{F}\rangle},
\label{ftildemin}
\end{equation}
then the integral is equal to $\mathcal{F}$ with zero variance. However,
we need to respect the constraints (\ref{constrftilde}), and even
worse, choosing a $\tilde{P}$ requires knowledge of $\langle\mathcal{F}\rangle$ prior
to evaluating the integral.

\subsubsection*{Baby Skyrme models}

The application of importance sampling to the partition function
(\ref{Zfunction}) is complicated, because the range of integration is
on a sphere of unit length. Therefore, we need to use a
$\tilde{P}(\ph(x_k))$ that is only non-zero  on the unit
sphere. Further, the maximum or most likely area of accepted values
depends on the present vector, and therefore $\tilde{P}(\ph_n(x_k))$
should not be restricted to a certain region of the sphere, but should
depend on the present vector $\ph_p(x_k)$. Looking at the results of a
uniform probability distribution function on the sphere, a Gaussian
distribution of the polar angle $\theta$ seems to be a good choice for
$\tilde{P}$, where $\theta=0$ is in the direction of the present
vector. This distribution is then rotated around the azimuths
$\phi$-axis and therefore the $\phi$-distribution is uniform. The
Gaussian-distributed $\theta$ and the uniformly distributed $\ph$
define the probability for the new vector. This vector is inserted
into $\mathcal{F}(C)$ as before, and the Metropolis algorithm accepts or
rejects this particular choice. The quantity $q$ is given by
\begin{equation}
q(C_2\mid C_1)=\frac{P(C_2)\tilde{P}(C_1)}{P(C_1)\tilde{P}(C_2)}.
\label{qImpSamp}
\end{equation}
To implement this method, a Gaussian is chosen centered along the $z$-axis in
the intrinsic frame. This Gaussian is, in terms of the $z$-coordinate,
\begin{equation}
 \tilde{f}=\exp[-A(1-\cos\theta)^{2}]=\exp[-A(1-z)^{2}].
\label{tildef}
\end{equation}
The integral
\begin{equation}
 \int_{-1}^{1}N\exp[-A(1-z)^{2}]\,\d z
\end{equation}
satisfies (\ref{constrftilde}) where $N$ is a normalization constant
and $A$ is an arbitrary parameter that changes the breadth of the
Gaussian and thereby alters the acceptance rate. This is sampled using
the Box-M\"{u}ller method, by choosing
\begin{equation}
 z=1-\frac{1}{\sqrt{A}}\sqrt{-\ln\xi_{1}}|\cos(2\pi\xi_{2})|,
\end{equation}
which is a shifted Gaussian so that the peak is at
$z=\cos\theta=1$. All values of $z<-1$ are rejected. The azimuth angle
$\phi$ is sampled uniformly along $[0,2\pi)$ by $\phi=2\pi\xi_{3}$.
The acceptance probability becomes
\begin{equation}
A(C_2\mid C_1)=\min\left(1,\exp \left[-\beta(\mathcal{V}_{n}-\mathcal{V}_{p})
+A(1-z)^{2}\right]\right),
\end{equation}
using (\ref{qImpSamp}) and (\ref{tildef}).

\subsubsection*{Nuclear Skyrme model}

We do importance sampling by prioritizing small $v$ values and
selecting $u$ and $\phi$ with uniform probability. The small $v$
region corresponds to the small $\chi$ region. Importance sampling is
used because the newly selected four dimensional unit vector should be
in the neighborhood of the previous vector. This new vector is
selected in the intrinsic frame, as discussed in Sec.~\ref{MCBabyImportance},
where the previous vector is pointing in the
$\chi=0$ direction. A good probability distribution is
\begin{equation}
 \tilde{f}=\e^{-A v}.
\label{Sktildef}
\end{equation}
The quantity $v$ is therefore selected using
\begin{equation}
 v=-\frac{1}{A}\log\left(\left[\exp\left(-\frac{\pi}{2}A\right)-1\right]\xi_{1}+1\right),
\end{equation}
where $\xi$ is a uniform random variable on $(0,1)$. The other angles
are sampled as in Eq.~(\ref{3dsampling}).
The acceptance probability now becomes
\begin{equation}
A(C_2\mid C_1)=\min\left(1,\exp\left[-\beta(\mathcal{V}_{n}-\mathcal{V}_{p})+A v\right] \right).
\end{equation}
A similar probability distribution function can be used for the baby
Skyrme model and is faster than the given gaussian.  In practice
importance sampling is not used, because it is computationally more
time-consuming than restricting the transition probability.

\bibliographystyle{prsty}
\bibliography{references}
\end{document}

sa.bbl
\begin{thebibliography}{10}

\bibitem{russell}
J.~S. Russell, {\em Rep. 14th Meet. BAAS, York} (John Murray, London, 1844),
  p.\ 311.

\bibitem{soliton:intro}
P.~G. Drazin and R.~S. Johnson, {\em Solitons: an introduction} (Cambridge
  University Press, Cambridge, UK, 1996).

\bibitem{ryder}
L.~H. Ryder, {\em Quantum {F}ield {T}heory} (Cambridge University Press,
  Cambridge, UK, 1994), chapter 10.

\bibitem{qhe}
S.~L. Sondhi, A. Karlhede, S.~A. Kivelson, and E.~H. Rezayi, Phys. Rev. B {\bf
  47},  16419  (1993).

\bibitem{holzwarth:review}
G. Holzwarth and B. Schwesinger, Rep. Prog. Phys. {\bf 49},  825  (1986).

\bibitem{cosmic:string}
A. Vilenkin and E.~P.~S. Shellard, {\em Cosmic strings and other topological
  defects} (Cambridge University Press, Cambridge, UK, 1994).

\bibitem{duff}
M. Duff, hep-th/9805177 (unpublished).

\bibitem{sc:second}
O. Schwindt, internal UMIST report (unpublished).

\bibitem{tom}
T. Weidig, Ph.D. thesis, University of Durham, 1999.

\bibitem{ingber}
L. Ingber, J. Math. Comp. Model. {\bf 18},  29  (1993).

\bibitem{num}
W.~H. Press, S.~A. Teukolsky, W.~T. Vetterling, and B.~P. Flannery, {\em
  Numerical recipes in C} (Cambridge University Press, Cambridge, UK, 1992).

\bibitem{baby_pots}
T. Weidig, Nonlinearity {\bf 12},  1489  (1999).

\bibitem{me:equation}
N. Metropolis, A.~W. Rosenbluth, M.~N. Rosenbluth, and A.~H. Teller, J. Chem.
  Phys. {\bf 21},  1087  (1953).

\bibitem{ka:monte}
M.~H. Kalos and P.~A. Whitlock, {\em Monte {C}arlo {M}ethods} (John Wiley, New
  York, 1986).

\bibitem{sim:science}
S. Kirkpatrick, C.~D. Gelatt, and M.~P. Vecchi, Science {\bf 220},  671
  (1983).

\bibitem{geman}
S. Geman and D. Geman, IEEE Trans. Patt. Anal. Mac. Int. {\bf 26},  721
  (1984).

\bibitem{raj}
R. Rajarman, {\em Solitons and instantons} (North-Holland, Amsterdam, 1996).

\bibitem{coleman:sine}
S. Coleman, Phys. Rev. D {\bf 11},  2088  (1975).

\bibitem{sa:cont}
D. Vanderbilt and S.~G. Louie, J. Comput. Phys. {\bf 56},  259  (1984).

\bibitem{derrick}
G.~H. Derrick, J. Math. Phys. {\bf 5},  1252  (1964).

\bibitem{pz:93}
B.~M. A.~G. Piette and W.~J. Zakrzewski, Chaos, Solitons and Fractals {\bf 5},
  2495  (1995).

\bibitem{psz:95a}
B.~M. A.~G. Piette, B.~J. Schoers, and W.~J. Zakrzewski, Z. Phys. C {\bf 65},
  165  (1995).

\bibitem{os:finite}
O. Schwindt and N. Walet (unpublished).

\bibitem{ma:statistical}
F. Mandel, {\em Statistical physics} (John Wiley, New York, 1988).

\bibitem{bi:monte}
K. Binder, {\em Monte {C}arlo {M}ethods} (Springer-Verlag, Berlin, 1979).

\bibitem{fr:understanding}
D. Frenkel and B. Smith, {\em Understanding {M}olecular {S}imulation} (Academic
  Press, London, 1996).

\bibitem{thooft:2}
G. 't~Hooft, Nucl. Phys. B {\bf 72},  461  (1973).

\bibitem{witten}
E. Witten, Nucl. Phys. B {\bf 160},  55  (1979).

\bibitem{anw}
G. Adkins, C. Nappi, and E. Witten, Nucl. Phys. B {\bf 228},  552  (1983).

\bibitem{braaten}
E. Braaten, S. Townend, and L. Carson, Phys. Lett. B {\bf 235},  147  (1990).

\bibitem{multi:num}
R.~A. Battye and P.~M. Sutcliffe, Phys. Rev. Lett. {\bf 79},  363  (1997).

\end{thebibliography}
