text:
# Bijective Proofs for “Enumerative Properties of Ferrers Graphs”

## 1 Introduction

In their paper , Richard Ehrenborg and Stephanie van Willigenburg introduce the concept of a “Ferrers graph”. Given a Ferrers diagram,111Recall that a Ferrers diagram is a set of lattice points in the upper-right quadrant $`\{(x,y):x,y\}`$ such that, if $`(x,y)`$ is in the set, so is $`(x^{},y^{})`$ for all $`x^{}`$ at most $`x`$ and $`y^{}`$ at most $`y`$. (We normally rotate the diagram by $`\frac{\pi }{2}`$ and draw boxes for the lattice points.) we construct its corresponding Ferrers graph by taking the vertices to be the rows and columns of the diagram, and the edges to be the squares of the diagram: row $`a`$ and column $`b`$ are connected by an edge if and only if there is a square in position $`(a,b)`$ of the diagram.

This bipartite graph is clearly equivalent to the original diagram, but in this new form we can raise new questions. Ehrenborg and Van Willigenburg examine several of these, including the number of spanning trees and Hamiltonian paths of the Ferrers graph; alas, the proofs they give for these are not bijective. In this paper, I give bijective proofs for both of these, which extend to weighted versions.

Ehrenborg and Jeffrey B. Remmel have independently found a bijection for the number of spanning trees; I have not yet seen it.

## 2 Hamiltonian paths

Let’s suppose our Ferrers diagram has the same number $`n`$ of rows and columns.222This is not much of a restriction. Since any Hamiltonian path in the Ferrers graph must alternate between rows and columns, there can be at most one more (or fewer) row than column, if there are to be any Hamiltonian paths at all. Ehrenborg and Van Willigenburg proved that the number of Hamiltonian paths for the corresponding Ferrers graph is the square of the number of $`n`$–rook placements on the Ferrers diagram. (By “number of $`n`$–rook placements”, we mean the number of ways to choose $`n`$ squares from the diagram so that no two are in the same row or column.) To prove this bijectively, we seek a procedure converting any Hamiltonian path to a pair of $`n`$–rook placements, and an inverse procedure converting any pair of $`n`$-rook placements back to a Hamiltonian path.

*Hamiltonian paths to pairs of $`n`$–rook placements*

One end of the Hamiltonian path is a row-vertex, and the other is a column-vertex. Suppose the vertices have labels $`a_1,b_1,a_2,b_2,\mathrm{},a_n,b_n`$, starting at the row end, so the “$`a`$” vertices are rows and the “$`b`$” vertices are columns.

1. Start with an empty Ferrers diagram which we’ll mark up with two distinct $`n`$–rook placements, labeled with A’s and B’s. Also set $`b_0:=1`$.
2. Repeat for each pair of vertices $`a_i,b_i`$:

1. Mark an A in the square $`(a_i,b_{i1})`$, unless there is already an A in column $`b_{i1}`$. If there is an A in that column, mark the A in the next available column without an A.
    2. Mark a B in the square $`(a_i,b_i)`$.
    3. Move on to the next pair of vertices $`a_{i+1},b_{i+1}`$.

We place one B in each row $`a_i`$, and one B in each column $`b_i`$, so the B’s constitute an $`n`$–rook placement. We place one A in each row $`a_i`$, and we never place more than one A in any column, so the A’s also constitute an $`n`$–rook placement. The only thing we need to check is that there always is a “next available column” in row $`a_i`$, in which to place our A.

This seems easy: after all, we’re about to mark a B in row $`a_i`$, column $`b_i`$, which guarantees that column $`b_i`$ is possible. Since we mark as many A’s as B’s at every stage, and every column we’ve marked so far with a B (that is, $`b_1`$ through $`b_{i1}`$) is guaranteed to have an A, we also know that $`b_i`$ doesn’t have an A and is thus available. So there is at least one available column.

The problem is that we need a *next* available column, that is, one that follows $`b_{i1}`$. The solution is to observe that the “next available column” is actually the *first* available column: every column above it has an A already. Suppose, at some pair of vertices $`a_i,b_i`$, that the column $`b_{i1}`$ we would have used already has an A. We’ve never been to $`b_i`$ before, so there are only two ways $`b_{i1}`$ could have acquired that A. Either it was put there in the first iteration of the loop, and $`b_{i1}=b_0`$ is the first column (in which case the next column without an A is obviously the first one without an A), or else it was put there as the “next available column” (so $`b_{i1}`$ was at one time the first available column, which means that all the columns before $`b_{i1}`$ contain an A). Hence the available column $`b_i`$ must follow $`b_{i1}`$.

*Pairs of $`n`$–rook placements to Hamiltonian paths*

1. Start with a Ferrers diagram marked with two different $`n`$–rook placements, labeled with A’s and B’s. Also set $`b_0:=1`$.
2. Repeat until we’ve identified each pair of vertices $`a_i,b_i`$:

1. Remove the A in column $`b_{i1}`$, or else the first column with an A, and set $`a_i`$ equal to the row it was in.
    2. Remove the B in row $`a_i`$, and set $`b_i`$ equal to the column it was in.
    3. Move on to identify the next pair of vertices $`a_{i+1},b_{i+1}`$.

This is clearly the inverse of the procedure given earlier. Note that “first column with an A” is equivalent to “next column with an A”, by our earlier argument.

One last comment: if we switch the roles of rows and columns in the above procedures, or if we switch the roles of A’s and B’s, we get another Hamiltonian path. It would be interesting to know if there is any deeper meaning to these two involutions.

## 3 Spanning trees

The number of spanning trees of a Ferrers graph is the product of the lengths of all but the first column, times the product of the lengths of all but the first row. We can prove this with a similar bijection.

We represent this product combinatorially in a similar fashion to the rook placements above: in each row but the first we place an R, and in each row but the first we place a C. (Note that now there may be more than one R in some column, or more than one C in some row; and we needn’t have equal number of rows and columns.) Also, place an X in the first row and first column, for a reason we’ll see later.

Given such a configuration of R’s and C’s, we will convert it to a spanning tree in two stages.

For the first stage, note that every row except the first has exactly one R, but it need not have any C’s; and likewise, each column but the first has one C and as few as zero R’s. If there is a row $`a`$ with an R in column $`b`$, but no other markings, we may construct a spanning tree by removing row $`a`$, constructing the spanning tree on the remaining rows and columns, and then joining row $`a`$ to the rest of the tree with a single edge to column $`b`$. Likewise, if there is a column $`b`$ with a C in column $`a`$, but no other markings, we may construct a spanning tree by removing column $`b`$, constructing the spanning tree on the remaining rows and columns, and then joining column $`b`$ to the rest of the tree with a single edge to row $`a`$. It clearly doesn’t matter in which order we prune these rows and columns. Conversely, given a spanning tree, we can remove a leaf node (say, row-vertex $`a`$), mark the square $`(a,b)`$ corresponding to the edge connecting it with the rest of the tree (we would mark with an R, since this is a row-vertex), and proceed to mark the remaining rows and columns of the spanning tree. Thus far, then, we have a bijection. Note that we cannot remove the first row or column in this stage.

For the second stage, we are thus left with only the irreducible objects of each type: the configurations of R’s and C’s for which every row but the first has at least one C and every column but the first has at least one R, against the paths from the first row-vertex to the first column-vertex. Remove the first row and column: then what we have left are simply pairs of $`n`$–rook placements (in the former case), and arbitrary Hamiltonian paths (in the latter), and we already have a bijection between these! If each row has a C then there must be at least as many C’s as R’s, and if each column has an R then there are as many R’s as C’s, and so there are exactly as many R’s as C’s: one per column, and also one per row. (So the R’s are an $`n`$–rook placement, and so are the C’s.) On the other hand, a Hamiltonian path starting at the first column and ending at the first row is simply an arbitrary Hamiltonian path on the second through $`n`$th rows and columns; just strip off the first-row and first-column endpoints, or put them back on again.

## 4 Weighted spanning trees

Our spanning-tree bijection also gives us a weighted formula for the total number of spanning trees. Specifically, if we assign weight $`x_ay_b`$ to the edge from row $`a`$ to column $`b`$, then the weight of any given spanning tree is equal to the weight of the corresponding configuration of R’s, C’s, and one X in the upper-left corner. This in turn lets us sum the weights of all spanning trees by summing the weights of all configurations, giving the product formula

$$(x_1x_2\mathrm{}x_m)(y_1y_2\mathrm{}y_n)\underset{2am}{}(y_1+\mathrm{}+y_{\lambda _a})\underset{2bn}{}(x_1+\mathrm{}+x_{\lambda _b^{}})$$

where $`\lambda _a`$ and $`\lambda _b^{}`$ are the lengths of row $`a`$ and column $`b`$ respectively.333This is not hard to see. Since each symbol in the configuration can be placed independently of the others, we simply calculate the weights each symbol can give, and multiply them together. The X must be placed in the upper left corner, contributing $`x_1y_1`$; the R in row $`a`$ can be placed in column $`1,2,\mathrm{},\lambda _a`$, for a factor of $`x_a(y_1+\mathrm{}+y_{\lambda _a})`$; and the C in row $`b`$ can be placed in rows $`1,2,\mathrm{},\lambda _b^{}`$, for a factor of $`y_b(x_1+\mathrm{}+x_{\lambda _b^{}})`$.

That our spanning-tree bijection is weight-preserving is evident in the first stage, because we remove an edge of weight $`x_ay_b`$ from both the configuration and the spanning tree; in the second stage it is only slightly less obvious. Each Hamiltonian path in the second stage starts in the first column and ends in the first row, but passes through all other vertices; hence it contributes $`x_1y_1(x_2\mathrm{}x_m)^2(y_2\mathrm{}y_n)^2`$ to the total. Likewise, there is a single X in the first row and the first column, but an R and a C in each other row and column, so these marks also contribute $`x_1y_1(x_2\mathrm{}x_m)^2(y_2\mathrm{}y_n)^2`$.

## 5 Future research

One natural extension444Suggested by Ehrenborg and Van Willenburg in their paper. of the definition of a Ferrers graphs is to skew Ferrers diagrams. In fact, the definition works with no further changes; unfortunately, the results above fail, and no reasonable fix is apparent. Nothing interesting is yet known in this case.

Let me suggest another extension which may prove more fruitful. One of the key invariants of a graph is its Tutte polynomial, which can be specialized to give such quantities as the number of spanning trees and the chromatic polynomial. (Thus, a combinatorial determination of the Tutte polynomial would encompass the result above.555Ehrenborg and Van Willenburg also suggested investigating the Tutte polynomial. ) The Tutte polynomial can be computed recursively by deleting and contracting edges, but if we try to calculate the Tutte polynomial of a Ferrers graph recursively by removing rows, these contractions yield multiple edges. Specifically, we get what we might call “plane-partition graphs”: the number of edges from row $`a`$ to column $`b`$ (corresponding to square $`(a,b)`$) does not increase as we increase $`a`$ or $`b`$ — or, stated another way, if we were to write in square $`(a,b)`$ the number of edges corresponding to that square, we would get a plane partition.


format: text

text:
# Cationic DMPC / DMTAP Lipid Bilayers: Molecular Dynamics Study

## I Introduction

Gene therapy based on the introduction of genetic material into cells is one of the most promising biomedical approaches to treat human diseases Langer, (1998); Mönkkönen and Urtti, (1998); de Gennes, (1999). The majority of the delivery vectors proposed are of viral nature. The viral vectors have been demonstrated to be very efficient but their use is restricted by accompanied toxicity Anderson, (1998). This has stimulated a search for non-viral delivery systems which should be characterized by greater safety and ease of manufacturing Langer, (1998). Numerous examples of non-viral delivery vectors include cationic liposomes, cationic polymers (such as polyamidoamine dendrimers, polyethylenimine, and spermine), and block copolymers Pitard et al., (1997); Mönkkönen and Urtti, (1998); Smedt et al., (2000); KukowskaLatallo et al., (1996); Astafieva et al., (1996); Gao and Huang, (1996); Fischer et al., (1999).

In the light of the above, it is surprising how little attention has been devoted to computational studies of membranes containing cationic lipids. Bandyopadhyay et al. Bandyopadhyay et al., (1999) performed an atomistic molecular dynamics (MD) study of a mixture of dimyristoylphosphatidylcholine (DMPC) and dimyristoyltrimethylammonium propane (DMTAP) in the presence of a short DNA fragment. Apart from the very elegant piece of work above there are, to the best of our knowledge, no published atomistic computational studies of systems containing cationic lipids – this is very much in contrast to the great number of computational studies of various neutral and anionic phospholipid bilayer systems Feller, (2000); Saiz and Klein, (2002); Saiz et al., (2002); Tieleman et al., (1997); Tobias, (2001). Another related example is the recent molecular dynamics study of Böckmann et al. Böckmann et al., (2003) who showed the importance of monovalent ions on the properties and organization of lipid membranes – ions are always present in cationic lipid systems. The above examples demonstrate that detailed molecular dynamics studies can provide valuable insight into the atomistic organization of systems containing cationic lipids and yield useful information for experimentalists about the underlying mechanisms on the atomic and molecular levels.

In the present work, our objective is to gain insight into the structural and electrostatic properties of cationic lipid bilayers through atomic classical molecular dynamics simulations. We concentrate on a bilayer mixture composed of two kinds of lipids: Neutral (zwitterionic) DMPC and cationic DMTAP (see Fig. 1 for their chemical structures). Since DMTAP is positively charged under physiological conditions, we have neutralized its positive charges by chloride counter-ions. From the computational point of view, this choice for a model system is motivated by the fact that DMPC and DMTAP have the same non-polar hydrocarbon chains and differ only by their headgroups. On the practical side, DMPC / DMTAP binary lipid mixtures have been widely studied in the presence of DNA by various experimental techniques Artzner et al., (1998); Zantl et al., 1999a ; Zantl et al., 1999b ; Pohle et al., (2000), and also by computational methods Bandyopadhyay et al., (1999).

We focus mostly on how the composition of the cationic lipid bilayer affects the structural and electrostatic properties of these lipid bilayer systems. To this end, we consider mixtures of DMPC and DMTAP with various different mole fractions of the cationic DMTAP component under conditions corresponding to the liquid-crystalline phase. We find that DMTAP plays a prominent role leading to considerable changes compared to the pure DMPC bilayer. As discussed in the present study, this is characterized by the strong interplay between electrostatics and structural changes in the vicinity of the membrane-water interface. In particular, we find that DMTAP gives rise to a non-monotonic dependence of the average area per lipid on the DMTAP concentration, substantial changes in the electrostatic profile of the membrane, and significant re-orientation of P – N dipole vectors in DMPC headgroups. The spatial re-arrangement of P – N dipoles is particularly interesting as it likely plays a significant role in the stability of DNA-membrane complexes.

## II System

### II.1 Model and simulation details

We have performed atomistic simulations of fully hydrated lipid bilayers consisting of a mixture of cationic DMTAP and zwitterionic DMPC lipids. In all simulations, the total number of lipid molecules was fixed to 128, and the number of water molecules ranged from 3527 (pure DMTAP) to 3655 (pure DMPC).

Force-field parameters for the lipids were taken from the recent united atom force-field Berger et al., (1997). This force-field has been previously validated Tieleman and Berendsen, (1996); Lindahl and Edholm, (2000) and is essentially based on the GROMOS forcefield for lipid head groups, the Ryckaert-Bellman’s potential Ryckaert and Bellemans, (1975, 1978) for hydrocarbon chains, and the OPLS parameters Jorgensen and Tirado-Rives, (1988) for the Lennard-Jones interactions between united CHn groups of acyl chains re-parameterized for long hydrocarbon chains to reproduce experimentally observed values of volume per lipid Nagle and Wiener, (1988). The parameters for this force field are available on-line at http://moose.bio.ucalgary.ca/Downloads/files/lipid.itp. Water was modeled using the SPC water model Berendsen et al., (1981). The unit positive charge carried by each DMTAP molecule is compensated by the introduction of the corresponding number of explicit Cl- counter-ions. While being aware of the effects of different models for chloride Patra and Karttunen, (2004), we decided to use the default set of chloride parameters supplied within the Gromacs force field Berendsen et al., (1995); Lindahl et al., (2001).

Following the original parameterization Berger et al., (1997), the Lennard-Jones interactions were cut off at 1 nm without shift or switch function. Since long-range electrostatic interactions are essential in the present study, and since truncation of these interactions has been shown to lead to artifacts in simulations of phospholipids bilayers Patra et al., (2003), we employ the particle-mesh Ewald (PME) method Darden et al., (1993). The long-range contribution to the electrostatics is updated every 10-th time step.

The simulations were performed in the $`NpT`$ ensemble. The temperature was kept constant using a Berendsen thermostat Berendsen et al., (1984) with a coupling time constant of 0.1 ps. Lipid molecules and water (including counter-ions) were separately coupled to a heat bath. Pressure was controlled by a Berendsen barostat Berendsen et al., (1984) with a coupling time constant of 1.0 ps. Pressure coupling was applied semi-isotropically: The extension of the simulation box in the $`z`$ direction (i. e., in the direction of the bilayer normal) and the cross-sectional area of the box in the $`x`$-$`y`$ plane were allowed to vary independently of each other.

We considered 11 DMPC / DMTAP mixtures ranging from pure DMPC to pure DMTAP. The molar fractions of the cationic DMTAP, $`\chi _{\text{T}AP}`$, were taken to be 0.0, 0.06, 0.16, 0.25, 0.31, 0.39, 0.50, 0.63, 0.75, 0.89, and 1.0.

The main transition temperature of a pure DMPC bilayer is $`T_m=24^{}`$Cevc and Marsh, (1987). For DMPC / DMTAP binary mixtures it has been found Zantl et al., 1999b  that the main transition temperature changes non-monotonically with the mole fraction of DMTAP, demonstrating a maximum of approximately $`37^{}`$C at $`\chi _{\mathrm{TAP}}0.45`$. All our simulations were done at a temperature of $`50^{}`$C, such that the bilayers are in the liquid-crystalline phase.

All bond lengths of the lipid molecules were constrained to their equilibrium values using the LINCS algorithm Hess et al., (1997) whereas the SETTLE algorithm Miyamoto and Kollman, (1992) was used for water molecules. The time step in all simulations was set to 2 fs. All simulations were performed using the GROMACS package Berendsen et al., (1995); Lindahl et al., (2001). The combined simulated time of all simulations amounts to 250 nanoseconds. Each simulation was run in parallel over 4 processors on an IBM eServer Cluster 1600 system. In total, the simulations took about 20,000 hours of CPU time.

### II.2 Simulation setup

Mixtures of DMPC and DMTAP were prepared and equilibrated in several steps as follows:

1. We used the equilibrated DPPC bilayer structure of Patra et al. Patra et al., (2003) as our initial configuration. The structure is available electronically at http://www.softsimu.org/downloads.shtml.
2. DPPC and DMPC molecules differ only by length of their tail. Thus, we created a DMPC bilayer by shortening the DPPC acyl chains by two hydrocarbons. This procedure does not disturb the acyl chain region or the water-lipid interface.
3. The next step was to compress the DMPC bilayer in order to eliminate the gap between leaflets created in the previous step. For this purpose a pre-equilibration run for 1 ns in the $`NpT`$ ensemble was performed. After this, the gap between leaflets had disappeared, and the obtained DMPC bilayer structure was used as initial configuration for all DMPC / DMTAP mixtures described in the following.
4. The chemical structures of DMPC and DMTAP differ only by their head groups, see Fig. 1. With that in mind, the following procedure was used to prepare mixtures. For each DMTAP concentration the corresponding number of randomly chosen PC headgroups in a pure DMPC bilayer were converted to TAP head groups Bandyopadhyay et al., (1999). To neutralize the unit charges in DMTAP headgroups, randomly chosen water molecules were replaced by chloride ions while ensuring a minimum separation of 0.5 nm between ions. To retain symmetry between the two leaflets, each of them contains the same number of cationic DMTAP.
5. Since TAP headgroups occupy a smaller volume than those of DMPCs, a short 10 ps run in the $`NVT`$ ensemble was performed to let the water molecules adjust at the lipid-water interface. This step completes pre-equilibration.
6. The actual equilibration runs were performed in the $`NpT`$ ensemble. The needed equilibration times prior to actual production runs ranged from 10 ns to 20 ns depending on the DMTAP mole fraction. We concluded that equilibration was completed when the average area per lipid had become stable and fluctuated around its mean with a standard deviation not exceeding the standard deviation for a pure DMPC bilayer.

After equilibration, for each DMPC / DMTAP mixture, we performed a production run of 10 ns in the $`NpT`$ ensemble to collect the data for analysis.

## III Results and Discussion

### III.1 Area per lipid

The average area per lipid, $`A`$, is one of the most fundamental characteristics of lipid bilayers Nagle and Tristram-Nagle, (2000). While being one of the rather few structural quantities that can be measured accurately from model membranes via experiments, its also plays a major role in a number of quantities, including the ordering of acyl chains and the dynamics of lipids in a bilayer. Further, from computational point of view, it is highly useful as a means of monitoring the equilibration process.

Due to the lack of experimental data for the average area per lipid in binary DMPC / DMTAP bilayer mixtures, or in a pure DMTAP bilayer \[apart from the low-temperature phase Lewis et al., (2001)\], reproduction of the experimental data available for pure DMPC membranes is essential to validate our approach. To this end, let us first consider the temporal behavior of the area per lipid, $`A(t)`$, presented in Fig. 2. It shows that the obtained average area per lipid for a pure DMPC system has a value of $`A=0.656\pm 0.008`$ nm2, in very good agreement with the experimentally observed value of 0.654 nm2 Petrache et al., (2000), thereby validating our model in this respect. As for the pure DMTAP bilayer, Lewis et al. Lewis et al., (2001) have been able to extract the area per lipid in the low-temperature phase, finding $`A=0.40`$ nm2 at 25 C. Studies of $`A`$ above the main transition temperature are lacking, however.

We find that the average area per lipid shows a non-monotonic dependence on DMTAP mole fraction $`\chi _{\mathrm{TAP}}`$, with a pronounced minimum roughly at $`\chi _{\mathrm{TAP}}=0.5`$, see Fig. 3. This behavior is not trivial, as modest amounts of the cationic DMTAP lead to a compression of the bilayer, while high concentrations lead to a major expansion of the membrane. More specifically we find that for $`0<\chi _{\mathrm{TAP}}0.8`$, the average area per lipid is smaller than the corresponding counterpart for any of the pure lipid systems. Such a behavior cannot be explained by steric interaction alone but most likely is rather of electrostatic origin.

While no published experimental data exists (to the best of our knowledge) for the area per lipid of DMPC / DMTAP bilayer mixtures, similar behavior has been observed in related systems.

For example, Zantl et al. Zantl et al., 1999b  considered DMPC / DMTAP monolayers using Langmuir-type film balance and found that, for pressures corresponding to the liquid-crystalline phase, the headgroup area decreased monotonically for small $`\chi _{\mathrm{TAP}}`$, then had a minimum at about the equimolar ratio ($`\chi _{\mathrm{TAP}}0.5`$), and increased for larger DMTAP mole fractions. This is in accord with our observations.

Another related work concerns Langmuir balance studies of mixed monolayers of zwitterionic palmitoyloleoylphosphatidylcholine (POPC) and cationic 2,3-dimethoxy-1,4-bis($`N`$-hexadecyl-$`N`$;$`N`$-dimethyl-ammonium)butane dibromide (SS–1) Säily et al., (2001). Although SS–1 is dicationic, one can qualitatively compare this system to our DMPC / DMTAP mixture. For the POPC / SS–1 system Säily et al. Säily et al., (2001) found that the average area per lipid has a non-monotonic behavior with a minimum at $`\chi _{SS1}0.38`$. They also found that this effect depends on the charge of the head group, and it disappeared when POPC (having a zwitterionic head group) was replaced by neutral dioleylglycerol (DOG).

Our results in Fig. 3 suggest local extrema in $`A`$ when $`\chi _{\mathrm{TAP}}`$ is between 0.16 and 0.5 (in addition to global minimum at $`\chi _{\mathrm{TAP}}=0.5`$). In addition to the above-mentioned POPC and SS–1 study Säily et al., (2001), similar and even more dramatic effects have been observed for mixtures of POPC and sphingosine Säily et al., (2003). The existence of critical concentrations in lipid membranes has also been theoretically postulated by Somerharju et al. Virtanen et al., (1998); Somerharju et al., (1999). In the present case, the local extrema for DMPC / DMTAP mixtures are within error bars, and therefore may be interpreted as fluctuations of $`A`$. To study such features in more detail one needs to decrease the fluctuations in the average area per lipid by, e.g., increasing the system size Lindahl and Edholm, (2000). This, however, is beyond the scope of the present study.

Nevertheless, we decided to approach this issue from a different perspective. To determine the average area per lipid separately for the two different components, we used the Voronoi tessellation technique in two dimensions Patra et al., (2003). In Voronoi tessellation, we first calculated the center of mass (CM) positions for the lipids and projected them onto the $`x`$-$`y`$ plane. A point in the plane is then considered to belong to a particular Voronoi cell, if it is closer to the projected CM of the lipid molecule associated with that cell than to any other CM position. As there is no unique definition for the area per molecule in a multi-component system, it is clear that the Voronoi results should be considered as suggestive rather than quantitative, providing insight mainly of the trends.

Figure 4 demonstrates that the areas occupied by DMPC and DMTAP are distinctly different. For small $`\chi _{\mathrm{TAP}}`$, the area per DMPC is considerably larger than that of DMTAP. For larger DMTAP mole fractions above $`\chi _{\mathrm{TAP}}=0.5`$, the situation is the opposite. This behavior is related to electrostatic effects and the ordering of acyl chains, and will be discussed in more detail in Sect. III.2. Here we only note that the fluctuations in $`A`$ (see Fig. 3) at $`0.1\chi _{\mathrm{TAP}}0.5`$ arise from fluctuations in the area occupied by DMTAP. Whether this is a true result due to, e.g., clustering of lipids in this region remains to be resolved.

### III.2 Ordering of lipid acyl chains

Ordering of non-polar hydrocarbon chains in lipid bilayers is typically characterized by the deuterium order parameter $`S_{CD}`$ measured through 2H NMR experiments. If $`\theta `$ is the angle between a CD-bond and the bilayer normal, the order parameter is defined as

$$S_{CD}=\frac{3}{2}\mathrm{cos}^2\theta \frac{1}{2}$$
(1)
separately for each hydrocarbon group. Since we employed a united-atom force field, the positions of the deuterium atoms are not directly available but have to be reconstructed from the coordinates of three successive nonpolar hydrocarbons, assuming an ideal tetrahedral geometry of the central CH2 group Chiu et al., (1995); Hofsäss et al., (2003). In practice, we calculated $`S_{CD}`$ following the standard approach described elsewhere Patra et al., (2003).

Figure 5 shows $`|S_{CD}|`$ averaged over the two similar atoms in the $`sn`$–1 and $`sn`$–2 chains, for both DMPC (top) and DMTAP (bottom) at different DMTAP concentrations. For the pure DMPC bilayer, we find $`|S_{CD}|0.18`$ close to the glycerol group of the molecule, in good agreement with recent experiments Petrache et al., (2000) and molecular dynamics simulation studies Róg and Pasenkiewicz-Gierula, (2001); Smondyrev and Berkowitz, (2001).

The order parameter profile for the first seven hydrocarbons (from C2 to C8) is a kind of plateau, and the average value of $`S_{CD}`$ in this region, denoted by $`S_{\mathrm{ave}}`$, is shown in Fig. 6. As expected, the results closely follow the change in the average area per lipid, i.e., a compression of the membrane is accompanied by enhanced ordering of nonpolar acyl chains. Results of similar nature have been observed, e.g., in bilayer mixtures of glycerophospholipids and cholesterol, where cholesterol both reduces the average area per molecule and enhances the ordering of lipid acyl chains Hofsäss et al., (2003).

For pure DMPC the plateau value $`S_{\mathrm{ave}}`$ equals to $`0.181\pm 0.009`$, which is in very good agreement with the experimentally measured value of $`0.184`$ Petrache et al., (2000). For a pure DMTAP bilayer, the plateau value of $`S_{CD}`$ is considerably smaller, being about $`0.147\pm 0.011`$. Therefore chains in a pure DMTAP bilayer are on average more disordered than in a DMPC system, in agreement with our findings for $`A`$.

### III.3 Orientation of phosphatidylcholine headgroups

Since the chemical structures of the acyl chains of DMPC and DMTAP are identical, it is obvious that the differences between their behavior are due to their headgroups. Since the headgroup of DMPC is zwitterionic (cf. Fig. 1), it possesses a dipole moment along the P – N vector. The electrostatic potential across a monolayer thus depends sensitively on the distribution of the angle $`\alpha `$ between the P – N vector and the interfacial normal $`\stackrel{}{n}`$ (where $`\stackrel{}{n}`$ has been chosen to point away from the bilayer center along the $`z`$ coordinate).

Figure 7 shows the probability distribution function $`P(\alpha )`$ for the angle in question. For a pure DMPC bilayer we find the distribution to be wide, thus allowing the P – N vector to fluctuate substantially, pointing at times in the direction of the membrane normal as well as into the bilayer interior. The average angle found in this case is about $`(80\pm 1)^{}`$ (see Fig. 8), i. e., the PC heads are on average almost parallel to the membrane surface. This is in agreement with experimental observations Hauser et al., (1981); Scherer and Seelig, (1989) as well as with recent computer simulations Gabdoulline et al., (1996); Smondyrev and Berkowitz, (1999); Pasenkiewicz-Gierula et al., (1999).

Figure 7 further reveals the major role played by DMTAP. The addition of even a small amount of DMTAP leads to a pronounced re-orientation of the P – N dipole vector. As $`\chi _{\mathrm{TAP}}`$ is increased, the profile of the distribution becomes considerably narrower and its maximum shifts to smaller angles. This trend continues up to the high-concentration limit $`\chi _{\mathrm{TAP}}0.75`$, beyond which the distribution is essentially similar with the case found for $`\chi _{\mathrm{TAP}}=0.75`$.

Results for the average angle $`\alpha `$ between the P – N vector and the membrane normal shown in Fig. 8 are consistent with this picture. On average, upon increasing $`\chi _{\mathrm{TAP}}`$, PC headgroups become more and more vertically oriented. Also this has been observed in several experiments Scherer and Seelig, (1989); Zantl et al., 1999b ; Säily et al., (2001, 2003). Moreover, our findings are in fairly good agreement with an atomistic MD study of a complex comprised of DNA and a mixture of DMPC and DMTAP Bandyopadhyay et al., (1999), in which the average angle between the P – N dipole vector and the bilayer normal was found to be $`(50\pm 8)^{}`$ at an almost equimolar mixture of DMPC and DMTAP. In the present case without DNA, we found $`\alpha =(42\pm 2)^{}`$ at $`\chi _{\mathrm{TAP}}=0.5`$.

Perhaps surprisingly, the correlation between the average area per lipid (Fig. 3) and the re-orientation of the P – N dipole is not complete. As Fig. 8 shows, the re-orientation extends by and large linearly up to $`\chi _{\mathrm{TAP}}=0.75`$, while the membrane compression completes at $`\chi _{\mathrm{TAP}}0.5`$. This contradicts the conclusions of Säily et al. Säily et al., (2001) who studied POPC / SS–1 cationic lipid mixtures using the Langmuir balance technique and suggested that the maximal average angle between the P – N dipole vector and membrane surface is achieved at the cationic lipid concentration that corresponds to the point where the membrane compression ends.

A closer inspection of Figs. 3 and 8 shows that the observed reduction of the average area per molecule is likely related to the re-orientation of PC dipoles, and this in turn is related to the role of electrostatic interactions between DMPC and DMTAP headgroups. To bridge the two issues, we propose the following schematic scenario. At small $`\chi _{\mathrm{TAP}}`$ where the DMTAP molecules are far apart and their mutual interaction is rather weak, we essentially suggest that the role of DMTAP is to re-orient the headgroups of those DMPC molecules that are beside a DMTAP molecule. This favors more dense packing at small $`\chi _{\mathrm{TAP}}`$, leading to a reduction in $`A`$, and consequently to a minimum in the area per molecule at intermediate concentrations as for large $`\chi _{\mathrm{TAP}}`$ the repulsive electrostatic interactions between TAP headgroups enforce $`A`$ to be expanded.

To validate this scenario, we complemented our results in Fig. 7 by calculating the probability distribution function $`P(\alpha )`$ for those DMPC molecules that are nearest neighbors to DMTAP.

As a criterion that a DMPC and a DMTAP form a pair, we monitored the distance between PC phosphorus and TAP nitrogen. For that, we first calculated the radial distribution functions (RDFs) between pairs of PPC and NTAP and determined the distance $`r_{nn}`$ at which the RDF had its first minimum after the main peak (see also Sect. III.6). The distance obtained in this fashion ($`r_{nn}0.665`$ nm) (and found not to depend on $`\chi _{\mathrm{TAP}}`$) was applied to identify the DMPCs residing next to a DMTAP. As shown in Fig. 8, the re-orientation of the P – N vector of these DMPCs is considerably stronger at small DMTAP concentrations as compared to that averaged over all DMPC lipids.

The above results imply that at small $`\chi _{\mathrm{TAP}}`$ the effect of DMTAP on the re-orientation is mainly local, i.e., the alternating PC and TAP headgroups pack more tightly than in a pure DMPC system. This idea is supported by the results for the radial distribution functions discussed in Sect. III.6. Beyond the small-$`\chi _{\mathrm{TAP}}`$ regime, for intermediate concentrations $`0.3\chi _{\mathrm{TAP}}0.5`$, further increase in the concentration of DMTAP continues to increase the number of units composed of PC and TAP heads, thus favoring a reduction in $`A`$. However, as repulsive electrostatic interactions between DMTAP molecules also become more and more important, the two effects compensate each other and $`A`$ is found to be approximately constant. Finally, for large $`\chi _{\mathrm{TAP}}`$, the repulsive electrostatic interactions between TAP groups dictate the case discussed here and lead to an enhancement of the average area per molecule.

Though this picture does not account for the explicit influence of counter-ions, it grasps essence of the process. The effect of counter-ions is discussed separately in Sect. III.6.

### III.4 Density profiles of lipid headgroups and chloride ions

To quantify the locations of charge groups and counter-ions, we computed the density profiles across the bilayer, separated into the different constituents of the system. The positions of all atoms in the system were determined with respect to the instantaneous center of mass position of the bilayer, exploiting mirror symmetry such that atoms with $`z<0`$ were folded to $`z>0`$ (the bilayer center being at $`z=0`$).

Figure 9 shows the scaled number densities $`\rho _N(z)`$ for a few selected cases. Additionally, we note that the essential information is given by the positions of the density maxima depicted in Fig. 10.

At small $`\chi _{\mathrm{TAP}}`$, the density maxima of the nitrogen and phosphorus atoms in the DMPC heads almost coincide (see Fig. 9). The density profile of nitrogen in DMPC is nevertheless broader and extends further out of the bilayer plane. For larger $`\chi _{\mathrm{TAP}}`$, the density profiles of phosphorus and nitrogen are distinctly separated, and nitrogen in particular extends rather deeply into the water phase. The TAP group represented by the nitrogen atom, however, is found to be deep in the bilayer. It seems obvious that these two issues are related, i. e., the density profiles of nitrogens in PC and TAP groups are well separated due to the electrostatic repulsion that essentially leads to the re-orientation of PC headgroups. These results are hence consistent with those in Sect. III.3 and reflect the dependence of DMPC headgroup orientation on $`\chi _{\mathrm{TAP}}`$.

Interestingly, while being attracted by the DMTAP headgroups, the chloride anions cannot penetrate the outer boundary of the bilayer formed by the DMPC choline groups. This is in a sense to be expected since the DMPC headgroup is longer than the TAP group and thus extends further outward from the bilayer. There is thus a significant amount of shielding of the chloride ions in the presence of DMPC. Only for an almost pure DMTAP layer are the chloride ions located in the vicinity of the TAP headgroups.

### III.5 Charge density, electrostatic potential, and orientation of water

The charge distribution shown in Fig. 11 was calculated in the same fashion as the density profiles. The results are clearly reminiscent of the density profiles in Fig. 9 and demonstrate the competition between charged PC and TAP groups, Cl anions, and water. The role of the TAP group is prominent.

From the charge densities, we also computed the electrostatic potential across the bilayer. The results are shown in Fig. 12 where the potential at the center of the bilayer has been set to zero. For a pure DMPC bilayer we obtain $`0.578`$ V, which is in agreement with previous MD simulation studies Chiu et al., (1995); Gabdoulline et al., (1996); Smondyrev and Berkowitz, (1999). Experimental data for phospholipid membranes ranges from $`0.2`$ V to $`0.6`$ V Hladky and Haydon, (1973); Pickard and Benz, (1978); Flewelling and Hubbel, (1986); Gawrisch et al., (1992); Clarke, (2001).

The charge of the DMTAP headgroup is mainly compensated by, depending on the value of $`\chi _{\mathrm{TAP}}`$, chloride ions or DMPC phosphate groups. Most of the electrostatic potential across the bilayer thus is not due to the DMTAP itself but rather due to the re-orientation of the DMPC headgroups. A clear indication of this is that the potential build-up saturates at $`\chi _{\mathrm{TAP}}0.75`$, i. e., at the same value at which the distribution of headgroup orientation saturates (see Fig. 8). The total potential of the bilayer increases with increasing DMTAP concentration, with a difference of 0.6 V between pure DMPC and pure DMTAP. This increase agrees well with the experimental data on cationic POPC / SS–1 monolayers Säily et al., (2001), and the authors offer the same explanation for their observations.

Many of the conclusions drawn from the charge density already follow from the number densities presented in Sect. III.4, since for charged particles number density and charge density are trivially related. Water, however, has an additional internal degree of freedom, and a quick discussion of the orientation of the water molecules seems appropriate. As seen from Fig. 13, the average direction of the water dipoles in the membrane-water interface region is inverted for $`\chi _{\mathrm{TAP}}=0.5\chi _{\mathrm{TAP}}=1.0`$. This is closely related with the familiar “hump” close to the interface, which is due to a subtle imbalance between the orientation of the water molecules and lipid headgroups Chiu et al., (1995). At higher DMTAP concentrations this “hump” disappears. A related issue concerns the pure DMTAP bilayer, in which case the density profile of water penetrates rather deep into the membrane (see Fig. 9), extending up to the interface region between the polar TAP group and the hydrophobic core. This is in accord with the interpretation of Fourier-transform infrared spectroscopic measurements by Lewis et al. Lewis et al., (2001).

### III.6 Radial distribution functions and coordination numbers

To characterize the structure of the membrane-water interface region in more detail, we computed various radial distribution functions (RDFs) between the atoms in the headgroups, Cl, and the oxygens. Here we briefly discuss the most relevant results.

The RDFs between the center of mass positions indicated that the leading (main) peak for DMPC – DMPC and DMTAP – DMTAP pairs was rather broad and at about 1.0 nm (data not shown). For DMPC – DMTAP pairs, however, the main peak of the RDF was much closer, at around 0.7 nm. This supports the conclusion made in Sect. III.3, i.e., DMPC and DMTAP form units that allow more dense packing than in a pure DMPC bilayer.

The NPC – NPC and NTAP – NTAP pairs were found to be rather far apart, the position of their main peak being at about 0.83 nm, while the NPC – NTAP pair was slightly closer (0.8 nm). The positions of the main peaks did not depend on $`\chi _{\mathrm{TAP}}`$. As for the RDFs of the phosphorus atoms in the PC headgroups, its main peak with respect to NPC and NTAP was found to be at a much closer distance, at 0.465 nm for NPC and 0.485 nm for NTAP. Again, the positions of these peaks did not depend on the DMTAP concentration.

We also calculated the coordination numbers for the phosphorus and nitrogen atoms at different DMTAP concentrations. These are shown in Fig. 14 (top). It turns out that in the range from $`\chi _{\mathrm{TAP}}=0`$ to 0.5 the PC nitrogens are to an increasing extent being replaced by NTAP in the vicinity of P. This has twofold consequences: First, the electrostatic attraction between N+ (TAP) and P- (PC) enhances the compression of the bilayer for $`0.0<\chi _{\mathrm{TAP}}0.5`$ (see Fig. 3). Second, the decreasing coordination number for P – NPC with $`\chi _{\mathrm{TAP}}`$ support the view that the DMPC nitrogens are pushed towards water, thereby PC headgroups are re-oriented to a more vertical alignment with respect to the membrane plane (Fig. 8).

We conclude this section with a discussion of the location of chloride counter-ions. The positions of the main peaks of the RDFs for both types of nitrogens (in PC and TAP groups) with respect to Cl ions are identical, about 0.475 nm, and do not depend on $`\chi _{\mathrm{TAP}}`$ (data not shown). In Fig. 14 (bottom) we plot the coordination numbers of chlorides in the vicinity of both types of nitrogens as a function of DMTAP concentration. The figure confirms that Cl ions are preferentially bound to PC nitrogens rather than to NTAP. This holds up to a DMTAP mole fraction of about 0.75. The explanation for this is straightforward: As $`\chi _{\mathrm{TAP}}`$ increases, the PC headgroups become more and more vertically oriented with respect to the bilayer plane. This, in turn, makes PC nitrogens more easily accessible for the Cl ions. In contrast, small TAP heads are located much deeper in the membrane surface region than the PC heads and therefore are able to attract fewer chlorides regardless of the fact that TAP heads carry a net positive charge. Interestingly, when the re-orientation of PC heads is accomplished (at $`\chi _{\mathrm{TAP}}0.75`$) the coordination number for NPC – Cl pairs seems to saturate, see Fig. 14 (bottom).

## IV Summary and conclusions

Drug delivery and gene therapy have attracted substantial interest due to their importance in treating human diseases. As far as experimental work is concerned, particular attention has been paid to non-viral delivery vectors such as cationic liposomes characterized by a number of desired properties such as high efficiency and lack of toxicity. Consequently, it is surprising how little is known about the atomic-level details of cationic lipid bilayers. Essentially, this is due to the lack of molecular simulations of these systems, the study by Bandyopadhyay et al. Bandyopadhyay et al., (1999) being the only exception, to the authors’ knowledge, in this regard.

As a first step toward a detailed understanding of cationic membrane–DNA complexes on atomic level, we have employed extensive molecular dynamics simulations of lipid bilayer mixtures composed of cationic DMTAP and neutral (zwitterionic) DMPC. Such binary DMPC / DMTAP mixtures have been studied widely through experiments, and have been shown to form stable complexes with DNA Artzner et al., (1998); Zantl et al., 1999a ; Zantl et al., 1999b ; Pohle et al., (2000). In the present work, we have focused on the influence of the composition of the cationic bilayer on its structural and electrostatic properties. For this purpose we studied numerous DMPC / DMTAP mixtures in the liquid-crystalline phase by varying the mole fraction of DMTAP, $`\chi _{\mathrm{TAP}}`$, from the pure DMPC to the pure DMTAP bilayer.

We have found that the properties of the DMPC / DMTAP bilayer mixture are largely dominated by the electrostatic properties of the headgroup region around the membrane-water interface. Most notably, our results indicate that there is a strong interplay between the PC and TAP groups together with the Cl counter-ions that concentrate in the vicinity of the bilayer-water interface.

The interplay between the PC and TAP groups leads to a number of intriguing observations. The key factor here is the re-orientation of PC groups due to an introduction of DMTAP in the bilayer. The re-orientation of the PC headgroups arises from electrostatic interactions that lead phosphate and choline groups to rearrange their positions with respect to the cationic TAP. This effect is enhanced as $`\chi _{\mathrm{TAP}}`$ is increased, and extends up to large molar fractions of approximately $`\chi _{\mathrm{TAP}}=0.75`$. Beyond this limit a further increase of DMTAP concentration has no additional effect on the orientation of PC headgroups. Interestingly, at small $`\chi _{\mathrm{TAP}}`$ the effect of the re-orientation is of local nature, i.e., the P – N dipoles of DMPCs beside DMTAP molecules re-orient considerably.

At small molar fractions of DMTAP, the re-orientation of PC dipoles leads to considerable compression of the bilayer as alternating PC and TAP groups are able to pack more tightly than in a pure DMPC bilayer. The minimum of the area per lipid at $`\chi _{\mathrm{TAP}}0.5`$ is about 12 % smaller than in the pure DMPC bilayer. A further increase of $`\chi _{\mathrm{TAP}}`$ leads to major expansion of the bilayer. This is essentially due to an increasing number of neighboring TAP groups whose cationic nature leads to repulsive electrostatic interactions that do not favor close packing. As expected, the ordering of acyl chains closely follows the change in the area per lipid. When these results are summarized, the present view of the average area per lipid coupled to the re-orientation of the headgroups can be summarized schematically as in Fig. 15.

In view of future studies of DNA-membrane systems, it is important to pay attention to the influence of DMTAP on the electrostatic properties of the membrane, including the increase in the electrostatic potential across the bilayer and the ordering of water in the vicinity of the membrane-water interface. Another often ignored aspect of electrostatics is that the ionic buffer liquids may affect membranes significantly Böckmann et al., (2003).

Perhaps the most significant observation in this study is the spatial re-arrangement of PC and TAP headgroups which is expected to play a significant role in the condensation of DNA onto the membrane surface. The cationic TAP and choline groups then play a key role as the anionic phosphate groups of DNA come into contact with the membrane. While the present study clarifies some of the underlying questions related to binary mixtures of cationic and neutral (zwitterionic) lipid membranes, further atomic-level studies are essential to resolve other important issues related to DNA-membrane systems, such as the influence of salt and its screening effects, and the stability and interface properties under those conditions. Work in this direction is under way.

## V Acknowledgments

We wish to thank M. T. Hyvönen, P. Niemelä, and X. Periole for fruitful discussions. This work has been supported by the Academy of Finland Grant Nos. 202598 (A. G.), 54113 (M. K.), and 80246 (I. V.), the Academy of Finland through its Center of Excellence Program (A. G., I. V.), and by the European Union through the Marie Curie fellowship HPMF–CT–2002–01794 (M. P.). We would also like to thank the Finnish IT Center for Science (CSC) and the HorseShoe (DCSC) supercluster computing facility at the University of Southern Denmark for computer resources.


text:
# Numerical study of metastable states in Ising spin glasses

(December 19, 2003)

## Abstract

We study numerically the structure of metastable states in the Sherrington-Kirkpatrick spin glass. We find that all non-paramagnetic stationary points of the free energy are organized into pairs, consisting in a minimum and a saddle of order one, which coalesce in the thermodynamic limit. Within the annealed approximation, the entropic contribution of these states, that is the complexity, is compatible with the supersymmetry-breaking calculation performed in \[A.J. Bray and M.A. Moore, J. Phys. C 13 L469 (1980)\]. This result indicates that the supersymmetry is spontaneously broken in the Sherrington-Kirkpatrick model.

In this last year there has been an outburst of new interest in the complexity of disordered systems, and in particular of spin-glasses zec ; brst1 ; brst2 ; brst3 ; leuzzi1 ; leuzzi2 ; abm ; plefka . The complexity is the entropic contribution due to the exponentially large number of metastable states, and in mean-field disordered models it can be computed by calculating the number of minima of the Thouless-Anderson-Palmer (TAP) free energy tap . Metastable states are crucial for the dynamical behaviour of disordered systems, and so is the complexity. Most notably, in glasses and some spin-glasses glass ; pspin  the complexity triggers a dynamical transition which is not associated to any static transition. Issues related to the analytic and numerical calculation of the complexity are therefore very relevant.

At the origin of the recent new studies there has been the observation done in brst1 , that the original calculation of the complexity for the Sherrington-Kirkpatrick (SK) model sk , performed by Bray and Moore in 1980 bm1 , breaks an intrinsic symmetry of the problem, namely a generalized form of the Becchi-Rouet-Stora-Tyutin (BRST) supersymmetry brs ; juanpe . Moreover, the supersymmetric complexity computed in brst1  was found to be much smaller that the supersymmetry-breaking complexity of bm1 . This discovery reopened, after more than twenty years, the problem of how to compute the complexity in the SK model, and, more in general, it raised the question of whether the BRST supersymmetry is spontaneously broken or not in disordered systems.

At the theoretical level, the situation is presently rather open. After the annealed calculation of brst1 , the supersymmetric (SS) complexity has been computed also at the quenched level brst2 ; brst3 , showing that it is consistent with the static solution of the SK model and that it coincides with the one obtained from the Legendre transform method monasson , at any level of replica symmetry breaking. However, it has been shown in leuzzi1 ; leuzzi2  that the SS complexity is stable only at the ground state free energy, where by definition it is zero. Thus, were the supersymmetry unbroken, this last result would lead to the conclusion that in the SK model the complexity of metastable states is in fact zero. On the other hand, the supersymmetry-breaking (SSB) complexity (which has only been computed at the annealed level bm1 ) is non-zero and stable on a finite range of free energy densities. Moreover, it has recently been proved in abm  that the SSB complexity describes pairs of solutions of the TAP equations which are either minima or saddles of order one of the free energy. According to abm , the two solutions of each pair coalesce in the thermodynamic limit. In conclusion, the SS and SSB predictions differ considerably.

Given the ambiguous theoretical situation, a way to discriminate between the two pictures above, and thus to understand whether or not the BRST supersymmetry is spontaneously broken in the SK model, is to perform a numerical experiment. Of course, one may question what would be the dynamical role of the marginally unstable states described by the SSB complexity, but we do not deal with this problem in the present work. Our goal here is just to understand what is the structure of TAP solutions in the SK model, and which theoretical framework correctly describes this structure.

Our numerical results show that: i) all solutions of the TAP equations are either minima or saddles of order one of the free energy; ii) the solutions are organized into minimum-saddle pairs, connected along a mode which is softer the larger the system size; iii) the free-energy difference of the paired solutions decreases with increasing system size; iv) the annealed complexity of TAP solutions as a function of their free energy density is consistent with the SSB complexity of bm1 . Our results therefore support the idea that the BRST supersymmetry is in fact broken in the SK model, and that the SSB complexity of bm1  is correct within the annealed approximation.

The TAP equations for the SK model are,

$$\mathrm{tanh}^1m_i+\beta ^2(1q)m_i\beta \underset{j}{}J_{ij}m_j=0,$$
(1)
where $`m_i`$ are the local magnetizations, with $`i=1,\mathrm{},N`$, and the random couplings $`J_{ij}`$ are drawn from a Gaussian distribution of mean zero and variance $`1/N`$. These equations are obtained by finding the stationary points of the TAP free energy,

$`F_{TAP}(m)={\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}J_{ij}m_im_j{\displaystyle \frac{N}{\beta }}\mathrm{log}2{\displaystyle \frac{N\beta }{4}}(1q)^2`$
$`+{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{i}{}}\mathrm{log}(1m_i^2)/2+m_i\mathrm{tanh}^1(m_i),`$
where $`q=(1/N)_im_i^2`$ is the self-overlap of a solution. We solve numerically equations (1) by using the C routine BROYDN of Numerical Recipes bro  starting with random initial conditions. We analyze $`20429`$ samples at $`N=20`$, $`5000`$ samples at $`N=30`$, $`1000`$ samples at $`N=40`$, $`347`$ samples at $`N=60`$ and $`77`$ samples at $`N=80`$. All data are for $`T=0.2T_c`$. In order to make our search of solutions as exhaustive as possible, we use two methods. First, all solutions we find in a sample are found at least $`5`$ times. Second, we monitor the Morse sum, $`_\alpha (1)^{k_\alpha }=1`$, where $`k_\alpha `$ is the number of negative eigenvalues of solution $`\alpha `$. For $`N=20,30,40,60`$ we discard samples where this identity is violated for more than $`\pm 2`$ (the fraction of discarded samples is always lower than $`1\%`$). For $`N=80`$, however, due to the low number of samples, we do not enforce this last check.

The first interesting result is that all solutions we find, at any value of $`N`$, are either minima, or saddles with one negative eigenvalue of the Hessian (the paramagnet is a minimum). As an example, the solutions of a sample at $`N=80`$ are reported in Fig.1. Moreover, all non-paramagnetic solutions are organized into minimum-saddle pairs: by making a gradient descent starting from a saddle and moving along the negative mode, we reach in one direction a very close minimum foot , and in the opposite direction the paramagnetic state. Thus, TAP pairs are not directly connected to each other, but they are all connected to the central paramagnetic minimum in a star-like structure. The self-overlap $`q`$ of a saddle is always smaller than that of the connected non-paramagnetic minimum. The important point is that the two solutions of each pair have very close free energy densities $`f`$: by plotting the average free energy density difference $`\mathrm{\Delta }f`$ as a function of $`N`$, we see that $`\mathrm{\Delta }f`$ decreases with $`N`$ (Fig. 2). A power law fit of the data gives $`\mathrm{\Delta }f1/N^{1.26}`$. The paired solutions are connected along the saddle negative mode, which becomes the smallest positive mode of the nearby minimum. The fact that the two solutions have very close free energy suggests that this direction is extremely flat. In Fig.2 we also plot the average modulus of the smallest eigenvalue $`\lambda _{\mathrm{min}}`$ as a function of $`N`$. A power law fit of the data gives $`|\lambda _{\mathrm{min}}|1/N^{1.40}`$. Thus the direction connecting the two solutions of a pair is flatter the larger $`N`$. The exponent $`1/2`$ proposed in abm  seem to fit the data less satisfyingly.

These data suggests that for $`N\mathrm{}`$ the two solutions of each pair coalesce, forming a single, marginally unstable TAP solution. This conclusion is in qualitative agreement with the theoretical results of abm , results found in connection with the SSB annealed complexity. It is then natural to ask whether also the annealed numerical complexity is compatible with the SSB result of bm1 . The annealed complexity is defined as,

$$\mathrm{\Sigma }(f)=\frac{1}{N}\mathrm{log}𝒩(f),$$
(2)
where $`𝒩(f)`$ is the number of solutions with free energy density $`f`$, averaged over the disorder $`J`$. The SSB annealed complexity is a smooth function of $`f`$, with a maximum at a free energy density $`f_{\mathrm{max}}`$. Let us define $`\mathrm{\Sigma }_{\mathrm{max}}=\mathrm{\Sigma }(f_{\mathrm{max}})`$ and $`\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime }=|\mathrm{\Sigma }^{\prime \prime }(f_{\mathrm{max}})|`$. For large $`N`$ the total number of TAP solutions $`𝒩`$ is given by,

$$𝒩=𝑑fe^{N\mathrm{\Sigma }(f)}=e^{N\mathrm{\Sigma }_{\mathrm{max}}}\sqrt{\frac{2\pi }{N\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime }}},$$
(3)
where we have expanded $`\mathrm{\Sigma }(f)`$ to second order close to $`f_{\mathrm{max}}`$. The probability $`p(f)`$ of finding a solution with free energy density $`f`$ is therefore given by,

$$p(f)=\frac{𝒩(f)}{𝒩}=\sqrt{\frac{N\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime }}{2\pi }}e^{\frac{1}{2}N\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime }(ff_{\mathrm{max}})^2},$$

which becomes a $`\delta `$-distribution for $`N\mathrm{}`$. For finite $`N`$ the probability distribution $`p(f)`$ is a numerically directly accessible quantity. In particular, its variance $`\sigma `$ can be easily checked without the need of any binning of the data. In Fig.3 we plot the variance $`\sigma `$ as a function of $`N`$. The full line is the prediction given by the SSB complexity, i.e. $`\sigma =1/\sqrt{N\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime }}`$, with $`\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime (\mathrm{SSB})}=8.9`$. Numerical data are compatible with the SSB prediction. Moreover, if we perform a fit to $`\sigma N^{1/2}`$ we find $`\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime }=10.3`$, and if we leave the exponent free, fitting to $`\sigma N^\gamma `$ we get $`\gamma =0.58`$ (not shown). We also show the best fit to $`\sigma N^{1/3}`$, which is the exponent reported in plefka . This value, however, does not seem compatible with our data.

By making a (arbitrary) binning of the data in free energy, we can compute the full probability $`p(f)`$ and thus our numerical estimate of the complexity, $`\mathrm{\Sigma }(f)=\mathrm{log}[p(f)𝒩]/N`$, where $`𝒩`$ is the total average number of solutions per sample we find numerically. Results are shown in Fig.4 and compared to the theoretical curve of the SSB complexity $`\mathrm{\Sigma }(f)`$. Even though the numerics is not excellent (especially for $`N=80`$), the data seem compatible with the SSB prediction and clearly indicate that the annealed complexity is nonzero. However, a more quantitative comparison, not dependent on the binning, is clearly desirable. To do this we study the position of the maximum $`f_{\mathrm{max}}`$ and the value of the complexity at this point $`\mathrm{\Sigma }_{\mathrm{max}}`$, at various $`N`$. In order to avoid a measurement depending on the binning, we study the integral of $`p(f)`$, i.e. $`g(f)=_0^f𝑑xp(x)`$, and interpolate $`g(f)`$ with a cubic function around its inflection point. In this way we have an estimate of both $`f_{\mathrm{max}}`$ and $`\mathrm{\Sigma }_{\mathrm{max}}`$, reported in Fig.5 as a function of $`1/N`$. The value of $`f_{\mathrm{max}}`$ is basically constant, and it agrees quite well with the theoretical SSB value $`f_{\mathrm{max}}^{(\mathrm{SSB})}=0.654`$. On the same scale we report the value of the ground state (equilibrium) free energy density $`f_0`$: our data do not seem compatible with the law $`f_{\mathrm{max}}f_0`$ for $`N\mathrm{}`$, proposed in plefka . The value of $`\mathrm{\Sigma }_{\mathrm{max}}`$ depends more strongly on $`N`$, as expected from equation (3), i.e. $`\mathrm{\Sigma }_{\mathrm{max}}^N=\mathrm{\Sigma }_{\mathrm{max}}^{\mathrm{}}+\mathrm{log}(N\mathrm{\Sigma }_{\mathrm{max}}^{\prime \prime })/(2N)+\mathrm{O}(1/N)`$. The extrapolation for $`N\mathrm{}`$ of our data for $`\mathrm{\Sigma }_{\mathrm{max}}`$ agrees very well with the SSB value $`\mathrm{\Sigma }_{\mathrm{max}}^{(\mathrm{SSB})}=0.052`$.

A last important consistency check between our numerical data and the SSB calculation concerns the behaviour of the overlap $`q`$ and the stability parameter $`x_p`$, defined as, $`x_p=1\beta ^2_i(1m_i^2)^2/N`$. In the thermodynamic limit the condition $`x_p0`$ must be satisfied in order to have a physically acceptable TAP state xp1 . In fact, it is precisely this last condition which is violated by the SS complexity in all points but the equilibrium free energy $`f_0`$. Both $`q`$ and $`x_p`$ depend on the free energy density $`f`$ of the solutions. In Fig.6 we plot $`x_p`$ as a function of $`q`$, parametrically in $`f`$, at various values of $`N`$. We see that $`x_p`$ is in fact positive for all solutions in the region where $`\mathrm{\Sigma }(f)0`$. The data agree quite well with the SSB prediction, even in the phase $`f<f_0`$, where $`\mathrm{\Sigma }(f)<0`$.

Our algorithm finds more easily solutions with lower free energy, as can be seen by checking the frequency each solution is found, and this effect is stronger the larger $`N`$. Therefore, at $`N=80`$ (where the Morse check is not enforced) we are probably missing a number of solutions with $`f>f_{\mathrm{max}}`$. We believe this is the reason why the values of $`\sigma `$ and $`f_{\mathrm{max}}`$ at $`N=80`$ in Figs. 3 and 5 are slightly smaller than expected, irrespective of error bars. The same problem arises in the reconstruction of $`\mathrm{\Sigma }(f;N=80)`$ for $`f>f_{\mathrm{max}}`$ in Fig.4. The problem of finding solutions at finite $`T`$ of the TAP equations is well known, and for this reason numerical investigations in the past have been done by using alternative methods: $`T=0`$ studies of one spin-flip stable states T0 , naive TAP equations naive , or modified TAP equations plefkamod . To the best of our knowledge, the present work is the first extensive study of the full TAP equations at finite $`T`$.

In summary, our results support the conclusion that at the annealed level the complexity of the SK model is given by the SSB solution of bm1 , and therefore that the BRST supersymmetry is spontaneously broken in this system. Considering the fact that the supersymmetry is not broken in the $`p`$-spin spherical model, it may be argued that systems characterized by full replica symmetry breaking (as the SK model), and systems solved by one step of replica symmetry breaking (as the $`p`$-spin), belong to two different classes also for what concerns the spontaneous breaking of the supersymmetry. The origin of this different behaviour must lie in the geometric structure of the states. In particular, the marginality of TAP saddle-minima pairs in the SK model is probably the cause of the spontaneous supersymmetry breaking. This fact is also likely to be connected to the different dynamical behaviour of these two classes of systems: the dynamics of the $`p`$-spin model is heavily influenced by the TAP states, while in the SK model, despite the non-zero SSB complexity, the dynamics asymptotically reproduces the static results, with no role played by the metastable states. The supersymmetry could then be an elegant way of distinguishing different dynamical classes.

We thank A. Bray, L. Leuzzi, M. Moore, T. Plefka, T. Rizzo and the Rome Complex Systems Group for discussions, and A. Annibale and E. Trevigne for the SSB data. We acknowledge support of the ESF-SPHINX program.


text:
# HI holes in galactic disks: Tracing the dark matter distribution

## 1  Introduction

The interstellar medium (ISM) shaped by the energy release from exploding massive stars provides the best sample for studying star formation processes on a “microscopic” level. The diversity of HI structures and associated H$`\alpha `$ distributions discovered in our own and nearby galaxies shows that star formation processes and dynamics of interaction of massive stars with the ISM are strongly environmentally dependent. Prominent examples are found in the Milky Way (Heiles 1984), M31 (Brinks Brinks (1981), Brinks & Bajaja BB (1986)), M33 (Deul & den Hartog Deul (1990)), Holmberg II (Puche et al. Puche (1992)), NGC 891 (Dettmar & Schulz Dettmar1 (1992)), M101 and NGC 6946 (Kamphuis Kam (1993)), NGC 253 (Sofue, Wakamatsu & Malin Sofue (1994)), SMC (Stavely-Smith et al. SS (1997)) and LMC (Kim et al. Kim (1999)), IC 10 (Wilcots & Miller Wilcots (1998)), IC 2574 (Walter & Brinks Walter (1999)), Holmberg I (Ott et al. Ott (2001)), and others. Edge-on (such as NGC 891, NGC 4631, NGC 253) and face-on (such as Holmberg I and II, IC 2574, SMC and LMC) galaxies show these “microscopic” structures in vertical and horizontal directions, respectively, and develop our understanding of not only the interrelation between the observed HI/H$`\alpha `$ structures and star formation activity (Dettmar Dettmar2 (1992), Rand Rand (1996), Sofue, Wakamatsu & Malin Sofue (1994), Rossa & Dettmar 2003, Stewart et al. Stewart (2000), Walter & Brinks Walter (1999)), but also the factors governing the star formation process itself.

Among these factors, the dark matter distribution in radial and vertical directions is of principle importance. At the same time, the influence of the dark matter distribution on the dynamics of HI and H$`\alpha `$ structures in edge-on and face-on galaxies was left aside general discussion. Only recently Ott et al. (2001) paid closer attention to the relevance of HI dynamics in the dwarf irregular galaxy Holmberg I to the vertical gas distribution. However, available estimates of the vertical scale heights are still based on qualitative arguments and assumptions on the contribution of the dark matter to the total gravitational potential. In this paper we address the question of whether the observations of morphologies and quantitative characteristics of the HI holes in face-on galaxies can provide us with firm conclusions about the vertical scale height of the HI distribution. For this purpose, we study numerically the characteristics of HI flows associated with multiple SN explosions. We determine the dependence of these characteristics on the vertical scale height of the HI distribution and, consequently, on the vertical shape of the gravitational potential and dark matter distribution.

The paper is organized as follows. In Sect. 2 the numerical model is formulated, and the quantitative characteristics of HI holes are numerically investigated for different vertical scale heights of gas distribution in face-on galaxies. In Sect. 3 the influence of the extragalactic background ionizing radiation on the vertical shape of the HI distribution is discussed. The main results are summarized in Sect. 4.

## 2 SNe driven supershells: dynamics and morphology

### 2.1 Model

The equilibrium HI distribution is obtained by solving the steady-state momentum equation in cylindrical $`(r,z)`$ coordinates in a fixed gravitational potential determined by the stellar and dark matter components. The self-gravity of gas is neglected. The dark matter density distribution ($`\rho _\mathrm{d}`$) is assumed to be that of a modified isothermal sphere (Binney & Tremaine Binney (1987))

$$\rho _\mathrm{d}=\frac{\rho _{\mathrm{d0}}}{(1+r/r_0)^2},$$
(1)
where $`r`$ is the radial distance from the galactic center, $`\rho _{\mathrm{d0}}`$ is the central density, and $`r_0`$ is the characteristic scale length. The stellar component is assumed to be vertically stratified (Spitzer Spitzer (1942))

$$\rho _{}=\rho _0\mathrm{sech}^2(z/h_{}),$$
(2)
where $`\rho _{}`$ is the stellar density, $`\rho _0`$ is the stellar midplane density, and $`h_{}`$ is the vertical scale height of stellar distribution.

The resulting equilibrium HI distribution is distinct for different dark matter and stellar distributions. Its vertical shape is neither Gaussian nor exponential (Celnik et al. Celnik (1979)), however, it is better fitted by a Gaussian than by an exponent. We consider three different models with three different gas vertical scale heights. The parameters of the stellar and dark matter distributions for each model are shown in Table 1. Fitting the equilibrium vertical HI distribution of each model by a Gaussian, we obtain the 1 $`\sigma `$ gas scale heights of $`h`$=210, 370, and 500 pc at the galactocentric radius of 2 kpc. For comparison, the corresponding exponential gas scale heights are $`h_{\mathrm{exp}}`$=230, 400, and 540 pc, respectively.

Further, we assume that the size of SN-driven shells is much smaller than that of a galaxy. The radial dependence of the dark matter gravitational potential within the shell can therefore be neglected and the initial equilibrium configuration of gas becomes a function of only the distance above the midplane of a galaxy. In all models the initial gas velocity dispersion was taken as $`v_\mathrm{g}=9`$ km s-1.

The energy of supernova explosions is released in a sphere with a radius of four zones. We use the constant wind approximation described in detail in Mac Low & Ferrara (MF (1999)). We convert the energy of each SN explosion ($`10^{51}`$ ergs) totally into the thermal energy, because in the present simulations we deal with large stellar clusters capable of producing a hundred supernovae. With such a number of SN explosions, the surrounding ISM will be quickly heated and diluted, making radiative cooling within the injection sphere ineffective. We choose the energy input phase to last for 30 Myr, which roughly corresponds to a difference in the lifetimes of the most and least massive stars capable of producing SNe in a cluster of simultaneously born stars.

SNe generate a supersonically expanding wind that compresses the gas, thus creating a bubble filled with the hot ejected gas surrounded by a shell of compressed cold material. The gas dynamics is followed by solving the usual set of hydrodynamical equations in cylindrical coordinates using the method of finite-differences with a time-explicit, operator split solution procedure of the ZEUS-2D numerical hydrodynamics code described in detail in Stone & Norman (Stone (1992)). We have implemented the cooling curve of Böhringer & Hensler (BH (1989)) for a metallicity of one tenth of solar. We use an empirical heating function tuned to balance the cooling in the background atmosphere so that it maintains the gas in hydrostatic equilibrium and may be thought of as a crude model for the stellar energy input. Cooling and heating are treated numerically using Newton-Raphson iterations, supplemented by a bisection algorithm for occasional zones where the Newton-Raphson method does not converge. In order to monitor accuracy, the total change in the internal energy density in one time step is kept below $`15\%`$. If this condition is not met, the time step is reduced and a solution is again sought.

### 2.2 Qualitative description

At the initial stages the remnant comes into a radiative phase relatively early, approximately at

$$t_r10^5\left(\frac{L_{38}}{n}\right)^{1/2}\mathrm{yr},$$
(3)
when its radius is only

$$R25\left(\frac{L_{38}}{n}\right)^{1/2}\mathrm{pc},$$
(4)
where $`L_{38}`$ is the mechanical luminosity in $`10^{38}`$ ergs s-1 released by multiple SN explosions, $`n`$ is the ambient volume density. This implies that the remnant is always radiative before it expands out to $`R2h`$ and enters a breakthrough phase (i.e. the shell breaking out of the disk and losing metal-enriched material to the intergalactic space). The breakthrough occurs at approximately

$$t10^3\left(\frac{n}{L_{38}}\right)^{3/2}h^4\mathrm{yr},$$
(5)
where $`h`$ is in parsecs. The further evolution shows a very fast drop of pressure inside the hot remnant and the “snowplow” expansion of a quasi-cylindrical shell of compressed cold material in the disk. Due to a sharp drop of pressure inside the remnant the expanding shell starts thickening and can no longer be treated as a thin shell. Figure 1a shows the temporal evolution of the inner and outer radii of the shell expanding in a face-on galaxy with a gas scale height of $`h=500`$ pc and HI surface density of $`\mathrm{\Sigma }_{\mathrm{HI}}=10.7M_{}`$ pc-2. The surface density of the swept-up gas is generally not constant across the shell, it reaches a maximum near the outer edge of the shell and falls off on both sides. We define the inner radius of the shell as a radial distance where $`\mathrm{\Sigma }_{\mathrm{HI}}`$ drops to its unperturbed value in the inner side of the shell. The outer radius is defined in the same manner. The mean radius is then obtained as an arithmetic average of the inner and outer radii. The mechanical luminosity of $`L_{38}=1`$ is assumed, which is equivalent to 100 SN explosions during the first 30 Myr. The cooling rate for $`T10^4`$ K is taken from Böhringer & Hensler (BH (1989)), for $`T<10^4`$ K the cooling function rapidly declines to zero. It is seen that the shell remains thinner than a quarter or less of its mean radius only in the SN-driven expansion phase during the first 30-35 Myr. At this evolutionary stage, the mean radius of the shell is well described by the $`R(L/\rho )^{1/5}t^{3/5}`$ law (Castor et al. Castor (1975)) plotted by the dotted line in Fig. 1a, where $`\rho =0.5\rho _0`$ and $`\rho _0`$ is the midplane mass density of the ambient gas. The shell grows very fast after $`t=35`$ Myr when a breakthrough is to occur. At the time of a breakthrough, the mean radius of the shell is $`R1.6h`$. After the breakthrough the shell enters the “snowplow” phase and its mean radius grows approximately as $`Rt^{1/3}`$ (as shown in Fig. 1a by the dotted line).

In case of a smaller gas scale height ($`h=210`$ pc), the breakthrough occurs at t=13.5 Myr when $`R1.8h`$. The mean radius of the shell is poorly described by both the spherical ($`Rt^{1/4}`$, McCray & Kafatos McCray (1987)) and cylindrical ($`Rt^{1/3}`$) thin-shell “snowplow” expansion laws as shown in Fig. 1b by the dashed and dotted lines, respectively. Thus, it is obvious that when an HI hole seen in face-on galaxies is fitted by the spherical thin-shell “snowplow” expansion law $`Rt^{1/4}`$, the resultant energy input into the hole from SN explosions may be largely overestimated.

### 2.3 Dynamics versus vertical gas scale height

At the early expansion stages, when the shells remain nearly spherical, their kinematical and dynamical characteristics scaled to the proper mechanical luminosity and local density are similar in different environments (see discussion in Castor et al. Castor (1975)). Therefore, one should explore the later expansion stages in order to determine the large scale properties of a parent galaxy. At the later stages, the expanding shells are radiative and their expansion velocities are only weakly supersonic so that the corresponding Mach number is $`M\stackrel{>}{}1`$ (in most face-on galaxies the expanding HI shells are marginally supersonic, or even subsonic as in the case of Holmberg II, see discussion in Shchekinov et al. YUS (2001)). In such conditions the compression factor of gas inside the shells (mass density contrast) is not high: for isothermal shock waves it is $`\rho /\rho _0=\gamma M^2`$, where $`\gamma `$ is the ratio of specific heats, and for weakly supersonic shocks at the later stages ($`t3040`$ Myr) it can be around $`\gamma `$.

A directly observed quantity is the HI column density. For an HI hole in a face-on galaxy created by consecutive SN explosions, the HI column density in a ring of compressed gas surrounding the hole is

$$N(\mathrm{HI})\frac{2\sqrt{2}}{3}\frac{R^{3/2}n_0}{\sqrt{\mathrm{\Delta }R}},$$
(6)
where $`R`$ is the ring radius, $`\mathrm{\Delta }R`$ is its thickness, $`n_0`$ is the midplane density of ambient gas $`n_0=N_0(\mathrm{HI})/h`$, $`N_0(\mathrm{HI})`$ is the column density of the background (unperturbed) HI gas. Substituting $`\mathrm{\Delta }R=R\rho _0/3\rho `$ and $`\rho /\rho _0\gamma `$ into Eq. (6), one obtains a linear dependence of the relative HI column density on its radius in a face-on ring

$$𝒦_{\mathrm{HI}}\frac{N(\mathrm{HI})}{N_0(\mathrm{HI})}2\frac{R}{h},$$
(7)
with the slope inversely proportional to the HI vertical scale height of a parent galaxy. The highest relative column density $`𝒦_{\mathrm{HI}}`$ is reached when the breakthrough occurs at $`R\alpha h`$, where $`\alpha `$ is a slowly varying function of the vertical gas scale height. Simulations in Sect. 2.2 indicate that $`\alpha 1.6`$ for $`h=500`$ pc, while for $`h=210`$ pc $`\alpha 1.8`$. An expanding ring reaches maximal relative column density $`𝒦_{\mathrm{HI}}4`$ only in galaxies with rather low vertical scale heights $`h200`$ pc, where the compressed shell always remains supersonic. After the breakthrough, the compressed gas in the shell starts to move into the cavity because of a sharp drop in pressure, and the column density in the ring decreases approximately as $`R^1`$.

When the gas scale height is sufficiently large, the shell expansion becomes sonic before breaking out of the disk, the associated perturbation propagates with the sound speed and simultaneously the shell itself starts swelling into the cavity, which results in a decrease of the gas column density. This occurs when

$$R230L_{38}^{1/2}h_{100}^{1/2}\mathrm{\Sigma }_{10}^{1/2}\mathrm{pc},$$
(8)
where $`h_{100}=h/100`$ pc, $`\mathrm{\Sigma }_{10}`$ is the unperturbed surface density of gas ($`\mathrm{\Sigma }_{g0}`$) in units of 10 $`M_{}`$ pc-2. For instance, when the gas scale height is $`h=500`$ pc and the unperturbed surface density $`\mathrm{\Sigma }_{g0}=10.5M_{}`$ pc-2, the shell becomes sonic when its radius is $`R500`$ pc. As a consequence, the shell starts swelling into the cavity and the corresponding relative column density of the ring decreases.

Figure 2 shows the relative column density in the ring, $`𝒦_{\mathrm{HI}}`$, obtained in our numerical simulations for three different vertical gas scale heights and two values of the undisturbed surface density $`\mathrm{\Sigma }_{g0}`$. The curves exhibit the characteristic behaviour described above, with the maximal values of $`𝒦_{\mathrm{HI}}`$ being higher for the lower vertical gas scale heights. Specifically, for $`h=210`$ pc the maximum value of $`𝒦_{\mathrm{HI}}`$ is approximately 3.4 at the time when the shell starts swelling because of a sharp pressure drop after the breakthrough. For larger $`h`$ the maximum value of $`𝒦_{\mathrm{HI}}`$ is smaller. Note that $`𝒦_{\mathrm{HI}}`$ is virtually independent of $`\mathrm{\Sigma }_{g0}`$. For the total period of energy injection $`30`$ Myr and for the vertical gas scale heights $`<600`$ pc, the maximum value of $`𝒦_{\mathrm{HI}}`$ and the corresponding radius of a shell do not depend on the luminosity $`L`$, while in gaseous disks with $`h600`$ pc the maximum value of $`𝒦_{\mathrm{HI}}`$ and the corresponding radius scale as $`L^{1/2}`$. However, since a widely accepted value for the mechanical luminosity $`L=10^{38}`$ ergs s-1 is typical for OB associations with a Salpeter IMF, one may expect that the curves depicted in Fig. 2 represent a universal relation and may be used to draw firm conclusions about the vertical scale heights of HI distributions in disk galaxies. It is obviously seen that the $`𝒦_{\mathrm{HI}}`$ versus $`R`$ relation is characteristic for each vertical gas scale height only in the early phases when the shell has not yet broken out of the disk, while after the breakthrough $`𝒦_{\mathrm{HI}}`$ can barely be distinguished among different $`h`$. Note in Fig. 2 a range in which the relative column density in the ring varies in galaxies with low vertical gas scale heights. It is 1.4-3.4 for $`h=210`$ pc, while $`𝒦_{\mathrm{HI}}`$ is restricted to a narrower range of 1.3-2 in galaxies with $`h=700`$ pc. Thus, a narrow spread of observed column densities may indirectly indicate that the vertical gas scale height is large. Moreover, considering a universal and single-valued dependence of $`𝒦_{\mathrm{HI}}(h,R)`$ for $`R600700`$ pc, the local gas scale height $`h`$ can be inferred directly from the measurements of $`𝒦_{\mathrm{HI}}`$. Further, the measurements of $`𝒦_{\mathrm{HI}}`$ for the individual HI holes located at different galactocentric distances allow for the determination of the radial variations in the gas scale height.

The increase in the mass of the dark matter generally makes a gas disk thinner, i.e. the gas scale height $`h`$ decreases as the dark matter mass increases. Assuming the local plane-parallel gas distribution, the total density $`\rho _{\mathrm{tot}}(0,r)`$ in the disk at the galactocentric radius $`r`$ and $`z=0`$ can be derived from $`h(r)=v_\mathrm{g}/\sqrt{4\pi G\rho _{\mathrm{tot}}(0,r)}`$ (van der Kruit Kruit (1981)). Provided that the radial distribution of the visible mass is known from independent measurements, one can obtain the radial distribution of the dark matter in the galactic plane. These are, of course, very approximate estimates; nevertheless, a similar approach applied to the dwarf irregular galaxy Holmberg I by Ott et al. (Ott (2001)) has yielded a dark matter mass within the HI content of the galaxy, the value of which is in agreement with that obtained in a more sophisticated numerical modeling by Vorobyov et al. (Vor (2004)). Note, however, that the flaring and warping of the gas disk may further complicate this analysis.

The velocity field of a shell expanding in a nearly face-on galaxy can additionally hint at the vertical structure of gas. At the initial stages far from breakthrough, the vertical component of gas velocity in an expanding shell, $`v_\mathrm{z}`$, scales as $`\sqrt{1(r/R)^2}`$, where $`r`$ is the projected distance from the geometrical center of a shell and $`R`$ is the shell radius. However, when the shell is near breakthrough, a substantial amount of gas constituting the shell becomes involved in strong vertical motion. In Fig. 3 we plot the mass-weighted $`v_\mathrm{z}`$ as a function of the projected distance from the geometrical center of a shell for different vertical gas scale heights and different fixed shell radii. The solid lines represent the shells with a radius of $`R=0.5`$ kpc, while the dotted lines correspond to the shells with a radius of $`R=0.7`$ kpc.

It is seen that the expanding gas is only weakly supersonic (the corresponding Mach number of 2-3) in the shells that are at the early phases of expansion, well before breakthrough ($`R=0.5`$ kpc, $`h=370`$ and 500 pc). On the other hand, the shell that approaches the breakthrough phase ($`R=0.5`$ kpc and $`h=210`$ pc) shows a steep velocity increase up to $`v_\mathrm{z}=60`$ km s-1 at the inner edge of the projected shell. The $`v_\mathrm{z}`$ radial profiles remain qualitatively similar at different phases of the shell expansion, namely at $`R=0.5`$ kpc and $`R=0.7`$ kpc.

At the breakthrough phase, the radius of a shell is almost twice as large as the gas scale height $`h`$. From this point of view, one may expect that the most extended HI holes in face-on galaxies (such as HI holes number 2, 5, 8, 10 and others in Holmberg II, Puche et al. Puche (1992) ) whose $`z`$-velocity reveals such a behaviour at the inner edges can provide us with a direct measure of the local gas scale height. Specifically, scanning the HI radial velocity profiles within a hole with a beam smaller than its size (for the Ho II case, the beam must be less than $`4^{\prime \prime }6^{\prime \prime }`$) could help to identify nearly breaking-through shells, the radii of which could then serve as a direct measure of the local gas scale height.

Out-of-plane SN explosions in nearly face-on galaxies can also provide quite a robust tool for the determination of the gas vertical scale height. In Fig. 4 we plot typical HI spectra of the shells produced by 100 successive SNe located at 100 pc above the midplane. The HI spectrum is obtained for a 1 kpc diameter beam centered on the expanding shell. Two model galaxies with different gas scale heights are considered: $`h=500`$ pc and $`h=210`$ pc. A one-dimensional velocity dispersion of $`\sigma _\mathrm{g}=3`$ km s-1 is assumed for the gas forming an expanding shell when constructing the model HI spectra, while $`\sigma _\mathrm{g}=9`$ km s-1 is adopted for the unperturbed gas. The HI spectrum of the shell expanding in a galaxy with $`h=500`$ pc is characterized by a typical two-humped structure. A smaller part of the accelerating shell expanding upwards ($`v_\mathrm{z}>0`$) and a more massive part of the shell expanding downwards ($`v_\mathrm{z}<0`$) with approximately constant velocity are clearly seen. The HI spectrum of the shell expanding in a galaxy with $`h=210`$ pc has a single peak located at $`v_\mathrm{z}0`$. Velocity centroids of gas expanding downwards are obviously smaller for lower gas scale heights $`h`$, because lower scale heights lead to a faster onset of breakthrough and, therefore, to a faster drop of pressure inside the SN-driven shell.

In Fig. 5 we plot the absolute values of the peak $`v_\mathrm{z}`$ velocities of gas expanding downwards as a function of the corresponding HI intensities (see Fig. 4) for three different gas scale heights: $`h=210`$, 370, and 500 pc. It is clearly seen that the smaller peak z-velocities are associated with the shells expanding in galaxies with lower gas scale heights. A distinguishing feature of this dependence is that the peak expansion velocity is virtually independent of intensity in a relatively wide range covering the latest stages of expansion before and after the breakthrough, and varies approximately as

$$v_z^{\mathrm{peak}}2\left(\frac{h}{100\mathrm{pc}}\right)^{0.9}\left(\frac{z_{}}{100\mathrm{pc}}\right)^{0.24}\mathrm{km}\mathrm{s}^1.$$
(9)
where $`z_{}`$ is the location of SNe above the midplane.

## 3 Effects of external ionizing radiation

One of the most striking facts related to HI holes in dwarf face-on galaxies is the lack of an underlying stellar population that could produce a strong energy input into the ISM. This tendency is obviously seen in Holmberg II (Ho II) where most extended HI holes are located outside the stellar disk (Rhode et al. Rhode (1999), Stewart et al. Stewart (2000)). Hence, the SN origin of those HI holes in Ho II is questionable. An alternative explanation of the HI hole formation connected with high-velocity HI clouds (HVCs) falling onto the disk seems to be excluded in the case of Ho II, because no HVCs are found nearby. In this case, Efremov et al. (Efremov (1998)) argue that a supernova in the merging event of a compact binary system may be responsible for HI hole creation. However, a distinguishing feature of multiple HI holes on the periphery of Ho II is that they form quite a regular structure resembling a spiral wave (see Stewart et al. Stewart (2000)), and seem unlikely to be produced by regularly placed supernovae explosions. In this connection it is worth mentioning the role of an external ionizing radiation field in shaping the HI distribution, particularly on the galactic periphery.

The position of the boundary between HI and HII layers, $`z_i`$, in an exponential vertical gas distribution ionized by a given flux of external UV photons $`I`$ is determined by the following equation $`2\pi I=\alpha _rn_0^2_{z_i}^{\mathrm{}}\mathrm{exp}(2z/h)𝑑z`$. The observed column density of atomic hydrogen, $`N(\mathrm{HI})=2n_0_0^{z_i}\mathrm{exp}(z/h)𝑑z`$, can then be found as

$$N(\mathrm{HI})=N_0(\mathrm{H})\sqrt{\frac{4\pi Ih}{\alpha _r}},$$
(10)
where $`\alpha _r`$ is the hydrogen recombination rate and $`N_0(\mathrm{H})=2n_0h`$ is the total column density of hydrogen (both neutral and ionized), $`n_0`$ is the total midplane density of hydrogen. It is readily seen that for a typical value of the background ionizing flux of $`I10^6`$ photons cm-2 s-1 sr-1 and for $`h300`$ pc the column density of ionized hydrogen can reach $`3\times 10^{20}`$ cm-2, nearly coincident with the observed in Ho II at the radii of $`6`$ arcmin (Puche et al. Puche (1992)). This implies that the observed $`N(\mathrm{HI})`$ may represent a small fraction of the total amount of hydrogen, and thus a relatively weak perturbation in $`n_0`$ and/or in $`N_0(\mathrm{H})`$ can cause a rather strong variation in the amount of neutral hydrogen. This effect increases for a flaring gaseous disk with the scale height $`h`$ growing outwards. In Fig. 6 we plot the total column density of hydrogen, $`N_0(\mathrm{H})`$, as a function of Ho II radius. $`N_0(\mathrm{H})`$ is derived from the observed column density of neutral hydrogen $`N(\mathrm{HI})`$ for two models of the vertical gas distribution in Ho II: a plane-parallel distribution with $`h`$=625 pc and a flaring disk with $`h=625+18.6(R1)^2`$ pc, where $`R`$ is given in arcminutes.

In the outer regions at $`R\stackrel{>}{}6`$ arcmin, the observed $`N(\mathrm{HI})`$ may represent a relatively small (less than 50 %) fraction of the total hydrogen content, and thus can be easily disturbed, ionized, and transformed into a hole of depleted HI column density.

In principle, HI holes created by multiple SN explosions can become radially asymmetric due to radial variations in the gas density, particularly at larger distances from the galactic center where the flaring of the gas disk may dominate. For an exponential disk with the density distribution $`\rho \mathrm{exp}(z/hR/R_0)`$, the expansion velocity of a shell at a given point varies as $`v_s\mathrm{exp}(z/2h+R/2R_0)`$. In this case, the geometrical center of an HI shell can be shifted outward with respect to the position of SN explosions (Kalberla, private communication). Typical scales in vertical and radial directions are of the order of hundreds and thousands of parsecs, respectively. For example, for the Ho II galaxy $`h600`$ pc, while the radial scale for the HI surface density $`R_03`$ kpc (Puche et al. 1992, Bureau & Carignan, 2002). For face-on galaxies, direct measurements of the radial distributions of the volume density or flaring are not available, but one can assume that the radial scale lengths for them are of the same order as for the HI surface density. The corresponding deviations of the expansion velocities in radial $`\delta v_R`$ and vertical $`\delta v_z`$ directions from a precise spherical shape can be estimated then as $`\delta v_R/\delta v_zh/R_01`$. In the case of Ho II, it becomes $`\delta v_R/\delta v_z0.2`$. Hence, one may conclude that the expected radial asymmetry of a shell, and the corresponding displacement of the geometrical center of a shell with respect to the location of SN explosions, is less than $`20\%`$ of the gas scale height.

## 4  Summary

Clustered supernova explosions in the disk of a face-on galaxy produce an expanding shell of compressed material, which is seen as an HI ring surrounding a central HI depression. In this paper we show that several physical characteristics of expanding shells are sensitive to the gas scale height and, hence, can be used for the determination of the dark matter content in face-on disk galaxies.

$``$ A functional dependence of the relative column density of gas in the HI ring ($`𝒦_{\mathrm{HI}}`$) on the radius of the ring is found to be specific for a given gas scale height $`h`$, which makes it possible to infer $`h`$ from observations of even a single HI ring. This is particularly true for the shells near breakthrough phase.

$``$ The vertical component of gas velocity in an expanding shell reveals a characteristic behaviour (particularly, at the stages close to a breakthrough), which can be used to trace the gas scale height $`h`$ in galaxies with a sufficiently thin HI layer.

$``$ Out-of-plane SN explosions produce asymmetric expansion, with a more massive part of the shell propagating towards the denser regions of the disk. The velocity of this part of the shell is a single-valued function of the scale height $`h`$ and the height above the midplane at which the SN explosions take place.

$``$ Some of HI holes, particularly in the outer regions of a galaxy, can be connected to the ionization of HI layer by extragalactic UV photons.

## Acknowledgments

We would like to thank the referee, Dr. P. M. W. Kalberla, for his suggestions and critical comments that substantially improved the paper. This work was supported by the RFBR (projects No 00-02-17689) and the INTAS grant YSF-2002-33. YS acknowledges financial support from *Deutsche Forschungsgemeinschaft, DFG* (project SFB N 591, TP A6).