Physica C 503, 123 (2014)

Vortex Transport and Pinning in Conformal Pinning Arrays

D. Ray1,2, C. Reichhardt2, C. J. Olson Reichhardt2, and B.Jankó1

1Department of Physics, University of Notre Dame, Notre Dame, Indiana 46556, USA
2Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

ARTICLE INFO
Article history:
Received 30 December 2013
Received in revised form 17 April 2014
Accepted 22 April 2014
Available online 5 May 2014
Keywords:
Vortex
Conformal pinning
Dynamics

Abstract
We examine the current driven dynamics for vortices interacting with conformal crystal pinning arrays and compare to the dynamics of vortices driven over random pinning arrays. We find that the pinning is enhanced in the conformal arrays for field densities less than 2.5 times the matching field. At higher fields, the effectiveness of the pinning in the moving vortex state is enhanced in the random arrays compared to the conformal arrays, leading to crossing of the velocity force curves.


1. Introduction
2. Simulations
3. Transport
4. Summary
References

1  Introduction

Many of the applications of type-II superconductors require that the system maintain a large critical current or the pinning of vortices in the presence of a magnetic field [1]. One of the key questions is, for a fixed number of pinning sites, what is the arrangement of the sites that maximizes the pinning. One approach to this question is the use of lithographic means to create various types of periodic arrays of artificial pinning sites [2,3,4,5,6] out of nanohole lattices [2,3,4,5] or arrays of magnetic dots [6]. In these periodic pinning systems, strong commensuration or matching effects can occur when the number of vortices equals an integer multiple n of the number of pinning sites [2,3,7,8]. At the matching conditions, there can be a peak in the critical current when the vortices form an ordered state [3,7,8]. For fields close to matching fields, interstitials or vacancies can appear in the ordered vortex structure and act as effective particles that are weakly pinned [9]. The enhancement of pinning at commensurate fields has also been observed in colloidal experiments [10] on periodic optical trap arrays. The colloids are repulsively interacting particles that have behavior similar to that of vortices, showing that understanding vortex dynamics on periodic or semi-periodic substrates is also useful for the general understanding of dynamics near commensurate-incommensurate transitions [11].
In superconducting samples at matching fields with two or more vortices per pinning site, there are several possibilities for how the vortices can be pinned. Multi-quanta vortices can be stabilized by each pin at the matching fields [8,12], multiple singly-quantized vortices can be stabilized at or near each pinning site [13,14], or there can be only a single vortex located at each pinning site with the remaining vortices trapped more weakly in the interstitial regions between pins [2,3,7,8]. It is also possible for combinations of these scenarios to occur, or for there to be a crossover from one regime to another as the field increases [2,8,15,16,17].
Another approach to pinning enhancement is the use of quasiperiodic pinning arrays such as Penrose tilings [18,19,20]. Commensuration effects still occur for such arrays; however, the strength of the pinning at nonmatching fields is generally greater than that found for nonmatching fields in periodic or random pinning arrays [18,19]. The random dilution of periodic pinning arrays produces peaks in the critical current not only at the matching fields but also at fields where the number of vortices matches the number of pinning sites in the original undiluted array [21]. Strong non-integer matching peaks in the critical current also appear in honeycomb pinning arrays, which are an example of an ordered diluted triangular pinning array [22]. Enhancement of the pinning at fractional fields can be achieved using artificial ice pinning array geometries [23], and these fractional matching peaks can be as strong as or even stronger than the integer matching peaks [24].
Fig1.png
Figure 1: The conformal pinning geometry used for simulated transport measurements. The drive is applied in the positive x-direction.
Recently, a new type of pinning geometry called a conformal crystal pinning array (illustrated in Fig. 1) was proposed [25]. This array is constructed by performing a conformal transformation [26] on a triangular lattice to create a gradient in the pinning density while preserving the local sixfold ordering of the original triangular array. Flux gradient simulations show that the overall critical current in the conformal pinning array (CPA) is enhanced over that of a uniform random pinning array [25] and is also higher than that of uniform periodic pinning arrays except for fields very close to integer matching, where the periodic pinning gives a marginally higher critical current. In the CPA, commensuration effects are absent and there are no peaks in the critical current. The gradient in pinning density present in the CPA enhances the pinning since it can match the gradient in the vortex density. Random pinning arrays with a density gradient equivalent to that of the CPA produce a small critical current enhancement compared to uniform random arrays; however, a CPA with the same number of pinning sites gives a substantially larger critical current, indicating that the preservation of the sixfold ordering of the pinning array is important for enhancing the pinning [25]. The simulation predictions were subsequently confirmed in experiments which compared CPAs to random and periodic arrays [27,28]. Other work on non-conformal pinning arrays containing gradients includes numerical studies of hyperbolic tesselations [29], as well as experiments on non-conformal pinning arrays with gradients in which the pinning was enhanced compared to uniform arrays [30].
The first numerical work on CPAs focused on flux-gradient driven simulations where the critical current is proportional to the width of the magnetization loop [25]. In such simulations, there is a gradient in the vortex density across the sample. One question is whether the CPA still produces enhanced pinning in systems driven with an applied current. Previous work indicated that the CPA produces a pinning enhancement compared to random arrays in this case as well [25]. In this work we further explore the current-driven system by varying the applied magnetic field and analyzing the vortex velocity as a function of external drive to produce a measurement that is proportional to an experimentally measurable voltage-current curve. We find that at very low vortex densities, the difference in critical current between random arrays and CPAs is small, and that as the field increases, the conformal arrays have stronger effective pinning, producing both a larger depinning force and a lower average vortex velocity in the moving state compared to random arrays. At higher fields, the CPA still has a high depinning threshold; however, once the vortices are in the moving state, the average vortex velocity for the random arrays can be lower than that for the CPA, indicating that the effectiveness of the pinning in the dynamic regime is suppressed for the CPA compared to the random pinning. We show that this arises due to the earlier onset of dynamical ordering [31,32] in CPAs compared to random pinning arrays at these higher magnetic fields.

2.  Simulations

We consider an effective 2D model of vortices where a single vortex i obeys the following equation of motion:
η d Ri

dt
= Fvvi + FPi + FDi.
(1)
Here η = ϕ02d/2πξ2ρN is the damping constant, d is the sample thickness, ϕ0 = h/2e is the elementary flux quantum, and ρN is the normal-state resistivity of the material. The vortex vortex interaction force is Fvvi = ∑NvjiF0K1(Rij/λ)Rij, where K1 is the modified Bessel function, λ is the London penetration depth, F0 = ϕ02/(2πμ0λ3), Rij = |RiRj| is the distance between vortex i and vortex j, and the unit vector Rij = (RRj)/Rij. The force from the pinning sites is given by FP. Various models for the pinning can be considered; here, we use parabolic attractive sites with
FPi = − Np

k=1 
(Fp/rp)(RiR(p)k)Θ[(rp − |RiR(p)k|)/λ],
(2)
where R(p)k is the location of pinning site k, Fp is the maximum pinning force, and Θ is the Heaviside step function. Other short-range pinning potentials should give essentially the same results as those reported here. An externally applied current produces a Lorentz force FDi=J ×B that is perpendicular to the applied current. Previous studies of this system predominantly focused on flux gradient driven geometries [25,33] in which vortices are added to or subtracted from a pin-free region outside of the sample. Here, we apply the driving force FD=Fdx in the x-direction and measure the vortex velocity in the x-direction. 〈Vx〉 = ∑Nvivi·x, where vi = dRi/dt. Using a spatially heterogeneous drive of the type that would arise in an experimental sample would not significantly alter the results; we assume that heating effects do not occur. The system geometry is illustrated in Fig. 1. There are Np pinning sites of radius rp = 0.12λ with a maximum force of Fp. The system size is 36λ×36λ, with periodic boundary conditions in the x and y-directions.

3.  Transport

Fig2.png
Figure 2: 〈Vx〉 versus Fd curves for the CPA (lower dark lines) and random arrays (upper light lines) in samples with np=1/λ2, rp=0.12λ, and Fp=0.55F0 for B/Bϕ= (a) 0.5, (b) 0.8, (c) 1.4, (d) 1.7, (e) 2.0, and (f) 0.2 (right lines) and 2.2 (left lines). In panels (a-e) the effectiveness of the pinning for the CPA is higher than for the random array.
Fig3.png
Figure 3: The difference between the velocity response in the random array and the CPA, ∆V = 〈Vxrand〉− 〈VxCPA〉, vs Fd for samples with B/Bϕ= 0.5, 1.4, and 2.2 (lower left to upper left). In this field range, ∆V is positive, indicating that the pinning is more effective in the CPA at low and intermediate values of Fd. At the highest drives, ∆V goes to zero.
Fig4.png
Figure 4: ∆V vs Fd for samples with B/Bϕ=2.5, 2.8, and 3.0 (upper right to lower right). For low Fd, the pinning is more effective in the CPA, as indicated by the positive value of ∆V; however, at intermediate Fd the pinning is more effective in the random arrays. At the highest drives, ∆V goes to zero.
In Fig. 2 we plot 〈Vx〉 versus Fd for samples with np = 1.0/λ2, rp=0.12λ, and Fp = 0.55F0 for uniform random arrays and CPAs at fields ranging from B/Bϕ=0.2 to 2.2. In Figs. 3 and 4 we plot the difference between the velocity response in the random array and the CPA, ∆V = 〈Vrandx〉− 〈VCPAx〉, as a function of Fd. At lower fields, shown in Fig. 3, the pinning in the CPA is more effective than in the random array, giving a positive ∆V, for all but the highest values of Fd. At B/Bϕ = 1.4, ∆V reaches a maximum value of 0.045 in Fig. 3, while for B/Bϕ=2.2 the maximum value of ∆V is only 0.015, indicating that as B/Bϕ increases, the difference in the pinning effectiveness of the CPA compared to the random pinning array is reduced. At higher fields, shown in Fig. 4, for low Fd above depinning, the pinning in the CPA is more effective than in the random array, but at higher Fd, the drag on the moving vortices created by the CPA is smaller than that of the random array. In each case, for large Fd, ∆V goes to 0 as the system enters the Ohmic regime where all the vortices are moving and the effects of the pinning are minimal. The transition to the Ohmic response regime occurs near Fd = 0.26 for B/Bϕ = 1.4 in Fig. 3, while for B/Bϕ = 2.2, the transition drops to a lower value of Fd = 0.07. Our results are not sensitive to a minor misalignment of the driving force with the x axis. Unlike a periodic pinning array, the CPA does not have easy flow channels, and we find almost no dependence of the 〈Vx〉 versus Fd curves on driving direction until quite large misalignments of order 40 degrees are applied.
Fig5.png
Figure 5: Lower solid curves: 〈Vx〉 vs Fd for a random array (light lines) and a CPA (dark lines). Upper symbols: P6, the fraction of sixfold coordinated vortices, vs Fd for a random array (light symbols) and a CPA (dark symbols). (a) B/Bϕ = 1.7. (b) B/Bϕ = 2.2. (c) B/Bϕ = 2.8. In (c), the vortices dynamically order at a lower drive for the CPA than for the random array, giving a lower value of 〈Vx〉 for the random array at the intermediate drives 0.1 < Fd < 0.2.
In Fig. 5 we plot 〈Vx〉 versus Fd for the same system in Fig. 2 at higher values of B/Bϕ = 1.7, 2.2, and 2.8. The critical depinning force Fc is higher in the CPA than in the random array for B/Bϕ=1.7 and 2.2, but at B/Bϕ=2.8 〈Vx〉 for intermediate drives is lower in the random array than in the CPA. This can be seen more clearly in the ∆V plot in Fig. 4, where for each of the fields ∆V is initially positive, but drops below zero as Fd increases, indicating that CPA pinning becomes less effective at intermediate values of Fd and that the average vortex velocity is higher in the CPA than in the random array. The reversal of the effectiveness of the pinning in the moving state produces as a crossing in the velocity vs force curves as shown in Fig. 5(c).
The reversal of the effectiveness in the pinning at intermediate drives occurs because the vortices dynamically order or partially crystallize at a lower drives in the CPA than in the random pinning array. It is known from current-driven simulation studies of random pinning arrays that a dynamical reordering transition can occur into a moving state that is partially crystalline or smectic-like [31,32]. In the dynamically ordered state, the vortex velocities are generally higher than in moving states with more random ordering, since the shear modulus of a random structure is much lower. A disordered vortex configuration has a higher probability of some vortices being temporarily pinned by the substrate, while in a moving crystal state, the vortices all move together and can not be individually trapped by pinning sites. For the random array, as the field increases, the drive FdOr at which the vortices begin to dynamically order decreases. FdOr is also a function of the pinning density np, and as np decreases, FOrd also decreases. In the CPA, the pinning density has a gradient, and as a result, there is a gradient in the value of FOrd across the sample. At the higher magnetic fields, the vortices can start to locally dynamically order in the lower pin density portions of the CPA sample. The partially ordered state forms in the low pin density regions during a transient time τo, and this state becomes disordered while passing through the high pin density regions during a transient time τd. As the field increases, these transient times change, and the vortices remain disordered if τd > τo, while for τo < τd the vortices can order. This means that in a random pinning array, the vortices are disordered when Fd < FOr; however, for a CPA at the same value of Fd, if τo < τd, an ordered moving vortex state will form and hence 〈Vx〉 for the CPA will be higher than for the random array. As B/Bϕ increases, τo decreases. This is consistent with the behavior in Fig. 2, where the extent of the range of Fd over which 〈VxCPA/〈Vxrand > 1 grows as B/Bϕ increases. It may be possible that at high enough B/Bϕ, the vortices in the CPA would immediately dynamically order as soon as they depin; in this case, the critical current for the random array would be higher than that of the CPA.
Fig6.png
Figure 6: The difference in the external drive Fd at which Vx=0.05 between the random and the conformal arrays, ∆F = Frandd(Vx = 0.05) − FCPAd(Vx = 0.05), vs B/Bϕ. At low fields, ∆Fd is small, at intermediate fields the CPA has stronger pinning (∆Fd > 0), and for B/Bϕ > 2.5 the random array has stronger pinning (∆Fd < 0).
In Fig. 5 we plot simultaneously 〈Vx〉 and the fraction of six-fold coordinated vortices P6 versus Fd for the random pinning and the CPA. In the dynamically ordered moving crystal state, P6 is close to 1 [31,32]. In Fig. 5(a) at B/Bϕ = 1.7, the pinning is more effective in the CPA over the entire window of Fd shown in the figure. At depinning, P6 drops for both types of pinning as the system enters a plastic flow regime. At higher Fd, P6 increases when the vortices begin to reorder, and in Fig. 5(a), P6 for the random array is higher than that for the CPA for Fd > 0.15. In Fig. 5(b) for B/Bϕ = 2.2, we find a similar trend; however, in Fig. 5(c) for B/Bϕ = 2.8, P6 is higher for the CPA than for the random array for Fd > 0.1. This also corresponds to the range of Fd over which 〈Vx〉 in the CPA is higher than in the random array. At Fd = 0.2, P6 reaches nearly the same value for both arrays, and the difference in 〈Vx〉 between the two arrays also vanishes. This result confirms that at high magnetic fields, the vortices dynamically order at a lower drive for the CPA than for a random pinning array.
We can roughly estimate the relative effectiveness of the pinning in the CPA and random pinning arrays by plotting the difference in the value Fd at which 〈Vx〉 = 0.05 for the two arrays, ∆Fd = Frandd(Vx = 0.05) − FCPA(Vx = 0.05). Figure 6 shows that at low B/Bϕ, ∆Fd is small and the difference between the random and conformal array is minimal due to the weak vortex-vortex interactions. In the limit of single vortex pinning, there is no difference between the two arrays and ∆Fd=0. Near B/Bϕ=1.4, ∆Fd reaches a maximum value; it then decreases and becomes negative for B/Bϕ > 2.5. The value of B/Bϕ at which ∆Fd drops below zero depends on the velocity value chosen for the measurement; however, Fig. 6 indicates that the enhanced effectiveness of the CPA is limited to B/Bϕ < 2.5. The CPA effectiveness may also depend on the size of the pinning sites and on Fp, both of which can be sample dependent.

4.  Summary

We investigated the current driven dynamics of vortices interacting with conformal pinning arrays and compared the effectiveness of the pinning to that of random pinning arrays with the same total number of pinning sites. The conformal pinning array is constructed by performing a conformal transformation of a triangular pinning array to create a new pinning array that has a density gradient but still conserves the local sixfold ordering of the original triangular array. We find that the for vortex densities less than twice the matching field, the critical depinning force for the conformal array is higher than that of the random array, and that in the moving vortex state that the velocity of the vortices in the conformal array is also lower than in the random array. At higher drives the difference between the two types of arrays is washed out due to dynamical reordering of the vortices. At higher fields, the critical depinning force for the conformal array remains higher than that of the random array, but at intermediate drives the average vortex velocity in the random arrays is lower than that in the conformal array, leading to a crossing of the velocity-force curves. This reversal of the pinning effectiveness arises because the vortices dynamically order at a lower drive in the conformal array than in the random array. There are still issues to consider in the conformal pinning array, such as performing conformal transformations on lattice structures other than a triangular array. It would also be interesting to investigate vortex ratchet effects of the type previously found in samples with random or periodic pinning arrays [34,35], as these effects may be enhanced in the conformal pinning arrays.
This work was carried out under the auspices of the NNSA of the U.S. DoE at LANL under Contract No. DE-AC52-06NA25396.

References

[1]
M. Tinkham, Introduction to Superconductivity (McGraw-Hill, New York, 1975); G. Blatter, M.V. Feigelman, V.B. Geshkenbein, A.I. Larkin, V.M. Vinokur, Rev. Mod. Phys. 66 (1994) 1125.
[2]
M. Baert, V.V. Metlushko, R. Jonckheere, V.V. Moshchalkov, Y. Bruynseraede, Phys. Rev. Lett. 74 (1995) 3269; Europhys. Lett. 29 (1995) 157.
[3]
K. Harada, O. Kamimura, H. Kasai, T. Matsuda, A. Tonomura, V.V. Moshchalkov, Science 274 (1996) 1167.
[4]
U. Welp, Z.L. Xiao, J.S. Jiang, V.K. Vlasko-Vlasov, S.D. Bader, G.W. Crabtree, J. Liang, H. Chik, J.M. Xu, Phys. Rev. B 66 (2002) 212507.
[5]
A.N. Grigorenko, S.J. Bending, M.J. Van Bael, M. Lange, V.V. Moshchalkov, H. Fangohr, P.A.J. de Groot, Phys. Rev. Lett. 90 (2003) 237001; A.D. Thakur, S. Ooi, S.P. Chockalingam, J. Jesudasan, P. Raychaudhuri, K. Hirata, Appl. Phys. Lett. 94 (2009) 262501; I. Swiecicki, C. Ulysse, T. Wolf, R. Bernard, N. Bergeal, J. Briatico, G. Faini, J. Lesueur, J.E. Villegas, Phys. Rev. B 85 (2012) 224502.
[6]
J.I. Martín, M. Vélez, J. Nogués, I.K. Schuller, Phys. Rev. Lett. 79 (1997) 1929; M.J. Van Bael, K. Temst, V.V. Moshchalkov, Y. Bruynseraede, Phys. Rev. B 59 (1999) 14674.
[7]
C. Reichhardt, C. J. Olson, F. Nori, Phys. Rev. B 57 (1998) 7937.
[8]
G.R. Berdiyorov, M.V. Milosevic, F.M. Peeters, Europhys. Lett. 74 (2006) 493; Phys. Rev. B 74 (2006) 174512.
[9]
C. Reichhardt, C.J. Olson, F. Nori, Phys. Rev. Lett. 78 (1997) 2648; C. Reichhardt, G.T. Zimányi, Niels Gr onbech-Jensen, Phys. Rev. B 64 (2001) 014501.
[10]
T. Bohlein, J. Mikhael, C. Bechinger, Nat. Mater. 11 (2012) 126.
[11]
C. Reichhardt, C.J. Olson Reichhardt, Phys. Rev. E 79 (2009) 061403; A. Vanossi, N. Manini, E. Tosatti, Proc. Natl. Acad. Sci. USA 109 (2012) 16429; D. McDermott, J. Amelang, C.J. Olson Reichhardt, C. Reichhardt, Phys. Rev. E 88 (2013) 062301.
[12]
C. Reichhardt, G.T. Zimányi, R.T. Scalettar, A. Hoffmann, I.K. Schuller, Phys. Rev. B 64 (2001) 052503.
[13]
A. Bezryadin, B. Pannetier, J. Low Temp. Phys. 102 (1996) 73.
[14]
C. Reichhardt, N. Gr onbech-Jensen, Phys. Rev. Lett. 85 (2000) 2372.
[15]
G. Karapetrov, J. Fedor, M. Iavarone, D. Rosenmann, W.K. Kwok, Phys. Rev. Lett. 95 (2005) 167002.
[16]
C.J. Olson Reichhardt, A. Libál, C. Reichhardt, Phys. Rev. B 73 (2006) 184519.
[17]
W.J. Zhang, S.K. He, H.F. Liu, G.M. Xue, H. Xiao, B.H. Li, Z.C. Wen, X.F. Han, S.P. Zhao, C.Z. Gu, X.G. Qiu, V.V. Moshchalkov, EPL 99 (2012) 37006.
[18]
V. Misko, S. Savel'ev, F. Nori, Phys. Rev. Lett. 95 (2005) 177007.
[19]
M. Kemmler, C. Gürlich, A. Sterck, H. Pöhler, M. Neuhaus, M. Siegel, R. Kleiner, and D. Koelle, Phys. Rev. Lett. 97 (2006) 147003; A.V. Silhanek, W. Gillijns, V.V. Moshchalkov, B.Y. Zhu, J. Moonens, and L.H.A. Leunissen, Appl. Phys. Lett. 89 (2006) 152507.
[20]
C. Reichhardt, C.J. Olson Reichhardt, Phys. Rev. Lett. 106 (2011) 060603.
[21]
C. Reichhardt, C. J. Olson Reichhardt, Phys. Rev. B 76 (2007) 094512. M. Kemmler, D. Bothner, K. Ilin, M. Siegel, R. Kleiner, D. Koelle, Phys. Rev. B 79 (2009) 184509.
[22]
C. Reichhardt, C.J. Olson Reichhardt, Phys. Rev. B 76 (2007) 064523; R. Cao, L. Horng, T.C. Wu, J.C. Wu, T.J. Yang, J. Phys.: Condens. Matter 21 (2009) 075705.
[23]
A. Libál, C.J. Olson Reichhardt, C. Reichhardt, Phys. Rev. Lett. 102 (2009) 237004; J. Trastoy, M. Malnou, C. Ulysse, R. Bernard, N. Bergeal, G. Faini, J. Lesueur, J. Briatico, J.E. Villegas, arXiv:1307.2881.
[24]
M.L. Latimer, G.R. Berdiyorov, Z.L. Xiao, F.M. Peeters, W.K. Kwok, Phys. Rev. Lett. 111 (2013) 067001.
[25]
D. Ray, C.J. Olson Reichhardt, B. Jankó, C. Reichhardt, Phys. Rev. Lett. 110 (2013) 267001.
[26]
F. Rothen, P. Pieranski, Phys. Rev. E 53 (1996) 2828.
[27]
Y.L. Wang, M.L. Latimer, Z.L. Xiao, R. Divan, L.E. Ocola, G.W. Crabtree, W.K. Kwok, Phys. Rev. B 87 (2013) 220501.
[28]
S. Guénon, Y.J. Rosen, A.C. Basaran, I.K. Schuller, Appl. Phys. Lett. 102 (2013) 252602.
[29]
V.R. Misko, F. Nori, Phys. Rev. B 85 (2012) 184506.
[30]
M. Motta, F. Colauto, W.A. Ortiz, J. Fritzsche, J. Cuppens, W. Gillijns, V.V. Moshchalkov, T.H. Johansen, A. Sanchez, A.V. Silhanek, Appl. Phys. Lett. 102 (2013) 212601.
[31]
A.E. Koshelev and V.M. Vinokur, Phys. Rev. Lett. 73 (1994) 3580.
[32]
C.J. Olson, C. Reichhardt, F. Nori, Phys. Rev. Lett. 81 (1998) 3757.
[33]
C. Reichhardt, C.J. Olson, J. Groth, S. Field, F. Nori, Phys. Rev. B 52 (1995) 10441; C. Reichhardt, J. Groth, C.J. Olson, S.B. Field, F. Nori, Phys. Rev. B 54 (1996) 16108.
[34]
C.J. Olson, C. Reichhardt, B. Jankó, F. Nori, Phys. Rev. Lett. 87 (2001) 177002.
[35]
W. Gillijns, A.V. Silhanek, V.V. Moshchalkov, C.J. Olson Reichhardt, C. Reichhardt, Phys. Rev. Lett. 99 (2007) 247002.



File translated from TEX by TTHgold, version 4.00.
Back to Home