text:
# Modelling and Measurement of Charge Transfer in Multiple GEM Structures

## 1 Introduction

A Time Projection Chamber (TPC) is foreseen as the main tracker of the detector for the TESLA linear collider . It should provide a momentum resolution of $`\delta (1/p_t)<210^4(\mathrm{GeV}/\mathrm{c})^1`$ to exploit especially the physics of final states with high energy leptons. Additionally the specific ionisation $`dE/dx`$ will be measured with a precision better than 5% to allow particle identification in hadronic jets.

To reach these goals readout planes with the finest possible granularity are required. They must provide a single point resolution of 100 – 150 $`\mu `$m and their systematic distortions must be controlled to better than 10 $`\mu `$m over the whole radius of 1.6 m. Enough charge amplification has to be provided by the system to keep the signals well above the noise level of modern readout electronics. At the same time the produced ions have to be suppressed intrinsically, because active gating seems impossible in between a TESLA bunch spacing of only 337 ns. Finally one aims for a minimum of material in front of the calorimeters, especially in the endcap region.

A conventional TPC using a multi wire plane for charge amplification is limited by the $`\stackrel{}{E}\times \stackrel{}{B}`$ effects in the region close to the wires. In strong magnetic fields these effects result in a broadening of the electron cloud and a worsening of the resolution. Additionally, as the wires define a preferred direction, the reconstructed hit location depends on the projected angle between track and wires. The readout is done via pads which detect the induced signals from the wires. This induced signal is broader than the cloud of arriving electrons and limits the granularity of the TPC as well as the double track resolution.

Using Gas Electron Multipliers (GEM)  as the charge amplifying system could solve some of the drawbacks of wire planes. A GEM consists of a polymer foil about 50 $`\mu `$m thick and metal coated on both sides. It is perforated by a high density of small holes with diameters of typically 50 – 100 $`\mu `$m. Voltages up to 500 V between the two conducting sides generate high electric fields in the order of 80 kV/cm inside the holes where the gas amplification of the drifting electrons occurs. When using GEM structures for the TPC end plate, the pads directly detect the amplified electron cloud which results in a fast and narrow charge signal. Also the slow ion tail is cut off since the ion cloud does not reach the induction region. A GEM foil shows no preferred direction, thus any $`\stackrel{}{E}\times \stackrel{}{B}`$ effects will be isotropic. And finally, using highly different electric fields in front and behind the GEM, the back drift of ions produced inside the GEM holes (ion feedback) can be largely suppressed.

To demonstrate the advantages of a TPC with GEM readout a prototype chamber will be built within the R&D activities of the linear collider TPC group . It will be used to investigate wether the momentum accuracy, the double track resolution and the quality of the $`dE/dx`$ measurement as demanded for a detector at TESLA can be reached. The end plates of this prototype will contain up to three planes of charge amplifying GEM foils. The operation conditions of this multi GEM structure should be optimised beforehand, because in the case of three GEM planes six electric fields have to be set. Therefore, a small test stand was set up at the RWTH Aachen.

The optimisation of a TPC readout plane should be done with respect to three important parameters:

* The electron transparency, $`T_{\mathrm{elec}}`$, i.e. the fraction of primary ionisation that experiences multiplication, should be near to 100%.
* The effective gain, $`G_{\mathrm{eff}}`$, i. e. the number of electrons reaching the anode pads per primary electron, should be sufficiently high.
* The ion feedback, $`B_{\mathrm{ion}}`$, i. e. the fraction of ions that reach the cathode plane, should be at a minimum.

To study the charge transfer mechanism of GEM structures we first introduce the following variables for a single GEM foil:

* The collection efficiency, $`C^{(+)}`$, i.e. the fraction of electrons (ions) collected into the GEM holes per number of electrons (ions) arriving.
* The gain, $`G`$, i.e. the factor by which the number of electrons is increased by gas multiplication inside the GEM holes
* The extraction efficiency, $`X^{(+)}`$, i.e. the fraction of electrons (ions) extracted from the GEM holes into the drift volume per number of electrons (ions) produced in the holes.
* The secondary extraction efficiency, $`X_{\mathrm{sec}}^+`$, i.e. the fraction of ions extracted from the GEM holes into the drift volume per number of ions which had been collected into the holes. This is different from the primary extraction, $`X^+`$, as will become clear in Section 4.

These transfer coefficients have been determined as a function of the electric and magnetic fields by measuring the various currents that appear on the electrodes of the GEM structure. Numerical simulations of the electric field fluxes allowed us to predict these variables and to compare the prediction with measurements. Additionally a parametrisation of the simulation has been established which can be used in future to optimize the working point of the GEM structure with respect to ion feedback and transparency.

## 2 The Test Chamber

### 2.1 Mechanical Setup

The test chamber is built on a PC board which is used as the anode plane. The sensitive area consists of a solid copper electrode of the same size as the GEM structures ($`10\times 10\mathrm{cm}^2`$). A bolt made from non conducting material is fixed in each of the four corners This allows the mounting of up to three GEM foils glued on frames with corresponding holes. The distances between the GEM foils (typically 2 mm) can be set by spacers. On top of the last GEM the cathode, consisting of a fine metal mesh, is fixed in the same way as the GEMs. The gas volume is closed by a composite frame to the sides. On top it is covered by a 24 $`\mu `$m thick mylar foil which allows the photons from an 55Fe source to penetrate into the chamber.

The chamber has been operated alternatively with Ar/CO2 82/18 for the measurements without a magnetic field (Section 4) or with Ar/CH4 95/5 for the magnet tests at a rate of 2-6 l/h (Section 5). It was irradiated from the top by an 55Fe source with an activity of the order of 100 MBq.

### 2.2 Electric Setup

Each electrode (GEM surfaces and cathode) is connected to an individual channel of a CAEN SY127 HV power supply, allowing the flexible setting of all electric fields and GEM voltages. Nano-Amperemeters , constructed in our electronics workshop and providing a resolution of approximately 0.1 nA, are inserted in the supply line of each HV channel. To measure the anode current, all pads and the surrounding copper are connected to ground via an individual nanoamperemeter. The high voltage control and current readout are handled by a custom application running on a Linux system.

### 2.3 Analysis

In Section 1, the charge transfer coefficients are defined as ratios of charges. We use the corresponding currents, as charges are produced and transferred continuously. These currents of drifting electrons and Ions ($`I_D`$, $`I_{T,n}`$, $`I_{Hole,n}`$ and $`I_I`$ in Figure 3) cannot be measured directly, but have to be calculated from the currents transferred on the electrodes ($`I_C`$, $`I_{GEM,n,C/A}`$, $`I_A`$). The fact that every current has an electron and an ion component must be considered when calculating transfer coefficients for electrons or ions separately. Depending on the individual setup and coefficient, this can become quite complex. Details are described in .

As an example, consider the electron extraction efficiency for GEM3. It is defined as the fraction of electrons extracted from the GEM holes into the induction region per number of electrons produced in the holes, and can be calculated from the measured currents as

$$X^{}=\frac{I_A}{I_A+I_{GEM3A}}.$$

This is a simple case, because in the region below GEM3 only electrons exist.

## 3 Numerical Simulations and Parametrisations

The numerical simulations are performed using the finite element program Maxwell 3D of Ansoft Co. . It allows to solve the Maxwell equations in three dimensional problems. For the simulation of the electric fields inside and between the GEM structures only the electrostatic problem solver is needed.

First the geometry of a unit cell of the GEM is constructed with the CAD part of the program. Material properties are assigned to the geometrical objects using the provided database. Above and below the GEM foil 200 $`\mu `$m of material-free drift field are added to the model. Then electric potentials are assigned to each conducting plane. The boundary conditions at the border of the simulated unit cell are chosen such that the electric field is parallel to the border planes (Neumann boundaries).

Maxwell generates the finite element mesh by dividing the model volume into small tetrahedrons. The mesh can be refined in critical regions like sharp edges for example. The electric potential is calculated at each knot and each centre of the connection lines of the mesh by solving the Poisson equation numerically. Inside the tetrahedrons the potential is interpolated by a second order polynomial.

After a numerical solution has been found, it can be visualized using the post processor. For example, field vectors or potential lines can be drawn, but also further calculations can be applied to the solution found. Especially mean electric fields or the electric flux through different areas is used extensively throughout this analysis. For details of the simulation studies, see .

### 3.1 Electric field inside the GEM hole

One of the most important parameters of a GEM structure is the mean electric field inside the GEM hole, $`E_{\mathrm{hole}}`$. We define this parameter from the integral of the electric field over the area of the hole centre, $`A`$ ($`\stackrel{}{n}`$ is its normal vector):

$$E_{\mathrm{hole}}=\frac{_A\stackrel{}{E}\stackrel{}{n}d^2r}{_Ad^2r}$$
(1)
The calculated field depends on the used GEM geometry, the voltage across the GEM foil and the electric fields of the drift regions above and below the GEM, as shown in Figure 5. The plot suggests the following linear parametrisation

$$E_{\mathrm{hole}}=aU_{\mathrm{GEM}}+b(E_D+E_I),$$
(2)
where $`a`$ and $`b`$ are parameters depending on the GEM geometry.

### 3.2 Electric field flux and charge transfer coefficients

Neglecting diffusion effects across the GEM structure and assuming zero magnetic field, all charges (electrons and ions) will follow the electric field lines. It is also assumed that during gas multiplication charges are spread homogeneously across the complete cross section of the GEM hole. In this approximation it is possible to calculate the charge transfer coefficients from the electric field flux.

Figure 6 shows a schematic plot of the electric flux through a GEM hole. The total flux coming from the anode $`\mathrm{\Phi }_{\mathrm{anode}}`$ and the flux going to the cathode $`\mathrm{\Phi }_{\mathrm{cathode}}`$ can be calculated from a surface integral over the electric field at the bottom and top of the GEM model, respectively. The total flux through the GEM hole is again the flux integral over the area of the hole centre.

Furthermore, the electric flux can be separated into contributions from distinct regions illustrated by different shades in Figure 6. There are the fluxes $`\mathrm{\Phi }_\mathrm{C}`$ and $`\mathrm{\Phi }_\mathrm{A}`$ which go directly from the electrodes to the copper planes of the GEM. The direction of the electric field changes between region $`\mathrm{\Phi }_\mathrm{C}`$ and $`R`$ and between $`\mathrm{\Phi }_\mathrm{A}`$ and $`Z`$. Therefore $`\mathrm{\Phi }_\mathrm{C}`$ and $`\mathrm{\Phi }_\mathrm{A}`$ can be extracted by a numerical integration over that part of the copper plane, where the electric field is pointing upwards. The electric flux $`D`$ which goes from one electrode through the GEM hole directly to the other electrode, the flux $`R`$ from the anode through the hole onto the top side of the GEM surface and the flux $`Z`$ going from the lower to the upper GEM surface cannot be calculated directly.

The collection efficiency $`C`$ is the fraction of the external flux that is collected into the GEM hole:

$$C_{\mathrm{top}}=\frac{D}{\mathrm{\Phi }_{\mathrm{cathode}}}=\frac{\mathrm{\Phi }_{\mathrm{cathode}}\mathrm{\Phi }_\mathrm{C}}{\mathrm{\Phi }_{\mathrm{cathode}}},C_{\mathrm{bottom}}=\frac{D+R}{\mathrm{\Phi }_{\mathrm{anode}}}=\frac{\mathrm{\Phi }_{\mathrm{anode}}\mathrm{\Phi }_\mathrm{A}}{\mathrm{\Phi }_{\mathrm{anode}}}$$
(3)
The extraction efficiency $`X`$ can be calculated from the fraction of the total electric flux originating from the hole and leaving the GEM structure towards cathode or anode:

$$X_{\mathrm{top}}=\frac{D}{\mathrm{\Phi }_{\mathrm{hole}}}=\frac{\mathrm{\Phi }_{\mathrm{cathode}}\mathrm{\Phi }_\mathrm{C}}{\mathrm{\Phi }_{\mathrm{hole}}},X_{\mathrm{bottom}}=\frac{D+R}{\mathrm{\Phi }_{\mathrm{hole}}}=\frac{\mathrm{\Phi }_{\mathrm{anode}}\mathrm{\Phi }_\mathrm{A}}{\mathrm{\Phi }_{\mathrm{hole}}}$$
(4)
### 3.3 Simulation results and parametrisation

Equations 3 and 4 suggest that the collection efficiency $`C`$ and extraction efficiency $`X`$ both vary with the ratio between the external flux and the flux through the GEM hole. This is shown in Figure 7, where the simulation results for both coefficients are plotted versus the ratio $`E_{\mathrm{ext}}/E_{\mathrm{hole}}`$. Different markers show simulation results with different GEM voltages ranging from 300 V to 500 V.

From the field flux model also a relation between the collection and the extraction efficiency can be derived. Let $`\mathrm{\Phi }_{\mathrm{ext}}`$ be the flux onto an external electrode (anode or cathode) of size $`A_{\mathrm{ext}}`$. As the electrode is far away from the GEM structure one can assume a homogeneous electric field and the external field strength is given by

$$E_{\mathrm{ext}}=\frac{\mathrm{\Phi }_{\mathrm{ext}}}{A_{\mathrm{ext}}}.$$
(5)
On the other hand, using the mean electric field inside the GEM hole $`E_{\mathrm{hole}}`$ as defined in Formula 1 one can write the total electric flux through the GEM hole as

$$\mathrm{\Phi }_{\mathrm{hole}}=E_{\mathrm{hole}}A_{\mathrm{hole}}.$$
(6)
Therefore, the ratio between extraction efficiency $`X`$ and collection efficiency $`C`$ is a linear function of the ratio between external and hole field,

$$\frac{X}{C}=\frac{\mathrm{\Phi }_{\mathrm{ext}}}{\mathrm{\Phi }_{\mathrm{hole}}}=\frac{E_{\mathrm{ext}}A_{\mathrm{ext}}}{E_{\mathrm{hole}}A_{\mathrm{hole}}}=\frac{1}{T_{\mathrm{opt}}}\frac{E_{\mathrm{ext}}}{E_{\mathrm{hole}}},$$
(7)
where the only geometrical parameter $`T_{\mathrm{opt}}`$ is the optical transparency of the GEM foil.

Figure 8 shows the predicted linear dependence between the ratio $`X/C`$ and the field ratio $`E_{\mathrm{ext}}/E_{\mathrm{hole}}`$ calculated for different GEM geometries. Six different geometries describing holes with different sizes and different shapes (conical and cylindrical) have been simulated. Geometries with equal optical transparency show the same slope.

The linear dependency between the logarithms of the collection efficiency $`C`$ and the field ratio $`E_{\mathrm{ext}}/E_{\mathrm{hole}}`$ leads to a parametrisation of the simulated results using the following ansatz ($`r,s\mathrm{}^+`$ are free parameters):

$$C=\{\begin{array}{cc}1\hfill & \text{for}E_{\mathrm{ext}}/E_{\mathrm{hole}}r^{1/s}\hfill \\ r(E_{\mathrm{ext}}/E_{\mathrm{hole}})^s\hfill & \text{for}E_{\mathrm{ext}}/E_{\mathrm{hole}}>r^{1/s}\hfill \end{array}$$
(8)
Following Formula 7 this gives in addition a parametrisation of the extraction efficiency $`X`$:

$$X=\{\begin{array}{cc}\frac{1}{T_{\mathrm{opt}}}(E_{\mathrm{ext}}/E_{\mathrm{hole}})\hfill & \text{for}E_{\mathrm{ext}}/E_{\mathrm{hole}}r^{1/s}\hfill \\ \frac{r}{T_{\mathrm{opt}}}(E_{\mathrm{ext}}/E_{\mathrm{hole}})^{1s}\hfill & \text{for}E_{\mathrm{ext}}/E_{\mathrm{hole}}>r^{1/s}\hfill \end{array}$$
(9)
Nearly perfect agreement between this parametrisation and the simulation results over the entire range of the field ratio can be seen from Figure 9.

## 4 Measurement Results and Comparison with Simulation

Figures 10 through 13 show the comparison between measurements made with the test chamber (Section 2) and the parametrisation obtained from the electrostatic simulations (Section 3).

All data except for the electron collection efficiency follow the simulation. That shows that those coefficients are governed mostly by the electrostatic conditions in the chamber, not by gas effects such as diffusion. This assumption is also supported by the fact that the extraction efficiencies for electrons and ions are equal within their errors.

For the measured collection efficiency of both electrons and ions, there is a steep decrease for field ratios approaching $`x=0`$ while the simulation remains constant at $`C=1`$. We believe that this effect is due to recombination of charge pairs at very low drift velocities caused by the small electric field in the drift region. In contrast, while showing a plateau of $`C=1`$ up to the same field ratio as the simulation, the collection efficiency for electrons decreases significantly steeper after that point. The origin of this effect is not yet understood.

The remaining difference between measurement and simulation for these variables is partly a consequence of the discrepancies between the ideal geometry of the GEM model in the simulations and the real GEMs. For example, changing the shape of the edges of the holes in the GEMs copper plating from perpendicular to the inclination of the conical holes in the polyimide causes a variation in the simulation results of the same order as the discussed differences.

The significant difference between primary and secondary ion extraction which does not occur for electrons, can be explained as follows: Through gas multiplication electrons are spread homogeneously across the cross section of every GEM hole they enter. Therefore a fraction of the electrons always occupies the radially outer regions of the hole penetrated by field lines ending on the top of the GEM and so decreasing the extraction efficiency. Primary ion extraction occurs through the same mechanism and shows the same values (compare Figures 10 and 12). However, ions subsequently collected into another GEM enter the hole in the middle, and because no multiplication and little diffusion occurs, they mostly remain on field lines continuing to the following GEM or electrode. This leads to the increased secondary extraction efficiency observed (Figure 12). Understanding these processes, calculations from electrostatic simulations describe even the effect of secondary extraction by taking the electric flux entering the GEM hole into account (Figure 6, 12).

## 5 Effect of a Magnetic Field

Our primary motivation for making tests in a magnet was to investigate whether there is a significant drop in collection efficiency for GEMs in high B-fields parallel to the electric fields as suggested by the Langevin Formula

$$\stackrel{}{v}_{Drift}\widehat{\stackrel{}{E}}+\omega \tau (\widehat{\stackrel{}{E}}\times \widehat{\stackrel{}{B}})+\omega ^2\tau ^2(\widehat{\stackrel{}{E}}\widehat{\stackrel{}{B}})\widehat{\stackrel{}{B}},$$
(10)
where $`\widehat{\stackrel{}{E}}`$ and $`\widehat{\stackrel{}{B}}`$ are unit vectors of the fields and $`\omega `$ is the cyclotron frequency

$$\omega =\frac{e}{m}B.$$

The last term proportional to $`B^2`$, which gives the contribution along the magnetic field lines, could cause a drop in collection efficiency for high magnetic fields. When this term dominates, most electrons will no longer be collected into a GEM hole but will stay on drift lines perpendicular to the GEM surface and eventually reach the GEMs copper coating. Those charges will be lost for the signal and consequently decrease the chamber’s $`dE/dx`$ capabilities due to the loss in primary ionisation statistics.

### 5.1 Simulation

To get a first qualitative view of the electrons behaviour, we calculated drift lines in the vicinity of a GEM foil pervaded by a 4 T magnetic field. This was done using the simulation program Garfield  which allows to take gas and magnetic field into account. We use a three dimensional electric field map of a GEM at 400 V with an external field of 200 V/cm obtained with Maxwell (see Section 3). The gas mixture selected was Ar/CH4 95/5 which is close to the mixture of Ar/CO2/CH4 93/2/5 proposed in the TESLA Technical Design Report . Each of the shown drift lines starts at a virtual straight ionisation track consisting of equally spaced electrons.

In Figure 14 the effect of the Lorentz term is well visible for electrons: they follow a track spiraling into the GEM hole. In contrast there is no contribution to the drift velocity along the magnetic field lines. Hence we have no indication for a significant drop in the collection efficiency for primary electrons.

### 5.2 Measurements

For a quantitative analysis, we have conducted a series of current measurements with our test chamber in a non superconducting magnet limited to a maximum field of 2 T. All tests were performed with a magnetic field parallel to the electric fields in the chamber.

For the test described here, we choose an electric field setup of 200 V/cm in the drift field, 5000 V/cm in the induction field and 2500 V/cm in the remaining fields. Each GEM is operated at a voltage of 310 V. The low drift field resembles the configuration in a TPC. The gas chosen is Ar/CH4 95/5, the same as in the Garfield simulation.

Opposite to the expected drop in collection efficiency, the signal current on the anode increases with $`B`$. At 2 T, it is approximately 6% higher than without a magnetic field. Though surprising, this result is still compatible with a decreasing collection as anticipated but being dominated by a rise in extraction efficiency or gas gain.

To find out about the mechanism for the observed rise in signal, another measurement was taken. The chamber was operated with all GEMs at 330 V and all fields at 1 kV/cm. This setting does not resemble the setting in a TPC (low drift field), but is chosen to allow simpler analysis of the charge transfer properties of the GEMs (same conditions at all GEMs, periodic setup): Because of the identical electrostatic conditions at all GEMs the corresponding transfer coefficients are equal and the signal current on the anode can be expressed as

$`I_A`$ $`=I_P`$ $`C_{GEM1}G_{GEM1}X_{GEM1}`$ (11)
$`C_{GEM2}G_{GEM2}X_{GEM2}`$
$`C_{GEM3}G_{GEM3}X_{GEM3}`$
$`=`$ $`I_PC_{GEM}^3G_{GEM}^3X_{GEM}^3`$
with $`I_P`$ being the current of primary ionisation by the photons from the 55Fe source.

Figure 16 shows the anode current normalised to its value at $`B=0`$. As before, it shows an increase with the magnetic field, this time by approximately 15% at 2 T. So obviously the influence of a magnetic field on the signal current is dependent on the electrostatic chamber setup.

From the currents in the chamber recorded during the measurement the extraction efficiency of GEM3, $`X_{GEM3}`$, was calculated. As discussed above, it is the same as for the two other GEMs, $`X_{GEM3}=X_{GEM}`$. The squares in Figure 16 show the values of $`X_{GEM}^3`$, again normalised to their value at $`B=0`$. The increase is even greater than that of the anode current (approximately 20%), which means that the product $`C_{GEM}^3G_{GEM}^3`$ must have decreased. Dividing the relative signal current by the calculated $`X_{GEM}^3`$ yields the value of $`C_{GEM}^3G_{GEM}^3`$, which is plotted as stars.

Obviously, the rise in the anode current is caused by an improved extraction efficiency for high magnetic fields. A drop in collection efficiency is hardly visible (note that for a single GEM, the effect is approximately $`\sqrt[3]{95\%}98\%`$ which is not significant taking the error of the measurement into account).

This result proves that the operation of triple GEM structures in a 2 T magnetic field is possible; for higher fields, corresponding measurements are scheduled.

## 6 Conclusions

It has been shown that the charge transfer coefficients of single GEM foils are predominantly determined by electrostatic conditions. Within the errors the results from measurements of the electric currents in the GEM structure are in agreement with simulations of the electric flux. The charge transfer coefficients can be described by a set of simple parametrisations which match the simulation almost identically.

A difference has been found between the extraction efficiencies of ions out of GEM holes where they were produced by gas multiplication (primary extraction) and out of GEM holes into which the ions were collected (secondary extraction). Electrons do not show this behavior, because they always experience multiplication in GEM holes (assuming a sufficiently large electric field) and are therefore distributed homogeneously across the volume of a GEM hole.

It has been demonstrated that for the working conditions in a TPC with low external field (the TPC drift field) and high electric field inside the GEM holes the ion feedback is intrinsically reduced to the percent level keeping full transparency for primary electrons. More sophisticated analysis of our parametrisation should allow to find an optimal field configuration of the GEM structure.

Applying a magnetic field up to 2 T perpendicular to the surface of a triple GEM structure does not result in a decrease of the effective gain. Further studies are needed to understand the quantitative behaviour of the electron collection efficiency in magnetic fields up to 4 T.

We thank F. Sauli and the CERN Gas Detector Development Group for supplying us with GEM foils and giving us a chance to look into the work of his team. We appreciate the support of T. Müller’s group at University of Karlsruhe, who manufactured the test chamber parts which helped us to get a quick start in conducting measurements. We are grateful to the IKP group of Forschungszentrum Jülich, especially O. Felden and T. Sagefka, for allowing us to use their test magnet. We acknowledge the support of T. Behnke, R.D. Heuer and R. Settles and the congenial atmosphere in the Linear Collider TPC group. Finally we thank S. Bachmann fruitful discussions and the careful reading of the manuscript.


format: text

text:
# New Detections of Optical Emission from Kiloparsec-scale Quasar Jets

## 1 Motivation

The number of active galaxies observed to have radio-bright kpc-scale relativistic jets emanating from their nuclei has grown to many 100’s over the last several decades (Bridle & Perley, 1984; Liu & Zhang, 2002), but only a handful of them were known to have X-ray counterparts prior to the launch of the Chandra X-ray Observatory. Chandra has in just the last three years netted X-ray detections of jet features in over 30 more radio sources (see Harris & Krawczynski, 2002, and the most recent updates in the associated website – http://hea-www.harvard.edu/XJET/), and many more are being discovered as a result of dedicated searches (see e.g. individual contributions by R. M. Sambruna, H. L. Marshall, and D. A. Schwartz). With the increased interest in the X-rays, studies in the optical/UV regime have also enjoyed a resurgence of activity. This is because the optical is the only other accessible window in the electromagnetic spectrum where comparable sub-arcsecond resolution can be readily achieved with the Hubble Space Telescope111Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555., and bridges the $``$6-7 orders of magnitude gap in frequency coverage between radio and X-ray studies of resolved features. It turns out that the optical fluxes and spectra can be a key diagnostic in distinguishing between different X-ray production mechanisms. Specifically, the faintness of the optical knots in the ‘first-light’ X-ray jet detection in the powerful quasar PKS 0637–752 (Chartas et al., 2000; Schwartz et al., 2000) helped to rule out canonical synchrotron and synchrotron self-Compton models, and established the importance of X-rays produced via inverse Compton scattering of CMB photons by synchrotron emitting electrons in the kpc-scale jets (Tavecchio et al., 2000; Celotti, Ghisellini, & Chiaberge, 2001). Also, if the emission mechanism responsible for the optical emission is synchrotron radiation, as is commonly accepted, observations in this waveband probe the physics of the most energetic and short-lived electrons, and are key to studying acceleration processes in the jet.

In the pre-Chandra era, optical jets were predominantly found in low power sources with the HST (e.g., Sparks, Biretta, & Macchetto, 1994; Scarpa & Urry, 2002), and optical counterparts to jets in quasars were difficult to detect with ground based facilities. Only a handful of examples existed, courtesy of the HST (e.g., Ridgway & Stockton, 1997), so optical jets in quasars were thought to be rare. However, Chandra’s successful X-ray detections of quasar jets has forced observers to examine new and archival HST data more carefully, resulting in many new detections in the optical. This resulted in the identifications of a number of bright optical knots coincident with peaks in the radio jets of about ten quasars in a recent Chandra/HST cycle-2 search for such emission (Sambruna et al., 2002, R. Scarpa et al., in prep.). We applied the same techniques employed in this survey and found a bright optical/UV knot in the kpc-scale radio jet of 3C 279 in archival HST data (Cheung, 2002) which encouraged us to search for more such emission in other available archival data.

## 2 New Detections and Search Strategy

Figures 1 & 2 show two of our first three optical detections in the quasars 3C 454.3 and PKS 1229–021 as a result of our search. The third detection is in the bright radio hotspot in the lobe-dominated quasar 3C 263 which was recently published by Hardcastle et al. (2002) along with the Chandra X-ray detection of the feature. These quasars were known to us to exhibit bright radio jets and were obvious candidates to search for optical counterparts. The optical emission in these and other known kpc-scale quasar jets (Ridgway & Stockton, 1997; Chartas et al., 2000; Sambruna et al., 2002; Cheung, 2002, R. Scarpa et al., in prep.) is compact, match peaks in the radio jets well, and is consistent with synchrotron emission from highly relativistic electrons in the jet.

With the growing number of optical jet detections in radio-loud sources, we can begin to make comparisons between the jet properties of different subclasses. One emerging trend is that the radio-to-optical spectral indices of the jet features in quasars (see above references) appear to be steeper by $`\mathrm{\Delta }\alpha 0.2`$ than those measured in lower power radio sources (see e.g., Sparks, Biretta, & Macchetto, 1994; Scarpa & Urry, 2002).

After these initial detections, we decided on a more systematic approach to mining through the available archival data. We are currently searching the HST archive database for images of radio bright ($`>`$ 0.5 Jy at 5 GHz) active galaxies obtained after $``$1994 (post-COSTAR), mostly concentrating on WFPC2 images which have $``$ 1 orbit of total exposure time. The resultant list of candidates are being correlated with the compilation of Liu & Zhang (2002) of sources known to exhibit radio jets. As a result of this ongoing effort, we have found optical counterparts to knots in an additional three quasar jets and have already obtained archival multi-frequency MERLIN and VLA data for our analysis. These results will be presented along with the complete analysis of the detections presented here in a forthcoming paper.


text:
# The TESLA Time Projection Chamber

## 1 Introduction

TESLA (TeV-Energy Superconducting Linear Accelerator) is an $`e^+e^{}`$ linear Collider designed to operate in the energy range beyond the mass of the Z0 and up to the TeV region. If approved, TESLA would allow a wide physics program, focusing among other things on precision tests of the established Standard Model, the Higgs mechanisms and the discovery of new particles.

To optimally exploit the physics potential of TESLA, the detector has to satisfy some stringent requirements. It has to have a very good track momentum resolution to measure the Z0 recoil mass, an excellent vertex resolution for heavy flavour identification, a high resolution in the reconstruction of electromagnetic and hadronic showers in the calorimeter, and has to be as hermetic as possible. The general concept of a detector fulfilling these criteria, discussed within the ECFA-DESY study, is described in details in .
In this contribution, we describe the current design of the Time Projection Chamber (TPC) and the different ongoing studies concerning the readout technologies.

## 2 The TESLA central Tracker

A large TPC is proposed as the main inner tracker at TESLA. The actual design, taking into account the different performance goals, is described in details in the TESLA TDR .

The TPC, as shown schematically in figure 1, is of cylindrical design with an inner radius of about 32 cm, an outer radius of about 170 cm and an overall length in the beam direction of about 2 $`\times `$ 273 cm. The internal radius is determined by the size of the mask system, whereas the external radius is constrained by the fact that the calorimeter is inside the coil. These large dimensions allow the tracks to be measured with about 200 $`(r,\varphi ,z)`$ space points, so that the global tracking efficiency remains close to 100 % in the overall volume and a momentum resolution $`\delta \left(1/p_t\right)<1.5\times 10^4`$ / GeV in the central region can be achieved. In addition the TPC should measure the specific energy loss of particles thus aiding in the identification of charged particles. With Ar-based gas mixtures, operating at atmospheric pressure, it is expected that a resolution of $`\delta (\mathrm{log}\left(\text{dE/dx}\right))4.5\%`$ can be reached, enough, to provide a separation between pions and kaons of better than two sigma for momenta between 2 and 20 GeV.

Operating a TPC in the environment at TESLA presents a number of challenges. Presently no dedicated hardware trigger is forseen within a bunch train, so that all detectors should be active throughout one full train about 1 ms. One readout cycle of the TPC will integrate over approximatively 150 bunch crossings at TESLA, and many readout cycles will follow continuously until a train has passed. This makes it important that the readout is granular enough to resolve the expected number of hits. In addition means have to be found to avoid the need for gating the TPC to remove the buildup of positive ions: during one bunch trains ions will travel significantly into the TPC drift volume, potentially distorting the drift field.

A possible solution to handle both the high granularity and to suppress the ions drift into the drift volume is the use of Micro-Pattern Gas Detectors as readout devices. Possible candidates are at the moment MicroMEGAS chambers  or Gas Electron Multipliers (GEM) foils .

Micro-Pattern Gas Chambers, due to their intrinsically smaller distances between amplification elements compared to the traditional wire chambers, offer the promise to significantly reduce systematic distortions in the electron drift close to the amplification region due to the so-called $`\stackrel{}{\text{E}}\times \stackrel{}{\text{B}}`$ effect.

## 3 Ongoing Studies

A challenging Research and Development program is underway to meet the current design goals for a TPC at a linear collider like TESLA. Institutes from a number of different countries are participating (Aachen, Carleton, DESY, IPN Orsay, Karlsruhe, Krakow, LAL, LBNL, MIT, MPI Munich, NIKHEF, Novosibirsk, Rostock, Saclay and Victoria). The main focus of the work is at the moment on the development of a valid amplification and ion-suppression scheme (see reference ). Over the next years, efforts to develop a compact readout scheme and to understand the design of the mechanical structure of such a large TPC are planned. In the following sections we present the different first results for these different key issues.

### 3.1 Gas Amplification Systems

Several studies focus on the TPC readout using Micro-Pattern Gas Detectors instead of the wire used conventionally to produce the gas amplification that is necessary to read the TPC signal. Gas Electron Multiplier (GEM) and MicroMEGAS, sketched in figure 2, are two attractive candidates for the readout planes.

A GEM foil is a sandwich consisting of thin conductive layers of copper separated by a thin polymer film of capton. The foil is perforated with small holes arranged in a hexagonal lattice. When operating in the detector, a potential applied to the two conductive layers generates a strong electric field, typically 80 kV/cm, in the small holes. The electric field lines guide charged particles through the GEM holes, with minimal losses due to collisions with the foil. Electrons are accelerated in the strong field lines within the GEM foil, producing secondary ionizations and charge multiplication which are then transfered by the electric field to the readout pads where they are collected. A unique property of the GEMs is their capability to operate in cascade, i.e. in a multi-GEM structure, offering several advantages like a higher gain and an effective ion feedback suppression. Most of the groups working with GEM readout use GEM towers with two or more GEMs.

An attractive alternative technology to GEM is that of MicroMEGAS. An uniform high-field, typically 30 kV/cm, is produced between a thin metallic mesh, few microns thick, stretched at a distance of 50-100 $`\mu `$m above the readout pad plane and held by insulating pillars. On top of the mesh, the electrons are transferred to the amplification gap, where they are multiplied in an avalanche process, resulting in a large number of electron-ion pairs. A MicroMEGAS has similar advantages to a GEM as far as simplicity of construction, low cost and efficiency of ion-feedback suppression are concerned. It has an excellent potential for dE/dx resolution due to the fact that the gain is nearly independent of the gap thickness. Moreover it is robust and can be built at least in parts using commercially available components.

R&D studies focus among other things on the gain stability and homogeneity to reach a dE/dx measurement with 5% precision. This is illustrated in figure 3 which shows the gain uniformity with time using a two-GEM tower structure.

### 3.2 Ion Feedback Suppression

A crucial issue at TESLA is to prevent the ions produced during the amplification to return to the TPC drift volume and cause field distortions. The use of GEMs or MicroMEGAS promises to reduce the ion feedback. The degree of reduction depends strongly on the exact choice of operational parameters and the Micro-Pattern Gas Detector used. Several ongoing studies, using GEMs and MicroMEGAS, demonstrated a suppression to the level of 10-2 (figure 4). However, a level of suppression to a level of 10-4 is needed to ensure that the positive-ion effects are at an acceptable level. In case this goal cannot be reached and to guarantee a stable and robust chamber operation, the insertion of a gating plane, gated in between bunch trains, might be needed in order to minimize the effects of ion feedback.

### 3.3 Simulations Studies

To optimally exploit GEMs and MicroMEGAS and in addition to the experimental program, extensive simulation studies are being performed. These concern mainly the study of the charge transfer mechanism in GEM or MicroMEGAS structures. These numerical simulations are of crucial importance, since they determine the optimal parameter choice needed to reach the different TESLA TPC milestones.

Figure 5 shows, as the validation of the simulation of a GEM hole, the comparison between measurement and simulation of the evolution of the primary and secondary extraction efficiencies respectively defined as the fraction of ions extracted from the GEM holes into the drift volume per number of ions produced in the holes and the fraction of ions extracted from the GEM hole into the drift volume per number of ions which have been collected in the GEM hole.

### 3.4 Resolution and Pad Geometry Studies

The charge produced during the avalanche in the GEMs is collected by readout pads located typically a few mm behind the last GEM. Since the number of pads is limited, the distance between pad centers is large compared to the size of the electron avalanche cloud (typically a few to several 100 $`\mu `$m). The charge from a track after amplification in the GEMs is therefore sometimes collected on a single pad. This effect is even more important with magnetic field at TESLA since the transverse diffusion is then smaller. The major consequence will be that the centre of gravity method cannot be applied anymore for position calculation. Ongoing simulation studies investigate more sophisticated pad geometries, like chevron-shaped pads which may lead to a better point resolution if the charge sharing between neighboring pads is more important.

Figure 6 shows the measured point resolution $`\sigma _x`$ for rectangular and chevron shaped pad with a size of 14 $`\times `$ 14 mm2 and for different drift lengths using GEMs and without magnetic field. For small drift distances the point resolution with chevrons is better than with rectangles. This effect of the pad geometry gets smaller with increasing drift distances since the influence of the diffusion of the electron cloud gets larger. The improvement is clearly visible but simulation studies concerning mainly the pad response function description, the optimal induction gap between the last GEM and the pad plane, are still required for a better description of the pad charge sharing mechanism.

Another approach to make the electron charge cloud broader, and get then a better point resolution, uses a resistive anode foil placed a few ten of $`\mu `$m between the readout plane and the last GEM. This work is illustrated by the plot on figure 7 showing the signal detected with a size of about 2 mm for an initial charge cluster size of 1 mm.

### 3.5 Induced Pulses in a GEM

The possibility to detect induced pulses in a GEM setup is also being studied. During the buildup of charges in the GEM hole a small induced pulse is expected in the pad opposite to this hole. This signal is also visible in the neighbor pads and can then be used to reconstruct the position of the avalanche and have a resolution significantly better. The price to pay for this improvement is that the readout electronics has to be able to detect the small pulses and that a very fast electronics is needed.

## 4 Conclusion and Outlook

Several institutes have started a joint R&D effort with the goal of developing the technologies needed for a large TPC at the next generation of linear colliders. The main challenge is to design a continuous-tracking, high performance TPC which should be substantially more performant than any existing TPC. Important progress has been made to better understand the different aspects of the TPC. Particular emphasis is placed on the amplification, the ion-suppression scheme and the large scale-integration electronics. While first results concerning these different points look very promising, some open questions like the readout behavior in a strong magnetic field, or the results with large TPC prototypes, still remain unanswered and need more investigations in the near future.

## 5 Acknowledgments

I wish to thank all the members of the ECFA-DESY TPC study for providing material for this contribution. I’m grateful to Ties Behnke for reading this manuscript.


text:
# Structure and Magnetization of Two-Dimensional Vortex Arrays in the Presence of Periodic Pinning

## I Introduction

In the mixed phase of type-II superconductors, magnetic flux penetrates the sample in the form of quantized vortex lines . The amount of flux carried by each vortex line is equal to the basic flux quantum $`\mathrm{\Phi }_0=hc/(2e)=2.07\times 10^7`$G.cm2. These vortex lines form a special physical system known as “vortex matter”. In the absence of any pinning sites in the material, the vortex lines form a triangular lattice known as the Abrikosov lattice .

Equilibrium and transport properties of the mixed phase of type-II superconductors are strongly affected by the presence of pinning centers, either intrinsic to the system or artificially generated. Understanding the effects of pinning in these systems is very important for practical applications because the presence of pinning strongly influences the value of the critical current in the mixed phase.

In recent years, a variety of nanofabrication techniques have been used to create periodic arrays of pinning centers in thin-film superconductors . Such arrays may consist of microholes (“antidots”) , defects produced by the bombardment of ion  or electron  beams, or magnetic dots . These pinning centers are “strong” in the sense that each pinning site can trap one or more vortices at low temperatures. The effects of periodic pinning depend strongly on the relative values of $`B_\varphi `$ and $`B`$, where $`B_\varphi =\rho _p\mathrm{\Phi }_0`$ ($`\rho _p`$ is the areal density of the pinning centers) is the so-called “matching field”, and $`B`$ is the magnetic induction that determines the areal density $`\rho _0`$ of vortices ($`\rho _0=B/\mathrm{\Phi }_0`$). The filling ratio, $`n`$, defined as $`nB/B_\varphi `$, measures the commensurability of the vortex system with the underlying pin lattice. The interplay between the lattice constant of the pin array (determined by $`B_\varphi `$) and the intervortex separation (determined by $`B`$) can lead to a variety of interesting effects in such systems.

Some of these effects have been observed in recent experiments. Imaging experiments using various techniques such as Lorentz microscopy  and scanning Hall probe microscopy  have shown the formation of ordered structures of the vortex system at low temperatures for commensurate values of $`n`$. Magnetization measurements  in the irreversible (vortex solid) regime have demonstrated the occurrence of anomalies at certain harmonics of $`B_\varphi `$. The effectiveness of pinning at integral values of $`n`$ has been found  to produce regularly spaced sharp minima in the resistivity versus field curve. A pinning-induced reconfiguration of the vortex lattice has been observed  in a thin-film superconductor with a rectangular array of magnetic dots. Some of these effects have been studied theoretically, using analytic  and numerical  methods. Experimental realizations of a system of interacting “particles” in the presence of an external periodic potential are also obtained in colloidal suspensions in interfering laser fields , and in periodic arrays of optical traps .

In this paper, we have used analytic and numerical methods to analyze the zero-temperature structure of vortex arrays in the presence of periodic pinning. We have also carried out a numerical study of the zero-temperature equilibrium magnetization of a superconducting film with a square array of pinning centers as a function of the applied field. In section II, we consider the ground states of a vortex system in a square array of pinning centers for fillings less than unity. We look at a class of structures that are Bravais lattices with one vortex per basis if the filling $`n`$ is of the form $`1/q`$, and with $`p`$ vortices per basis if $`n=p/q`$ ($`p`$ and $`q`$ are integers greater than unity, with $`p<q`$). The structure with the lowest energy in this class can be obtained rather easily. We find that the “ground-state” structure obtained this way matches those obtained from experiments  and simulations  for a large number of simple fractional values of $`n`$. The results obtained in this section can also be used to predict the ground-state structures for $`1<n<2`$. In section III, we consider the ground-state structures for fillings greater than two. In these calculations, we use simple geometric arguments to arrive at the ground states. This analysis is performed under the assumption that the range of the intervortex interaction, which is set by the penetration depth, is much smaller than the spacing between the pinning sites. We show that the ground-state structures obtained from this simple analysis match the ones obtained from simulated annealing. This analysis is extended to rectangular and triangular pin lattices in section IV. In section V, the zero-temperature equilibrium magnetization of a vortex system in a square array of pinning centers is obtained by first calculating the ground-state energy as a function of the magnetic induction and then finding the applied field from a numerical differentiation of the data. The ground-state energies for different values of the magnetic induction are obtained using a simulated annealing procedure. The calculated $`BH`$ curve exhibits interesting commensurability effects, manifested as plateaus occurring at simple rational values of the filling fraction $`n`$. The main results of our study are summarized in section VI.

## II Ground states for a square pin array with filling ratio less than one

We consider a superconducting film that has a square array of pinning sites with lattice constant $`d`$. The magnetic field is assumed to be perpendicular to the surface of the film. The “matching field” $`B_\mathrm{\Phi }`$ is then given by $`B_\varphi \mathrm{\Phi }_0/d^2`$, and the filling fraction $`n`$ is given by $`n=B/B_\varphi =Bd^2/\mathrm{\Phi }_0`$ where $`B`$ is the magnetic induction. We assume that the pinning potential is much stronger than the intervortex interaction, but is of extremely short range. The large strength of the pinning potential implies that the vortices must occupy pinning sites as long as the number of vortices does not exceed the number of pinning sites. We also assume that a pinning site can not accommodate more than one vortex. If the pinning centers in the film are microholes, then this assumption amounts to the requirement that the radius of each hole is close to two times the coherence length $`\xi `$ of the superconductor. These assumptions ensure that interstitial vortices appear only when the filling fraction $`n`$ is greater than unity. The assumed short range of the pinning potential can be justified if the defect diameter is small compared to the defect spacing $`d`$. Another assumption that we will make in most of our calculations is that the intervortex interaction falls off rapidly with distance. This is ensured if the penetration depth $`\lambda `$ is much smaller that the pin-lattice spacing $`d`$. In our calculations, we take the ratio $`\lambda /d`$ to be $`10`$. This value is appropriate for the pin lattice of Ref.. We consider temperatures that are low enough to neglect effects of depinning and vortex-lattice melting. The problem of finding the structure of the vortex system then reduces to locating the ground state in the presence of the pinning potential.

Consider now fillings of the form $`n=1/q`$, $`q`$ being an integer greater than one. Let us look at Bravais lattices that can be formed for a specific $`n`$ by distributing the vortices on the square pin lattice with one vortex per basis. The motivation for considering such lattices is that this will automatically ensure that there is no shear of the vortices with respect to the pin lattice, since the forces on a vortex due to other vortices will add up exactly to zero. The unit cell area of these structures has to be $`q.d^2`$. So the possible unit cells can be obtained by factorizing $`q`$ into products of the form $`r.s`$ ($`r`$ and $`s`$ are integers), arranging the vortices at the corners of rectangles of dimension $`rd\times sd`$, and then sliding the parallel sides relative to each other. This procedure produces a large number of structures depending on the value of $`n`$ and we have to pick the one that minimizes the energy. For small values of $`q`$, this can be done by hand, but as $`q`$ becomes large and highly factorizable, the number of possible structures increases rapidly. In such cases, we have resorted to computers to generate these structures and compare their energies.

The structures so obtained for fillings $`1/2`$ and $`1/4`$ match those found in the imaging experiment . Also for fillings $`1/2,1/3,1/4,1/5,1/8,1/9,1/10`$ and $`1/15`$, we find the same structures as those obtained by solving the “greedy lattice gas model”  exactly. This is understandable because when the intervortex interaction falls off rapidly as the distance is increased above the defect spacing $`d`$, the ground state can be attained by finding the lattice that maximizes the shortest distance between vortex pairs. If two structures have same value and number of shortest distances, then the next shortest distance should be maximized, and so on. For fillings $`1/2,1/3,1/4`$, and $`1/5`$, our analysis also yields the same structures as those found in the large $`U_e`$ ($`U_e`$ is the energy of on-site repulsion between two electrons) limit of the neutral Falicov-Kimball model . In Fig. 1(a,b,c), we have shown the structures so obtained for a few fillings of the form $`n=1/q`$. The ground-state structure shown in Fig. 1(b) for $`n=1/5`$ is different from that found in Ref. from a simulated annealing calculation. This difference is probably due to the use of a different (logarithmic) intervortex potential in Ref..

Let us now compare the energies of the nearest and next-nearest neighbors in one of these lattices. The interaction energy between two vortices separated by a distance $`r`$ is given by the expression,

$`U(r)={\displaystyle \frac{\mathrm{\Phi }_0^2t}{8\pi ^2\lambda ^2}}K_0\left({\displaystyle \frac{r}{\lambda }}\right)`$ (1)
where $`K_0`$ is the zeroth order Hankel function of imaginary argument and $`t`$ is the film thickness. For $`n=1/2`$, the nearest neighbor distance is $`\sqrt{2}d`$ and the next-nearest neighbor distance is $`2d`$. So the interaction energies are, to within a constant prefactor, given by

$`U_nK_0\left({\displaystyle \frac{\sqrt{2}d}{\lambda }}\right)=0.2\times 10^6,`$ (2)
$`U_{nn}K_0\left({\displaystyle \frac{2d}{\lambda }}\right)=0.6\times 10^{10}.`$ (3)
One can see here that there is orders of magnitude difference in these energies which cannot be compensated by differences arising from interactions with more distant neighbors. This difference is going to be more prominent at lower densities. This tells us that the maximization of the shortest intervortex distance

in a lattice for a given filling would lead to the ground states, provided the lattice spacing is large compared to the penetration depth of the film. This, in fact, is exactly the definition of the greedy lattice gas. However, one has to be cautious about this method because, as noted in Ref. , the structures can be strongly dependent on the form of the potential in certain ranges of $`n`$ and we can even have aperiodic structures as ground states. The ground-state structures shown here have been cross checked with simulations to ensure that they are indeed the lowest energy configurations. To give an example of a case where this treatment does not lead to the true ground state, we found that for filling $`1/16`$, the energy per vortex for the structure with minimum energy obtained this way was greater than that for filling $`1/15`$, implying that the structure obtained for $`n=1/16`$ was not the ground state.

When the filling fraction is of the form $`p/q`$ with $`p`$ not equal to one, one can look for ground states in a subset of structures where the unit cell has size $`qd^2`$ with $`p`$ vortices in a basis. We have shown in Fig. 1(d,e,f) some of the ground-state structures obtained this way. These structures match those obtained from our simulated annealing calculation. These ground states show the “stripe” structure predicted by Watson  and Kennedy  in appropriate density ranges.

## III Ground states for a square pin array with filling ratio greater than two

If $`n`$ is greater that $`1`$ but smaller than $`2`$, then the ground-state structures are similar to the ones for the case of $`n`$ less than one. The only difference is that the pinning sites are all occupied and the centers of the square unit cells of the pin lattice now act as new pinning centers for interstitial vortices. But things look different when the filling goes above the value $`n=2`$. For such values of $`n`$, we can no longer place the interstitial vortices at the centers of the squares and look for simple structures obtained this way. Also, we now have to start looking into the stability of the structures since the square symmetry would not be present.

### A The ground state for $`n=5/2`$

Here we are faced with the task of placing more than one vortex in a square. Before going to the problem of finding the ground state for $`n=5/2`$, let us ask a more basic question: given a single unit cell of the square pin lattice with each corner occupied by a vortex. how can we arrange two more vortices inside this square so as to minimize the energy. Since “greedy lattice gas” has been a good approximation for the previous cases, we try to tackle this problem by “maximizing the shortest distance” method. In order to stabilize an interstitial vortex by maximizing the shortest distance, its distance from at least three nearest vortices must be the shortest distance. It is also required that these nearest set of vortices must be spread in such a way that if we draw straight lines from the vortex in question to these neighbors, the angles formed by adjacent lines must be less than $`180^o`$. The proof of this statement is given in Appendix A.

It can be seen from the symmetry of the problem that we have to place the two vortices on the lines joining the centers of the sides to meet the condition mentioned above. This leaves us with only two possible ways of doing it, which are shown in Fig. 2.

In (a) the shortest distance $`D_{ma}`$ can be obtained by solving the equation

$`{\displaystyle \frac{1}{4}}(dD_{ma})^2+{\displaystyle \frac{d^2}{4}}=D_{ma}^2.`$ (4)
On solving this equation, we get $`AB=BC=BE=D_{ma}=(\sqrt{7}1)d/3`$. In (b) the vortices $`A^{},B^{}`$ and $`C^{}`$ form an equilateral triangle. Thus the nearest neighbor distance $`A^{}B^{}=B^{}C^{}=B^{}E^{}`$ is

$`D_{mb}=\mathrm{sec}(15^o){\displaystyle \frac{d}{2}}`$ (5)
The angle $`\varphi `$ in (a) is $`24.49^o`$, and the angle $`\theta `$ in (b) is $`15^o`$. The interaction energies corresponding to these two distances for $`d/\lambda =10`$ are

$`U_{AB}K_0\left({\displaystyle \frac{\sqrt{7}1}{3\lambda }}d\right)2.2\times 10^3`$ (6)
$`U_{A^{}B^{}}K_0\left({\displaystyle \frac{\mathrm{sec}(15^o)}{2\lambda }}d\right)3.0\times 10^3.`$ (7)
From comparing these two energies it is clear that (a) is the global minimum, whereas configuration (b) can at best be a local minimum.

Coming back to the $`n=5/2`$ case, we now have to build up the lattice with equal number of two types of squares - one with two interstitial vortices and the other with one interstitial vortex. Note that here we have neglected structures that have three or more interstitial vortices inside a square unit cell because such structures would drastically bring down the nearest neighbor distance. Let us now look at the possible units cells of size $`2d\times 2d`$ that can be made out of these two types of squares. These are shown in Fig. 3.

If one constructs the lattice with these unit cells, configuration (b) offers the least number of next-nearest neighbors, the number of nearest neighbors being the same in all the cases. So one can expect (b) to be the ground state unit cell, at least for large values of $`d/\lambda `$. However, unit cell (a) is preferred if $`d/\lambda `$ is not very large. This can be understood in the following way: the advantage that (b) has over (a) is that it has only half the number of next-nearest neighbors (interactions like that between vortices $`m`$ and $`n`$) compared to (a). But this is done at the cost of bringing in interactions like those between vortex pairs ($`p,q`$) and ($`p,t`$) for every “gain” of a next-nearest neighbor interaction. The energies of these two interactions for $`d/\lambda =10`$ are found to be quite close. These energies are

$`U_{mn}K_0(r_{mn}/\lambda )3.2\times 10^4,`$ (8)
$`U_{pq}K_0(r_{pq}/\lambda )1.9\times 10^4,`$ (9)
where $`r_{mn}`$ and $`r_{pq}`$ are the distances between vortices $`m`$ and $`n`$, and $`p`$ and $`q`$ in Fig.3, respectively. It is clear from this comparison that the unit cell (a) would be preferred for $`d/\lambda =10`$.

The ratio of interaction energies of the next-nearest and the nearest neighbors is 0.07 for this lattice when $`d/\lambda =10`$. This energy difference is appreciable here, so that we can expect that the unit cell we arrived at is the correct one. Note that any net force that might be present on one of the interstitial vortices due to the asymmetry in the structure can be compensated by extremely small displacements from the positions obtained from the “maximization of the shortest distance” method.

### B The ground state for $`n=3`$

When the filling fraction equals $`3`$, we have to build up the lattice using blocks of type (a) in Fig 2. Again, looking at unit cells of size $`2d\times 2d`$ or smaller, we have the configurations shown in Fig. 4 to consider.

Here it is easy to see that unit cell (b) is preferred over the others. This is because it is the configuration that maximizes the minimum distance between any two vortices in different squares, the distances between vortices within one square being the same in all the configurations. Again comparing the nearest interaction and the next-nearest one, we have

$`U_nK_0\left(r_{sp}/\lambda \right)2.2\times 10^3,`$ (10)
$`U_{nn}K_0\left(r_{pq}/\lambda \right)1.9\times 10^4.`$ (11)
There is appreciable difference between these two values, and hence, the ground state we have arrived at is reasonable. When the filling lies between $`2`$ and $`3`$, one can safely assume that the ground state structure can be built up using squares of type (a) in Fig. 2, and squares that have one vortex at the center. In fact we make use of this in our simulations to arrive at the ground states, as described in section V.

### C The ground state for $`n=4`$

Here we have to place three interstitial vortices in one square. This is nontrivial since even if we ensure that the shortest distance is maximized in one square, two vortices in nearby squares may be closer to each other than the shortest distance within a square when we create the lattice. Note that we did not come across this problem in the $`n=5/2`$ or $`n=3`$ fillings. To illustrate this problem, we show in Fig. 5(a) a single cell of a square lattice with three interstitial vortices, which may very well be a local minimum configuration.

But if we try to build the lattice using this cell, we can not do it without bringing the vortices in nearby squares closer than the minimum distance in an individual cell. So what we need to do is to look for a pattern that will include vortices in different squares while doing the minimization of the shortest distance. We can solve this problem trigonometrically. Consider the figure shown in Fig. 6. The solution we are looking for can be obtained by solving the equations

$`AB=BC=CA=D_s,`$ (12)
$`CB^{}=AC^{\prime \prime }=AP=D_s.`$ (13)
Note that here we have assumed a unit cell size $`d\times d`$. On solving these equations, we obtain the unit cell shown in Fig. 5(b). In the figure, the angle $`\theta `$ equals $`15^o`$ and the angle $`\phi `$ is $`60^o`$. In this lattice, the shortest distance is $`D_s=\mathrm{sec}(15^o)d/2`$ and the next shortest distance is $`D_{ns}=3d\mathrm{sec}(15^o)/(4\sqrt{2})`$.

This simple solution may not be the correct one if the lattice spacing is not large enough. For example, if $`d/\lambda =10`$, as we have been assuming when comparing energies, then the nearest and next-nearest neighbor interaction energies turn out to be really close. Hence we can not rule out the possibility of the lattice arranging in such a way that the shortest distance is reduced so as to decrease the number of nearest or next-nearest neighbors. The relevant energies for $`d/\lambda =10`$ are

$`K_0\left({\displaystyle \frac{D_s}{\lambda }}\right)3.1\times 10^3,`$ (14)
$`K_0\left({\displaystyle \frac{D_{ns}}{\lambda }}\right)2.1\times 10^3.`$ (15)
So the ground state obtained above is guaranteed to be the correct one only for much larger values of $`d/\lambda `$.

The ground states we have obtained for $`n=5/2`$ and $`n=3`$ match with the images from experiments  and also the results of simulations . But for $`n=4`$, the structures found in experiments and simulations are different from the one shown in Figs 5(b) and 6. This is expected, since in the experiment the value of the ratio $`d/\lambda `$ was close to ten. Our simulation with $`d/\lambda =10`$ gives the ground-state structure show in Fig. 7 (a), which is similar to the one obtained in the experiment . In our simulations with very large values of $`d/\lambda `$, we get structures similar to that in Fig. 6. The simulation result for $`d/\lambda =50`$ is shown in Fig. 7 (b), which matches well with the predicted structure. It is to be noted that the simulation result was obtained by starting the system near the expected ground state. So the claim is that it offers at least a local minimum of the interaction energy. The simulations were carried out for different system sizes form $`2d\times 2d`$ to $`10d\times 10d`$ to rule out any dependence of the results on the boundary condition.

## IV Ground states for rectangular and triangular pin arrays

One can extend this type of analysis to pinning arrays with other symmetries for finding the least energy structures for simple filling fractions. Let us first consider the case of a rectangular array of pinning sites with a pinning unit cell of dimensions $`l\times b`$, where we take $`l`$ to be the longer side of the rectangle. We shall consider here only the cases where the filling is greater than one. In the absence of the square symmetry, it is obvious that the ground state structure will depend not only on the penetration depth $`\lambda `$, but also on the ratio $`l/b`$. In the following analysis, we shall always assume that $`\lambda `$ is much smaller than $`b`$, the shorter side of the basic rectangular pinning cell. When the filling is $`2`$, for values of $`l/b`$ less than $`\sqrt{3}`$, the ground state is one where each interstitial vortex is at the center of the rectangle, since this ensures that the shortest distance is maximized (see Fig. 8(b)).

But when the aspect ratio exceeds $`\sqrt{3}`$ and the interstitial vortices are placed at the centers of the rectangles, the distance between two interstitial vortices in neighboring cells would be shorter than that between an interstitial vortex and the closest pinned vortex. This would lead to a displacement of the interstitial vortices sideways from the center, along the bisector of the shorter sides of the rectangle, to maximize the shortest distance. The resulting structure is shown in Fig. 8 (a). The displacement of the vortex from the center is given by

$`D_s={\displaystyle \frac{l+\sqrt{4l^29b^2}}{6}}.`$ (16)
It is worth noting that since the vortex in the center would be moving towards two of the pinned vortices, and away from only one interstitial vortex per unit cell, the displacement will approach the value given above only when the ratio $`l/b`$ is appreciably large and in the limit of small penetration depth compared to the sides of the rectangle. For example we have found in our simulations that even when the ratio $`l/b`$ is $`2`$, the ground state for $`b/\lambda =15`$ is one where the interstitial vortex is very close to the center, whereas for the $`l/b=3`$ and $`b/\lambda =15`$, the ground-state structure is quite close to the one obtained from maximizing the shortest distance. Some of our simulation results for filling equal to two are shown in Fig. 12 (a,b). Also, if the ratio $`l/b`$ becomes too large, the analysis will have to include more than two of the interstitial vortices, since now the solution like that shown in Fig. 8 (a) can lead to two vortices being closer in the next nearest cells or ones even further apart.

In trying to arrive at the lowest energy structures for fillings $`5/2`$ and $`3`$, it is important to determine how one can accommodate two vortices in a rectangular cell with the shortest distance being maximized. There are two possible minima that one has to consider: one where the vortices are arranged along the line dividing the shorter sides, and one where they are arranged along the line dividing the longer sides, as shown in Figs 9 (a) and 9 (b), respectively. The shortest distance in each case is given by,

$`D_{s1}={\displaystyle \frac{\sqrt{4l^2+3b^2}l}{3}},`$ (17)
$`D_{s2}={\displaystyle \frac{\sqrt{4b^2+3l^2}b}{3}}.`$ (18)
If one considers the distances within the cell, configuration (b) in Fig. 9 gives the lowest energy. But for large values of the ratio $`l/b`$, this configuration is disfavored since it allows the vortices in one rectangle to get close to those in a neighboring one. Also, for $`l/b>2`$, the interstitial vortices “spill over” into the next cell, since the distance $`D_{s2}`$ becomes greater than $`b`$. So, one has to work out the structures for fillings $`5/2`$ and $`3`$ case by case. For fillings between $`2`$ and $`3`$, one has to choose appropriate number of two-vortex rectangles of the right kind and single-vortex ones and arrange them so as to maximize the shortest distance appearing in the structure. We have looked at the $`2\times 2`$ unit cells possible for filling $`5/2`$ for two values of the aspect ratio, $`l/b=4`$ and $`l/b=5/4`$. The unit cells that provide the largest minimum distance are shown in Fig. 10. Note that when

$`l/b=4`$, the vortices are arranged parallel to the longer side and in the other case, parallel to the shorter side. Also when $`l/b=5/4`$, the vortex in the single-vortex rectangular cell is not located at the center, but slightly displaced sideways along the bisector of the shorter sides to facilitate the maximization of the second shortest distance involved. One should again keep in mind that this sort of analysis will go wrong if the aspect ratio is too large, since then the distances between vortices in next-nearest or further neighbor cells will become important. In Fig. 12 (c) we have shown the ground state structure obtained from simulations for aspect ratio $`4`$, with $`b/\lambda =5`$ and in Fig. 12 (d), the ground state for aspect ratio $`5/4`$, with $`b/\lambda =24`$. In Fig. 12 (d), one can notice the slight shift of the vortex from the center in rectangular units having a single interstitial vortex. The shift is not as large as expected from the above analysis. This is due to the large but finite value of $`b/\lambda `$.

For filling equal to $`3`$, the lowest energy structures obtained by considering $`2\times 2`$ unit cells for two values of the aspect ratio, $`2`$ and $`5/4`$, are shown in Fig. 11. Here too, for large values of the aspect ratio, the structure is composed of rectangular cells where the interstitial vortices are aligned parallel to the longer side (Fig. 11(a)), whereas when the aspect ratio is smaller, the structure is made up of an alternating arrangement of rectangular cells of both types (Fig. 11(b)). The simulated annealing results for similar values of the aspect ratio (Fig. 12 (e,f)) yield the structures obtained from the above analysis.

Ground state structures obtained by simulated annealing for a rectangular pin array with $`l/b=2`$ and integral values of $`n`$ are reported in Ref.. In that study, the intervortex interaction was assumed to depend logarithmically on the intervortex distance. The ground state structure found in Ref. for $`n=2`$ is similar to that shown in Fig. 12(b), but the structure found there for $`n=3`$ is quite different from those of Fig. 12(e,f). This is another example of the importance of the detailed nature of the intervortex interaction in determining the structure of the ground state.

For a triangular array of pinning sites, it is easy to see that when the filling is greater than $`1`$ and less than 3, the interstitial vortices will be placed on the centroids of the triangles in the limit where one can safely apply the method of maximization of the shortest distance. So the ground states when $`n`$ is between $`1`$ and $`3`$ will be made up of parallelogram cells of the form shown in Fig. 13. These unit cells match well with the results of molecular dynamics simulations .

## V Equilibrium magnetization of the ground states

In this section, we describe a calculation of the zero-temperature, equilibrium magnetization of a thin-film superconductor in the presence of a square array of pinning sites. The region in the $`BH`$ plain we are interested in is that just above $`H_{c1}`$, when the flux tubes start entering the sample. The idea is to find the free energy $`F`$ of the ground state as a function of the magnetic induction $`B`$, and then to obtain the applied magnetic field $`H`$ by taking a derivative of the free energy with respect to $`B`$. Since we are considering the zero-temperature case, the free energy is just the internal energy of the flux lattice. Since we are looking for a nearly continuous variation of the internal energy for taking the derivative, we need to locate the ground states for filling fractions separated by small intervals. This would be difficult to do analytically, since the unit cells for some filling fractions can be arbitrarily large. Also, as $`n`$ becomes large, the simple procedure of maximization of the shortest distance is not going to yield the correct ground-state structures. So we have resorted to simulations to determine the ground states. In particular, we have used the simulated annealing technique to locate the global minima (or at least low-lying local minima close in energy to the global ones) of the part of the internal energy associated with intervortex interactions.

The Helmholtz free energy per unit volume of the superconductor at zero temperature in the presence of the pinning sites is

$`F_s(n)={\displaystyle \frac{nϵ_l}{d^2}}+E_n{\displaystyle \frac{nϵ_p}{d^2}}`$ (19)
where the first term is the line energy, the second term is the interaction energy and the third term is the pinning energy. Here, $`ϵ_l`$ is the line energy per unit length, $`ϵ_p`$ is the pinning energy per unit length and $`E_n`$ is the interaction energy per unit volume for filling fraction $`n`$. We note here that the pinning energy increases linearly with $`n`$ till $`n`$ becomes one and then remains constant, since multiple occupation of a pinning center is not allowed. Further, for simplicity, we express the pinning energy as

$`ϵ_p=mϵ_l`$ (20)
where $`m`$ is a positive number whose magnitude depends on the nature of pinning. The interaction part of the free energy, $`E_n`$, is the computational input. Once we know the free energy, we can compute the applied magnetic field $`H`$ using the relation

$`{\displaystyle \frac{F_s}{B}}={\displaystyle \frac{H}{4\pi }}.`$ (21)
Using the standard expression  for $`ϵ_l`$ and taking the logarithm of the Ginzburg-Landau number to be equal to 2, we get the following expression for the applied field as a function of the filling fraction:

$`H={\displaystyle \frac{\mathrm{\Phi }_0}{2\pi \lambda ^2}}[1m\{1\mathrm{\Theta }(n1)\}]+{\displaystyle \frac{E_n^{}}{n}}.`$ (22)
Here $`E_n^{}`$ is given by the expression

$`E_n^{}={\displaystyle \frac{\mathrm{\Phi }_0}{2\pi \lambda ^2N}}{\displaystyle \underset{i>j}{}}K_0(r_{ij}/\lambda ),`$ (23)
where $`N`$ is the number of basic pinning squares in the system and $`r_{ij}`$ is the separation between vortices $`i`$ and $`j`$ in the ground state for the filling fraction $`n`$.

The size of the systems we simulated varies from $`2d\times 2d`$ to $`8d\times 8d`$. In all cases, we used periodic boundary conditions to minimize surface effects. So the minimum difference between two consecutive filling fractions was $`\mathrm{\Delta }n=1/64`$. The ratio $`d/\lambda `$ was taken to be $`10`$, as in our previous analysis. In order to save computation time, the vortices were allowed to stay only at the pinning sites when the filling was less than one. For fillings between one and two, every pinning site was occupied by a vortex which was never moved and the extra ones were allowed to move near the centers of the basic pinning squares. When the filling was greater that two and less than three, the vortex configurations were constructed using basic units of squares containing one vortex at its center and squares containing two vortices placed such that the shortest distance within a square is maximized (as in Fig. 2(a)). These units were then moved around and twisted while cooling to arrive at the minimum energy states. This procedure helped us to track low-lying minima faster than if we allowed vortices to move freely. Once the basic structure was thus obtained, the vortices were allowed to move freely during a second cooling schedule starting from a lower temperature to obtain the lowest-energy structure. In Fig. 14 we have shown some of the ground state structures we have obtained this way for fillings between $`2`$ and $`3`$. For fillings $`5/2`$ and $`3`$, we find that the structures match those obtained in experiments , as well as in our analysis using maximization of the shortest distance. The structures for $`n=9/4`$ and $`n=11/4`$ may not be the actual ground states, either due to the smallness of the unit cell of our simulation or due to the presence of many nearly degenerate local minima.

In Fig. 15 we have plotted the ground state energies obtained from the simulation for different fillings. The simulation unit cell was $`8d\times 8d`$ and the energies were computed for fillings $`1/64`$ to $`3`$. The upper curve shows the results obtained in the presence of the pinning sites and the lower curve is the energy of the triangular lattice for the same density of vortices. Note that we have not included the pinning energy in the plot. This would bring down the upper curve below the curve for the pin-free case.

From the energy versus filling fraction data, one can find the applied field using Eq.(11) and then compute the magnetization $`M`$ using the relation

$`B=H+4\pi M`$ (24)
In Fig. 16 we have plotted $`B`$ versus $`H`$ in the entire range of filling for which simulations were carried out, from $`n=0`$ to $`n=3`$. Figs 17 and 18 show magnified versions of this plot in the regions between filling fractions $`0`$ and $`1`$, and between $`2`$ and $`3`$, respectively. Note that we have not explicitly included the pinning energy term in our analysis. This term would just add a constant contribution to $`H`$ for fillings up to $`1`$. The features of the curve from $`n=1`$ to $`n=2`$ are the same as those in the interval between $`n=0`$ to $`n=1`$. This is due to the fact that the ground state structures are similar in the two regions (see section III). The $`BH`$ plot shows flat regions at values of $`B`$ corresponding to fillings $`1/8`$, $`1/5`$, $`1/4`$, $`1/2`$, $`3/4`$, $`4/5`$, $`7/8`$, $`1`$, $`9/8`$, $`6/5`$, $`5/4`$, $`3/2`$, $`7/4`$, $`9/5`$, $`2`$ in the filling fraction range between $`0`$ and $`2`$. Also, in the range of $`n`$ between $`2`$ and $`3`$, there are roughly two plateaus, appearing near $`n=2.3`$ and $`n=2.6`$.

The observed values of the filling fractions between $`n=0`$ and $`n=1`$ at which the plateaus occur indicate that these values of $`n`$ correspond to fillings where the introduction of a new vortex into the system leads to the appearance of a shorter distance than those existing in the lattice or in its dual (i.e. the lattice obtained by replacing particles by holes, and vice versa) . This makes sense because the introduction of this shorter distance brings in a larger energy scale, leading to a discontinuous change in the derivative of the energy with respect to the filling $`n`$. In the simulations, we have not scanned very small intervals of $`n`$. Also, for some of the fillings, the ground state may not have been obtained in our simulated annealing calculation. For these reasons, we can not say anything definite about the true nature of the $`BH`$ curve. It is possible that this curve has plateaus and discontinuities occurring at all scales (e.g. at all rational values of $`n`$). The noisy nature of the $`BH`$ plot in the range where $`n`$ lies between $`2`$ and $`3`$ (see Fig. 18) is also due to these difficulties. However, this plot shows clear signatures of two plateaus that appear near $`n=2.3`$ and $`n=2.6`$. These can be understood as happening when first there is an occurrence of squares containing two vortices coming next to each other diagonally (as in Fig 3(a)), and again when they have to be next to each other with a common side (as in Figs 3(b) and 3(c)), as the value of $`n`$ increases from $`2`$ to $`3`$.

## VI Summary and discussions

In this paper, we have reported the results of an analytic calculation of the lowest-energy states of a vortex system in the presence of a periodic array of strong pinning centers, each of which can trap only one vortex. We have considered several different lattice structures of the pin array and a large number of filling fractions in the range between zero and four. The analytic calculations are based on the principle of maximization of the shortest intervortex distance. We have argued that this principle leads to the exact ground states when the spacing of the defect lattice is large compared to the range of the intervortex interaction set by the value of the penetration depth. This principle has been used, in combination with simple geometric considerations, to obtain the ground states for several values of the filling fraction $`n`$. The ground-state structures so obtained are found to be identical to those found in imaging experiments  and in earlier simulations . We have also carried out simulated annealing calculations of the ground states in order to test some of the predictions of the analytic approach. In all cases, we found that the analytic results agree with those of our numerical calculation.

We have also described the results of a numerical calculation of the equilibrium magnetic induction $`B`$ and magnetization $`M`$ of a planar superconductor with a square array of pinning centers as functions of the externally applied field $`H`$. We show that the interplay between the lattice spacing of the pin array and the intervortex separation set by the value of $`B`$ leads to interesting commensurability effects, showing up as plateaus and discontinuities in the $`B`$ vs $`H`$ plot at simple rational values of the filling fraction $`n`$. Anomalies in the irreversible magnetization of thin-film superconductors with periodic arrays of pinning centers have been observed at certain integral values of $`n`$ in experiments  and simulations . The presence of a periodic array of pins is also expected  to produce anomalies in the equilibrium magnetization of the high-temperature vortex liquid at small integral values of $`n`$. Our results show that these commensurability effects are more pronounced in the field-dependence of the equilibrium magnetization and magnetic induction of such systems in the low-temperature vortex-solid regime. Experimental investigations of these effects would be most welcome.

## A Condition for maximizing the shortest distance

Given a particle fixed at some point in a plane, the problem is how to place three other particles around the first one in such a way that if we try to move the first particle from its position, it will get closer to at least one of the three particles. The solution is as follows. If we draw straight lines from the particle that we want to move to the other particles, then each angle between adjacent lines must be less than $`180^o`$. In other words, it should not be possible to draw a straight line through the particle in question in such a way that all the other three particles lie on one side of the line.

Let the particle that we want to move be at $`P_0`$ (see Fig. 19). Let us place two particles at $`P_1`$ and $`P_2`$, anywhere on the plane.

This can be done since one of the angles between any two lines will always be less than or equal to $`180^o`$. In Fig. 19, $`AE`$ and $`CF`$ are the tangents to the circles centered at $`P_1`$ and $`P_2`$ and passing through the point $`P_0`$. The presence of these particles restricts to within the angle $`\xi `$ the direction in which the particle at $`P_0`$ can move without decreasing its distance to $`P_1`$ and $`P_2`$. Now let us check where we can fix the third particle such that the aforementioned condition is met, that is the particle at $`P_0`$ cannot be moved without bringing it closer to one of the particles at $`P_1`$, $`P_2`$ and $`P_3`$. We will consider two possible cases - the third particle in the region $`BP_0D`$ of the plane and outside this region.

When the third particle is in the region $`BP_0D`$ (point $`P_3`$ in the figure), it is not possible for the particle $`P_0`$ to move from its position without decreasing the distance to any one of the particles. This can be seen if one drops the perpendiculars to the lines $`AE`$ and $`CF`$ from the point $`P_3`$. If the point $`P_0`$ is moved along either of these lines, it will be getting closer to the point $`P_3`$. But that in turn implies that it cannot be moved into the region $`FP_0E`$ at all without decreasing the distance from the point $`P_3`$.

If the particle is outside the region $`BP_0D`$, for example a point like $`P_3^{}`$, then it is easy to see that by moving along one of the tangents to the circles at the point $`P_0`$, the particle at $`P_0`$ can move away from all the three other points. Thus we find that the earlier statement we made is proved. This gives us a nice way to maximize the shortest distance since all one has to do is to make the three distances involved equal, so that if the central particles tries to move, then one of the distances has to decrease.