text:
# Noise in neurons is message-dependent (neuronal reliability/temporal coding/information theory/error-correcting)

## I Introduction

Brains represent signals by sequences of identical action potentials or spikes . Upon presentation of a stimulus, a given neuron will produce a certain pattern of spikes; upon repeating the stimulus, the pattern may repeat spike for spike in a highly reliable fashion, or may be similar only ”rate-wise”, or some spikes may be repeatable and others not. If individual spikes can be recognized and tracked across trials, then their timing accuracy can be ascertained unambiguously; this is not always the case. The existence of repeatable spike patterns and the reliability of their timing changes not only from neuron to neuron, but even for the same neuron under varying circumstances. Bryant and Segundo  first noticed (in the mollusk Aplysia) that spike timing accuracy depends on the particulars of the input driving the neuron. This intriguing property has received renewed attention, notably among those searching for experimental evidence supporting different theories of neural coding . In a study of the response of pyramidal neurons to injected current  the temporal pattern of firing was found to be unreliable when the injected current was constant, but highly reliable when the input current was “noisy” and contained high frequency components. This study showed explicitly the difference between the irregularity of the spike pattern (which was irregular in both cases) as opposed to the reliability or accuracy of spike timing, and also highlighted the connection that “natural” stimuli are noisy and contain sharp transitions. Similar results have been obtained from in vivo recordings of the H1 neuron in the visual system of the fly : constant stimuli produced unreliable or irreproducible spike firing patterns, but noisy input signals deemed to be more “natural” and representative of a fly in flight yielded much more reproducible spike patterns, in terms of both timing and counting precision. Although some aspects of the last study have been challenged recently in , the phenomenon of different reliabilities for varying inputs was reaffirmed: see e.g. Fig. 1 in . Thus one concludes that similar observations of the timing precision of spikes can be made in very different types of neurons under vastly different conditions, and that this must be a very universal, almost basic characteristic of neuronal function. We shall now show how these experimental observations follow from rather general and elementary considerations, and discuss the possible implications for brain function and architecture.

## II Stochastic spiking models

The essence of this phenomenon can be demonstrated easily in toy models. A simple example is the leaky integrate-and-fire (I-F) model , which assumes the neuron is a (leaky) capacitor driven by a current which simulates the actual synaptic inputs. We add membrane fluctuations representing several internal sources of noise (cluttering of ion channels, synaptic noise, etc.) to obtain a system described by the following Langevin  equation:

$$C\dot{V}=gV+I(t)+\xi (t)$$
(1)
where $`C`$ is the cell membrane capacitance, $`V`$ the membrane potential, $`gV`$ the leakage term ($`g`$ is a conductance), $`I(t)`$ is the input current, and $`\xi (t)`$ is Gaussian noise, $`\xi (t)=0`$ with autocorrelation $`\xi (t)\xi (t^{})=\sigma \delta (tt^{})`$; when the potential reaches the threshold $`V_0`$ an action potential is generated and the system returns to the equilibrium potential $`V_e`$, here set arbitrarily to zero.

Fig.1 shows the results of a numerical simulation of Eq. (1) in response to two different signals , in both cases in the presence of identical “internal” noise. Following refs. , the left panels show a time independent input and the neuronal responses to repeated stimulation, while the right panels ensue from a “noisy” input and its responses. When $`I(t)`$ contains high frequency components (right side panels) spikes are clustered more tightly to the upward strokes of the input and overall there is less trial-to-trial variability than in the time-independent case. The phenomenology described in Fig. 1 captures the main feature that renewed interest in this subject; i.e., that apparently more “noisy” (in the sense of “severely fluctuating”) inputs can generate less variable responses.

By concentrating on these two extreme stimuli, these studies have somewhat obscured an issue which is central to the variability phenomenon: the relationship of the output to the timescale or timescales over which the input itself varies. The constant stimulus can be thought to have no timescale, or at best, just a timescale since onset; while in principle the noise process has indefinitely fast timescales. Thus the two extreme cases under discussion have timescales outside the range of interaction with physiological processes.

A simple dimensional argument shows that timing variability is dependent on the characteristic timescale of the input. As a consequence of noise, the neuron’s membrane potential will have fluctuations of a characteristic size (or root mean square) $`\mathrm{\Delta }V`$. In order to translate this $`\mathrm{\Delta }V`$ into a timing noise $`\mathrm{\Delta }t`$, we have to scale the potential noise by something having units of $`\mathrm{\Delta }t/\mathrm{\Delta }V`$; the most natural such quantity at hand is the inverse of the potential’s time derivative.

While the distribution of the derivative of the potential is a reasonable measure of the degree of reliability of a particular input, one can see from a geometrical argument that the derivative of the potential at the particular time in which the potential crosses the threshold is a better estimator of stimulus reliability. The uncertainty in firing time is essentially the time during which the membrane potential is at a distance from the threshold smaller than the value of the voltage noise, and so is scaled by the derivative at this particular moment. Thus, the faster the voltage approaches the threshold (in the vecinity of the threshold), the more reliable the timing of the spike will be. The statement captures much of the phenomenology in , and in our own data, as we will show.

We proceed now to apply this method of estimation to the situation shown in Fig. 1. A rapidly fluctuating signal was constructed  and the times at which spikes occurred were determined first for the deterministic condition (i.e., $`\sigma =0`$), and subsequently for several hundreds stochastic realizations (i.e., with identical $`I(t)`$ but with different stochastic realizations of the noise term in Eq. 1). At the same time the voltage derivative preceding each spike was computed over the prior fifty time steps. The fluctuating character of $`I(t)`$ (as illustrated in the top right panel of Fig. 1) offers the opportunity to explore a wide range of d$`V`$/dt at threshold crossings. To estimate the precision of spike timing we define an index of “temporal jitter” as follow: it is the absolute difference between the time at which a spike is generated in the noise-free simulation and that of the corresponding nearest spike for the stochastic realization, averaged over realizations, and over all spikes which occurr within a d$`V`$/dt interval of $`0.1`$. This quantity is plotted in Fig. 2. The data points in Fig. 2 are scattered around the expected (dotted lines) inverse relationship between the temporal jitter of the spike and the speed at which the voltage crossed the threshold. The trajectories plotted in both insets help to visualize the geometrical argument already discussed: inputs that rise more rapidly to the membrane potential threshold produce less variable spike timing.

The same basic phenomenon will affect any model of spike generation, in which the noise in the voltage-like variable is translated into jitter of its threshold crossings. To provide for a specific example, we choose a widely used model of excitable media with continuous dynamics, the FitzHugh-Nagumo model (FHN) with additive stochastic forcing , described by the following Langevin system:

$`\dot{V}`$ $`=`$ $`VV^3/3W+I(t)+\xi (t)`$ (2)
$`\dot{W}`$ $`=`$ $`\varphi (V+abW)`$ (3)
The variable $`V`$ is the fast voltage-like variable, $`W`$ is the slow recovery variable; $`a`$, $`b`$ and $`\varphi `$ are constant parameters, and $`\xi (t)`$ is a zero-mean Gaussian white noise of intensity $`D`$. According to these equations, after a spike the recovery produces an absolute refractory period during which a second spike can not occur, followed by a longer relative refractory period during which firing requires stronger perturbations.

As previously argued, it is biologically reasonable to assume that noise affects the voltage-like variable. Thus, according to our argument, a higher steady input will result in a faster approach to the threshold crossing and therefore will reduce the probability of an untimely crossing. This is confirmed in Fig. 3, where the distributions of periods are plotted for two inputs amplitudes. Note the large change in variance (1.5 to 2.34, 56$`\%`$) for comparably small change in the average period (36.4 to 39.1, 7$`\%`$) in going from high to low input, also note that the low input distribution has a supra-exponential tail. For comparison, the distribution of inter-spike intervals corresponding to the I-F model with constant input is shown. In both cases, the distribution for high input tends to a Gaussian (parabola on semi-log plot), while the low input has an exponential tail (linear on semi-log plot) as expected for Poisson-like statistics. It is interesting to notice that even though this phenomenon will affect more strongly systems with zero-frequency quiescent-to-firing bifurcations, given the longer time to integrate fluctuations at low input (i.e. the system can evolve near a homoclinic orbit), the example presented here shows that it also affects non-zero-frequency transitions like the Hopf bifurcation present in the FHN model.

## III Experimental results

Thus, in neural data, one should expect to find that the output noise generically depends on the structure of the input. Experiments discussed next demonstrate the relevance of this observation for brain function. We recall that an assumption usually made in the application of information theory to neurons is that the noise introduced by the communication channel is independent of the message being transmitted, and now explore to what extent this assumption is valid in real neural data. The result is presented in Fig. 4, where different visual stimuli (i.e., the four messages) were presented to an anesthetized cat while electrophysiological recordings were made in the lateral geniculate nucleus (LGN), which is the second stage of processing in the visual pathway (see  for methods).

The raster plot in Panel A shows the response to a moving bar with high contrast and high speed, and Panel B shows the responses to a low contrast bar moving at low speed, both results are shown in a window of 200 msec centered at the peak of the ON response. The responses in Panel A clearly display a higher temporal precision than those for the condition in Panel B. Two further conditions were also presented, for high contrast and low speed, and for low contrast and high speed. In the inset of Fig. 4 we quantify the noise for the four messages as the “jitter” of the response onset, which we define as the average of the absolute value of the position with respect to the mean of the first spike in each trial. The messages include two conditions for the velocity and two for the contrast. For each velocity condition, the high contrast is more reliable than the low one: 8.4 msec vs. 22.1 msec, and 16.8 msec vs. 23.8 msec. Note however that the high contrast/high speed and low contrast/low speed stimuli, although they have a very similar spike response (4.17 vs. 3.86) have the largest difference in jitter (8.4 msec vs. 23.8 msec). Alternatively, we can measure the entropy difference between the corresponding first spike probability distributions (1 msec bins), which varies from 2.0 bits for low contrast/high speed vs. high contrast/high speed, to 0.6 bits for high contrast/low speed vs. high contrast/high speed. Similar results showing contrast-dependent temporal precision in LGN have been reported previously in , where it is also proposed that the variability is due to intrinsic noise. We also show that messages with indistinguishable mean rate responses can show large differences in their variability, in agreement with the hypothesis of the threshold crossing derivative as the relevant parameter, which can be independent of the firing rate for particular inputs.

## IV Network architecture and reliability

The observations discussed in the previous section are restricted to the behavior of individual neurons. How can they be extended to networks of spiking elements? In principle, we can assume that input-dependent noise will be present in early stages of processing, given the relative simplicity of the neural architecture. The relevance of this phenomenon in higher order areas remains to be explored; nevertheless, we may hypothesize that unless dedicated architectures are implemented to eliminate it, significant input-dependent noise will be ubiquitous in neural processing.

Indeed, a recurring question in brain theory and computer theory is whether there are reliable ways to compute which use unreliable components. We shall now show explicitly the existence of a network topology for which the overall time reliability can be made arbitrarily better than that found in its individual neurons.

Assume a single input source $`I(t)`$, which “fans out” to $`N`$ noisy neurons $`X_i`$; these in turn “fan in”, or input into a single (also noisy) neuron $`Y`$. Let’s imagine an $`I(t)`$ constructed so that the $`X_i`$ fire once each at times $`t_i`$. If the size of the connections $`XY`$ is set so that $`Y`$ fires after seeing half of the expected number of spikes coming from the $`X`$s, then clearly neuron $`Y`$ will be fire near the median time of the spike times $`t_i`$. The median, as a descriptor of the $`t_i`$, has a number of exceedingly nice properties in this context: (a) its variance decreases as $`1/\sqrt{N}`$, and thus the timing accuracy of $`Y`$ can be made better than the individual $`X`$s simply by choosing an appropriately large $`N`$; (b) the median, unlike the mean, is exceedingly robust against outliars and heavy-tailed noise; (c) the median does not require that it see the entire set $`t_i`$ but only half of it, and (d) the median is expected to lie at the time of highest concentration of the $`t_i`$, thus maximizing the timing accuracy of $`Y`$. This can be fully apreciated in Fig. 5. Thus, we see that different topologies will propagate spike timing accuracy in different fashions, and that one probably should expect architectural correlates in place in the brain.

To gain further insight into this issue and to study the implications of our observations on the synchronization of a neuronal network, we now concentrate on random but fixed time delays across pathways, rather than on individual timing accuracy. We modify the previous model: all neurons of the middle layer are now noise-free and in principle identical, and so are their spike trains; the neurons are integrate-and-fire units as in Eq. 1, with the further addition of a refractory period. We next model a dispersion of the transmission delay times by rigidly shifting the spike train of each individual neuron by a random value uniformly distributed between $`0`$ and $`\delta `$. The total charge that the integrating neuron receives does not depend on this shift, but when all neurons are synchronized, the derivative of the charge will be high and then we expect more reliable responses. Fig. 6 shows how the timing reliability (open circles), measured by the method of spike metrics , decreases with the decorrelation of the input layer. Interestingly, when the derivatives are high the response is more reliable, but by the same token, the current supplied during the refractory period will be lost. This means that smooth distributions will result in an increased effective charge and thus the rate of response of the integrating neuron is bigger when the input layer is decorrelated, as shown in Fig. 6 (open squares) ***There is a second regime (not shown), if the decorrelation is large compared to the decay time of the neuron, in which rate decreases with increasing $`\delta `$ due to current loses..

We think that this very simple concept might be of physiological relevance, in that it implies a trade-off between temporal reliability and the dynamic range of the signal: when the dynamic range is large, it is possible to afford the cost of losing part of the charge to achieve reliability, while in a situation in which the dynamic range is limited, reliability has to be sacrificed in order to preserve and propagate the signal. This suggests the possibility of integration pathways which multiplex the signal, as found in retinal adaptation.

## V Conclusion

The relevance of the message-dependent nature of noise in spiking elements is, to our understanding, two-fold. One aspect is its consequences for the use of Shannon’s Information Theory as a framework to measure the information content present in the output of a neuron. This is the subject of much work recently . In this regard, a basic assumption commonly made is that the noise introduced by the communication channel is independent of the message being transmitted, which allows the modeling of a neuron as a Gaussian channel . This is indeed desirable given that departures from the Gaussian channel have proved to be rather cumbersome; for instance, no general theory for the computation of channel capacity exists . Our results render average measures of information per spike less meaningful than usually thought, and speak for the necessity of concentrating on particular individual spikes. This is specially critical when the message consists of only a few spikes (for instance, in the auditory pathway), or in the event of fast conductance modulations , which can dramatically affect the reliability. Thus, perhaps the most fundamental consequence of Mainen $`\&`$ Sejnowski’s result  is the implicit demonstration that cortical neurons are not always classical Gaussian channels.

A second relevant aspect is the possible anatomical correlate of this phenomenon. In this regard, we have shown that unless dedicated architecture is implemented to reduce the multiplication of noise along the processing pathways, the encoding of behaviorally relevant (not only “artificial” stimuli) can be highly degraded. Finally, we demonstrated that it is possible to design network topologies with arbitrarily large temporal reliability, and consequently one should expect an evolutionary pressure to implement them in specific areas of the brain where time accuracy is essential.

Supported in part by the Mathers Foundation (MM), NIH (JMA), Human Frontiers (LM), Rosita & Norman Winston Foundation (GAC) and Burroughs Wellcome (MS). Supercomputer resources are supported by NSF ARI Grants. Discussions with T. Gardner, S. Ribeiro and R. Crist are appreciated. We also acknowledge the valuable input from D. Reich and B. Knight.


format: text

text:
# Weak ferromagnetism of quasi-one-dimensional 𝑆 = 1/2 antiferromagnet BaCu2Ge2O7

## Abstract

Weak ferromagnetism of quasi-one-dimensional $`S`$ = 1/2 antiferromagnet BaCu2Ge2O7 is studied by the magnetization measurement. The spontaneous magnetization appears along the $`b`$ axis below 8.8 K. The local symmetry between the intra-chain nearest neighbor spins allows the presence of Dzyaloshinskii-Moriya interaction, and the only possible spatial configuration of the weak ferromagnetic moment per spin uniquely determines the sign of the inter-chain interaction. A weak $`a`$-axis magnetic field can change the direction of the magnetization to the $`a`$-axis direction, which shows that the spin chain forms a weakly coupled weak-ferromagnetic chain system.

One-dimensional (1D) quantum antiferromagnet (AF) is of particular interest in both theory and experiment of quantum magnetism. It is widely accepted that a pure 1D S=1/2 antiferromagnet has no long-range order at $`T`$ = 0 K.  On the other hand, almost all the actual spin-chain systems show magnetic three-dimensional long-range order (3D-LRO) at their ground states due to a finite inter-chain interaction. While in general there are many kinds of magnetic 3D-LRO states, the transition to a normal 3D AF-LRO state has been studied mostly.  One remarkable exception is a weak ferromagnetic transition that was recently found in 1D compound Sr0.73CuO2.  However, this compound cannot be treated as a uniform spin-chain system but rather is well explained as an alternating-chain system. In this paper, we report weak ferromagnetism (WF) of a uniform $`S`$ = 1/2 spin chain in BaCu2Ge2O7, which is isomorphous to antiferromagnet BaCu2Si2O7.  In BaCu2Ge2O7, intra-chain Dzyaloshinskii-Moriya (DM) interaction and inter-chain Heisenberg interaction cooperatively add a spontaneous magnetization in contrast to AF BaCu2Si2O7.

The single-crystal sample of BaCu2Ge2O7 was grown by a floating-zone method, and was cut into a rectangular shape ($`a`$ $`\times `$ $`b`$ $`\times `$ $`c`$ = 1.7 $`\times `$ 1.6 $`\times `$ 3.1 mm3). Since the $`a`$\- and the $`c`$-axis lengths are almost the same, particular attention was paid to determine the crystal orientation by x-ray diffraction. We took the diffraction at randomly selected three different points at every surface, and found no trace of a twinned structure. However, the possibility of inclusion of misoriented grains inside the crystal was not eliminated. Magnetization measurement was carried out with a commercial SQUID magnetometer (MPMS, Quantum Design) in the field range from 0 to 5 T. We also measured transverse magnetization with this system.

Figure 1 shows the magnetic susceptibility measured at $`H`$ = 1000 Oe. A broad peak is observed around 300 K, and the data can be fitted with a calculation by Bonner and Fisher (BF),  which indicates that the system is well described as an 1D Heisenberg AF. $`\chi _b`$ and $`\chi _c`$ are well fitted assuming different $`g`$ values, and $`J`$ = 540 K is obtained as an intra-chain coupling showing a good agreement with that for the polycrystalline sample.  As the temperature decreases to the WF transition temperature at $`T_N`$ = 8.8 K, deviation from the BF curves becomes significant. Below $`T_N`$, spontaneous magnetization appears as shown in the inset of Fig. 1. The magnetization along the $`b`$ axis ($`M_b`$) shows a typical evolution for ferromagnetic transition and is saturated at approximately 9$`\times `$10-2 emu/g for $`T`$ $``$ 0, which corresponds only to 0.42% of what we may observe when all the $`S`$ = 1/2 spins on Cu2+ site are aligned toward the same direction. Thus we eliminate the possibility of normal ferromagnetic transition. Ferrimagnetism is impossible, because BaCu2Ge2O7 has no magnetic ions other than Cu2+ that occupies an equivalent site. Since the sample is highly insulating, itinerant or band ferromagnetism is impossible.

To see the anisotropic magnetization, the field dependence of magnetization was measured at the WF state. Figure 2 shows low-field magnetization along three principal axes. A characteristic steep increase of the magnetization is found along the $`a`$ and $`c`$ axes up to $`H`$ = 280 Oe ($`H_c`$). $`M_a`$ reaches more than 1.6 $`\times `$ 10-1 emu/g, while $`M_b`$ remains approximately 9$`\times `$10-2 emu/g. Once magnetic field exceeds $`H_c`$, the slopes become similar to one another as shown in the inset of Fig. 2. The magnetization continues to increase along all the axes, which is a typical behavior of WF-LRO state. We did not observe a spin-flop anomaly in this field region, which indicates the lack of a strong easy-axis anisotropy in this system.

Let us first discuss the crystal structure to reveal this WF-LRO state. BaCu2Ge2O7 has a structure isomorphous to BaCu2Si2O7 (see Fig. 1 of Ref. 4) with slightly longer lattice parameters: $`a`$ = 7.028 Å, $`b`$ = 13.403 Å, and $`c`$ = 7.044 Å.  Two neighboring Cu2+ ions are connected by one O2- ion and they form a chain along the $`c`$ axis. However, the bond of Cu-O-Cu is not straight as shown in Fig. 3 (a). The angle of $`\mathrm{}`$ Cu-O-Cu (135) is larger than that of BaCu2Si2O7 (124), which is consistent with the larger intra-chain interaction. Two mechanisms of WF are known up to now: one is a single-ion anisotropy for NiF2,  and the other is a DM antisymmetric interaction for $`\alpha `$-Fe2O3,  where the antisymmetric exchange interaction $`𝐃[𝐒_i\times 𝐒_{i+1}]`$ gives spin canting.  There cannot be a single-ion anisotropy in BaCu2Ge2O7, because a magnetic moment in BaCu2Ge2O7 is attributed to a $`S`$ = 1/2 spin on Cu2+ ions. Then it necessarily follows that the DM interaction is the origin of WF in this compound. The local symmetry of nearest neighbor (nn) Cu2+ ions is actually low enough to allow the DM interaction. According to the rule given by Moriya,  the Dzyaloshinskii vectors $`𝐃_i`$ at bond $`i`$ are approximately defined as schematically drawn in Fig. 3 (a). The positions of only two neighboring Cu2+ ions and one intermediate O2- ion are taken into account. According to the room-temperature structural data,  $`𝐃_i`$ is written as $`𝐃_i=D`$ ((-1)i0.8332, (-1)i0.5462, 0.0860). Here, the components perpendicular and parallel to the chain have different effects on the spin order; the former ($`𝐃_i^{}`$ = (-1)$`{}_{}{}^{i}D`$(0.8332, 0.5462, 0)) changes its sign from one bond to the next, while the latter ($`𝐃_i^{}`$ = $`D`$(0, 0, 0.0860)) keeps its sign along the chain. For the rest of the paper we shall ignore the effect of $`𝐃_i^{}`$ and simply consider only $`𝐃_i^{}`$ from the following two reasons: 1) the magnitude of $`𝐃_i^{}`$ is far smaller than that of $`𝐃_i^{}`$, and 2) $`𝐃_i^{}`$ will stabilize a spiral order that is against the experimental observation.

We must know the arrangement of D vectors in the whole magnetic unit cell, because the bulk spontaneous magnetization is closely related to how D vectors are distributed. It is noted that D vectors and spin operators are both an axial vector when we consider the symmetry operations. For the $`a`$-axis directions, one chain is transformed to the neighboring one by the successive operations of a) reflection in the ac plane and b) rotation around the axis located at the middle of two adjacent Cu2+ ions along the  direction and parallel to the $`b`$ axis. Through these operations, the $`a`$ component of $`𝐃_i^{}`$ changes its sign twice, the $`b`$ component is unchanged, and the order of the concerning spin operators are exchanged. As a result, $`𝐃_i^{}`$ at the $`a`$-axis neighboring bonds points an opposite direction. We can also transform a chain to the $`b`$-axis neighboring one by the operation of reflection in the $`ac`$ plane and following parallel shift to the $`b`$ direction. This process changes the sign of the $`a`$ component only. In this way we determined all the D vectors in the magnetic unit cell as shown in Fig. 3 (b).

Next we discuss how these D vectors tilt spins to create the $`b`$-axis spontaneous magnetization. For this purpose, we must determine the principal direction of the spins. The crystal symmetry is high enough that it must be along the $`a`$, $`b`$, or $`c`$ direction. If we take only DM interaction into account, the fact that D vectors are not collinear along the $`b`$ direction can differentiate the three cases because there is a finite inter-chain interaction along the $`b`$ axis. However, the recent theoretical studies proved that DM interaction is always accompanied by an additional term quadratic to $`D`$, which recovers $`O(3)`$ symmetry of the spin-spin interaction.  This term is called as KSEA interaction, and has been established also experimentally.  Therefore we cannot predict the principal spin direction from DM and KSEA interactions even in the presence of inter-chain interaction.

On the other hand, experimental results indicates that the principal spin direction is along the $`c`$ axis. This is explained as follows; 1) Each chain must have weak ferromagnetic moment along the $`b`$ direction, while $`𝐃_i^{}`$ has the $`a`$ and $`b`$ components, 2) if the principal direction is along the $`a`$ or $`b`$ direction, spins are canted only toward the $`c`$ direction and we will have no $`b`$-axis magnetization, 3) only when the principal direction is along the $`c`$ axis, we will observe weak ferromagnetic moment along both the $`a`$ and $`b`$ directions. Thus, we safely concluded that the $`c`$ axis is the easy-axis direction in the actual system. Under this condition, spins on a single chain are confined in the plane perpendicular to D vectors (we shall call it as easy plane for convenience) as shown in Fig. 3 (c). It is convenient to define the WF moment per chain as the sum of the two sub-lattice magnetizations: $`𝐌_{WF,k}`$ = ($`𝐌_{1,k}+𝐌_{2,k}`$)/2, where $`𝐌_{i,k}`$ denotes the magnetic moment at the $`i`$-th sub-lattice in the $`k`$-th chain. Obviously the direction of the weak ferromagnetic moment is parallel to the easy-plane and simultaneously perpendicular to the $`c`$ axis. Once we obtain the direction of $`𝐌_{WF,k}`$ on each spin chain, the spatial configuration of $`𝐌_{WF,k}`$ determines the bulk magnetization. Relative directions of $`𝐌_{WF,k}`$ between the neighboring chain depend on the inter-chain interaction. In Fig. 3 (d), four different cases according to the signs of $`J_a`$ and $`J_b`$ are schematically drawn. Among the four cases, a spontaneous magnetization along the $`b`$ axis is realized only when $`J_a>`$ 0 and $`J_b>`$ 0, and therefore BaCu2Ge2O7 is considered to be the case. Consequently, the spin configurations in the magnetic unit cell os determined as shown in Fig. 3 (e).

So far, we have assumed that the DM interaction is the dominant perturbation to the spin system and the inter-chain interaction gives only a secondary effect. This assumption is probably valid by the following discussion. The magnitude of DM interaction can be estimated from the magnetization data. Using the angle between the easy plane and the $`ac`$ plane being 57.3, the canting angle of each spin from the $`c`$-axis direction is estimated as $`\theta `$ = 0.95.  Then the coarse relation, $`\mathrm{tan}2\theta D/J`$, gives $`D`$ as 18 K. On the other hand, we can also estimate the inter-chain interaction as an averaged value according to the chain mean-field theory:  $`J`$ = 540 K and $`T_N`$ = 8.8 K give the averaged inter-chain interaction as $`J^{}`$ = 2.9 K, which is sufficiently smaller than $`D`$. Of course, the direct estimation of $`J_a`$ and $`J_b`$ is necessary for further quantitative discussion, but $`|D||J_a|`$ and $`|D||J_b|`$ are probably valid in BaCu2Ge2O7. The above discussion means that BaCu2Ge2O7 is a good system to study the magnetism of weak-ferromagnetic spin chains that are weakly coupled with each other, and the sign of the inter-chain interaction can easily change the magnetic ground state demonstrated as the difference between BaCu2Ge2O7 and BaCu2Si2O7.

Finally we discuss the low-field anisotropic magnetization of BaCu2Ge2O7 upon this model structure. $`M_a`$ steeply increases up to $`H_c`$ first and suddenly changes its slope. This is qualitatively explained by the rotation of WF moment as schematically shown in Fig. 4 (a). When the field is applied parallel to the $`a`$ axis, $`𝐌_{WF,1}`$ and $`𝐌_{WF,2}`$ hold their directions, while $`𝐌_{WF,3}`$ and $`𝐌_{WF,4}`$ are rotated almost by 180, which can produce a large magnetization along the $`a`$ axis. We consider that the spins on the chains 3 and 4 are cooperatively rotated without changing their relative angle so much, because both $`J_b`$ and the Zeeman energy are too weak to change the relative angle of the spins against the dominant DM interaction. After this rotation has been finished at 280 Oe, large magnetization appears along the $`a`$ axis, and simultaneously net magnetization along the $`b`$ axis should disappear. In order to confirm this rotation, we carried out the transverse magnetization measurement, and found that the $`b`$-axis magnetization is actually suppressed with increasing field along the $`a`$ axis as shown in Fig. 4 (b). Since the geometry of the pick-up coil for our transverse magnetization measurement was the second differential type, finite $`M_b`$ gives an even-function response to the direct signal. Actually we got such a signal around 0 Oe, and observed the evolution of $`M_b`$. However, as the field approach 280 Oe, the direct signal changes its shape and finally turns to an odd-function response, which means that the $`M_a`$ induced by the field dominates the response of the pick-up coil. $`M_b`$ is now negligibly small, and we conclude that the spin structure shown in Fig. 4 (b) is realized. Such spin rotation is one of the characteristics of weakly-coupled spin chain, where DM interaction roughly keeps WF moment per chain, while the direction of WF moment is determined as a result of the competition between Zeeman energy and inter-chain coupling.

The difference of the magnetic ground states in BaCu2Ge2O7 and BaCu2Si2O7 provides us useful spin system to study the effect of randomness in the inter-chain interaction. Both compounds are the end materials of the solid-solution system of BaCu2(Ge1-xSix)2O7. The mixture of Ge and Si will introduce randomness to the sign of the inter-chain coupling along the $`a`$ axis. This randomness reduces the evolution of 3D long-range order as is reported in a polycrystalline samples by Yamada and Hiroi.

To summarize, we discovered the weak-ferromagnetic state in the 1D Heisenberg antiferromagnet BaCu2Ge2O7, which is to our knowledge the first uniform $`S`$ = 1/2 spin chain system that shows weak-ferromagnetic long-rang order. The Dzyaloshinskii-Moriya interaction adds a weak-ferromagnetic moment to each spin chain, and the relatively weak inter-chain interaction allows us to treat this system as weakly-coupled weak-ferromagnetic chains. The spin rotation by weak external field becomes possible by the combination of weak ferromagnetic moment per each spin and relatively weak inter-chain interaction, which is expected only in such a weakly coupled 1D spin system.

We wish to acknowledge valuable discussions with A. Zheludev, T. Yamada, and Z. Hiroi. I.T. would like to thank also Y. Ando for renting me a SQUID magnetometer. Work at the University of Tokyo is supported in part by the Grant-in-Aid for COE Research by the Ministry of Education, Science, Sports, and Culture of Japan.


text:
# Teleportation Criteria: Form and Significance

## Abstract

Our criteria for continuous variable quantum teleportation \[T.C.Ralph and P.K.Lam, Phys.Rev.Lett. 81, 5668 (1998)\] take the form of sums, rather than products, of conjugate quadrature measurements of the signal transfer coefficients and the covariances between the input and output states. We discuss why they have this form. We also discuss the physical significance of the covariance inequality.

Recently we have proposed criteria for the characterization of continuous variable quantum teleportation . It was shown that for any classical teleportation scheme (i.e. where no entanglement is shared) the following inequalities cannot be violated for minimum uncertainty Gaussian input states:

$`T_t`$ $`=`$ $`T_s^++T_s^{}1`$ (1)
$`V_t`$ $`=`$ $`{\displaystyle \frac{1}{2}}(V_{cv}^++V_{cv}^{})1`$ (2)
Here $`T_s^+`$ ($`T_s^{}`$) are the amplitude (phase) quadrature signal transfer functions from the input to output fields of the teleporter defined by

$$T_s^\pm =\frac{\mathrm{SNR}_{\mathrm{out}}^\pm }{\mathrm{SNR}_{\mathrm{in}}^\pm }$$
(3)
where the $`\mathrm{SNR}`$ are the signal to noise ratios of small classical test signals. The RF frequency test signals are placed on each quadrature of the input beam. The conditional variances between the input and output for each of the quadrature amplitudes are defined by

$$V_{cv}^\pm (\omega )=V_{\mathrm{out}}^\pm (\omega )(1C^\pm (\omega ))$$
(4)
where $`V_{\mathrm{out}}^\pm (\omega )=|\delta X_{\mathrm{out}}^\pm (\omega )|^2`$ are the standard spectral variances for the amplitude ($`\delta X^+`$) and phase ($`\delta X^{}`$) quadrature fluctuations of the output state assessed at some RF frequency $`\omega `$. The correlation function, $`C`$, is defined by

$$C^\pm =\frac{|\delta X_{\mathrm{in}}^\pm \delta X_{\mathrm{out}}^\pm |^2}{V_{\mathrm{in}}^\pm V_{\mathrm{out}}^\pm }$$
(5)
and is directly related to the $`\mathrm{SNR}`$’s via $`C^\pm =T_s^\pm `$  when cross-coupling between the quadratures can be ignored . The criteria in Eqs. (1) and (2) are semi-independent. Violation of either inequality indicates entanglement is present. A strong test of teleportation would require both inequalities be violated simultaneously.

Although there is no question of the validity of the inequalities in Eqs. (1) and (2) there has been some discussion about their significance . In particular, it has been queried as to why sums are formed of the conjugate quadrature measurements instead of the “usual” procedure of forming products (cf: Heisenberg uncertainty principle). In the following we justify the form of $`T_t`$ and $`V_t`$ by deriving them from more fundamental expressions.

Suppose the quantum fluctuations of our input field are described in the usual way by the zero-mean annihilation operator $`\delta a_{\mathrm{in}}`$ whilst our output field is described by $`\delta a_{\mathrm{out}}`$. The amplitude and phase quadrature fluctuations of the input field are defined respectively by $`\delta X_{\mathrm{in}}^+=\delta a_{in}+\delta a_{in}^{}`$ and $`\delta X_{\mathrm{in}}^{}=i(\delta a_{in}^{}\delta a_{in})`$. The quadrature fluctuations of the output field can be described by

$$\delta X_{\mathrm{out}}^\pm =Y^\pm \delta X_{\mathrm{in}}^\pm +Z^\pm \delta X_N^\pm $$
(6)
where $`Y^\pm `$ and $`Z^\pm `$ are c-numbers and $`\delta X_N^\pm `$ includes all added noise sources. Here again we will assume there is no cross coupling between the amplitude and phase quadratures of the input and output ,. Then

$$V_{\mathrm{out}}^\pm =|Y^\pm |^2V_{\mathrm{in}}^\pm +|Z^\pm |^2V_N^\pm $$
(7)
and

$$T_s^\pm =\frac{V_{\mathrm{in}}^\pm |Y^\pm |^2}{|Y^\pm |^2V_{\mathrm{in}}^\pm +|Z^\pm |^2V_N^\pm }$$
(8)
Because the information used to produce the output field traveled through a classical channel the added noise terms must be sufficient to ensure the generalized uncertainty principle is maintained for the photo-currents in the classical channels ,, this implies

$$\frac{|Z^+|^2}{|Y^+|^2}V_N^+\frac{|Z^{}|^2}{|Y^{}|^2}V_N^{}1$$
(9)
For a minimum uncertainty input state ($`V_{\mathrm{in}}^+V_{\mathrm{in}}^{}=1`$) and using Eq. (8) the inequality in Eq. (9) reduces to that of Eq. (1). Thus we see that the origin of the signal transfer inequality is, in fact, the standard product uncertainty. In quantum teleportation the shared entanglement acts as a “quantum key” that enables the inaccessible (because of the added noise) quantum information on the classical channel to be retrieved . The noise on the output field can then be less than the minimum required for the classical channel thus allowing the inequality to be violated.

We now consider how one might generalize the conditional variance used in QND  to teleportation. In QND one is only interested in how well the properties of one quadrature are preserved. In teleportation we wish to quantify how well the entire state is preserved. As we are working in the Heisenberg picture, in which it is the operators which evolve in time, this is equivalent to quantifying how well the field operator, $`a`$, is preserved. Caves  defined the variance of the field operator as

$$V_f=\frac{1}{2}\delta a^{}\delta a+\delta a\delta a^{}$$
(10)
Following the standard procedure for constructing a conditional variance we are thus motivated to consider the field correlation

$$C_f=\frac{|\frac{1}{2}\delta a_{\mathrm{out}}^{}\delta a_{\mathrm{in}}+\delta a_{\mathrm{in}}\delta a_{\mathrm{out}}^{}+\delta a_{\mathrm{out}}\delta a_{\mathrm{in}}^{}+\delta a_{\mathrm{in}}^{}\delta a_{\mathrm{out}}|^2}{\delta a_{\mathrm{out}}^{}\delta a_{\mathrm{out}}+\delta a_{\mathrm{out}}\delta a_{\mathrm{out}}^{}\delta a_{\mathrm{in}}^{}\delta a_{\mathrm{in}}+\delta a_{\mathrm{in}}\delta a_{\mathrm{in}}^{}}$$
(11)
If the fields are identical, i.e. $`\delta a_{\mathrm{in}}=\delta a_{\mathrm{out}}`$, then $`C_f=1`$. If the fields are completely independent, i.e. $`[a_{\mathrm{out}},a_{\mathrm{in}}^{}]=[a_{\mathrm{in}},a_{\mathrm{out}}^{}]=0`$, then $`C_f=0`$. We then construct the field conditional variance as

$$V_{cvf}=\delta a_{\mathrm{out}}^{}\delta a_{\mathrm{out}}+\delta a_{\mathrm{out}}\delta a_{\mathrm{out}}^{}(1C_f)$$
(12)
Using $`\delta a=\frac{1}{2}(\delta X^++i\delta X^{})`$, we can rewrite Eq. (12) as

$`V_{cvf}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(V_{\mathrm{out}}^++V_{\mathrm{out}}^{}{\displaystyle \frac{|\delta X_{\mathrm{out}}^+\delta X_{\mathrm{in}}^++\delta X_{\mathrm{out}}^{}\delta X_{\mathrm{in}}^{}|^2}{V_{\mathrm{in}}^++V_{\mathrm{in}}^{}}})`$ (13)
$`=`$ $`{\displaystyle \frac{1}{2}}(V_{\mathrm{out}}^++V_{\mathrm{out}}^{}{\displaystyle \frac{(\sqrt{T_s^+V_{\mathrm{out}}^+V_{\mathrm{in}}^+}+\sqrt{T_s^{}V_{\mathrm{out}}^{}V_{\mathrm{in}}^{}})^2}{V_{\mathrm{in}}^++V_{\mathrm{in}}^{}}})`$ (14)
For independent fields $`V_{cvf}1`$. If we assume that $`Y^+=Y^{}`$, i.e. the teleporter acts symmetrically on the two quadratures of the input field, then it is straightforward to show that in fact

$$V_{cvf}=V_t$$
(15)
Thus the inability of classical teleportation schemes to violate the inequality of Eq. (2) can be seen as showing that in a certain sense the input and output remain independent fields. The sum form of the inequality is seen to have its origin in the structure of the field operator as a sum of the two quadrature components. For asymmetric manipulations of the input quadratures ($`Y^+Y^{}`$) we find $`V_{cvf}V_t`$. Our analysis suggests that $`V_{cvf}`$ may be the more appropriate measure when quantifying such asymmetric schemes.

We now discuss the physical significance of violating the field covariance inequality Eq. (2). It has been argued by some that our criteria are too stringent because the presence of entanglement in the teleporter can be demonstrated without exceeding our inequalities. For example the field covariance of a lossless, symmetric teleportation scheme which has an entanglement resource $`V_{ent}`$ ($`V_{ent}=1`$ represents no entanglement whilst $`V_{ent}0`$ represents maximal entanglement) and a gain of $`\lambda `$ (assumed real) is given by

$$V_{cvf}=\frac{1}{2}(1+\lambda )^2V_{ent}+\frac{1}{2}(1\lambda )^2\frac{1}{V_{ent}}$$
(16)
Consider the case of unity gain ($`\lambda =1`$). Without entanglement ($`V_{ent}=1`$) we find $`V_{cvf}=2`$. Only by introducing entanglement ($`V_{ent}<1`$) can we obtain $`V_{cvf}<2`$. However it is not till we have introduced more than 50% entanglement ($`V_{ent}<.5`$) that we can obtain $`V_{cvf}<1`$ as per Eq. (2). Indeed it is in this intermediate region, where $`2>V_{cvf}>1`$, that the only experimental demonstration of continuous variable teleportation presently lies . One may ask if there is any qualitative difference between the types of correlation that can be observed when $`2>V_{cvf}>1`$ and those that can be observed when $`1>V_{cvf}>0`$. If qualitative differences exists then the more stringent definition of teleportation may be justified. We now give brief examples which illustrate that such differences do exist.

Entanglement Requirements. Firstly, it should be pointed out that true EPR entanglement is not required to reach the intermediate region. A single-mode squeezed beam, split on a 50:50 beamsplitter, is a sufficient resource . Although entangled , the beams so produced will exhibit non-classical correlations on only one quadrature with classical correlations on the conjugate quadrature. The signature of true EPR entanglement is non-classical correlations on both quadratures. The field covariance for a beam teleported with such a resource is given by

$$V_{cvf}=\frac{1}{4}(1+\lambda )^2(1+V_s)+\frac{1}{4}(1\lambda )^2(1+\frac{1}{V_s})$$
(17)
where $`V_s`$ is the level of single-mode squeezing in the resource. The value of Eq.(17) can never fall below 1 for any level of the squeezing or gain. Only when non-classical correlations can be observed on both quadratures can the field covariance fall below 1. Thus only with true EPR entanglement can we observe $`V_{cvf}<1`$.

Preservation of Non-classical Statistics. It is clearly of interest to ask under what conditions can teleportation of states with a non-classical character result in output states which retain that non-classical character. Suppose our input state is squeezed along the amplitude quadrature. We may ask under what conditions will the output still be squeezed. For symmetric, lossless teleportation the spectral variance of the output amplitude quadrature is given by

$`V_{out}^+`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+\lambda )^2V_{ent}+{\displaystyle \frac{1}{2}}(1\lambda )^2{\displaystyle \frac{1}{V_{ent}}}+\lambda ^2V_{in}^+`$ (18)
$`=`$ $`V_{cvf}+\lambda ^2V_{in}^+`$ (19)
It is immediately clear from Eq.(19) that it is impossible for squeezing to appear on the output unless $`V_{cvf}<1`$ is satisfied. In the limit of very strong squeezing ($`V_{sq}^+0`$) squeezing on the output is guaranteed when the field covariance is less than 1, however in general a field covariance less than 1 is a neccessary but not sufficient condition for squeezing to be preserved.

Preservation of Entanglement. Of even greater importance for quantum information applications is to ask under what conditions entanglement between two systems, of which one has been teleported, is preserved. We have recently quantified this question by looking at the violation of a Bell-type inequality . In particular we looked at the value of the Clauser-Horne variable , $`S`$, between entangled photon beams when one of the beams is teleported using a continuous variable method. Local realistic hidden variable theories place the following restriction on the value of $`S`$; $`S1`$. Quantum mechanical states allow $`S`$ to violate this inequality. The maximum violation occurs for non-maximally entangled states  and has the value; $`S=1.5`$ . Our result (in the limit of no loss) can be written

$$S=\frac{\frac{1}{2}(V_{cvf}1)+\lambda ^2(S_i+\frac{1}{2})}{(V_{cvf}1)+2\lambda ^2}$$
(20)
where $`S_i`$ is the value of $`S`$ which would be obtained between the beams before teleportation. Eq.(20) shows that if $`V_{cvf}<1`$ and $`S_i=1.5`$ then $`S>1`$. That is, provided that the entangled beams show a maximum violation of $`S`$ before teleportation, then some violation of local realism is guaranteed after teleportation if the field covariance falls below 1. As for squeezing this condition is neccessary, but not sufficient if $`S_i`$ does not have its maximum value.

We have derived the criteria for continuous variable teleportation from more fundamental arguments. We have shown by example that the properties that can be exhibited by the teleported system when the strong inequality $`V_{cvf}<1`$ is satisfied are qualitatively more “quantum mechanical” than in the region $`2<V_{cvf}<1`$. We have also noted that the strong inequality can only be satisfied by using true EPR entanglement in the teleporter.

This work was supported by the Australian Research Council.


text:
# Collective Multi-Vortex States in Periodic Arrays of Traps

## Abstract

We examine the vortex states in a 2D superconductor interacting with a square array of pinning sites. As a function of increasing pinning size or strength we find a series of novel phases including multi-vortex and composite superlattice states such as aligned dimer and trimer configurations at individual pinning sites. Interactions of the vortices give rise to an orientational ordering of the internal vortex structures in each pinning site. We also show that these vortex states can give rise to a multi-stage melting behavior.

Vortex lattice states in periodic pinning arrays have recently been attracting considerable attention as they represent ideal examples of an elastic lattice interacting with a periodic substrate. Experiments with periodic pinning arrays using microholes , blind holes,  and magnetic dot arrays  as well as simulations  have found numerous commensurability effects in which the pinning is enhanced at applied fields where the number of vortices equals an integer or rational fraction of the number of pinning sites. Direct imaging using Lorentz microscopy, of these commensurate vortex lattice structures in periodic pinning arrays, in which only one vortex is trapped at a single pinning site, has been conducted by Harada et al.  where a remarkable variety of different kinds of ordered vortex crystals were found. These experiments also show that for vortex densities greater than the first matching field, vortices can be positioned in the interstitial regions between the vortices trapped at the pinning sites.

Another possibility at the matching fields is that multiple vortices can occupy individual pinning sites. Bitter decoration experiments by Bezryadin et al.  have shown multi-vortex states in which up to nine vortices can be trapped at a single site. The individual vortices in the pinning sites are observed to form various types of structures and position themselves as far apart as possible to reduce their repulsive interaction energy. This work also demonstrates that for smaller pinning sites there can be a coexistence between multi-vortex states at the pinning sites with interstitial vortices between the pinning sites. As the pinning size is increased the number of vortices trapped at a single site also increases. Several other experiments have interpreted magnetization and creep measurements in terms of multi-vortex and composite vortex states

Vortex states in single large pinning sites and mesoscopic disks have also been studied . Kalfin and Shapiro  have predicted the vortex patterns for an increasing number of fluxons in large pinning sites which are similar to those observed in experiments . Multi-vortex states and configurations in thin superconducting disks can have preferred structures for a specific number of vortices that are captured . Similar systems of repulsive particles in a confining potential include electron crystals in Coulomb islands , and colloids in parabolic traps .

Despite the work done on vortex configurations in a single defect or disk the case of multi-vortex configurations in coupled arrays of defects has not been addressed. In this work we show through numerical simulations of logarithmically interacting vortices in 2D superconductors, that as a function of pinning size and strength, a series of novel collective multi-vortex and composite vortex lattice states can be realized. The long range interactions of the vortices causes the vortex structures internal to the pinning sites to give rise to an orientational ordering with respect to the vortex structures in the other pinning sites. Such states include orientationaly ordered dimer, trimer and composite states. Transitions between different vortex states can be observed as jumps in the critical depinning force as a function of pinning size. We also show that in these systems multi-stage melting processes can occur where the orientational degrees of freedom of the internal vortex lattice structure in the pinning sites melt at a much lower temperature than the whole vortex lattice does. Particular systems in which these states can be realized include thin-film superconductors with arrays of blind holes, disks or weak magnetic dots. Other possible physical realizations include charged colloidal particles or Wigner crystals in 2D trap arrays.

We consider a 2D system with periodic boundary conditions We numerically integrate the overdamped equation of motion for a vortex $`i`$

$$\eta \frac{d𝐫_i}{dt}=𝐟_i=𝐟_i^{vv}+𝐟_i^{vp}.$$
(1)
The vortex-vortex interaction potential is chosen to be logarithmic $`U_v=\mathrm{ln}(r)`$. The force on vortex $`i`$ from the other vortices is $`𝐟_i^{vv}=_{ji}^{N_v}_iU_v(r_{ij})`$ where $`r_{ij}=|𝐫_i𝐫_j|`$ is the distance between vortices $`i,j`$, and $`\eta `$ is the Bardeen-Stephen friction. We evaluate the periodic long-range logarithmic interaction with the resummations given in Ref . The pinning is modeled as a square array of attractive parabolic wells with $`f_i^{vp}=_{k=1}^{N_p}(f_p/r_p)\mathrm{\Theta }(|𝐫_i𝐫_k^{(p)}|)\widehat{𝐫}_{ik}^{(p)}`$. Here $`\mathrm{\Theta }`$ is the step function, $`𝐫_k^{(p)}`$ is the location of pinning site $`k`$, $`f_p`$ is the maximum pinning force, and $`r_p`$ is the radius of the pinning site. To obtain vortex configurations we start from a high temperature where the vortices are in a molten state and gradually cool until $`T=0`$. We have verified that our cooling rate is sufficiently slow so that the final state no longer depends on the cooling rate. In this work we will focus on the case of $`B/B_\varphi =`$ 2, 3, and 4, where $`B_\varphi `$ is the field at which there is one vortex per pinning site. This is adequate to capture the general features of the vortex states for higher matching fields. Results for higher $`B/B_\varphi `$ and incommensurate filling fractions will be presented elsewhere. The results for this work are for $`8\times 8`$ pinning arrays. We have also conducted simulations with larger systems and for different pinning lattice constants and have found the same features as seen for the $`8\times 8`$ systems. When the pinning sites are sufficiently large, multiple vortices can be captured per pinning site. We note that it is possible that vortices in the pinning sites can for certain parameters form individual giant vortex states  which we do not consider here.

In Fig. 1 we show the four vortex states that are possible for $`B/B_\varphi =3`$ for varying pinning strength and size. In Fig. 1(a) for the weak $`f_p=0.25`$ and small pinning $`r_p=0.1`$ the vortex-vortex interactions dominate and the vortex lattice forms a nearly triangular lattice. The vortex lattice still takes advantage of the pinning; however, only half the pinning sites can be occupied in order for the vortex lattice to have triangular ordering. We label this phase the commensurate elastic lattice. In Fig. 1(b), for stronger pinning, the pinning sites each capture one vortex and the vortex lattice is no longer triangular. The overall vortex lattice is still ordered with pairs of interstitial vortices alternating in positions. The state in Fig. 1(b) is identical to the state observed experimentally by Harada et al. , for the third matching field. For increased pinning radius $`r_p=0.5`$, and $`f_p=1.25`$ in Fig. 1(c) each pinning site can capture two vortices giving rise to a composite lattice of multi-vortices at the pinning sites and interstitial vortices. The vortices in the pinning sites repel one another and move to the edges of the pinning site. The interactions between the vortices in the pinning sites with the vortices in the other pinning sites gives rise to an orientational ordering between the dimers with a rotation of 45 degrees from one another. The interstitial vortices also form a periodic distorted square sub-lattice with unit cell for the overall vortex lattice, outlined in Fig. 1(c), consisting of two dimers and two interstitial vortices. In Fig. 1(d) for larger $`f_p=1.25`$ and $`r_p=0.8`$, each pinning site captures three vortices with the vortices forming triangles inside the pinning wells, producing a trimer state. The vortices in the pinning site again show nontrivial ordering with respect to each other. The vortex structures in alternating rows of pinning sites are rotated about 40 degrees from one another. Additionally the vortex structures in every other pinning site in the same row show a smaller rotation by about 5 degrees. The unit cell as outlined in Fig. 1(d) consists of four pinning sites and 12 vortices.

In Fig. 2 we show the evolution of the vortex states for $`B/B_\varphi =4.0`$ for increasing pinning size and strength. We again observe four vortex states. For weak and small pinning sites \[Fig. 2(a)\], $`f_p=0.25,r_p=0.15`$, only one vortex is captured per pinning site with the overall vortex lattice having a slightly distorted triangular ordering. We do not find a state where only a fraction of the vortices are filled, as in region I for $`B/B_\varphi =3.0`$, which can be understood by considering that in Fig. 2(a) the vortex lattice is already triangular. The state in Fig. 2(a) was also observed in experiments . In Fig. 2(b) where two vortices can be captured we do not find a completely ordered overall lattice. Here the lattice breaks up into domains with two separate orientations. This can be a finite size effect where the $`8\times 8`$ system is incommensurate with the unit cell of order. In Fig. 2(c) three vortices are captured per pinning site where the trimers are oriented with respect to one another and the interstitial vortices form a square sub-lattice. In Fig. 2(d) where each pinning site captures four vortices, the vortices in the pinning sites form a square lattice with the same orientation as the pinning lattice.

We have also conducted simulations for $`B/B_\varphi >4.0`$ and observe the same general features of the vortex states as outlined above, in particular the orientational ordered multi-vortex lattice states and ordered interstitial sublattices. The number of different kinds of states increases with the field.

In Fig. 3(a) we show the evolution of the phases for the $`B/B_\varphi =3.0`$ case for varied $`r_p`$ and $`f_p`$ with region I through IV corresponding to the phases in Fig. 1(a-d). As $`r_p`$ and $`f_p`$ are increased region II slowly decreases, region III maintains a roughly constant width, and region IV grows. Region I disappears for $`f_p>1.0`$. In Fig. 4(b) we plot the evolution of the phases for $`B/B_\varphi =4.0`$ with regions I$`^{^{}}`$ through IV$`^{^{}}`$ corresponding to the phases in Fig. 2(a-d). Here region I$`^{^{}}`$ is considerably larger than the other phases. For increasing $`f_p`$ region IV$`^{^{}}`$ grows while regions II and III keep roughly the same width.

The onset of the different vortex states as a function of $`r_p`$ at constant $`f_p`$ can also be observed as discrete jumps in the critical depinning force. The depinning is determined by adding an increasing driving force term to Eq. (1) in the symmetry direction of the pinning lattice and monitoring when the vortex velocities become non-zero.

In Fig. 4 we show that the collective multi-vortex states can give rise to a novel multi-stage melting behavior. We focus here on the simplest collective multi-vortex state at $`B/B_\varphi =2.0`$, when there are two vortex states. In the first state, every pinning site captures one vortex while interstial vortices form a square sub-lattice as previously observed in simulations  and experiments . The second state, shown in Fig. 4(a), is an aligned dimer state with all the vortex dimers being aligned at 45 degrees. We apply a temperature by adding Gaussian noise to Eq. (1). The dimers stay aligned until $`T=0.004`$ at which point they begin to freely rotate inside the wells as shown in the vortex trajectories in Fig 4(c). This destroys the orientational order, as seen in Fig. 4(b), and results in a liquid dimer state. As $`T`$ is increased the vortices remain confined in the pins until $`T=0.01`$ when the overall lattice melts with vortices diffusing randomly in the sample as shown in Fig. 4(d). To measure the melting transitions quantitatively in Fig. 4(e) we plot the angular correlation of the dimers $`C_\theta =_{nm}\mathrm{exp}(2i(\theta _n\theta _m))/N`$. We also plot a measure of the vortex displacements from their initial position at $`T=0`$ and compared to the displacements at a higher $`T`$: $`d(T)=(r(\tau )r(0))^2`$. Here $`\tau `$ is the time interval between temperature increments. For low $`T`$, $`C_\theta `$ is near unity and $`d(T)0`$ as the dimers remain aligned. At $`T=0.004`$, the dimers begin to freely rotate, as seen in the drop in $`C_\theta `$ to $`0.3`$ (notice that $`C_\theta `$ would be close to zero for a large system simulated over a long time) as well as the finite jump in in $`d(T)`$. Since the dimers are still confined to the pinning wells $`d(T)`$ stays at a constant value in the molten dimer state. For $`T>0.09`$ $`d(T)`$ rises rapidly as the vortices begin to jump out of the wells and diffuse randomly in the liquid state. The melting temperature of the dimer states is reduced as the pinning lattice constant is increased or the pinning radius is reduced. For higher matching fields a similar multi-stage melting behavior is observed where the loss of orientational ordering of the vortices in the pinning sites occurs before the loss of order in the overall lattice. More work is needed to determine the nature of the melting of the dimer or trimer states, such as whether it is continuous and similar to the melting in XY type models. It may also be possible that additional melting stages occur when the vortices are still inside the pinning site similar to the melting behaviors of vortices  or colloids  inside individual traps.

We point out that in addition to simulations with logarithmically interacting vortices we have also conducted 2D simulations with a finite range Bessel function interaction and find (provided the interaction range is sufficiently large) similar features in the vortex structures and multi-stage melting indicating that many of the phases we observed here are general features of systems of 2D repulsive particles in periodic arrays.

In conclusion we have studied the vortex states in thin-film superconductors interacting with periodic pinning arrays in which multiple vortices can be trapped at individual sites. We find that a rich variety of novel vortex states are possible as a function of pinning strength and pinning size. These states include collective dimer, trimer and composite states in which the vortex structures in the pinning sites exhibit an orientational ordering with each other. Transitions between the different states can be observed as a series of discrete jumps in the critical depinning force for varied pin radius. We also show that these systems exhibit a multi-stage melting where the structures internal to the vortex lattice melt before the overall vortex lattice melts. Besides vortices in superconductors these states may be observable for charged colliods in multi-trap arrays.

We thank C.J. Olson for critical reading of this manuscript and G.T. Zimanyi, R.T. Scalettar and F. Nori for usefull discussions. This work was supported by the Director, Office of Advanced Scientific Computing Research, Division of Mathematical, Information and Computational Sciences of the U.S. Department of Energy under contract number DE-AC03-76SF00098 as well as CLC and CULAR (Los Alamos National Laboratory).