# Simulation of the Zero Temperature Behavior of a 3-Dimensional Elastic Medium

## Abstract

We have performed numerical simulation of a 3-dimensional elastic medium, with scalar displacements, subject to quenched disorder. In the absence of topological defects this system is equivalent to a $`(3+1)`$-dimensional interface subject to a periodic pinning potential. We have applied an efficient combinatorial optimization algorithm to generate exact ground states for this interface representation. Our results indicate that this Bragg glass is characterized by power law divergences in the structure factor $`S(k)Ak^3`$. We have found numerically consistent values of the coefficient $`A`$ for two lattice discretizations of the medium, supporting universality for $`A`$ in the isotropic systems considered here. We also examine the response of the ground state to the change in boundary conditions that corresponds to introducing a single dislocation loop encircling the system. The rearrangement of the ground state caused by this change is equivalent to the domain wall of elastic deformations which span the dislocation loop. Our results indicate that these domain walls are highly convoluted, with a fractal dimension $`d_f=2.60(5)`$. We also discuss the implications of the domain wall energetics for the stability of the Bragg glass phase. Elastic excitations similar to these domain walls arise when the pinning potential is slightly perturbed. As in other disordered systems, perturbations of relative strength $`\delta `$ introduce a new length scale $`L^{}\delta ^{1/\zeta }`$ beyond which the perturbed ground state becomes uncorrelated with the reference (unperturbed) ground state. We have performed scaling analysis of the response of the ground state to the perturbations and obtain $`\zeta =0.385(40)`$. This value is consistent with the scaling relation $`\zeta =d_f/2\theta `$, where $`\theta `$ characterizes the scaling of the energy fluctuations of low energy excitations.

Observation of glassy behavior in flux line arrays in high $`T_c`$ superconductors  calls for a thorough theoretical description of such behavior. In this system the collective pinning of the flux line array, rather than the interactions of a single flux line with the disorder, can dominate the physics . For weak pinning, where dislocations are believed to be unimportant at large length scales , the entire flux line array can be modeled as a single medium subject to a pinning potential. Analytic calculations carried out using the approximation of linear elasticity and including the effects of the short range order in this system indicate that quasi-long range order exists in 3 dimensions . The elastic medium assumption was justified a posteriori and is further supported by an approximate domain wall renormalization calculation . The structure factor of a topologically ordered system was predicted by these calculations to have power law divergences of the form $`S(k)k^3`$. We will consider only the case of scalar displacements, which also models a charge density wave pinned by charge impurities .

We have numerically generated ground states for an elastic medium subject to quenched point disorder in the topologically ordered phase. Our results for the coefficient of the divergence of $`S(k)`$ lie between the renormalization group and Gaussian variational method results obtained by Giamarchi and Le Doussal . In addition to supporting their analysis, we are able to examine the response of the system to changes in boundary conditions and pinning potential. By a suitable choice of the boundary conditions we can simulate the domain wall of elastic deformations induced by a dislocation loop . The energy of the domain wall dominates the random part of the energy cost of introducing a single topological defect . Our results on the energetics of the domain walls thereby indirectly support the analysis carried out by Fisher , which indicated that the this system is marginally stable with respect to the introduction of dislocations. The numerically generated domain walls were found to have a fractal dimension $`d_f=2.60(5)`$. At large length scales the ground state is highly sensitive to small perturbations in the disorder potential, as in spin glasses and other disordered systems . Perturbations of relative strength $`\delta `$ in the disorder decorrelate the ground state on length scales $`L^{}\delta ^{1/\zeta }`$, with $`\zeta =0.385(40)`$. We are able to relate this response to disorder perturbations to the properties of the domain walls.

We have generated exact ground states for a discrete model whose energy in the continuum limit is given by

$$H=d^3x\frac{c}{2}[u(\stackrel{}{x})]^2+V(u(\stackrel{}{x}),\stackrel{}{x})$$
(1)   
with distortions of the medium represented by $`u(\stackrel{}{x})`$, which is assumed to be slowly varying over the system. The coefficient $`c`$ is the elastic constant. The potential felt by the medium due to the randomly placed impurities is represented by $`V(u(\stackrel{}{x}),\stackrel{}{x})`$. In microscopic descriptions of an elastic system subject to weak disorder there is a length scale $`\xi _p`$ below which the elastic energy dominates and the medium is ordered. This short range order manifests itself here as correlations in the disorder potential of the form $`V(u,\stackrel{}{x})=V(u+a,\stackrel{}{x})`$, with $`a`$ the intrinsic period of the medium. The period of the potential, $`a`$, is the lattice spacing in a flux line array or the wavelength of a charge density wave. Although this Hamiltonian is insufficient to describe the core of a dislocation loop, it can serve to describe the sheet of elastic deformations which span the loop, since the approximation of linear elasticity breaks down only in the region near the dislocation core. Comparing the ground state of Eq. (1) in a system of size $`L`$ subject to periodic boundary conditions to that for the same disorder realization with twisted boundary conditions in one direction (i.e., $`u(x,y,0)=u(x,y,L)+a`$) allows for the identification of the domain wall which would be caused by a single dislocation loop encircling the system .

The elastic Hamiltonian, Eq. (1), describes a $`(3+1)`$-dimensional interface subject to a disorder potential. In this picture, the displacement variable $`u(\stackrel{}{x})`$ maps to the height of the directed interface. This interface model has a natural discrete representation in which the configuration of the interface is specified by the set of bonds it cuts in a 4-dimensional lattice . The bonds of this lattice are assigned weights which directly correspond to the disorder potential. The sum of the weights of the bonds which the interface cuts give the energy of the configuration corresponding to the disorder energy $`d^3xV(u,\stackrel{}{x})`$. In these discrete models an effective elastic constant arises from a dependence of the number of configurations on the average gradient of the interface. Maximal flow algorithms , a subclass of combinatorial optimization algorithms, allow for the generation of ground states of this discrete representation of the interface .

The lattices numerically studied here are composed of $`L^3\times U`$ nodes, where $`L`$ is the linear size of the elastic medium, and $`U`$ is the extent of the lattice in the $`\widehat{u}`$ direction in which the displacement variable fluctuates. Unlike the simulation of elastic manifolds , where the bond weights are non-periodic, there are long range correlations in the disorder in the $`\widehat{u}`$ direction. We generated random integer weights, chosen from a uniform distribution over $`[0,V_{max}]`$, independently for each of the forward bonds in a layer of unit cells at constant height. The data we will present were obtained with $`V_{max}=5000`$, but we have verified that our results on the structure of the interface are not significantly altered for $`V_{max}`$ as low as $`100`$. Throughout the following discussion, we have normalized the energy of the system, so that the effective range of the bond weights is in $`[0,1]`$. This set of bond weights on one layer is sufficient to fix the value of the disorder on all of the bonds because we require that bonds which differ only by translations along $`\widehat{u}`$ have the same weight. In order to study the universality of the coefficient of the divergence in $`S(\stackrel{}{k})`$, we have simulated interfaces on the simple hypercubic lattice (SHC) as well as the Z-centered hypercubic lattice (ZHC) . For the SHC lattice $`\widehat{u}`$ was chosen to be along the (1111) crystallographic axis. The elementary bonds in this lattice are: $`(x,y,z,u)=`$ $`\pm (0,\sqrt{2}/2,1/2,1/2)`$, $`\pm (0,\sqrt{2}/2,1/2,1/2)`$, $`\pm (\sqrt{2}/2,0,1/2,1/2)`$, $`\pm (\sqrt{2}/2,0,1/2,1/2)`$. In the ZHC lattice we considered the following 12 bonds extending from each node: $`(x,y,z,u)=`$ $`(\pm \sqrt{2}/2,0,0,\pm \sqrt{2}/2)`$, $`(0,\pm \sqrt{2}/2,0,\pm \sqrt{2}/2)`$, $`(0,0,\pm \sqrt{2}/2,\pm \sqrt{2}/2)`$ . These lattices are the natural extensions of the two types of lattices used in simulations of the ground state of a 2-dimensional elastic medium . Both types of lattices were simulated using periodic boundary conditions in the transverse directions. In addition, in the ZHC lattice the ground state was computed for each realization of disorder with twisted boundary conditions in one direction.

We have developed a custom implementation of the push-relabel maximal flow algorithm  optimized for application to the regular lattices considered here. Our modifications to the Push-Relabel algorithm reduced the memory requirements by nearly a factor of ten, allowing for simulation of systems composed of up to approximately $`6\times 10^6`$ nodes using less than $`512`$MB. The primary modification involves computing nearest neighbor relations as needed rather than storing this information. Overhangs in the interface are precluded by assigning a large weight to backwards arcs . Since these backwards arcs have an effectively infinite weight, they cannot be be part of the minimal cut. Thus our algorithm can operate without storing their weight, provided that flow is always allowed to move along the backwards arcs. These modifications increased the running time of the algorithm, but obtaining the ground state for each realization of disorder still took less than one hour of processor time on a single 400 MHz Pentium II CPU for the largest system sizes studied. The memory requirement is linear in the number of nodes $`N`$; the processor time was found to scale approximately as $`N^{1.3}`$, compared with the worst case bound of $`N^2`$ .

The modest computational requirements of this algorithm have allowed us to average the properties of the ground states for a variety of system sizes over a large number of disorder realizations. In addition to generating the value of the minimal energy, the algorithm produces the configuration of the interface. The interface can be then represented by $`u(\stackrel{}{x})`$, which is defined on the 3-dimensional lattice formed by projecting the interface along the $`\widehat{u}`$ direction . Due to the periodicity of the disorder, the energy is invariant under global translations of integer multiples of $`a`$ in the displacement variable $`u(\stackrel{}{x})u(\stackrel{}{x})+na`$. Considering the set of the forward bonds cut at each location provides a characterization of the configuration equivalent to measuring the gradient of the interface. This representation of the interface is useful when comparing different ground states since it insensitive to global shifts of $`u`$. For the SHC lattice, we have generated at least $`10^3`$ realizations of disorder for systems of size $`L=8,16,24,32,40,48,60,80`$. We chose the extent of the lattice in the displacement direction $`\widehat{u}`$ to ensure that the boundaries of the system do not affect the ground state. For the SHC lattice $`U=20`$ was sufficient. Our simulations for the ZHC lattice were more extensive, with at least $`10^4`$ realizations for systems, subject to both periodic and twisted boundary conditions, of size $`L=8,16,32,48`$, and at least $`10^3`$ realizations for systems of size $`L=64,80`$. The largest systems here required $`U=12`$ to prevent the configuration from being affected by the boundaries of the lattice in the $`\widehat{u}`$ direction.

We have examined the displacement correlations of the minimal energy configurations by computing the disorder averaged structure factor $`\overline{S(\stackrel{}{k})}`$ of the displacement variables. This allows us to more clearly distinguish the large length scale behavior; direct measurement of the width is more difficult to analyze due to finite size effects. The orientationally averaged structure factor $`\overline{S(k)}`$ has been obtained by averaging the value of $`\overline{S(\stackrel{}{k})}`$ over radial bins of size $`\mathrm{\Delta }k=0.025`$, and is presented in Fig. 1. The error bars represent the fluctuations of $`\overline{S(\stackrel{}{k})}`$ within each spherical shell. These fluctuations, which measure the anisotropy of the structure factor, generally decrease with decreasing $`k`$; for $`k<0.75`$, these fluctuations saturate at a value comparable to the statistical fluctuations in $`S(\stackrel{}{k})`$ indicating the range in $`k`$ where the system is isotropic within the statistical fluctuations. To extract the coefficient of the leading order divergence we have fit $`k^3\overline{S(k)}`$ with the form

$$k^3\overline{S(k)}=A+Bk$$
(2)   
over the region $`k<0.5`$. The leading order term of $`\overline{S(k)}k^3`$ indicates the quasi-long range order of the ground state. For the ZHC lattice we obtain $`k^3\overline{S(k)}=1.08(5)+1.02(6)k`$; for the SHC lattice we obtain $`k^3\overline{S(k)}=1.01(4)+0.46(4)k`$. The error estimates on the parameters in these fits represent the statistical uncertainty over the given fit range and systematic errors arising from the choice of the cutoff in $`k`$. Giamarchi and Le Doussal have applied a renormalization group technique to order $`ϵ=4d`$ to this system and obtained $`A=1.0`$ . They have also carried out a Gaussian variational calculation and obtained $`A=1.1`$. In both of these approximations the value of the coefficient is universal.

The form $`Bk^2`$ of the leading order corrections to $`\overline{S(k)}`$ can be obtained by a renormalization group calculation. Upon renormalization, the periodic pinning potential $`V(u(\stackrel{}{x}),\stackrel{}{x})`$ introduces a new term into the Hamiltonian of the form $`\stackrel{}{\mu }(\stackrel{}{x})\text{.}u(\stackrel{}{x})`$ as when $`d=2`$ . This random tilting field is short range correlated with $`\overline{\stackrel{}{\mu }(\stackrel{}{x})\text{.}\stackrel{}{\mu }(\stackrel{}{y})}=g^2\delta ^3(\stackrel{}{x}\stackrel{}{y})`$. Unlike the case where $`d=2`$, the strength of this field $`g^2`$ undergoes only a finite renormalization for $`d=3`$ . The effects of this term can be determined from the effective small length scale Hamiltonian

$$H_{eff}=d^3x\frac{c}{2}[u(\stackrel{}{x})]^2+\stackrel{}{\mu }(\stackrel{}{x})\text{.}u(\stackrel{}{x}),$$
(3)   
which ignores the periodic pinning potential. Solving the equations of motion and averaging over realizations of the tilting field predicts corrections of the form $`Bk^2`$ with $`B=g^2/c^2`$ at $`T=0`$. In order to consider the effects of the pinning potential and the tilting field separately we relied on the fact that the renormalization group flow of $`V(u(\stackrel{}{x}),\stackrel{}{x})`$ is unaffected by the presence of a non-zero $`g`$ . The results of this analysis can be confirmed by examining the stability of the functional renormalization group fixed point obtained by Giamarchi and Le Doussal in an $`ϵ=4d`$ expansion . The most slowly decaying perturbation to the fixed point decays $`L^ϵ`$, implying corrections to the structure factor of the form $`k^{ϵ3}`$. This order $`ϵ`$ calculation also predicts the form of the corrections which are observed in our simulation data.

Naturally, real-space measurements of the width must be consistent with the structure factor. We measured the disorder averaged squared width

$$w^2=\overline{<u^2><u>^2},$$
(4)   
where $`<>`$ denotes a spatial average over the system. These data, shown in Fig. 2, have been fit using the real space version of Eq. (2): $`w^2=a+bln(L)+c/L`$. The constant term arises from the short wavelength fluctuations, while the second and third terms arise from the $`k^3`$ and $`k^2`$ terms in the structure factor respectively. The real space coefficient $`b`$ is related to the leading order behavior of $`S(k)`$ by $`b=A/4\pi ^2`$. This three parameter fit gives estimates for the values of $`A`$ and $`B`$, the coefficients describing the long wavelength form of the structure factor, which are consistent with that obtained by fitting $`S(k)`$ at a single system size. If the form of $`S(k)`$ depended on $`L`$, then this consistency would not be maintained. Direct comparison of the structure factor for various system sizes also demonstrates that $`\overline{S(k)}`$ exhibits negligible system size effects, except for the change in $`k_{min}=L/2\pi `$.

Other measures of the displacements provide us with additional information on the structure of the ground state. The disorder averaged extremal displacement difference, $`\mathrm{\Delta }H=\overline{u_{max}u_{min}}`$, Fig. 3, was found to grow logarithmically with system size for both lattice types. We computed least squares fits of the form $`\mathrm{\Delta }H=\stackrel{~}{a}+\stackrel{~}{b}ln(L)`$ to obtain $`\stackrel{~}{b}_{ZHC}=0.76(1)`$ and $`\stackrel{~}{b}_{SHC}=0.70(1)`$. The coefficients of the logarithmic term differ by less than $`10`$% for the two lattices studied here, suggesting that this measure of the system is weakly, if at all, dependent on the lattice discretization of the medium. This logarithmic growth is consistent with the following picture of the ground state structure developed by Fisher . At each length scale $`R=b,b^2,b^3,\mathrm{}`$ the displacement undergoes one shift of amount $`\pm a`$. Furthermore the sign of the displacement shift is random at each scale, leading to the logarithmic growth of the squared width. When traversing from the minimum to the maximum the signs of the displacements are strongly correlated, leading to a coherent sum, and a logarithmic dependence on the system size for the extremal differences results.

We have also determined the effect of coarse-graining the displacement variable. The coarse grained displacement is defined as the average of $`u`$,

$$u_R(\stackrel{}{y})=\frac{1}{R^3}_{\mathrm{\Omega }_R(\stackrel{}{y})}d^3xu(\stackrel{}{x}),$$
(5)   
over $`\mathrm{\Omega }_R(\stackrel{}{y})`$, a cube of size R centered at the point $`\stackrel{}{y}`$. We measured the fluctuations in these coarse grained height variables $`|\mathrm{\Delta }u_R|^2=\overline{(u_R(\stackrel{}{y})u_R(\stackrel{}{y}+\stackrel{}{b}))^2}`$, with $`\stackrel{}{b}`$ the vector between the centers of the cubes which touch at one corner. This spatial averaging procedure is similar to a real space renormalization transformation. Villain and Fernandez have explicitly carried out a real space renormalization calculation for a 3-dimensional elastic medium with cubic symmetry . Their calculations indicate that $`|\mathrm{\Delta }u_R|^2`$ has a finite limit as $`R,L\mathrm{}`$. We have directly measured $`|\mathrm{\Delta }u_R|^2`$ for the ZHC lattice (Fig. 4). The coarse grained height fluctuations are related to the structure factor by:

$$|\mathrm{\Delta }u_R|^2=_{BZ}\frac{d^3k}{(2\pi )^3}|G(\stackrel{}{k})|^2e^{i\stackrel{}{k}\stackrel{}{b}}S(\stackrel{}{k})$$
(6)   
with $`G(\stackrel{}{k})=_{\mathrm{\Omega }_R}d^3xe^{i\stackrel{}{k}\stackrel{}{x}}`$. In the infinite volume limit $`|\mathrm{\Delta }u_R|^2`$ depends only on the leading order behavior, $`Ak^3`$, of the structure factor and the limit of the ratio $`R/L`$. We have numerically evaluated the right hand side of Eq. (6) in this limit to obtain the infinite size limit presented in Fig. 4. For finite sized systems the sub-leading order corrections to $`S(k)`$ contribute to the coarse grained height fluctuations. These terms lead to the divergence of $`|\mathrm{\Delta }u_R|^2`$ as $`R/L0`$. The dominant corrections arise from the $`Bk^2`$ corrections to $`S(k)`$ seen in the structure factor data, but decay as $`L^1`$. The data are consistent with convergence to a finite limit as $`L\mathrm{}`$.

For the ZHC lattice, we have also investigated the behavior of the system subject to the twisted boundary conditions defined previously. For each realization of disorder we have compared the ground state energy with periodic boundary conditions, $`E_p`$, to that obtained with twisted boundary conditions along one of the lattice directions, $`E_t`$. This allows us to investigate the properties of the excitations induced by the change in boundary conditions. We identify the energy difference $`E_{DW}=E_tE_p`$ with the energy of the domain wall. The domain wall is identified by the set of bonds that the interface intersects in one set of boundary conditions but not the other. Even though the domain wall could be identified by examining the values of $`u_p(\stackrel{}{x})`$ and $`u_t(\stackrel{}{x})`$, it is more efficient to identify the domain wall by examining the sets of cut bonds for each boundary conditions. The cut bonds are projected along the $`\widehat{u}`$ direction for each boundary condition. The domain wall is then the symmetric difference between these two sets of bonds. Without directly simulating a dislocation loop itself, we are able to investigate the properties of the domain wall induced by the introduction of a single loop encircling the system.

The energetics of these domain walls dominate the random part of the energy cost of introducing a dislocation loop into an elastic medium . The mean energy difference $`\overline{E}_{DW}`$ grows linearly with the (linear) size $`L`$ of the domain wall (Fig. 5), a result consistent with the scaling of the elastic contribution to the energy and the statistical symmetry of the disorder potential of the continuum model. We have also analyzed the variance $`\sigma ^2(E_{DW})`$ of the distribution of domain wall energies (Fig. 6). No single power law fit of the form $`\sigma ^2(E_{DW})L^{2\theta }`$ can adequately fit our data over the range of sizes simulated. However the data are well described by the empirically determined form $`\sigma ^2(E_{DW})=0.031L^2+0.24L`$ displayed as the solid line in Fig. 6. The finite size correction leads to a size dependent effective value of the exponent, $`\theta _{eff}`$, characterizing the scaling of the sample to sample fluctuations in the energy of low energy excitations. Depending on the lower limit imposed on the fit, single power-law fits $`\sigma ^2(E_{DW})L^{2\theta }`$ give $`0.85<\theta _{eff}<0.92`$ over this range of sizes. The distribution of domain wall energies depends on $`\overline{E}_{DW}`$ and $`\sigma (E_{DW})`$ only through the combination $`ϵ=(E_{DW}\overline{E}_{DW})/\sigma (E_{DW})`$, as can be seen in Fig. 7. This collapsed distribution has more highly weighted tails than a unit normal distribution. The frequency with which negative energy domain walls were observed in the simulations is significantly higher that what one expects for a Gaussian distribution with the measured mean and standard deviation at each system size (Fig. 8). This behavior also occurs for high energy domain walls; to within the statistical uncertainty the distribution is symmetric about its mean value. Our data are consistent with both $`\overline{E}_{DW}`$, and $`\sigma (E_{DW})`$ increasing linearly with $`L`$ for large systems. Fisher’s argument assumed this behavior of the domain wall energetics in his domain wall renormalization calculation indicating the marginal stability of the Bragg glass phase .

Because of the balance between the elastic energy scale and the scale of the energy fluctuations due to the disorder potential, the domain walls are highly convoluted and expected to have a fractal dimension $`d_f`$ between 2 and 3. Similar to the approach used for a 2-dimensional elastic medium , we have measured the size of the domain wall by counting the number of bonds $`N_b`$ in the wall. The data for the area of the wall as a function of system size, averaged over disorder, can be fit by a simple power law, $`N_bL^{d_f}`$, with the fractal dimension of the domain wall $`d_f2.60`$, shown in Fig. 9(a). This fit has been taken after excluding the smallest system size $`L=8`$. However, the form of the residuals, see Fig. 9(b), indicate the presence of sub-leading order corrections, suggesting that this may underestimate the value of $`d_f`$. In order to verify this, we have also fit the whole range of data using the form $`N_b=aL^{d_f}+bL^2`$. The second term arises from the effectively two dimensional nature of the domain walls at small length scales. This three parameter fit gives $`d_f=2.65`$; thus we conclude that our systematic errors in estimating the fractal dimension are approximately $`0.05`$.

In order to investigate other low energy excitations of this system, we have examined the sensitivity of the ground state to perturbations in the disorder potential. Similar studies have been carried out for disordered systems such as spin glasses , $`(1+1)`$-dimensional directed polymers in random media , and 2-dimensional elastic media . In all of these disorder dominated systems, perturbations of relative strength $`\delta `$ in the disorder potential introduce a length scale $`L^{}\delta ^{1/\zeta }`$, $`\zeta =d_f/2\theta `$, beyond which the ground state becomes uncorrelated with the reference ground state. The exponent $`\zeta `$ characterizing the sensitivity of the ground state is referred to as the chaos exponent. These studies have been done by comparing the ground state for two correlated choices of the disorder potential. In our simulations, we have obtained the ground states for the two pinning potentials $`V^\pm (u(\stackrel{}{x}),\stackrel{}{x})=b(u(\stackrel{}{x}),\stackrel{}{x})\pm d(u(\stackrel{}{x}),\stackrel{}{x})`$, with both terms periodic in the $`\widehat{u}`$ direction. The constant part of the potential $`b(u(\stackrel{}{x}),\stackrel{}{x})`$ was an integer chosen uniformly from $`[1000,2000]`$. The term which generates differences between the realizations, $`d(u(\stackrel{}{x}),\stackrel{}{x})`$, was chosen uniformly from $`[d_{max}/2,d_{max}/2]`$. The parameter $`\delta =d_{max}/2000`$ characterizes the relative strength of the perturbations. This prescription was chosen to ensure that for a fixed value of $`\delta `$ the distribution of the bond weights is the same for both realizations of disorder. Our simulations include values of $`\delta `$ ranging from 0.01 to 0.75, with at least $`500`$ independent realizations of disorder at each $`\delta `$ and $`L`$. By performing scaling analysis of both the energetic and structural correlations between the ground states for these realizations of disorder we can extract the value of the chaos exponent.

We have found that both the energetic and structural correlations are governed by the same length scale. First we calculated the domain wall energy $`E_{DW}^\pm `$ for the two disorder realizations $`V^\pm `$, which can then be used to compute the domain wall energy correlation function

$$G=\frac{(\overline{E_{DW}^+\overline{E_{DW}^+}})\overline{(E_{DW}^{}\overline{E_{DW}^{}})}}{\sigma (E_{DW}^+)\sigma (E_{DW}^{})}.$$
(7)   
The simple scaling form $`G=f(\delta L^\zeta )`$ describes our data well (Fig. 10). We found reasonable data collapse for $`\zeta =0.38(4)`$, taking into account the statistical errors. The value of the chaos exponent can be related to the domain wall fractal dimension and the energy fluctuation exponent by a simple scaling argument as in the case of spin glasses . The perturbations introduce a random change in the energy of the domain wall of order $`\delta L^{d_f/2}`$ because the perturbations are uncorrelated with the location of the domain wall. The typical fluctuations in the domain wall energy scale as $`L^\theta `$. When these energy scales become comparable, at a length scale $`L^{}\delta ^{1/\zeta }`$, the domain wall energies become uncorrelated. When using the effective value of the energy fluctuation exponent, $`\theta _{eff}0.9`$, this scaling relation holds to within 5% accuracy.

We can understand the structural deformations induced by the bond perturbations by reasoning similar to that for domain wall correlations. Here we consider the differences in the ground states of the system (with periodic boundary conditions) due to the changes in the pinning potential. Again, the ratio of the energy change due to the random perturbations, and that of the fluctuations in the energy landscape, $`\delta L^{d_f/2}/L^\theta `$, determines the behavior of the system. In response to Zhang’s simulation of directed polymers in random media, Feigel’man and Vinokur had argued that the probability of a positional excitation grows linearly with $`\delta L^\zeta `$ for small values of the perturbation strength . Following their argument, we expect the probability of a change in the displacement variable at a single location to grow linearly with $`\delta L^\zeta `$. Unlike the case for the directed polymer, the magnitude of the differences is bounded for periodic pinning since excitations with $`|\mathrm{\Delta }u|>a`$ probe the same energy landscape as those with $`|\mathrm{\Delta }u|<a`$. Thus, the spatially and disorder averaged mean squared displacement difference $`\chi =\overline{<(u^+(\stackrel{}{x})u^{}(\stackrel{}{x}))^2>}\delta L^\zeta `$ for small perturbations. For large values of $`\delta `$ the ground states are completely decorrelated, and we expect that $`\chi ln(L)`$. We have found that the data for $`\chi `$ collapse according to the scaling form $`\chi =f(\delta L^\zeta )`$, with $`\zeta =0.39(2)`$ (Fig. 11). Our data are consistent with the results obtained by the scaling argument in both limiting cases. Before computing $`\chi `$, we made the transformation $`u^+u^++na`$, where $`n`$ is the integer which maximizes the number of locations at which $`u^+(\stackrel{}{x})u^{}(\stackrel{}{x})=0`$ for each realization. This transformation minimizes $`\chi `$ over the discrete set of global translations which leave the energy of the medium invariant. Our scaling ansatz is significantly different from that proposed for the 2-dimensional elastic medium $`\chi =\delta ^{1/\zeta }f(L)`$ . This form cannot adequately collapse our data over the range of parameters we have simulated.

Implicit in this discussion is the assumption that both the perturbation induced deformations, and the boundary condition induced domain walls are characterized by the same fractal dimension. Even for relatively small perturbations in the disorder, the deformations are typically composed of a set of disconnected clusters. We have directly measured these clusters’ fractal dimension. The size of a cluster $`R`$ is defined as the average of the sides of the bounding box which encloses the cluster and is measured in units of the lattice spacing. Our algorithm identifies the sets of nodes on which $`u^+(\stackrel{}{x})u^{}(\stackrel{}{x})0`$ after performing the translation which minimizes $`\chi `$. The surface area of a cluster $`s`$ is the number of nodes with neighbors not in the cluster. We have collapsed the data for $`\delta =0.05`$ using the finite size scaling form $`s=L^{d_f}f(R/L)`$ (Fig. 12). We expect that the scaling function $`f(R/L)`$ should have the form $`f(R/L)^{d_f}`$ in the region $`R/L<<1`$, $`R>>1`$, but this regime is not clearly visible in Fig. 12 due to lattice and finite size effects. Despite this, the best collapse of the data, for $`0.5<R/L<1.0`$, provides an estimate $`d_f=2.65(10)`$ which is consistent with the estimate for the domain wall fractal dimension. Similar analysis at other values of $`\delta `$ provide equivalent values for the cluster fractal dimension. The anomalous data at $`R/L1`$ arise from the rare clusters which span the system in all directions. This scaling also breaks down for clusters with $`R<8`$, where lattice effects make the surface effectively 2-dimensional (Fig. 13). The equality of the fractal dimension of the boundary condition induced domain walls and bond perturbation induced deformations can be justified by a simple argument. For small values of the disorder perturbation parameter, the cluster boundaries lie in regions where there is a small energy cost to deforming the medium. If one considers only a small volume containing a portion of the cluster boundary, the structural difference is the same as would be caused by the change from periodic to twisted boundary conditions on that volume. Thus both the deformations induced by small changes in the disorder potential and those caused by a change in boundary conditions should be characterized by the same fractal dimension. Despite the fact that the scaling regime is inaccessible due to the limits on the size of systems studied, our data are consistent with the conclusion that the fractal dimension of the clusters which compose the deformations is the same as that of the boundary condition induced domain walls.

We have performed extensive numerical simulation of a model 3-dimensional elastic medium subject to quenched disorder with scalar discrete displacements. Our results for the structure in the Bragg glass phase indicate that the structure factor has divergences of the form $`S(k)Ak^3`$. Our results for the coefficient A fall between the approximate values $`A=1.0`$, and $`A=1.1`$, obtained via a renormalization group and a replica approach . The observed energetics of the boundary condition induced domain walls indirectly support arguments for the stability of the Bragg glass phase. These domain walls correspond to the elastic deformations due to the introduction of a single dislocation loop winding around the system. Our data are consistent with the hypothesis that the mean energy and the energy fluctuations of a section of domain wall both scale linearly with the linear size of the section for large sizes. This balance is a crucial element of the analysis carried out by Fisher indicating the marginal stability of the Bragg glass phase to the introduction of dislocations . We are also able to measure the spatial structure of these domain walls and obtain their fractal dimension $`d_f=2.60(5)`$. We have observed that random changes in the disorder potential of relative strength $`\delta `$ decorrelate the ground state on length scales larger than $`L^{}\delta ^{1/\zeta }`$ with $`\zeta =0.385(40)`$. The properties of the domain walls and this sensitivity to disorder perturbations can be related to each other by the scaling relation $`\zeta =d_f/2\theta `$, where $`\theta `$ characterizes the fluctuations in the low energy excitations.