# Evidential strength of categorical expert witness statements

R.J.F.Ypma

In many jurisdictions, forensic evidence is presented in the form of
categorical statements by forensic experts. In recent years, the scientific
validity of forensic science practices have been questioned. In response,
several large-scale performance studies have been performed, reporting various
error rates for different forensic fields. This important work has elucidated
the uncertainty associated with categorical statements.

There is growing scientific consensus that the likelihood ratio (LR) framework
is the logically correct form of presentation for forensic evidence
evaluation. Yet, the very relevant results from the large-scale performance
studies have not been cast in this framework.

Here, I show how to straightforwardly calculate an LR for any given
categorical statement using data from the performance studies. This number
quantifies how much more we should believe the hypothesis of same source vs
different source, when provided a particular expert witness statement. LRs are
reported for categorical statements resulting from the analysis of latent
fingerprints, bloodstain patterns, handwriting, footwear and firearms. The
highest LR found for statements of identification was 376 (fingerprints), the
lowest found for statements of exclusion was 1/28 (handwriting).

The LRs found may be more insightful for those used to this framework than the
various error rates reported previously. An additional advantage of using the
LR in this way is the relative simplicity; there are no decisions necessary on
what error rate to focus on or how to handle inconclusive statements. The
values found are closer to 1 than many would have expected. One possible
explanation for this mismatch is that we undervalue numerical LRs. Finally, a
note of caution: the LR values reported here come from a simple calculation
that does not do justice to the nuances of the large-scale studies and their
differences to casework, and should be treated as ball-park figures rather
than definitive statements on the evidential value of whole forensic
scientific fields.

##  1 Introduction

Forensic evidence plays a large role in criminal law. In many jurisdictions,
expert witnesses give categorical statements, i.e. claiming identification of
a person or object from forensic traces. In recent decades forensic sciences
have been criticized for ignoring or being unable to quantify the uncertainty
associated with such statements [1, 2]. In response, several large-scale
studies have been published that investigate this uncertainty for a range of
forensic fields [3, 4, 5, 6, 7]. In these studies, uncertainty is quantified
as error rates. This is however not a trivial computation, as several relevant
error rates exist (e.g. false negative rate, false positive rate, false
discovery rate) and there is no single correct way of handling statements of
the form ‘unsuitable for analysis’ or ‘inconclusive’ [8].

The likelihood ratio (LR) framework provides an alternative form of reporting
evidence, and is widely recommended as logically correct [9, 10, 11, 8]. The
likelihood ratio is defined as the likelihood of the evidence $E$ under one
hypothesis ($H_{1}$) divided by the likelihood under a competing hypothesis
($H_{2}$):

$LR=P(E|H_{1})/P(E|H_{2})$

Adoption of numerical LRs is largest in the field of forensic DNA analysis,
which boasts well-developed statistical models. However, an increasing number
of forensic fields is following suit in developing numerical LR systems [12].
The framework is equally applicable to subjective judgements, e.g. taking the
form ‘the evidence provides strong support for the first proposition relative
to the alternative’.

No trier-of-fact calculates a probability of guilt. Rather, there is a
psychological process that incorporates the evidence and results in a
conviction of proven or not-proven. It is unclear whether this process takes
into account uncertainty, expressed by error rates or likelihood ratios, in
the way that the forensic scientist intended. Most likely, it does not [13,
14, 15, 16, 17]. It seems reasonable that training in and exposure to
particular forms of expressing evidence will improve interpretation [18].
Thus, for those used to the likelihood ratio framework, it may be of use to
quantify the uncertainty measured in the performance studies in the form of an
LR.

Below, I will illustrate how likelihood ratios for expert witness statements
can be straightforwardly computed from performance studies. In fact, if only
the two statements of ‘same source’ or ‘different source’ are used, the LR is
simply (1 - false negative rate) / false positive rate. The interpretation of
these calculated LRs is the support for same source relative to different
source that is given by a particular statement, from an expert and sample pair
comparable to those in the performance study. None of these numbers should be
treated as a thorough assessment of the evidential strength of a statement in
the context of a specific case.

##  2 Methods

We assume we have access to many evaluations made by a set of experts for
which ground truth is known. Scott et al. [4] provides such data for firearm
examiners. Possible conclusions for the examiner to reach were
individualisation (‘ID’), three levels of inconclusive (‘Inconcl. A-C’),
‘Elimination’ and ‘Other’. For example, for the statement of individualisation
the LR would be:

$LR=\frac{P(\textrm{expert concludes bullet fired from firearm}|\textrm{bullet
was fired from firearm})}{P(\textrm{expert concludes bullet fired from
firearm}|\textrm{bullet was not fired from firearm})}$

Table 1 shows how often conclusions were reached for matching and non-matching
pairs.

The LR can be calculated in a straightforward manner. The probability in the
numerator is equal to the number of times that an expert said ‘ID’ when
$H_{1}$ was true, i.e. the bullet was fired from the firearm, divided by the
total number of times $H_{1}$ was true. The first is 1076, the second is
1076+127+125+36+41+24 = 1429. This gives us a probability of

$\begin{split}&P(\textrm{expert concludes bullet fired from
firearm}|\textrm{bullet was fired from firearm})\\\ &=1076/1429\\\ &\approx
0.75\end{split}$

Likewise for the denominator, we get

$\begin{split}&P(\textrm{expert concludes bullet fired from
firearm}|\textrm{bullet was not fired from firearm})\\\ &=20/2891\\\ &\approx
0.007\end{split}$

The ratio of these probabilities is the likelihood ratio: $LR\approx
0.75/0.007\approx 109$. This calculation is exactly the same for all studies
we looked at. Sometimes an additional tallying step is needed as the data are
not given in aggregated form. Code is available from
github111github.com/NetherlandsForensicInstitute/evidential_value_of_expert_witness_statements.

| ID | Inconcl.-A | Inconcl.-B | Inconcl.-C | Elimination | Other
---|---|---|---|---|---|---
Matching | 1,076 | 127 | 125 | 36 | 41 | 24
Non-matching | 20 | 268 | 848 | 745 | 961 | 49
Table 1: Statistics on bullet evaluations by expert witnesses, from [4].

##  3 Results

Table 2 gives the LRs for all categories considered for bullet evaluations.
LRs for all different statement categories considered for the other
performance studies (footwear, cartridge, handwriting, fingerprints,
bloodstain patterns) are given in the tables in the appendix. Table 3
summarises these by giving for every study the LR calculated for the
identification and exclusion statements.

| ID | Inconcl.-A | Inconcl.-B | Inconcl.-C | Elimination | Other
---|---|---|---|---|---|---
LR | 109 | 1 | 1 / 3 | 1 / 10 | 1 / 12 | 1
Table 2: Likelihood ratios for bullet evaluation conclusions by expert witnesses, calculated from data in [4]. | LR (identification) | LR (exclusion)
---|---|---
pattern type from bloodstain pattern [3] | 6 | 1/4
writer from handwriting sample [7] | 17 | 1/28
firearm from cartridge [4] | 81 | 1/28
firearm from bullet [4] | 109 | 1/12
footwear from print [6] | 113 | 1/7
person from latent fingerprints [5] | 376 | 1/11
Table 3: Likelihood ratios for expert witness from various fields, for
statements of identification and exclusion.

##  4 Discussion

This study presents the likelihood ratios for expert witness statements for
five forensic fields from large performance studies. The LRs are an
alternative to error rates for quantifying the uncertainty associated with
these statements. This alternative may be useful to those who are more used to
the LR framework than error rates, and may help bridge the gap between the two
frameworks. Incidentally, the LR framework simplifies the uncertainty analysis
as we don’t have to consider various different error rates and ways of
handling ‘inconclusive’ statements.

The likelihood ratios reported should be treated as ‘ballpark’ figures, giving
some insight into what LRs to associate with particular statements, such as
identification. We should be very hesitant in interpreting the LRs reported
here as the evidential strength of an expert witness statement in a given
case. The performance studies measure the performance of a set of experts,
which may not be representative for the case, on a specific set of samples,
which may not be representative for the case. For example, several studies
explicitly select different-source samples that are very challenging to
evaluate [7, 5, 4]. As an illustration, considering only the top 1% hardest
different-source comparisons would result in an LR 100 times lower (assuming
no errors would be made on the remaining 99%). Furthermore, it is common
practice in casework to have forensic witness evaluations reviewed by an
independent expert - this extra quality check is absent in the performance
studies. Lastly, we stress that evidential strength is only one aspect of the
value of forensic evidence; moderately strong evidence that a suspect held the
murder weapon may be much more relevant than very strong evidence that he was
near the crime scene.

The actual values of LRs found are interesting. As expected, statements of
increasing certainty result in higher LRs and ‘inconclusive’ statements
generally result in LRs close to 1. More surprising is that even the highest
LRs found for inclusion were below 400, and the lowest found for exclusion
above 1/50. This range is small. For example, both experts and the public
expect a statement of fingerprint identification to be stronger than an LR of
100,000 [19]. The caveats mentioned above only partly explain the difference
of several orders of magnitude between expected and calculated LRs. It seems
we either overvalue categorical statements, undervalue numerical LRs, or both.
This undervaluation is something we anecdotally encounter in practice, and is
perhaps partly caused by the forensic science community. After all, an LR of
1000 is strong enough to convert a weak prior probability of 0.10 to a
posterior of 0.99. Indeed, an LR of 1000 is referred to as ‘decisive’ in
general scientist’s nomenclature [20, 21]. Yet in forensic science we call it
‘moderately strong’ [9, 10], and in DNA analysis this number may even be
considered too low to report [22]. Clearly, more research is needed in how
‘mid range’ LRs are actually interpreted, and how we can improve this.

## Acknowledgements

My sincere thanks to Marjan Sjerps, Wauter Bosma, David vd Vloed and Charles
Berger for critical reading of this manuscript.

## References

  * [1] President’s Council of Advisors on Science and Technology.  Forensic science in criminal courts: Ensuring scientific validity of feature-comparison methods.  Technical report, 2016.
  * [2] National Research Council, Division on Engineering and Physical Sciences, Committee on Applied and Theoretical Statistics, Policy and Global Affairs, Committee on Science, Technology, and Law, and Committee on Identifying the Needs of the Forensic Sciences Community.  Strengthening Forensic Science in the United States: A Path Forward.  National Academies Press, August 2009.
  * [3] R Austin Hicklin, Kevin R Winer, Paul E Kish, Connie L Parks, William Chapman, Kensley Dunagan, Nicole Richetelli, Eric G Epstein, Madeline A Ausdemore, and Thomas A Busey.  Accuracy and reproducibility of conclusions by forensic bloodstain pattern analysts.  Forensic Sci. Int., 325:110856, August 2021.
  * [4] L Scott Chumbley, Max D Morris, Stanley J Bajic, Daniel Zamzow, Erich Smith, Keith Monson, and Gene Peters.  Accuracy, repeatability, and reproducibility of firearm comparisons part 1: Accuracy.  July 2021.
  * [5] Bradford T Ulery, R Austin Hicklin, Joann Buscaglia, and Maria Antonia Roberts.  Accuracy and reliability of forensic latent fingerprint decisions.  Proc. Natl. Acad. Sci. U. S. A., 108(19):7733–7738, May 2011.
  * [6] R. Austin Hicklin, Brian C. McVicker, Connie Parks, Jan LeMay, Nicole Richetelli, Michael Smith, JoAnn Buscaglia, Rebecca Schwartz Perlman, Eugene M. Peters, and Brian A. Eckenrode.  Accuracy, reproducibility, and repeatability of forensic footwear examiner decisions.  Forensic Science International, 339:111418, 2022.
  * [7] R Austin Hicklin, Linda Eisenhart, Nicole Richetelli, Meredith D Miller, Peter Belcastro, Ted M Burkes, Connie L Parks, Michael A Smith, Joann Buscaglia, Eugene M Peters, Rebecca Schwartz Perlman, Jocelyn V Abonamah, and Brian A Eckenrode.  Accuracy and reliability of forensic handwriting comparisons.  Proc. Natl. Acad. Sci. U. S. A., 119(32):e2119944119, August 2022\.
  * [8] Jonathan J Koehler, Jennifer L Mnookin, and Michael J Saks.  The scientific reinvention of forensic science.  Proc. Natl. Acad. Sci. U. S. A., 120(41):e2301840120, October 2023\.
  * [9] Association of Forensic Science Providers.  Standards for the formulation of evaluative forensic science expert opinion.  Sci. Justice, 49(3):161–164, September 2009.
  * [10] ENFSI (European network of forensic science institutes).  ENFSI guideline for evaluative reporting in forensic science, 2015.
  * [11] American Statistical Association.  American statistical association position on statistical statements for forensic evidence, January 2019.
  * [12] R van Lierop, Ramos D, Sjerps M, and Ypma RJF.  An overview of log likelihood ratio cost in forensic sciences – where is it used and what values can we expect?  submitted.
  * [13] Raymond Marquis, Alex Biedermann, Liv Cadola, Christophe Champod, Line Gueissaz, Geneviève Massonnet, Williams David Mazzella, Franco Taroni, and Tacha Hicks.  Discussion on how to implement a verbal scale in a forensic laboratory: Benefits, pitfalls and suggestions to avoid misunderstandings.  Sci. Justice, 56(5):364–370, September 2016.
  * [14] William C Thompson and Eryn J Newman.  Lay understanding of forensic statistics: Evaluation of random match probabilities, likelihood ratios, and verbal equivalents.  Law Hum. Behav., 39(4):332–349, August 2015.
  * [15] Marjan Sjerps and Dirk B Biesheuvel.  The interpretation of conventional and ’bayesian’ verbal scales for expressing expert opinion: a small experiment among jurists.  Int. J. Speech Lang. Law, 6(2):214–227, August 1999.
  * [16] Heidi Eldridge.  Juror comprehension of forensic expert testimony: A literature review and gap analysis.  Forensic Sci Int Synerg, 1:24–34, March 2019.
  * [17] Jan de Keijser and Henk Elffers.  Understanding of forensic expert reports by judges, defense lawyers and forensic professionals.  Psychol. Crime Law, 18(2):191–207, February 2012.
  * [18] Eric Rassin, Nurul Arbiyah, Irena Boskovic, Henry Otgaar, and Harald Merckelbach.  Likelihood ratios in psychological expert opinion, and their reception by professional judges.  The International Journal of Evidence & Proof, 26(4):325–341, October 2022.
  * [19] Thomas Busey and Morgan Klutzke.  Calibrating the perceived strength of evidence of forensic testimony statements.  Sci. Justice, 63(1):38–53, January 2023.
  * [20] Jeffreys H.  Theory of Probability (3rd Ed.).  Oxford University Press, 1961.
  * [21] Andrew F Jarosz and Jennifer Wiley.  What are the odds? a practical guide to computing and reporting bayes factors.  The Journal of Problem Solving, 7(1):2, 2014.
  * [22] Meinhard Hahn, Katja Anslinger, Martin Eckert, Rolf Fimmers, Stefanie Grethe, Carsten Hohoff, Sebastian Kranz, Christoph Leuker, Claus Oppelt, Sven Razbin, Thomas Rothämel, Harald Schneider, Michael Templin, Marielle Vennemann, Andrea Wächter, Volker Weirich, Peter Zimmermann, and Peter M Schneider.  Gemeinsame empfehlungen der projektgruppe ‘biostatistische dna-berechnungen’ und der spurenkommission zur biostatistischen bewertung forensischer DNA-analytischer befunde mit vollkontinuierlichen modellen (VKM).  Rechtsmedizin, 33(1):3–12, February 2023.

certainty of bloodstain pattern type classification | LR
---|---
identification (‘definitive’) | 6
possible (‘included’) | 1
excluded | 1 / 4
Table 4: Likelihood ratios for bloodstain pattern evaluation conclusions by expert witnesses.[3] | LR
---|---
The questioned sample was written by the known writer | 17
The questioned sample was probably written by the known writer | 7
No conclusion | 1 / 2
The questioned sample was probably not written by the known writer | 1 / 22
The questioned sample was not written by the known writer | 1 / 28
Table 5: Likelihood ratios for handwriting evaluation conclusions by expert witnesses.[7] | LR
---|---
identification | 113
high association | 13
association | 1
limited association | 1
no association | 1 / 6
exclusion | 1 / 7
inconclusive | 1 / 2
not suitable | 1 / 2
Table 6: Likelihood ratios for footwear print evaluation conclusions by expert witnesses.[6] | LR
---|---
identification | 81
inconclusive A | 2
inconclusive B | 1 / 2
inconclusive C | 1 / 14
elimination | 1 / 28
other | 1
Table 7: Likelihood ratios for cartridge case evaluation conclusions by expert witnesses.[4] | LR
---|---
individualisation | 376
inconclusive | 2
exclusion | 1 / 11
Table 8: Likelihood ratios for latent fingerprint evaluation conclusions by
expert witnesses.[5]


040005 2012 M. C. Barbosa 040005

With the aim of establishing a criterion for identifying when a lipid bilayer
has reached steady state using the molecular dynamics simulation technique,
lipid bilayers of different composition in their liquid crystalline phase were
simulated in aqueous solution in presence of CaCl2 as electrolyte, at
different concentration levels. In this regard, we used two different lipid
bilayer systems: one composed by 288 DPPC (DiPalmitoylPhosphatidylCholine) and
another constituted by 288 DPPS (DiPalmitoylPhosphatidylSerine). In this
sense, for both type of lipid bilayers, we have studied the temporal evolution
of some lipids properties, such as the surface area per lipid, the deuterium
order parameter, the lipid hydration and the lipid-calcium coordination. From
their analysis, it became evident how each property has a different time to
achieve equilibrium. The following order was found, from faster property to
slower property: coordination of ions $\approx$ deuterium order parameter $>$
area per lipid $\approx$ hydration. Consequently, when the hydration of lipids
or the mean area per lipid are stable, we can ensure that the lipid membrane
has reached the steady state.

# A criterion to identify the equilibration time in lipid bilayer simulations

Rodolfo D. Porasso  [inst1] J.J. López Cascales E-mail<EMAIL_ADDRESS>[inst2]

(3 September 2012; 24 October 2012)

††volume: 4

99 inst1 Instituto de Matemática Aplicada San Luis (IMASL) - Departamento de
Física, Universidad Nacional de San Luis/CONICET, D5700HHW, San Luis,
Argentina. inst2 Universidad Politécnica de Cartagena, Grupo de Bioinformática
y Macromoléculas (BioMac) Aulario II, Campus de Alfonso XIII, 30203 Cartagena,
Murcia, Spain.

##  1 Introduction

Over the last few decades, different computational techniques have emerged in
different fields of science, some of them being extensively implemented and
used by a great number of scientists around the globe. Among others, the
Molecular Dynamics (MD) simulation is a very popular computational technique,
which is widely used to obtain insight with atomic detail of steady and
dynamic properties in the fields of biology, physics and chemistry. In this
regard, a critical aspect that must be identified in all the MD simulations is
related to the required equilibration time to achieve a steady state. This
point is crucial in order to avoid simulation artifacts that could lead to
wrong conclusions. Currently, with the increment of the computing power
accessible to different investigation groups, much longer simulation
trajectories are being carried out to obtain reliable information about the
systems, with the purpose of approaching the time scale of the experimental
phenomena. However, even when this fact is objectively desirable without
further objections, nowadays, much longer equilibration times are arbitrarily
being required by certain reviewers during the revision process. From our
viewpoint, this should be thoroughly revised due to the following two main
reasons: first, because it results in a limiting factor in the use of this
technique by other research groups which cannot access to very expensive
computing centers (assuming that authors provide enough evidence of the
equilibration of the system). Second, to avoid wasting expensive computing
time in the study of certain properties which do not require such long
equilibration times, once the steady state of the system has been properly
identified.

Phospholipid bilayers are of a high biological relevance, due to the fact that
they play a crucial role in the control of the diffusion of small molecules,
cell recognition, and signal transduction, among others. In our case, we have
chosen the PhosphatidylCholine (PC) bilayer because it has been very well
studied by MD simulations [1, 2, 3, 4, 5, 6, 7] and experimentally as well [8,
9, 10, 11, 12, 13, 14]. Furthermore, studies of the effects of different types
of electrolytes on a PC bilayer have also been studied, experimentally [15,
16, 17, 18, 19, 20, 21, 22, 23, 24] and by simulation [25, 26, 27, 28, 29, 30,
31, 32, 33].

As mentioned above, the Molecular Dynamics (MD) simulations have emerged
during the last decades as a powerful tool to obtain insight with atomic
detail of the structure and dynamics of lipid bilayers [34, 35, 36]. Several
MD simulations of membranes under the influence of different salt
concentrations have been carried out. One of the main obstacles related to
these studies has been the time scale associated to the binding process of
ions to the lipid bilayer. Considering the literature, a vast dispersion of
equilibration times associated to the binding of ions to the membrane has been
reported, where values ranging from 5 to 100 ns have been suggested for
monovalent and divalent cations [25, 27, 28, 29, 32, 37, 38]. In this regard,
we carried out four independent simulations of a lipid bilayer formed by 288
DPPC in aqueous solutions, for different concentrations of CaCl2 to provide an
overview of their equilibration times. Among other properties, the surface
area per lipid, the deuterium order parameters, lipid hydration and lipid-
calcium coordination were studied.

Finally, in order to generalize our results, a bilayer formed by 288 DPPS in
its liquid crystalline phase, in presence of CaCl2 at 0.25N, was simulated as
well .

##  2 Methodology

Different molecular Dynamics (MD) simulations of lipid bilayer formed by 288
DPPC were carried out in aqueous solutions for different concentrations of
CaCl2, from 0, up to 0.50 N. Furthermore, with the aim of generalizing our
results, a bilayer of 288 DPPS in presence of CaCl2 at 0.25 N was simulated as
well. Note that the concentration of CaCl2 in terms of normality is defined
as:

$\textrm{normality}=\frac{n_{\textrm{equivalent
grams}}}{l_{\textrm{solution}}}$ (1)

where $n_{\textrm{equivalent
grams}}=\frac{\textrm{gr(solute)}}{\textrm{equivalent weight}}$ and
$\textrm{equivalent weight}=\frac{\textrm{Molecular weight}}{n}$, being $n$
the charge of the ions in the solution.

In Table 1, the number of molecules that constitute each system, applying Eq.
(1), is summarized.

Type of Lipid | [CaCl2] N | Ca2+ | Cl- | Water
---|---|---|---|---
DPPC | 0 | 0 | 0 | 10068
DPPC | 0.06 | 5 | 10 | 10053
DPPC | 0.13 | 12 | 24 | 10032
DPPC | 0.25 | 23 | 46 | 9999
DPPC | 0.50 | 46 | 92 | 9930
DPPS | 0.25 | 204 | 120 | 26932
Table 1: The simulated bilayer systems. Note that the salt concentration is
given in normal units. The numerals describe the number of molecules contained
in the simulation box.

To build up the original system, a single DPPC lipid molecule, or DPPS lipid
(Fig. 1), was placed with its molecular axis perpendicular to the membrane
surface ($xy$ plane). Next, each DPPC, or DPPS, was randomly rotated and
copied 144 times on each leaflet of the bilayer. Finally, the gaps existing in
the computational box (above and below the phospholipid bilayer) were filled
using an equilibrated box containing 216 water molecules of the extended
simple point charge (SPC/E) [39] water model.

Figure 1: Structure and atom numbers for DPPC and DPPS used in this work.

Thus, the starting point of the first system of Table 1 was formed by 288 DPPC
in absence of CaCl2. Once this first system was generated, the whole system
was subjected to the steepest descent minimization process to remove any
excess of strain associated with overlaps between neighboring atoms of the
system. Thereby, the following DPPC systems in presence of CaCl2 were
generated as follows: to obtain a [CaCl${}_{2}]$=0.06 N, 15 water molecules
were randomly substituted by 5 Ca2+ and 10 Cl-. An analogous procedure was
applied to the rest of the systems, where 36, 69 and 138 water molecules were
substituted by 12, 23 and 46 Ca2+ and 24, 46 and 92 Cl-, to obtain a
[CaCl${}_{2}]$ concentration of 0.13 N, 0.25 N and 0.50 N, respectively.
Finally, the DPPS bilayer was generated following the same procedure described
above for the DPPC, starting from a single DPPS molecule and once the lipid
bilayer in presence of water passed the minimization process, 324 water
molecules were substituted by 204 Ca2+ and 120 Cl- (note that 144 of the 204
calcium ions were added to balance the negative charge associated with the
DPPS).

GROMACS 3.3.3 package [40, 41] was used in the simulations, and the properties
showed in this work were obtained using our own code. The force field proposed
by Egberts et al. [2] was used for the lipids, and a time step of 2 fs was
used as integration time in all the simulations. A cut-off of 1.0 nm was used
for calculating the Lennard-Jones interactions. The electrostatic interaction
was evaluated using the particle mesh Ewald method [42, 43]. The real space
interaction was evaluated using a 0.9 nm cut-off, and the reciprocal space
interaction using a 0.12 nm grid with a fourth-order spline interpolation. A
semi-isotropic coupling pressure was used for the coupling pressure bath, with
a reference pressure of 1 atm which allowed the fluctuation of each axis of
the computer box independently. For the DPPC bilayer, each component of the
system (i.e., lipids, ions and water) was coupled to an external temperature
coupling bath at 330 K, which is well above the transition temperature of 314
K [44, 45]. For DPPS bilayer, each component of the system was coupled to an
external temperature coupling bath at 350 K, which is above the transition
temperature [46, 47]. All the MD simulations were carried out using periodic
boundary conditions. The total trajectory length of each simulated system was
of 80 ns of MD simulation, where the coordinates of the system were recorded
every 5 ps for their appropriate analysis.

Finally, in order to study the effect of the temperature, only the case
corresponding to 0.25 N CaCl2 was investigated at two additional temperatures,
340 K and 350 K.

##  3 Results and discussion

###  3.1 Effect of the CaCl2 concentration

####  3.1.1 Surface area per lipid

Surface area per lipid $\langle A\rangle$ is a property of lipid bilayers
which has been accurately measured from experiments [48]. The calculation of
mean area per lipid can be determined from the MD simulation as:

$\langle A\rangle=\frac{x\cdot y}{N}$ (2)

where $x$ and $y$ represent the box sizes in the direction $x$ and $y$
(perpendicular to the membrane surface) over the simulation, and $N$ is the
number of lipids contained in one leaflet, in our case $N=144$.

Figure 2: Running area per lipid at T = 330 K in presence of [CaCl2] at (A)
0.06 N, (B) 0.13 N, (C) 0.25 N, (D) 0.50 N and (E) 0.25 N (in this case, T =
350 K). Solid lines represent the mean area obtained from the last 70 ns of
the simulated trajectories (see text for further explanation). The type of
lipid is indicated in the legends.

Focusing on the study of the time evolution of the area per lipid, Figure 2
depicts the running surface area per lipid for different concentrations of
CaCl2 and type of lipid. In general, for the 5 bilayers formed by DPPC or
DPPS, the area per lipid achieved a steady state after 10 ns of simulation,
being this equilibration time almost independent of the concentration of CaCl2
and type of lipid which composed the membrane.

In absence of salt, an average area per lipid of $\langle A\rangle=0.663\pm
0.008$ nm2 was calculated from the last 70 ns of the simulated trajectory,
discarding the first 10 ns corresponding to the equilibration time. This value
agrees with experimental data, where values in a range from 0.55 to 0.72 nm2
have been measured [10, 11, 48, 49, 50, 51]. Table 2 shows the mean surface
area per lipid (again, after discarding the equilibration time of 10 ns) with
their corresponding error bar.

From the simulation results, a shrinking in the surface area per lipid with
the increment of the ionic strength of the solution is observed. This
shrinking is expected and attributed to the complexation of lipid molecules by
calcium, such as it has been pointed out in previous studies [28, 29, 52].

####  3.1.2 Deuterium order parameter

The deuterium order parameter, $S_{CD}$, is measured from 2H-NMR experiments.
This parameter provides relevant information related to the disorder of the
hydrocarbon region in the interior of the lipid bilayers by measuring the
orientation of the hydrogen dipole of the methylene groups with respect to the
perpendicular axis to the lipid bilayer. Due to the fact that hydrogens of the
lipid methylene groups (CH2) have not been taken into account (in an explicit
way) in our simulations, the order parameter $-S_{CD}$ on the $i+1$ methylene
group was defined as the normal unitary vector to the vector defined from the
$i$ to the $i+2$ CH2 group and contained in the plane formed by the methylene
groups $i$, $i+1$ and $i+2$. Thus, the deuterium order parameter $-S_{CD}$ on
the $i-th$ of the CH2 group can be estimated by Molecular Dynamics simulations
as follows:

$-S_{CD}=\frac{1}{2}\langle{3\cos^{2}(\theta)-1}\rangle$ (3)

where $\theta$ is the angle formed between the unitary vector defined above
and the $z$ axis. The expression in brackets $\langle\dots\rangle$ denotes an
average over all the lipids and time. Hence, note that the $-S_{CD}$ can adopt
any value between -0.5 (corresponding to a parallel orientation to the
lipid/water interface) and 1 (oriented along the normal axis to the lipid
bilayer).

Type of | [CaCl2] N |  $\langle A\rangle$ (nm2) | Hydration
---|---|---|---
Lipid |  |  | Number
DPPC | 0 | 0.663 $\pm$0.008 | 1.758 $\pm$0.009
DPPC | 0.06 | 0.658 $\pm$0.008 | 1.740 $\pm$0.009
DPPC | 0.13 | 0.651 $\pm$0.007 | 1.719 $\pm$0.010
DPPC | 0.25 | 0.641 $\pm$0.009 | 1.680 $\pm$0.015
DPPC | 0.50 | 0.628 $\pm$0.010 | 1.610 $\pm$0.015
DPPS | 0.25 | 0.522 $\pm$0.007 | 2.552 $\pm$0.010
Table 2: Area per lipid and lipid hydration number as a function of salt
concentration (see text for further explanation). Note that the salt
concentration is given in normal units. Error bars were calculated for each
system separately from subtrajectories of 10 ns length. Simulation temperature
= 330 K.

Figure 3 shows the running $-S_{CD}$ for different carbons of the DPPC and
DPPS tails and salt concentrations. Only the carbons which correspond to the
initial (hydrocarbons 2 and 6), the middle (hydrocarbon 10) and final
(hydrocarbons 13 and 15) methylene groups of the lipid tails were depicted in
this figure. Each point of the figure represents the average values of
$-S_{CD}$ over 5 ns of subtrajectory length, and the lines represent the mean
values calculated from the last 70 ns of the trajectories simulated. From this
figure, it is observed how in all the cases, the required equilibration time
is less than 10 ns of simulation time, independently of the salt concentration
and the type of lipid. Finally, it is noted that Figure 3 exhibits an increase
in the deuterium order parameters with the salt concentration, consistent with
the shrinking of the area per lipid described above.

Figure 3: Running deuterium order parameter, $-S_{CD}$, in presence of [CaCl2]
at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N, (D) 0.50 N and (E) 0.25 N. DPPC
simulations were performed at 330 K and DPPS simulation temperature was 350 K.
Solid lines represent mean values $-S_{CD}$ obtained from the last 70 ns of
the simulated trajectories. The type of lipid is indicated in the legends.
Symbols: $\circ$ hydrocarbon 2; $\diamond$ hydrocarbon 6; $\triangleleft$
hydrocarbon 10; $+$ hydrocarbon 13 and $\times$ hydrocarbon 15. Note that the
error bars have the same size as symbol.

####  3.1.3 Lipid hydration

To analyze the lipid hydration, the radial distribution function $g(r)$ of
water around one of the oxygens of the phosphate group (atom number 10 in Fig.
1 for DPPC and DPPS) was calculated. The radial distribution function
$g\left(r\right)$ is defined as follows:

$g(r)=\frac{N(r)}{4{\pi}r^{2}{\rho}\delta r}$ (4)

where $N\left(r\right)$ is the number of atoms in a spherical shell at
distance $r$ and thickness $\delta r$ from a reference atom. $\rho$ is the
density number taken as the ratio of atoms to the volume of the total
computing box.

From numerical integration of the first peak of the radial distribution
function, the hydration numbers can be estimated for different atoms of the
DPPC or DPPS. Figure 4 depicts the hydration number of phosphate oxygen (atom
10 in Fig. 1 for DPPC and DPPS) in presence of CaCl2, where each point
represents the average of 5 ns subtrajectory length. These results show how
this property reached a steady state for the cases (A), (B) and (E), after 10
ns of simulation. However, for the cases (C) and (D), 5 ns of extra simulation
trajectory were required to reach a steady state. Table 2 shows the hydration
numbers for the last 70 ns of the trajectory length, corresponding to four
concentrations of CaCl2 and both types of lipids, DPPC and DPPS. In this
regard, from Fig. 4, the significant lipid dehydration with the increment of
the ionic strength of the solution is evident, in good accordance with
previous results [52].

Figure 4: Hydration number of the phosphate oxygen (atom 10 in Fig. 1) along
the simulated trajectories in presence of [CaCl2] at (A) 0.06 N, (B) 0.13 N,
(C) 0.25 N, (D) 0.50 N (for DPPC T = 330 K) and (E) 0.25 N. In this case, T =
350 K. Solid lines represent the mean value of the hydration number calculated
from the last 70 ns of the simulated trajectories. The type of lipid is
indicated in the legends. Note that the error bars have the same size as
symbol.

####  3.1.4 Phospholipid-calcium coordination

Some authors have reported how the lipid coordination by divalent cations
widely varies . Thus, on the one hand, some authors [25] have reported that
this is a very slow process, which requires about 85 ns of simulation time,
but, on the other hand, other authors [26] have suggested that this process
results much more rapid, taking less than 1 ns. In this sense, the
coordination of DPPC-Ca2+ was studied by monitoring the oxygen-calcium
coordination of the carbonyl oxygens (atoms 16 and 35 in Fig. 1) and phosphate
oxygens (atoms 9 and 10 of DPPC in Fig. 1), as a function of time. The left
column in Fig. 5 represents the oxygen coordination number, while the right
one depicts the percentage of calcium ions involved in the coordination
process with respect to the total number of calcium ions present in the
aqueous solution. Figure 5 shows how the DPPC coordination by calcium is a
quick process, taking less than 5 ns of simulation time to achieve a steady
state. The kinetic of this process appears to be related to the ratio between
calcium/lipid. Thus, after the first 5 ns of simulation time, the Ca–lipid
coordination presents some fluctuation along the rest of the simulated
trajectory. However, in Fig. 5 (A) and (B) (for the cases of lower
concentration), it is observed how the percentage of coordination fluctuates
between a 60% and a 100%. We consider that this broad fluctuation is related
to the limited sample size of our simulations that introduces a certain noise
in our results.

Figure 5: Left column represents the number of Ca2+ coordinated to lipids in
presence of [CaCl2] at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N and (D) 0.50 N, for
T = 330 K. Right column shows the quantity of calcium ions coordinated to
lipids expressed in percentage, along the simulated trajectory.

###  3.2 Effect of temperature

This section focuses on the study of the role played by temperature on the
equilibration process. In this regard, only the system corresponding to a
concentration of 0.25 N in CaCl2 was studied, for a range of temperatures from
330 K to 350 K (all of them above the transition temperature of 314 K [44, 45]
for the DPPC).

Figure 6 shows the running area along the trajectory. In this case, it was
noticed how the systems achieve a steady state after a trajectory length of
roughly 10 ns, where Table 3 shows the mean area per lipid calculated from the
last 70 ns of simulation time. Figure 7 shows the deuterium order parameter of
the methylene groups along the lipid tails, calculated from Eq. (3). Figure 7,
on the one hand, clearly shows that for the three temperatures the systems
have reached the steady state before the first 10 ns of simulation. On the
other hand, it shows an increase in the disorder of the lipid tails with
temperature, which is closely related with the increase of the area per lipid,
such as it was pointed out above. Figure 8 depicts the results of the
hydration numbers of DPPC for the three temperatures studied, where the
equilibrated state was achieved after 10 ns of simulation time. Table 3
provides the calculated hydration numbers in the equilibrium, showing how the
lipid hydration remained invariable with the rising of the temperature.
Concerning the lipid-calcium coordination, Fig. 9 represents the lipid-calcium
coordination number, and the right column represents the calcium that
participates in the coordination expressed in percentage respect the total of
calcium ions in solution. From simulation, it becomes evident how calcium ions
required less than 5 ns to achieve an equilibrated state for the three
temperatures studied. In summary, for all the properties studied in this
section, a slight decrease in the equilibration time with the increasing
temperature was appreciated.

T (K) |  $\langle A\rangle$ (nm2) | Hydration Number
---|---|---
330 | 0.642 $\pm$0.009 | 1.680 $\pm$ 0.010
340 | 0.650 $\pm$0.007 | 1.683 $\pm$ 0.020
350 | 0.666 $\pm$0.008 | 1.689 $\pm$ 0.015
Table 3: Area per lipid and the lipid hydration number as a function of
temperature. Error bars were calculated from subtrajectories of 10 ns length.
DPPC bilayer in presence of [CaCl2] = 0.25 N. Figure 6: Running area per lipid
for [CaCl2] = 0.25 N at different temperatures, (A) T = 330 K, (B) T = 340 K
and (C) T = 350 K. Solid lines represent the mean values obtained from the
last 70 ns of simulation. Figure 7: Deuterium Order Parameter, $-S_{CD}$,
along the simulated trajectory for a concentration of [CaCl2] = 0.25 N, for
the following temperatures: (A) T = 330 K, (B) T = 340 K and (C) T = 350 K.
Solid lines represent the average order parameter for the last 70 ns of
simulation. Note that the error bars have the same size as symbol. Figure 8:
Hydration number of phosphate oxygen (atom 10 in Fig. 1) along the simulated
trajectories for a [CaCl2] = 0.25 N at different temperatures: (A) T = 330 K,
(B) T = 340 K and (C) T = 350 K. Solid lines represent the average hydration
number for the last 70 ns of simulation. Note that the error bars have the
same size as symbol. Figure 9: Left column represents the number of calcium
ions involved in the lipid coordination along time for a concentration of
[CaCl2]= 0.25 N at different temperatures: (A) T = 330 K, (B) T = 340 K and
(C) T = 350 K. The right column shows the same information expressed as a
percentage of the total number of calcium ions in solution.

##  4 Conclusions

The present work deals with the simulation time required to achieve the steady
state for a lipid bilayer system in presence of CaCl2. In this regard, we
studied two different systems: one with DPPC and another one with DPPS
bilayer; both systems in presence of CaCl2 (at different level concentration).
The salt free case was also studied, as control. The analysis of various lipid
properties studied here indicates that some properties reach the steady state
more quickly than others. In this sense, we found that the area per lipid and
the hydration number are slower than the deuterium order parameter and the
coordination of cations. Consequently, to ensure that a system composed by a
lipid bilayer has reached a steady state, the criterion that we propose is to
show that the area per lipid or the hydration number have reached the
equilibrium.

From our results, two important aspects should be remarked:

  1. 1.

The equilibration time is strongly dependent on the starting conformation of
the system. Wrong starting conformations will require much longer
equilibration times, even of one order of magnitude higher than the requested
from a more refined starting conformation.

  2. 2.

Temperature is a critical parameter for reducing the equilibration time in our
simulations, due to the fact that higher temperatures increase the kinetic
processes, i.e., the sampling of the configurational space of the system.

###### Acknowledgements.

Authors wish to thank the assistance of the Computing Center of the
Universidad Politécnica de Cartagena (SAIT), Spain. RDP is member of ‘Carrera
del Investigador’, CONICET, Argentine.

## References

  * [1] O Berger, O Edholm, F Jahnig, Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature, Biophys. J. 72, 2002 (1997).
  * [2] E Egberts, S J Marrink, H J C Berendsen, Molecular dynamics simulation of a phospholipid membrane, Eur. Biophys. J. 22, 423 (1994).
  * [3] U Essmann, L Perera, M L Berkowitz, The origin of the hydration interaction of lipid bilayers from MD simulation of dipalmitoylphosphatidylcholine membranes in gel and crystalline phases, Langmuir 11, 4519 (1995).
  * [4] S E Feller, Y Zhang, R W Pastor, R B Brooks, Constant pressure molecular dynamics simulations: The Langevin piston method, J. Chem. Phys. 103, 4613 (1995).
  * [5] W Shinoda, T Fukada, S Okazaki, I Okada, Molecular dynamics simulation of the dipalmitoylphosphatidylcholine (DPPC) lipid bilayer in the fluid phase using the Nosr-Parrinello-Rahman NPT ensemble, Chem. Phys. Lett. 232, 308 (1995).
  * [6] D P Tieleman, H J C Berendsen, Molecular dynamics simulations of a fully hydrated dipalmitoylphosphatidylcholine bilayer with different macroscopic boundary conditions and parameters, J. Chem. Phys. 105, 4871 (1996).
  * [7] K Tu, D J Tobias, M L Klein, Constant pressure and temperature molecular dynamics simulation of a fully hydrated liquid crystal phase dipalmitoylphosphatidylcholine bilayer, Biophys. J. 69, 2558 (1995).
  * [8] M F Brown, Theory of spin-lattice relaxation in lipid bilayers and biological membranes. Dipolar relaxation, J. Chem. Phys. 80, 2808 (1984).
  * [9] M F Brown, Theory of spin-lattice relaxation in lipid bilayers and biological membranes. ${}^{2}H$ and ${}^{14}N$ quadrupolar relaxation, J. Phys. Chem. 77, 1576 (1982).
  * [10] J F Nagle, R Zang, S Tristam-Nagle, W S Sun, H I Petrache, R M Suter, X-ray structure determination of fully hydrated L. phase dipalmitoylphosphatidylcholine bilayers, Biophys. J. 70, 1419 (1996).
  * [11] R P Rand, V A Parsegian, Hydration forces between phospholipid bilayers, Biochim. Biophys. Acta 988, 351 (1989).
  * [12] J Seelig, Deuterium magnetic resonance: Theory and application to lipid membranes, Q. Rev. Biophys. 10, 353 (1977).
  * [13] J Seelig, A Seelig, Lipid conformation in model membranes and biological systems, Q. Rev. Biophys. 13, 19 (1980).
  * [14] W J Sun, R M Suter, M A Knewtson, C R Worthington, S Tristram-Nagle, R Zhang, J F Nagle, Order and disorder in fully hydrated unoriented bilayers of gel phase dipalmitoylphosphatidylcholine, Phys. Rev. E. 49, 4665 (1994).
  * [15] H Akutsu, J Seelig, Interaction of metal ions with phosphatidylcholine bilayer membranes, Biochemistry 20, 7366 (1981).
  * [16] M G Ganesan, D L Schwinke, N Weiner, Effect of Ca2+ on thermotropic properties of saturated phosphatidylcholine liposomes, Biochim. Biophys. Acta 686, 245 (1982).
  * [17] L Herbette, C A Napolitano, R V McDaniel, Direct determination of the calcium profile structure for dipalmitoyllecithin multilayers using neutron diffraction, Biophys. J. 46, 677 (1984).
  * [18] D Huster, K Arnold, K Gawrisch, Strength of Ca2+ binding to retinal lipid membrane: Consequences for lipid organization, Biophys. J. 78, 3011 (2000).
  * [19] Y Inoko, T Yamaguchi, K Furuya, T Mitsui, Effects of cations on dipalmitoyl phosphatidylcholine/cholesterol/water systems, Biochim. Biophys. Acta 413, 24 (1975).
  * [20] R Lehrmann, J J Seelig, Adsorption of Ca2+ and La3+ to bilayer membranes: Measurement of the adsorption enthalpy and binding constant with titration calorimetry, Biochim. Biophys. Acta 1189, 89 (1994).
  * [21] L J Lis, W T Lis, V A Parsegian, R P Rand, Adsorption of divalent cations to a variety of phosphatidylcholine bilayers, Biochemistry 20, 1771 (1981).
  * [22] L J Lis, V A Parsegian, R P Rand, Binding of divalent cations to dipalmitoylphosphatidylcholine bilayers and its effect on bilayer interaction, Biochemistry 20, 1761 (1981).
  * [23] T Shibata, Pulse NMR study of the interaction of calcium ion with dipalmitoylphosphatidylcholine lamellae, Chem. Phys. Lipids. 53, 47 (1990).
  * [24] S A Tatulian, V I Gordeliy, A E Sokolova, A G Syrykh, A neutron diffraction study of the influence of ions on phospholipid membrane interactions, Biochim. Biophys. Acta 1070, 143 (1991).
  * [25] R A Böckmann, H Grubmüller, Multistep binding of divalent cations to phospholipid bilayers: A molecular dynamics study, Angewandte Chemie 43, 1021 (2004).
  * [26] J Faraudo, A Travesset, Phosphatidic acid domains in membranes: Effect of divalent counterions, Biophys. J. 92, 2806 (2007).
  * [27] A A Gurtovenko, Asymmetry of lipid bilayers induced by monovalent salt: Atomistic molecular-dynamics study, J. Chem. Phys. 122, 244902 (2005).
  * [28] P Mukhopadhyay, L Monticelli, D P Tieleman, Molecular dynamics simulation of a palmitoyl-oleoyl phosphatidylserine bilayer with Na+ counterions and NaCl, Biophys. J. 86, 1601 (2004).
  * [29] S A Pandit, D Bostick, M L Berkowitz, Molecular dynamics simulation of a dipalmitoylphosphatidylcholine bilayer with NaCl, Biophys. J. 84, 3743 (2003).
  * [30] U R Pedersen, C Laidy, P Westh, G H Peters, The effect of calcium on the properties of charged phospholipid bilayers, Biochim. Biophys. Acta 1758, 573 (2006).
  * [31] J N Sachs, H Nanda, H I Petrache, T B Woolf, Changes in phosphatidylcholine headgroup tilt and water order induced by monovalent salts: Molecular dynamics simulations, Biophys. J. 86, 3772 (2004).
  * [32] K Shinoda, W Shinoda, M Mikami, Molecular dynamics simulation of an archeal lipid bilayer whit sodium chloride, Phys. Chem. Chem. Phys. 9, 643 (2007).
  * [33] N L Yamada, H Seto, T Takeda, M Nagao, Y Kawabata, K Inoue, SAXS, SANS and NSE studies on “unbound state” in DPPC/water/CaCl2 system, J. Phys. Soc. Jpn. 74, 2853 (2005).
  * [34] D Frenkel, B Smit, Understanding molecular simulations, Academic Press, New York (2002).
  * [35] J J López Cascales, J García de la Torre, S J Marrink, H J C Berendsen, Molecular dynamics simulation of a charged biological membrane, J. Chem. Phys. 104, 2713 (1996).
  * [36] W F van Gunsteren, H J C Berendsen, Computer simulations of molecular dynamics: Methodology, applications and perspectives in chemistry, Angew. Chem Int. Ed. Engl. 29, 992 (1990).
  * [37] R A Böckmann, A Hac, T Heimburg, H Grubmüller, Effect of sodium chloride on a lipid bilayer, Biophys. J. 85, 1647 (2003).
  * [38] A A Gurtovenko, I Vattulainen, Effect of NaCl and KCl on phosphatidylcholine and phosphatidylethanolamine lipid membranes: Insight from atomic-scale simulations for understanding salt-induced effects in the plasma membrane, J. Phys. Chem. B. 112, 1953 (2008).
  * [39] H J C Berendsen, J R Grigera, T P Straatsma, The missing term in effective pair potentials, J. Phys. Chem. 91, 6269 (1987).
  * [40] H J C Berendsen, D van der Spoel, R van Drunen, A message-passing parallel molecular dynamics implementation, Comp. Phys. Comm. 91, 43 (1995).
  * [41] E Lindahl, B Hess, D van der Spoel, GROMACS 3.0: A package for molecular simulation and trajectory analysis, J. Mol. Mod. 7, 306 (2001).
  * [42] T Darden, D York, L Pedersen, Particle mesh Ewald: An N.log(N) method for Ewald sums in large systems, J. Chem. Phys. 98, 10089 (1993).
  * [43] U Essmann, L Perea, M L Berkowitz, T Darden, H Lee, L G Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys. 103, 8577 (1995).

  * [44] L R De Young, K A Dill, Solute partitioning into lipid bilayer-membranes, Biochemistry 27, 5281 (1988).
  * [45] A Seelig, J Seelig, The dynamic structure of fatty acyl chains in a phospholipid bilayer measured by deuterium magnetic resonance, Biochemistry 13, 4839 (1974).
  * [46] G Cevc, A Watts, D Marsh, Titration of the phase transition of phosphatidilserine bilayer membranes. Effect of pH, surface electrostatics, ion binding and head-group hydration, Biochemistry 20, 4955 (1981).
  * [47] H Hauser, F Paltauf, G G Shipley, Structure and thermotropic behavior of phosphatidylserine bilayer membranes, Biochemistry 21, 1061 (1982).
  * [48] J F Nagle, S Tristam-Nagle, Structure of lipid bilayers, Biochim. Biophy. Acta 1469, 159 (2000).
  * [49] B A Lewis, D M Engelman, Lipid bilayer thickness varies linearly with acyl chain length in fluid phosphatidylcholine vesicles, J. Mol. Biol. 166, 211 (1983).
  * [50] R J Pace, S I Cham, Molecular motions in lipid bilayer. I. Statistical mechanical model of acyl chains motion, J. Chem. Phys. 76, 4217 (1982).
  * [51] R L Thurmond, S W Dodd, M F Brown, Molecular areas of phospholipids as determined by 2H NMR spectroscopy, Biphys. J. 59, 108 (1991).
  * [52] R D Porasso, J J López Cascales, Study of the effect of Na+ and Ca2+ ion concentration on the structure of an asymmetric DPPC/DPPS + DPPS lipid bilayer by molecular dynamics simulation, Coll. and Surf. B. Bioint. 73, 42 (2009).


June 29, 2022

# Nematic Tomonaga-Luttinger Liquid Phase in an $S\\!=\\!1/2$ Ferromagnetic-
Antiferromagnetic Bond-Alternating Chain

Takashi Tonegawa1,2,3 Kiyomi Okamoto2 Kiyohide Nomura4 and Tôru Sakai2,5 1
Professor Emeritus1 Professor Emeritus Kobe University Kobe University Kobe
657-8501 Kobe 657-8501 Japan
2 Graduate School of Science Japan
2 Graduate School of Science University of Hyogo University of Hyogo Hyogo
678-1297 Hyogo 678-1297 Japan
3 Department of Physics Japan
3 Department of Physics Osaka Metropolitan University Osaka Metropolitan
University Sakai 599-8531 Sakai 599-8531 Japan
4 Department of Physics Japan
4 Department of Physics Kyushu University Kyushu University Fukuoka 812-8581
Fukuoka 812-8581 Japan
5 National Institute for Quantum Science and Technology (QST) Japan
5 National Institute for Quantum Science and Technology (QST) SPring-8
SPring-8 Hyogo 679-5148 Hyogo 679-5148 Japan Japan<EMAIL_ADDRESS>
###### Abstract

We numerically investigate the ground-state phase diagram of the $S\\!=\\!1/2$
ferromagnetic-antiferromagnetic bond-alternating chain, in which the
ferromagnetic interactions are stronger than the antiferromagnetic ones, and
the anisotropies of the former and latter interactions are of the Ising-type
and the $XY$-type, respectively. We use various numerical methods, such as the
level spectroscopy and phenomenological renormalization-group analyses of the
numerical data obtained by the exact diagonalization method, and so on. The
resultant phase diagrams contain the ferromagnetic, $XY$1, singlet-dimer, and
up-up-down-down phases as well as the nematic Tomonaga-Luttinger liquid (nTLL)
phase which appears in a wide region of the interaction parameters.
Perturbation calculations from the strong limit of the ferromagnetic
interactions reproduce fairly well the numerical results of the phase boundary
lines associated with the nTLL phase in the phase diagrams.

$S\\!=\\!1/2$ ferromagnetic-antiferromagnetic bond-alternating chain, ground-
state phase diagram, nematic Tomonaga-Luttinger liquid phase, numerical
calculation, perturbation theory

##  1 Introduction

Over the past several decades, a great deal of numerical, theoretical, and
experimental work has been devoted to clarifying the quantum phase transition
in one-dimensional $S\\!=\\!1/2$ systems. As is well known, due to the strong
quantum fluctuation, a variety of exotic quantum phases appear in the ground-
state phase diagrams of these systems. A typical example attracted recently
much attention is the nematic Tomonaga-Luttinger liquid (nTLL) phase [1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12], which is characterized not only by the
formation of two-magnon bound pairs but also by the dominant nematic four-spin
correlation function.

Recently, we[12] have numerically explored the ground-state phase diagram of
an $S\\!=\\!1/2$ anisotropic two-leg ladder with different leg interactions in
the absence of external magnetic field. This system has frustration when the
signs of the two kinds of leg interactions are different from each other.
Then, we have found that, when the ferromagnetic rung interactions with the
Ising-type anisotropy is much stronger than the antiferromagnetic leg
interactions with the $XY$-type anisotropy, the nTLL phase appears in the
unfrustrated region as well as in the frustrated region. For this phase in the
former region, the asymptotic form of the nematic correlation function shows
the power-law decay with the uniform character; on the other hand, for this
phase in the latter region, that shows the power-law decay with the staggered
character. Thus, both nTLL phases are different phases. As far as we know,
this is the first report of the realization of the nTLL phase in an
$S\\!=\\!1/2$ unfrustrated one-dimensional spin system under no external
magnetic field.

According to this result, it is reasonably expected that the nTLL state
appears as the zero-field ground state in general $S\\!=\\!1/2$ unfrustrated
one-dimensional systems in which pairs of $S\\!=\\!1/2$ spins coupled strongly
with the Ising-type ferromagnetic interaction are connected by the weak
antiferromagnetic interactions. As an example of such systems, we investigate
in this paper an $S\\!=\\!1/2$ ferromagnetic-antiferromagnetic bond-
alternating chain. We express the Hamiltonian describing this system as

${\cal H}=-J_{{\rm F}}\sum_{j=1}^{N/2}\bigl{\\{}\Gamma_{{\rm
F}}\bigl{(}S_{2j-1}^{x}S_{2j}^{x}+S_{2j-1}^{y}S_{2j}^{y}\bigr{)}+S_{2j-1}^{z}S_{2j}^{z}\bigr{\\}}+J_{\rm
AF}\sum_{j=1}^{N/2}\bigl{\\{}S_{2j}^{x}S_{2j+1}^{x}+S_{2j}^{y}S_{2j+1}^{y}+\Delta_{{\rm
AF}}\,S_{2j}^{z}S_{2j+1}^{z}\bigr{\\}}\,,$ (1)

where $S_{j}^{\mu}$ ($\mu\\!=\\!x,y,z$) is the $\mu$-component of the
$S\\!=\\!1/2$ operator $S$j at the $j$th site; $J_{{\rm F}}$ and $J_{{\rm
AF}}$ denote, respectively, the magnitudes of the ferromagnetic and
antiferromagnetic interactions; $\Gamma_{{\rm F}}$ and $\Delta_{{\rm AF}}$
are, respectively, the parameters representing the $XXZ$-type anisotropies of
the former and latter interactions; $N$ is the total number of spins in the
system, which is assumed to be a multiple of four. We assume that $J_{{\rm
F}}\\!>\\!J_{{\rm AF}}\\!>\\!0.0$, $1.0\geq\Gamma_{{\rm F}}\geq 0.0$, and
$1.0\geq|\Delta_{{\rm AF}}|$, that is, the ferromagnetic interactions are
stronger than the antiferromagnetic ones, and the anisotropies of the former
and latter interactions are of the Ising-type and the $XY$-type, respectively.

##  2 Ground-State Phase Diagram

First of all, we show in Fig. 1 our numerical results of the ground-state
phase diagrams obtained in the following cases: (a) the phase diagram on the
$\Delta_{\rm AF}$ versus $\Gamma_{\rm F}$ plane in the case of $J_{{\rm
F}}\\!=\\!1.0$ and $J_{{\rm AF}}\\!=\\!0.1$, (b) that on the $\Delta_{\rm AF}$
versus $J_{\rm AF}$ plane in the case of $J_{{\rm F}}\\!=\\!1.0$ and
$\Gamma_{{\rm F}}\\!=\\!0.8$, and (c) that on the $\Gamma_{\rm F}$ versus
$J_{\rm AF}$ plane in the case of $J_{{\rm F}}\\!=\\!1.0$ and $\Delta_{{\rm
AF}}\\!=\\!-0.12$. These diagrams contain the ferromagnetic (F), $XY$1,
singlet-dimer (SD), and up-up-down-down (UUDD) phases in addition to the nTLL
phase which appears in wide regions. The physical pictures of the SD and UUDD
states are sketched in Fig. 2. The SD state is a unique and symmetry protected
topological gapped state and the UUDD state is a doubly degenerate gapped
state.

Figure 1: Ground-state phase diagrams (a) on the $\Delta_{\rm AF}$ versus
$\Gamma_{\rm F}$ plane in the case of $J_{{\rm F}}\\!=\\!1.0$ and $J_{{\rm
AF}}\\!=\\!0.1$, (b) on the $\Delta_{\rm AF}$ versus $J_{\rm AF}$ plane in the
case of $J_{{\rm F}}\\!=\\!1.0$ and $\Gamma_{{\rm F}}\\!=\\!0.8$, and (c) on
the $\Gamma_{\rm F}$ versus $J_{\rm AF}$ plane in the case of $J_{{\rm
F}}\\!=\\!1.0$ and $\Delta_{{\rm AF}}\\!=\\!-0.12$, determined numerically in
the present work. In (a’), (b’), and (c’), parts of (a), (b), and (c) are
enlarged, respectively. The results of the perturbation calculations are shown
in (a”), (b”), and (c”); by these calculations the $XY$1-SD, SD-UUDD, and
F-$XY1$ phase transition lines cannot be obtained.
Figure 2: Physical pictures of the SD (left) and UUDD (right) states. Open
circles denote $S\\!=\\!1/2$ spins, and solid and dotted lines denote
ferromagnetic and antiferromagnetic bonds, respectively. Ellipses represent
singlet pairs of two $S\\!=\\!1/2$ spins, while arrows denote fixed projection
values of the $S\\!=\\!1/2$ spins. As can be seen from these figures, the SD
phase is a unique gapped phase, while the UUDD phase is a doubly degenerate
gapped phase.

Let us now discuss how to determine numerically the phase boundary lines in
the phase diagram shown in Fig. 1. We denote, respectively, by $E_{0}^{\rm
P}(N,M)$ and $E_{1}^{\rm P}(N,M)$, the lowest and second-lowest energy
eigenvalues of the Hamiltonian ${\cal H}$ under periodic boundary conditions
within the subspace of fixed $N$ and $M$, where $M(=\\!0,\pm 1,\cdots,\pm
N/2)$ is the total magnetization. Furthermore, we denote by $E_{0}^{\rm
T}(N,M,P)$ the lowest eigenvalue of ${\cal H}$ under twisted boundary
conditions within the subspace of fixed $N$, $M$, and $P$, where $P(=\\!\pm
1)$ is the eigenvalue of the space inversion operator with respect to the
twisted bond. We have numerically calculated these energies for finite-size
systems with up to $N\\!=\\!28$ spins by means of the exact-diagonalization
method. The ground-state energy of the finite-$N$ system is given by
$E_{0}^{\rm P}(N,N/2)$ in the F phase and by $E_{0}^{\rm P}(N,0)$ in the other
phases. In the following way, we have estimated the finite-size critical
values of the interaction parameters for each phase transition. Then, the
phase boundary line for the transition has been obtained by connecting the
results for the $N\\!\to\\!\infty$ extrapolation of the finite-size critical
values.

First, the phase transition between the $XY$1 and SD phases is the
Berezinskii-Kosterlitz-Thouless (BKT) transition[13, 14]. It is known that for
this transition, the level spectroscopy method developed by Nomura and
Kitazawa[16] is very powerful for calculating the phase transition line. This
method implies that the finite-size critical values are estimated from

$\Delta E_{0}^{\rm P}(N,2)=\Delta E_{0}^{\rm T}(N,0,-1),$ (2)

under the condition that $\Delta E_{0}^{\rm T}(N,0,+1)$ is larger than these
excitations, where $\Delta E_{0}^{\rm P}(N,M)=E_{0}^{\rm P}(N,M)-E_{0}^{\rm
P}(N,0)$ and $\Delta E_{0}^{\rm T}(N,M,P)=E_{0}^{\rm T}(N,M,P)-E_{0}^{\rm
P}(N,0)$. This equation is also applicable to estimate the nTLL-UUDD phase
transition line (see Appendix). Secondly, the phase transition between the SD
and UUDD phases is the 2D Ising-type transition. Therefore, as is well known,
the phase transition line is determined by the phenomenological
renormalization-group (PRG) method[17]. Then, to estimate the finite-size
critical values, we solve the PRG equation,

$N\,\Delta E_{1}^{\rm P}(N,0)=(N+4)\,\Delta E_{1}^{\rm P}(N+4,0)\,,$ (3)

where $\Delta E_{1}^{\rm P}(N,M)\equiv E_{1}^{\rm P}(N,M)-E_{0}^{\rm P}(N,M)$.
Thirdly, as for the phase transition between the nTLL and $XY$1 phases, the
nTLL state accompanies two-magnon bound-states, while the $XY$1 state does
not. Therefore, in the ground-state magnetization curve for the finite-size
system, the magnetization increases from $M\\!=\\!0$ to $M\\!=\\!2$ in the
former state and from $M\\!=\\!0$ to $M\\!=\\!1$ in the latter state. Thus,
the finite-size critical values are estimated from

$\Delta E_{0}^{\rm P}(N,2)=2\,\Delta E_{0}^{\rm P}(N,1)\,.$ (4)

We note that the binding energy of two magnons[12] is defined by $E_{\rm
bind}(N)\\!\equiv\\!\Delta E_{0}^{\rm P}(N,2)\\!-\\!2\,\Delta E_{0}^{\rm
P}(N,1)$. Accordingly, Eq.(4) is also the condition of $E_{\rm
bind}(N)\\!=\\!0$. Lastly, it is apparent that the finite-size critical values
for the phase transitions between the F phase and one of the $XY$1 and nTLL
phases are estimated from

$E_{0}^{\rm P}(N,N/2)\\!=\\!E_{0}^{\rm P}(N,0)\,.$ (5)

##  3 Perturbation Theory

In the nTLL phase, the important states of two $S\\!=\\!1/2$ spins connected
by the Ising-like ferromagnetic coupling $J_{\rm F}$ are
$|\uparrow\uparrow\rangle$ and $|\downarrow\downarrow\rangle$, and the
remaining two states
$|(1/\sqrt{2})(\uparrow\downarrow\pm\downarrow\uparrow)\rangle$ have higher
energies. To describe the nTLL state, we introduce a pseudo-spin operator $T$
with $T\\!=\\!1/2$, where $|T^{z}=1/2\rangle=|\uparrow\uparrow\rangle$ and
$|T^{z}=-1/2\rangle=|\downarrow\downarrow\rangle$. We perturbationally derive
the effective Hamiltonian $\cal H_{\rm eff}$ described by $T$. We take the
first term of the right-hand side of Eq.(1) as the unperturbed Hamiltonian,
and the second term as the perturbation. Up to the second order perturbation
calculation, we obtain the following $XXZ$ model, apart from the energy shift,

${\cal H}_{\rm eff}=\sum_{j=1}^{N/2}\left\\{J_{\rm
eff}^{\perp}(T_{j}^{x}T_{j+1}^{x}+T_{j}^{y}T_{j+1}^{y})+J_{\rm
eff}^{z}T_{j}^{z}T_{j+1}^{z}\right\\}$ (6)

with

$\displaystyle J_{\rm eff}^{\perp}=c_{1}-2c_{2}+c_{3},~{}~{}~{}~{}~{}J_{\rm
eff}^{z}=b+c_{1}+2c_{2}+c_{3},$ (7) $\displaystyle b=J_{\rm AF}\Delta_{\rm
AF},~{}~{}~{}~{}~{}c_{1}={J_{\rm AF}^{2}\over 8J_{\rm F}(1-\Gamma_{\rm
F})},~{}~{}~{}~{}~{}c_{2}={J_{\rm AF}^{2}\over 8J_{\rm
F}},~{}~{}~{}~{}~{}c_{3}={J_{\rm AF}^{2}\over 8J_{\rm F}(1+\Gamma_{\rm F})},$
(8)

where $J_{\rm eff}^{\perp}\geq 0$. The very well known exact solution[14] for
Eq. (6) shows that the ground state phase diagram of ${\cal H}_{\rm eff}$
consists of the F phase, the TLL phase, and the Néel phase. These phases
correspond, respectively, to the F phase, the nTLL phase and the UUDD phase of
the original model (1). Thus, the F-nTLL boundary line of the original model
(1) is determined by $J_{\rm eff}^{z}=-J_{\rm eff}^{\perp}$, while the nTLL-
UUDD boundary line by $J_{\rm eff}^{z}=J_{\rm eff}^{\perp}$. We note that the
same effective Hamiltonian with the periodic boundary condition is obtained
irrespective of the periodic boundary condition or the twisted boundary
condition of the original Hamiltonian (1).

The nTLL-$XY1$ boundary can be estimated by considering the energy cost of
replacing a $|T^{z}=\pm 1/2\rangle$ pseudo-spin in the ${\cal H}_{\rm eff}$
picture by $|(1/\sqrt{2})(\uparrow\downarrow+\downarrow\uparrow)\rangle$. If
this energy cost is negative, a macroscopic number of spin pair states with
$|(1/\sqrt{2})(\uparrow\downarrow+\downarrow\uparrow)\rangle$ are generated,
which brings about the breakdown of the ${\cal H}_{\rm eff}$ picture. On the
other hand, if this energy cost is positive, the
$|(1/\sqrt{2})(\uparrow\downarrow+\downarrow\uparrow)\rangle$ spin pair state
are scarcely yielded, resulting in the stability of the ${\cal H}_{\rm eff}$
picture. The replacement of a pseudo-spin
$|T^{z}=-1/2\rangle\Rightarrow|(1/\sqrt{2})(\uparrow\downarrow+\downarrow\uparrow)\rangle$
in the $M=0$ ground state leads to the $M=1$ state. Thus, if this cost is
positive, the $M=1$ excitation is gapped, which is consistent with the picture
of the nTLL state. We note that the replacement of a $|T^{z}=-1/2\rangle$ spin
state by a $|T^{z}=+1/2\rangle$ one brings about the $M=2$ state. After some
calculations, we obtain this energy cost as $(1/2)\\{J_{\rm F}(1-\Gamma_{\rm
F})-J_{\rm AF}\\}$. Thus, the nTLL-$XY1$ boundary line is estimated by $J_{\rm
AF}=J_{\rm F}(1-\Gamma_{\rm F})$.

The above three boundary lines are depicted in Figs. 1(a”), (b”), and (c”).

##  4 Concluding Remarks

We have determined the ground-state phase diagram of the $S\\!=\\!1/2$
ferromagnetic-antiferromagnetic bond-alternating chain described by the
Hamiltonian of Eq. (1) by using mainly the numerical methods. Discussing the
case where the ferromagnetic interactions with the Ising-type anisotropy are
much stronger than the antiferromagnetic ones with the $XY$-type anisotropy,
we have shown that in the phase diagrams (see Fig. 1), there appear the F,
$XY$1, SD, and UUDD phases as well as the nTLL phase. We have also developed
the perturbation theory, which qualitatively explains the numerically obtained
phase boundary lines associated with the nTLL phase. This paper gives, in a
succession of our previous work[12], a report of the appearance of the nTLL
phase in an unfrustrated $S\\!=\\!1/2$ chain under no external magnetic field.

It should be noted that the SDW2[3] (’SDW’ is an abbreviation for ’spin-
density-wave’) does not appear in the ground-state phase diagram of the
present model. The reason for this is as follows. In both of the nTLL state
and the SDW2 state, both of the nematic four-spin correlation
$C_{2}(j)\equiv\langle S_{1}^{+}S_{2}^{+}S_{1+j}^{-}S_{2+j}^{-}\rangle$ and
the longitudinal two-spin correlation $C_{z}(j)\equiv\langle
S_{1}^{z}S_{1+j}^{z}\rangle$, where $j$ is even, decay in power laws. The
former correlation is dominant in the nTLL state, while the latter is dominant
in the SDW2 state. Namely, $\eta_{2}<\eta_{z}$ in the nTLL state, whereas
$\eta_{2}>\eta_{z}$ in the SDW2 state, where $\eta_{2}$ and $\eta_{z}$ are the
power-decay exponents of $C_{2}(j)$ and $C_{z}(j)$, respectively. Since there
is the TLL relation $\eta_{2}\eta_{z}=1$, it holds $\eta_{2}>1$ in the SDW2
state. However, this condition $\eta_{2}>1$ is nothing but the necessary
condition for the UUDD state (see Appendix). In the zero magnetization case,
the operator coming from the Umklapp process becomes relevant due to
$\eta_{2}>1$, which leads to the UUDD state. Thus, the SDW2 region does not
exist in the ground state phase diagram. On the other hand, in the finite
magnetization case, since the operator induced by the Umklapp process does not
exist, the SDW2 region appears in the phase diagram.

One can see several tricritical points in the phase diagrams given in Fig. 1;
these are the F-$XY$1-nTLL, $XY$1-SD-UUDD, and $XY$1-UUDD-nTLL tricritical
points in the case (a), the F-$XY$1-nTLL, $XY$1-SD-nTLL, and SD-UUDD-nTLL
tricritical points in the case (b), and the SD-UUDD-nTLL and $XY$1-SD-nTLL
tricritical points in the case (c). According to the discussion on the
university class (see Appendix), on the other hand, two tricritical points
associated with the $XY$1, SD, UUDD, and nTLL phases should merge into one
$XY$1-SD-UUDD-nTLL tetracritical point. In order to obtain numerically this
tetracritical point, it is indispensable to develop a new numerical method
which is applicable commonly to the phase transition between the $XY1$ and
nTLL phases and to that between the SD and UUDD phase. Since both of these
transitions are of the 2D Ising-type, it seems not to be impossible to find
this method. However, it is has not yet been succeeded at present, and thus
this problem is beyond the scope of the present study. Furthermore, it may be
expected that the nTLL phase appear even in the case where the
antiferromagnetic interactions are stronger than the ferromagnetic
interactions in the ground-state phase diagram of the present system, if
$1.0\gg\Gamma_{{\rm F}}\\!\geq\\!0.0.$ and $1.0\gg|\Delta_{{\rm AF}}|$. We are
now planning to perform this calculation, and the results will be reported in
the near future.

## Acknowledgments

This work has been partly supported by JSPS KAKENHI, Grant Numbers 16K05419,
16H01080 (J-Physics), 18H04330 (J-Physics), JP20K03866, JP20H05274, and
21H05021. We also thank the Supercomputer Center, Institute for Solid State
Physics, University of Tokyo and the Computer Room, Yukawa Institute for
Theoretical Physics, Kyoto University for computational facilities.

##  Appendix A The Level Spectroscopy Method for the nTLL-UUDD Transition

Let us consider the transverse two-spin correlation $C_{1}(j)\equiv\langle
S_{1}^{+}S_{1+j}^{-}\rangle$ and the nematic four-spin correlation
$C_{2}(j)\equiv\langle S_{1}^{+}S_{2}^{+}S_{1+j}^{-}S_{2+j}^{-}\rangle$, where
$j$ is even. Both of them decay in power laws in the $XY1$ phase with the
relation $\eta_{2}=4\eta_{1}$ [12], where $\eta_{1}$ and $\eta_{2}$ are the
power-decay exponents of $C_{1}(j)$ and $C_{2}(j)$, respectively. As is well
known, the BKT transition from the $XY1$ phase to the unique gapped phase (the
SD phase in the present case) occurs at $\eta_{1}=1/4$ (hence $\eta_{2}=1$),
whereas to the doubly degenerate gapped phase (the UUDD phase in the present
case) at $\eta_{1}=1$ (accordingly $\eta_{2}=4$) [15, 16]. In the level
spectroscopy method by Nomura and Kitazawa[16] for the BKT transition between
the $XY1$ phase and the SD phase, we search the $\eta_{2}=1$ line. This is
clear from the fact that the excitation $\Delta E_{0}^{\rm P}(N,2)$ is used in
Eq.(2), since $\Delta E_{0}^{\rm P}(N,2)$ is closely related to $\eta_{2}$ as
$\eta_{2}=\lim_{N\to\infty}[N\Delta E_{0}^{\rm P}(N,2)/\pi v_{\rm s}]$, where
$v_{\rm s}$ is the spin wave velocity [15, 16, 18, 19].

On the other hand, in the nTLL phase, $C_{1}(j)$ exhibits the exponential
decay, while $C_{2}(j)$ the power decay. Thus, the role of $\eta_{1}$ in the
BKT transitions from the $XY1$ phase is superseded by $\eta_{2}$ in those from
the nTLL phase. Namely, the BKT transition from the nTLL phase to the unique
gapped phase occurs at $\eta_{2}=1/4$, while that to the doubly degenerate
gapped state at $\eta_{2}=1$. Hence, the BKT transition from the nTLL phase to
the doubly degenerate gapped phase (the UUDD phase in the present case) can
also be determined by Eq.(2) in which the line $\eta_{2}=1$ is searched.

Neither the BKT transition from the $XY1$ phase to the UUDD phase nor that
from the nTLL phase to the SD phase occurs on the $\eta_{2}=1$ line. This is
because the former transition should occur at $\eta_{1}=1$ ($\eta_{2}=4$) even
if it occurs. Similarly the latter should occur at $\eta_{2}=1/4$ even if it
occurs. Accordingly there should exist the $XY$1-SD-UUDD-nTLL tetracritical
point on the $\eta_{2}=1$ line which is determined by Eq.(2). A similar
discussion can also be applied to the $S\\!=\\!1$ chain with the XXZ and on-
site anisotropies which has been discussed, for example, by Chen et al. [20].
Namely, the $XY$1-Haldane-Néel-nTLL tetracritical point should exist in this
$S\\!=\\!1$ chain.

## References

  * [1] A. V. Chubukov, Phys. Rev. B 44, 4693 (1991).
  * [2] T. Vekua, A. Honecker, H.-J. Mikeska and F. Heidrich-Meisner, Phys. Rev. B 76, 174420 (2007).
  * [3] T. Hikihara, L. Kecke, T. Momoi and A. Furusaki, Phys. Rev. B 78, 144404 (2008).
  * [4] J. Sudan, A. Lüscher and A. M. Läuchli, Phys. Rev. B 80, 140402(R) (2009).
  * [5] T. Sakai, T. Tonegawa and K. Okamoto, Phys. Status Solidi B 247, 583 (2010).
  * [6] M. Sato, T. Hikihara and T. Momoi, Phys. Rev. Lett. 110, 077206 (2013).
  * [7] O. A. Starykh and L. Balents, Phys. Rev. B 89, 104407 (2014).
  * [8] N. Büttgen, K. Nawa, T. Fujita, M. Hagiwara, P. Kuhns, A. Prokofiev, A. P. Reyes, L. E. Svistov, K. Yoshimura and M. Takigawa, Phys. Rev. B 90, 134401 (2014).
  * [9] A. Orlova, E. L. Green, J. M. Law, D. I. Gorbunov, G. Chanda, S. Krämer, M. Horvatić, R. K. Kremer, J. Wosnitza and G. L. J. A. Rikken, Phys. Rev. Lett. 118, 247201 (2017).
  * [10] A. Parvej and M. Kumar, Phys. Rev. B 96, 054413 (2017).
  * [11] M. Bosioc̆ić, F. Bert, S. E. Dutton, R. J. Cava, P. J. Baker, M. Poz̆ek and P. Mendels, Phys. Rev. B 96, 224424 (2017).
  * [12] T. Tonegawa, T. Hikihara, K. Okamoto, S. C. Furuya and T. Sakai, J. Phys. Soc. Jpn. 87, 104002 (2018).
  * [13] Z. L. Berezinskii, Sov. Phys. JETP 34, 610 (1971); J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973).
  * [14] T. Giamarchi, Quantum Physics in One Dimension (Clarendon Press, Oxford, 2003).
  * [15] A. Kitazawa, K. Nomura and K. Okamoto, Phys. Rev. Lett. 76, 4038 (1996).
  * [16] K. Nomura and A. Kitazawa, J. Phys. A 31, 7341 (1998).
  * [17] M. P. Nightingale, Physica A 83, 561 (1976).
  * [18] J. L. Cardy, J. Phys. A: Math. Gen. 17, L385 (1984).
  * [19] T. Sakai, Phys. Rev. B 58, 6268 (1998).
  * [20] W. Chen, K. Hida and B. C. Sanctuary, Phys. Rev.B 67, 104401 (2003).