# Contrasting global approaches for identifying and managing cybersecurity
risks in supply chains

Colin Topping
Bristol Cyber Security group
University of Bristol, UK

Ola Michalec
Bristol Cyber Security group
University of Bristol, UK

Awais Rashid
Bristol Cyber Security group
University of Bristol, UK

###### Abstract

Supply chains are increasingly targeted by threat actors. Using a recent
taxonomy, we contrast the diverse levels of detail given by national
authorities. The threat is commonly acknowledged, but guidance is disjointed.
NIST SP 800-161 aligns closely with the taxonomy and offers a potential
pathway towards a common set of principles.

Cybersecurity Supply Chain Risk Management (C-SCRM) has gained strategic
importance for organizations of all sizes and sectors globally. As
organizations matured in their cybersecurity posture, threat actors needed new
avenues to promote their nefarious activities. The evolution of digitization
and inter-connectivity of the Information and Communications Technology (ICT)
product and service supply chains has introduced a vector of attack that has
seen rich pickings of late [1].

Our previous research produced a C-SCRM taxonomy that is divided into four
categories; Ownership, Risk, Service, and End-2-End. [5]. The taxonomy allows
for a rapid international comparison of the current threat landscape and the
associated risk management guidance to understand whether there is commonality
towards establishing C-SCRM processes, procedures, and practices. Such an
approach is important as supply chains become increasingly global in nature
and transcend national borders. It would also allow suppliers to comply with
broadly common criteria, as opposed to a smorgasbord of national and sectorial
conditions that may create additional complexities and potentially introduce
conflicting requirements.

The NotPetya cyber attack by Russia on Ukraine in 2017 and the more recent
SolarWinds exploit [4] were deployed through the compromise of legitimate
software at source, with the client/supplier trust then used to deploy the
malware.

Incidents like these and the recent Apache Log4Shell/Log4j vulnerability [6]
highlight the increased use of the supply chain in general, and the software
supply chain in particular, as becoming a critical new frontier within
cybersecurity. Organizations want to be able to understand the supply chain
risk and to receive common guidance on how to mitigate it. Our review of North
American, European, and Asian Pacific (APAC) national organizations found a
large discrepancy in the level of detail provided at a national level, ranging
from Highly Detailed to Minimal levels of advice and guidance for both the
threat landscape and risk management.

For our comparison, we looked at major regions and focused on countries with
highly developed cyber security national approaches. This included the USA and
Canada for North America, whilst APAC was represented by Australia, New
Zealand, and Singapore. For Europe, we included the UK, Germany, France, The
Netherlands and, more widely, the European Union Agency for Cybersecurity
(ENISA).

We also identified the US National Institute of Standards and Technology
(NIST) as being a possible resource that authorities and organizations could
gravitate towards for both detailed risk management NIST SP 800-161 and a
supporting cybersecurity framework. This would support the goal of enabling a
set of commonly defined principles that could then be adopted within specific
national and sectoral operational contexts.

## Supply Chain Threat Landscape

### Regional and National Threat Assessments

Accepting the trend of increased targeting of the cyber supply chain, we
looked for progressive approaches taken to identify and communicate the
threats and how they influence their guidance.

The ENISA, in their recent report on the Threat Landscape for Supply Chain
Attacks, mapped and studied supply chain attacks between January 2020 and July
2021. This period included the SolarWinds attack [4] and, more recently, the
Kaseya ransomware attack [3]. The report develops a taxonomy (Figure 1) to
classify supply chain attacks for improved analysis. Half of the 24 supply
chain attacks studied were attributed to well-known Advanced Persistent Threat
(APT) groups. 62 percent of attacks took advantage of the customers’ trust in
their suppliers; 66 percent focused on the suppliers’ code as the vector of
compromise, with data (customer, personal, and intellectual property) being
the most popular target at 58 percent.

Figure 1: ENISA: Proposed taxonomy for supply chain attacks [2]

The UK NCSC predicted the increased vulnerability of the cyber supply chain to
compromise in its 2018 Annual Review. The report included reference to “the
role the supply chain plays in leaving organizations vulnerable to
compromise.”. It evidenced supply chain compromises for both Managed Service
Providers (MSPs) and legitimate software vendors. Compromise of the former
would potentially give access to commercially sensitive data of the MSP’s
clients or access to their clients’ networks via the trusted relationship,
whilst targeting the latter enables a compromise of a product used by multiple
organizations. Similar to many of the other organizations we reviewed, it
maintains this focus on supply chain threats in subsequent annual reviews,
whilst also providing both high-level and detailed advisories relating to
critical supply chain vulnerabilities or compromises.

The APAC output is limited to contents within the respective annual cyber
threat reports. Australia provides the greatest details of the cohort,
focusing on the types of threat actors that target Australian organizations
whilst also highlighting specific campaigns that were reliant upon supply
chain attacks, such as SolarWinds. Expecting similar supply chain attacks to
become prevalent, it focuses on software supply chain compromises and outlines
the difficulty in detecting such attacks. It recommends organizations focus on
incident response plans, almost with the assumption that compromise at some
level is inevitable.

The US Cybersecurity and Infrastructure Security Agency (CISA) launched the
ICT SCRM Task Force in October 2018. It was formed as a public-private
partnership to provide advice and guidance on a number of working groups (WG).
WG2 focuses on supplier, product, and service threat evaluation. The latest
version included nine threat categories:

  * •

Counterfeit Parts

  * •

External Attacks on Operations and Capabilities

  * •

Internal Security Operations and Controls

  * •

System Development Life Cycle (SDLC) Processes and tools

  * •

Insider Threat

  * •

Economic Risks

  * •

Inherited Risk (Extended Supply Chain)

  * •

Legal Risks

  * •

External End-to-End Supply Chain Risks (Natural Disasters, Geo-Political
Issues)

The report builds out a detailed threat list below these nine common threat
groups and identifies associated threat categories and sources. It then
creates threat scenarios against each threat group. This granular level of
detail is ideal for a specialist team, but a higher-level document is required
for the business senior stakeholders engaged in C-SCRM strategic planning.

There is a broad agreement amongst the organizations compared that the risks
encompass third party products and services throughout the life of the
contract. They routinely identify and provide tactical guidance on critical
vulnerabilities that threat actors may target (such as Log4J) or active
campaigns (such as SolarWinds). The level of detail differs for authorities,
but the main crux is the appreciation of the threat that the cyber supply
chain poses, and this should drive national authorities to promote best
practice C-SCRM guidance that organizations large and small in all sectors can
adopt.

Figure 2: Regional and National C-SCRM Guidance

## Cybersecurity Supply Chain Risk Management

The C-SCRM Taxonomy categories in Figure 2 allowed us to compare the different
regional and national levels of guidance to identify the sub-categories
commonly addressed against those receiving less focus. Those sub-categories
comprise the following:

  * •

Ownership: Foreign Ownership, Global Supply Chain, Private Entity.

  * •

Risk: Complex Supply Chain, Risk Management, Attack Surface, Information
Sharing.

  * •

Service:Service Provider, Sub-Contractor, Original Manufacturer, Authorised
Reseller, Hardware and Software.

  * •

End-2-End: Contracts and Procurement, Excluding Supplier, Removing Supplier,
Supply Chain Life Cycle.

We determined the levels of detail offered using the following five-point
scale:

  * •

Highly Detailed: A clear level of detail across a broad spectrum of the sub-
categories.

  * •

Detailed: A broad spectrum was covered, but to a lesser degree of detail.

  * •

Average: Less coverage, but to an adequate level of detail.

  * •

Low Detail: Limited subset of categories/sub-categories covered, but with
little detail.

  * •

Minimal: Although a sub-category is referenced, there is no detail or
definition.

  * •

None.

#### Ownership

All national authorities reviewed acknowledged that the supply chain is a
global entity and that any compromise can transcend national boundaries. The
UK is unique in the European cohort in addressing foreign ownership within the
cyber supply chain and the associated risk that this may create. It extends
that to also considering if data or services are offshored and the
implications if that were the case.

The Australia Cyber Security Centre (ACSC) also focuses heavily on foreign
control, influence, and interference and notes that it can be difficult to
identify the nationality of suppliers, particularly for multinational
corporations. Similarly, the New Zealand NCSC points to regulatory and
compliance requirements being clearly defined when employing international
suppliers. It recommends a security review if a critical supplier has a change
of ownership or major shareholding, or undergoes a merger.

Although the Task Force WG4 has created a highly detailed vendor SCRM
template, CISA doesn’t cover foreign ownership or investment in the
recommended controls.

#### Risk

The Dutch NCSC’s basic cybersecurity guidance doesn’t reference supply chain
risk management. This is despite identifying that small and medium-sized
enterprises (SMEs) often fail to have adequate security controls, yet can form
part of the supply chains of critical processes. It argues that a lack of
transparency means that third-party risk assessments offer no guarantees of
misuse. France provides tactical advice on critical vulnerabilities and
compromises of hardware and software platforms, which is in line with all
other national authorities assessed, but there is no evidence of C-SCRM
advice.

The UK provides good detail on securing the supply chain with a number of
guidance documents. It covers similar ground to others around the attack
surface that is created through a complex supply chain, allowing
vulnerabilities to be inherent. This leads to risk management guidance offered
at a high level, but starts with understanding one’s supply chain and critical
suppliers.

ACSC offers guidance on most of the risk management attributes, albeit with
less detail than some. It calls for cybersecurity incident reporting to be
open and transparent to customers and appropriate authorities in a timely
manner. New Zealand noted that traditional IT is increasingly circumvented as
business units deal directly with cloud and SaaS providers, thus introducing
risk through the side door.

CISA is aligned with most of those reviewed in acknowledging the complexity of
the global digital supply chain and associated attack surface. There is a gap
in guidance on information sharing, but increased detail regarding risk
management advice, especially relating to business resilience strategies, and
towards assurance and compliance.

#### Service

Both the Dutch and UK NCSCs highlight outsourcing services such as cloud and
SaaS as being a risk, although the Dutch do not currently offer any mitigation
guidance. It does, however, note the misuse of Industrial Control System (ICS)
supply chains and the increasing exploitation of vulnerabilities in hardware
and software.

The German BSI focuses on the common practice of software projects, including
software libraries that are reliant on others to maintain and that, if
compromised, can impact countless software platforms. It is promoting the use
of Software Bill of Materials (SBOMs) as a risk control measure. The IT
Security Act 2.0 subjects original manufacturers of critical components (both
hardware and software) to certain obligations that will flow down the supply
chain to protect critical infrastructure.

APAC agencies cover similar ground. ACSC defines the supply chain as
comprising the supplier, manufacturer, distributor, and retailer, as well as
the service provider. It calls for a commitment to secure-by-design products
using secure coding practices. Access control from vendors and service
providers is also covered.

CISA provides a broad and detailed question set within this category, with
most of the sub-categories and attributes covered multiple times. The
requirement for products to be secure by design is particularly detailed, as
is the life cycle maintenance and support.

#### End-2-End

Germany supports European certification schemes, such as consumer and public
Internet of Things (IoT) devices and the need for an entire life cycle across
its supply chain, whilst the UK calls for security controls and assurance to
be contracted and to include the removal of the supplier, be that at contract
end or through the termination of the contract. It was unique in the region to
reference the exclusion of suppliers. It recently published procurement
guidance for business leaders, practitioners, and suppliers. It mirrors a
similar approach taken by New Zealand, which also signposts the need for
contract cancellation provisions that ensure the return or disposal of assets
and data from the supplier.

Cybersecurity expectations should be clearly documented within contracts and
be justifiable, achievable, and proportionate to the risk exposed by the
supplier, advises ACSC. Some suppliers may be excluded from receiving
contracts if the government deems them a national security concern,
particularly for critical infrastructure providers. Consideration for the
supply chain life cycle is referenced repeatedly.

Singapore does not currently offer any specific C-SCRM guidance and didn’t
have sufficient output to populate the individual categories. However, it is
advanced in implementing security ratings for consumer IoT devices.

Assurance that products are not procured from excluded suppliers is included
in CISA guidance, although it doesn’t address the removal of a supplier, be
that at the natural end or termination of a contract. The supply chain
lifecycle does have multiple control questions within the template.

The only C-SCRM guidance from the Canadian Centre for Cyber Security (CCCS) is
limited and provides little detail.

### Threat Landscape and Risk Management comparison

This rapid review reveals that one cannot assume that the level of detail an
authority gives to the threat landscape will be similarly reflected in their
risk management output. Most authorities are closely aligned in their product
output, but there are exceptions.

Whilst ENISA has produced the detailed Threat Landscape for Supply Chain
Attacks, we found no advice on how to manage that risk. It may be delegated
down to individual member states, although, currently, the guidance coming
from those authorities does not reflect this. The supply chain is covered
within the NIS Directive. This would indeed then defer to member states, but
only for operators of essential services and digital service providers.

Conversely, New Zealand covers the threat landscape at a minimal level of
detail, whilst their risk management guidance is highly detailed.

The US provided the highest level of detail across the two areas. In general,
there is slightly more detail in risk management than in threat awareness
amongst national authorities. This would be expected as threat awareness
offers general awareness, whilst risk management must provide specifics on how
to address the threat.

All organizations reviewed offered some level of risk management advice, even
if it was not particularly detailed. UK guidance comes from several documents.
New Zealand and the US give very detailed advice, albeit with different
approaches. New Zealand goes for a more principles-based assurance approach,
whilst CISA outlines each possible control. This may be an issue if the
template issued by CISA doesn’t cover some of the taxonomy sub-categories.

#### Coverage of the C-SCRM Taxonomy

Nine national organizations were included in our review of their guidance
against the C-SCRM taxonomy. Any sub-categories that were referenced seven or
more times were deemed to be popular and are in bold font in Figure 2. The
first thing to note is that there was no reference of the Private Entity sub-
category within Ownership. This attribute relates to shareholder and investor
priorities. There are different business priorities between public and private
organizations. This was acknowledged by the US and UK water sectors and by
most of the frameworks in our previous paper [5], but ignored at an authority
level, and that has been replicated here.

All categories received at least one popular sub-category, with the exception
of End-to-End. Embedding C-SCRM into current and future contracts is key to
transparency and accountability. This is still in its infancy for many
organizations, although guidance and coverage are maturing. This will drive
improvements in the supply chain lifecycle that will include the off-boarding
of suppliers and service providers at contract end or termination.

Coverage of the Global and Complex supply chain are closely aligned and widely
reported. This lends itself to the appreciation of the Attack Surface that
threat actors target. This often involves the Hardware, Software, and Service
Providers, so it is understandable that these areas receive good coverage.

Just because a national authority covers many of the sub-categories doesn’t
mean that they are providing quality advice. The final column provides the
level of detail that we identified from the documents reviewed. It is the
substance of the guidance that defines the level of detail provided.

#### The search for clearly defined principles

There needs to be synergy in C-SCRM to allow businesses to be exposed to the
best and most appropriate authoritative advice and guidance, but also to allow
supply chains and service providers to offer consistent assurances, regardless
of the country or sector they support.

National authorities are increasingly identifying NIST as the institution to
support businesses in understanding and defining their cybersecurity posture,
and this includes C-SCRM guidance. Although a US government agency, the
guidance is freely available globally and transcends international borders.
Their products also often map to complementary standards.

### NIST Standards

#### NIST SP 800-53 Rev 5

Titled "Security and Privacy Controls for Information Systems and
Organizations"; Revision 5 was enhanced to include a number of controls and
enhancements that relate to C-SCRM, including:

  * •

SCRM Control Family (SR)

  * •

SCRM Strategy (PM-30)

  * •

Suppliers of Critical or Mission Essential Items (PM-30(1))

  * •

Integrated Situational Awareness (SI-4(17))

  * •

External System Awareness

  * •

Acquisition Process

  * •

Supply Chain Risk Assessment

  * •

Criticality Analysis

  * •

Incident Handling - Supply Chain Coordination (IR-4(10))

  * •

Incident Reporting - Supply Chain Coordination (IR-6(3))

  * •

Controlled Maintenance (MA-2)

  * •

Tamper Protection (PE-3(5))

It forms a key component of the Risk Management Framework (RMF) that also
incorporates NIST SP 800-37.

#### NIST SP 800-161 Rev 1

Whilst NIST publications have elements of C-SCRM embedded within them, the
newly published NIST SP 800-161 Rev 1 is the go-to document for "C-SCRM
Practices for Systems and Organizations". It integrates C-SCRM into risk
management activities through a multi-level approach that includes the
development of strategic implementation plans, C-SCRM policies, and risk
assessments for products and services that apply to both IT and Operational
Technology (OT) environments, and inclusive of IoT. It promotes the inclusion
of a wide array of stakeholder groups and covers activities that span the
entire System Development Life Cycle (SDLC).

It references NIST publications and legislative developments that have
influenced the content. It also references international standards,
guidelines, and best practice documents, the majority of which are provided by
the International Organization for Standardization and the International
Electrotechnical Commission (ISO/IEC).

When compared against the taxonomy, NIST SP 800-161 is broadly aligned, even
down to the attribute level. Some gaps may be covered by other NIST documents.
Exclusion of suppliers is not directly considered, but it does indicate
adherence to legislative requirements, where this would be covered. The
removal of a supplier is only briefly covered. These are key components of
C-SCRM and would benefit from being more obviously signposted.

There is no consideration that service providers have multiple customers,
which may include direct competitors and require contractual separation of
duties. Additionally, the business priorities of public and private
organizations are often very different. Promoting that awareness would afford
a greater appreciation of that potential dichotomy.

Despite these limitations, coupled with the observation that NIST is a US
government agency focused on promoting US innovation and industrial
competitiveness, the content is highly detailed and inclusive for an
international audience. This lends itself to becoming the pathway towards a
common set of principles. It aligns closely with the C-SCRM taxonomy,
addressing most of the gaps identified in this rapid review of global
guidance. This allows for a complementary approach whereby the taxonomy
promotes an appreciation of what is required, whilst NIST can provide the
systemic approach to addressing those needs. This will allow regional
authorities to deliver similar guidance at a level of detail that satisfies
both their mission statement and their audience needs.

## A glimmer of hope!

The challenge remains for national authorities to come together to speak with
one authoritative voice. On May 11, 2022, Australian, British, Canadian, New
Zealand, and United States agencies jointly issued a cybersecurity advisory to
protect MSPs and customers. Although restricted to traditional political
"5-Eyes" partners, it hopefully announces the start of more collaborative
approaches towards a global common set of C-SCRM principles.

## References

  * [1] Crowdstrike  “Crowdstrike 2022 Global Threat Report”, 2022  URL: https://go.crowdstrike.com/rs/281-OBQ-266/images/Report2022GTR.pdf
  * [2] European Union Agency for Cybersecurity.  “ENISA threat landscape for supply chain attacks.”  LU: Publications Office, 2021  URL: https://data.europa.eu/doi/10.2824/168593
  * [3] “Kaseya Ransomware Attack: Guidance for Affected MSPs and their Customers | CISA”  URL: https://us-cert.cisa.gov/kaseya-ransomware-attack
  * [4] Lavi Lazarovitz  “Deconstructing the SolarWinds breach”  In _Computer Fraud & Security_ 2021.6, 2021, pp. 17–19  DOI: 10.1016/S1361-3723(21)00065-8
  * [5] Colin Topping et al.  “Beware suppliers bearing gifts!: Analysing coverage of supply chain cyber security in critical national infrastructure sectorial and cross-sectorial frameworks”  In _Computers & Security_ 108, 2021, pp. 102324  DOI: 10.1016/j.cose.2021.102324
  * [6] Martin Zugec  “Technical Advisory: Zero-day critical vulnerability in Log4j2 exploited in the wild”  URL: https://businessinsights.bitdefender.com/technical-advisory-zero-day-critical-vulnerability-in-log4j2-exploited-in-the-wild

##### Colin Topping

is the cyber incident director at Rolls-Royce plc. He is also undertaking a
part-time Ph.D. at the Bristol Cyber Security Research Group, University of
Bristol, United Kingdom. This is funded by the National Cyber Security Centre.
His principal research interest is focused on cyber security within the supply
chain in an ever increasing global, technical, and interdependent environment.

##### Ola Michalec

is a Senior Research Associate at the Bristol Cyber Security Research Group,
University of Bristol, United Kingdom. Her interests span the design and
implementation of infrastructure policies regarding climate mitigation and
cyber security.

##### Awais Rashid

is professor of cyber security with the University of Bristol, United Kingdom.
His research interests are in cyber security risk assessment, security of
cyber-physical systems, and human and organizational aspects of security. He
is Director of the EPSRC Centre for Doctoral Training in Trust, Identity,
Privacy and Security in Large-scale Infrastructures, is a fellow of the Alan
Turing Institute and leads projects as part of two UK Research Institutes: the
Research Institute in Trustworthy Industrial Control Systems (RITICS), the
Research Institute on Science of Cyber Security (RISCS) and the National
Centre of Excellence on Cyber Security of IoT (PETRAS). He also heads a major
international effort to develop a Cyber Security Body of Knowledge (CyBOK).


# Where is the Quantum Spin Nematic?: Supplemental Material

Shengtao Jiang (gbsn蒋晟韬)  Department of Physics and Astronomy, University of
California, Irvine, California 92697, USA  Judit Romhányi  Department of
Physics and Astronomy, University of California, Irvine, California 92697, USA
Steven R. White  Department of Physics and Astronomy, University of
California, Irvine, California 92697, USA  M. E. Zhitomirsky  Université
Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France  Institut Laue-Langevin,
71 Avenue des Martyrs, CS 20156, 38042 Grenoble Cedex 9, France  A. L.
Chernyshev  Department of Physics and Astronomy, University of California,
Irvine, California 92697, USA

##  I Multipolar order parameters

We are interested in quantum spin orders that break full rotational symmetry
without displaying dipolar ordering. These phases are characterized by
multipole order parameters such as quadrupole, octupole, etc. Higher quantum
spins have larger local Hilbert spaces and naturally exhibit higher-multipole
degrees of freedom. In general, spins $S$ allows for on-site order parameters,
which transform as rank-$k$ tensor operators, with $k$ taking values of
$0,1,\ldots,2S$.

A rank-$k$ tensor operator $\bm{T}^{(k)}$ has $2k+1$ components that satisfy
the commutation relations [1]

$\displaystyle\big{[}S^{z},T^{(k)}_{q}\big{]}=qT^{(k)}_{q}\qquad\text{and}\qquad\big{[}S^{\pm},T^{(k)}_{q}\big{]}=\sqrt{k(k+1)-q(q\pm
1)}\,\,T^{(k)}_{q\pm 1}\;,$ (1)

which enable a systematic construction of the multipole order parameters
$T^{(k)}_{q}$ with $q\in\left[-k,k\right]$ for each $k$.

For example, the dipole operators correspond to the $k=1$ multiplet with
$T_{1}^{(1)}=S^{+}$, $T_{0}^{(1)}=S^{z}$, and $T_{-1}^{(1)}=S^{-}$, while the
quadrupole operators form a rank-2 tensor with

$\displaystyle T_{2}^{(2)}=S^{+}S^{+},\quad
T_{1}^{(2)}=-S^{+}S^{z}-S^{z}S^{+},\quad
T_{0}^{(2)}=\sqrt{\frac{2}{3}}\big{(}3S^{z}S^{z}-S(S+1)\big{)},\quad
T_{-1}^{(2)}={T_{1}^{(2)}}^{\dagger},\ \ \text{and}\ \
T_{-2}^{(2)}={T_{2}^{(2)}}^{\dagger}.$ (2)

Octupolar operators can be generated from $T_{3}^{(3)}=S^{+}S^{+}S^{+}$ by the
repeated application of $\big{[}S^{\pm},T^{(k)}_{q}\big{]}$, see Eq. (1).

For the $S=\frac{1}{2}$ systems, quadrupole (nematic) orders require a
consideration of the two-spin operators in the form of
$S^{\alpha}_{1}S^{\beta}_{2}$ defined on bonds instead of the on-site spin-
operators in case of the higher spins. The decomposition of such bond
operators consists of the trace, that is the scalar operator
$\bm{S}_{1}\cdot\bm{S}_{2}$, which is SU(2) invariant and commutes with the
total spin, the traceless antisymmetric part $\bm{S}_{1}\times\bm{S}_{2}$,
which is a vector (rank-1 tensor), and the traceless symmetric part
corresponding to spin-quadrupole operators (rank-2 tensor)
$Q^{\alpha\beta}_{12}=S^{\alpha}_{1}S^{\beta}_{2}+S^{\beta}_{1}S^{\alpha}_{2}-\delta_{\alpha\beta}\frac{2}{3}\big{(}\bm{S}_{1}\cdot\bm{S}_{2}\big{)}$.

The quadrupole bond operators can be expressed in the time-reversal invariant
basis

$\displaystyle\left|s\right>=\frac{1}{\sqrt{2}}\big{(}\left|\uparrow\downarrow\right>-\left|\downarrow\uparrow\right>\big{)},\qquad\left|t_{x}\right>=\frac{i}{\sqrt{2}}\big{(}\left|\uparrow\uparrow\right>-\left|\downarrow\downarrow\right>\big{)},\qquad\left|t_{y}\right>=\frac{1}{\sqrt{2}}\big{(}\left|\uparrow\uparrow\right>+\left|\downarrow\downarrow\right>\big{)},\qquad\left|t_{z}\right>=\frac{-i}{\sqrt{2}}\big{(}\left|\uparrow\downarrow\right>+\left|\downarrow\uparrow\right>\big{)},$
(3)

as $Q^{\alpha\beta}_{12}=-\frac{1}{2}\big{(}|t_{\alpha}\rangle\langle
t_{\beta}|+|t_{\beta}\rangle\langle
t_{\alpha}|\big{)}+\frac{1}{3}\delta_{\alpha\beta}\sum_{\gamma}|t_{\gamma}\rangle\langle
t_{\gamma}|$. A bond nematic state can be written as a linear combination of
the triplet states
$\left|\bm{t}\right>=\sum_{\alpha}\bar{t}_{\alpha}\left|t_{\alpha}\right>$
with real coefficients $\bar{t}_{\alpha}$. Such a state is time-reversal
symmetric, and has no dipole order,
$\left<\bm{t}\right|S_{i}^{\alpha}\left|\bm{t}\right>\equiv 0$, as only the
quadrupole matrix elements are finite.

While the expectation values of dipolar operators are zero in a purely
quadrupolar state, the reverse is not true; quadrupolar order parameters can
be finite in a magnetic state. For example, the magnetic state
$\left|\uparrow\uparrow\right>\\!=\\!\frac{1}{\sqrt{2}}\big{(}\left|t_{y}\right>-i\left|t_{x}\right>\big{)}$
gives $\left<S^{z}_{\rm{tot}}\right>=1$ and at the same time
$\left<Q^{zz}\right>=\frac{1}{3}$.

Bond-nematic phases in spin-half frustrated magnets can arise from a
condensation of bound magnon pairs in the field-polarized paramagnet. The two-
magnon state is directly related to the quadrupole operator
${T_{-2}^{(2)}}_{ij}=-S^{-}_{i}S^{-}_{j}$

$\sum_{i,j}\psi_{ij}{T_{-2}^{(2)}}_{ij}|\uparrow\uparrow\ldots\uparrow\rangle=\sum_{i,j}\psi_{ij}\big{(}Q^{x^{2}-y^{2}}_{ij}-iQ^{xy}_{ij}\big{)}|\uparrow\uparrow\ldots\uparrow\rangle\;,$
(4)

with the bound state in the $Q^{x^{2}-y^{2}}_{ij}$ and $Q^{xy}_{ij}$ channels.
The orbital (lattice) symmetry of the bound state is encoded in the
$\psi_{ij}$ coefficients that are determined from the solution of the two-
magnon Schroedinger equation (SE) considered in the next Section. While the
nomenclature of the $Q^{\alpha\beta}_{ij}$ operators is reminiscent of that of
the $L=2$ orbital states, it is unrelated to the $d$-wave nature of the bound
states discussed below.

The ordering of the dipole and quadrupole moments can be expressed via the
spin and quadrupole structure factors, ${\cal S}^{\alpha\beta}({\bf
k})=\langle{S}^{\alpha}_{\bf k}\,{S}^{\beta}_{\bf-k}\rangle$ and ${\cal
Q}^{\alpha\beta}({\bf k})=\langle{Q}^{\alpha\beta}_{\bf
k}\,{Q}^{\alpha\beta}_{\bf-k}\rangle$, respectively. One expects finite Bragg
peaks in both ${\cal S}({\bf k})$ and ${\cal Q}({\bf k})$ in the conventional
magnetic phase, whereas in the nematic state Bragg peaks are present only in
${\cal Q}({\bf k})$.

A general quadrupole-quadrupole correlation function has the form $\langle
Q^{\alpha\beta}_{i,j}Q^{\alpha\beta}_{i+r,j+r^{\prime}}\rangle$, which can be
measured numerically in finite clusters. However, for the symmetry-broken
states, the expectation value of the nematic order parameter $\langle
Q^{\alpha\beta}_{ij}\rangle$ can be also detected directly, similarly to the
ordered moment $\langle S^{z}_{i}\rangle$ in the dipolar magnetic states. This
is used in our DMRG approach, see Secs. VI below and the main text.

##  II Magnon bound states

###  II.1 Magnon interaction

We restrict ourselves to the bound-state problem in a fully polarized FM phase
of an isotropic spin-$S$ Heisenberg model on a Bravais lattice in a field

$\displaystyle{\cal H}=\frac{1}{2}\sum_{i,{\bm{\delta}}}J_{\bm{\delta}}\,{\bf
S}_{i}\cdot{\bf S}_{j}-H\sum_{i}S^{z}_{i},$ (5)

where ${\bm{\delta}}\\!=\\!{\bf r}_{j}-{\bf r}_{i}$ spans all non-zero $n$th-
neighbor couplings that are present in the model.

The kinetic energy of a magnon, or a single-spin-flip eigenstate
$|\psi_{1,{\bf q}}\rangle\\!=\\!1/\sqrt{N}\sum_{i}e^{-i{\bf q}{\bf
r}_{i}}S^{-}_{i}|0\rangle$, relative to the ground-state energy $E_{0}$ of
$|0\rangle\\!=\\!|\uparrow\uparrow\uparrow\dots\rangle$, can be written as [2]

$\displaystyle\varepsilon_{{\bf q}}=H+S\big{(}{\cal J}_{\bf q}-{\cal
J}_{0}\big{)},\ \ \ \ \mbox{with}\ \ \ \ {\cal J}_{\bf
q}=\sum_{\bm{\delta}}J_{\bm{\delta}}e^{i{\bf q}{\bm{\delta}}},$ (6)

where ${\cal J}_{\bf q}\\!=\\!{\cal J}_{-\bf q}$ as
$J_{\bm{\delta}}\\!=\\!J_{-\bm{\delta}}$ on a Bravais lattice. Assuming that
$\varepsilon_{{\bf q}}$ has an absolute minimum at the momentum ${\bf Q}$ and
given the simple Zeeman-like dependence of energy on the field, one can find
the critical field for a transition to a single-spin-flip BEC from the
condition $\varepsilon_{{\bf Q}}\\!=\\!0$, which yields

$\displaystyle H_{s1}=S\big{(}{\cal J}_{0}-{\cal J}_{\bf Q}\big{)}.$ (7)

In the absence of an attraction between magnons, the expression (7) is valid
for any value of $S$. Thus, one refers to $H_{s1}$ and to the single-magnon
BEC transitions itself as the the “classical” ones in this case.

Magnon interaction requires a consideration of the two-spin-flip wavefunction,
whose generic form is given by

$\displaystyle|\psi_{2}\rangle=\frac{1}{2}\sum_{i,j}\psi_{i,j}S^{-}_{i}S^{-}_{j}|0\rangle.$
(8)

The reason why this function cannot be built directly from the basis of
single-spin-flip operators, which are complete and orthogonal in the single-
particle sector, is the non-commutativity of spin operators, the problem
recognized early on [3, 4, 5]. The commonly followed resolution [6, 7] is to
operate directly in the basis of the $|\psi_{2}\rangle$ wavefunctions (8),
which is complete, but not orthogonal. Despite this latter inconvenience, the
exact two-particle Schroedinger equation (SE) can be obtained by standard, if
tedious, manipulations [2] to yield

$\displaystyle(E-\varepsilon_{\frac{\bf K}{2}+{\bf q}}-\varepsilon_{\frac{\bf
K}{2}-{\bf q}})\psi_{\bf K}({\bf q})=\frac{1}{2N}\sum_{\bf p}V_{\bf K}({\bf
q},{\bf p})\psi_{\bf K}({\bf p}),\ \ \ \mbox{with}\ \ \ \psi_{\bf K}({\bf
q})=\frac{1}{N}\sum_{i,j}e^{-i\frac{\bf K}{2}({\bf r}_{i}+{\bf
r}_{j})}e^{-i{\bf q}({\bf r}_{j}-{\bf r}_{i})}\psi_{i,j},$ (9)

and the magnon interaction potential given by

$\displaystyle V_{\bf K}({\bf q},{\bf p})={\cal J}_{{\bf q}-{\bf p}}+{\cal
J}_{{\bf q}+{\bf p}}-{\cal J}_{\frac{{\bf K}}{2}-{\bf q}}-{\cal J}_{\frac{\bf
K}{2}+{\bf q}}.$ (10)

The first two terms in $V_{\bf K}({\bf q},{\bf p})$ in (10) adhere to the
conventional “potential-like” form, as they depend on the momentum transfer
between initial and final particles in the scattering, with the momenta ${\bf
K}/2\pm{\bf q}$ and ${\bf K}/2\pm{\bf p}$, respectively. The last two terms
depend on the individual particle momenta and originate from the “hard-core”
repulsion of spin flips.

Given the single-magnon energy in (6) and magnon interaction potential in
(10), it is made explicit by SE (9) that the magnon interaction is manifestly
a $1/S$ effect, with the role of pairing diminishing in the large-$S$ limit
and single-magnon BEC retaining its classical value.

We also would like to remark that an expression identical to (10) is trivially
obtained in the $O(S^{0})$ order of expansion of (5) within the perfectly
commutative and orthogonal bosonic basis of the standard $1/S$-approximation,
obviously leading to a form of the two-magnon SE that is identical to (9).

###  II.2 Simple case, Cooper problem analogy, and higher harmonics

It is instructive to have a general intuition on solving SE in (9) for a
simple case. Consider a continuum limit for particles of mass $m$ and energies
$\varepsilon_{\bf k}\\!=\\!{\bf k}^{2}/2m$ that are attracted via a
$\delta$-functional potential in real space, $V({\bf r}_{1},{\bf
r}_{2})\\!=\\!-\alpha\delta({\bf r}_{1}-{\bf r}_{2})$, which corresponds to a
constant in the momentum space, $V_{\bf K}({\bf q},{\bf p})\\!=\\!-\alpha$.
This reduces (9) to an algebraic equation for $\psi$ that is trivially solved
as

$\displaystyle\psi_{\bf K}({\bf q})=-\frac{\alpha}{2}\cdot\frac{C_{\bf
K}}{E-E_{\bf K}-{\bf q}^{2}/m},\ \ \ \ \mbox{with}\ \ \ \ C_{\bf
K}=\frac{1}{N}\sum_{\bf p}\psi_{\bf K}({\bf p}),$ (11)

where $E_{\bf K}\\!=\\!{\bf K}^{2}/4m$ is the bottom of the two-particle
continuum for the total momentum ${\bf K}$. The second equation in (11)
provides a self-consistency condition, which is also an implicit equation on
the pair’s energy $E$

$\displaystyle 1=-\frac{\alpha}{2N}\sum_{\bf p}\frac{1}{E-E_{\bf K}-{\bf
p}^{2}/m}.$ (12)

Finding $E\\!<\\!E_{\bf K}$ from (12) constitutes a solution of the bound-
state problem. This consideration is analogous to a textbook problem of the
bound state for a single particle in a free space with a $\delta$-functional
potential well [8], also demonstrating a mapping of the two- to one-body
problem. A direct integration in Eq. (12) in the $D$-dimensional continuum
demonstrates the existence of a bound state for an arbitrarily weak attraction
in the dimensions $D\\!=\\!1$ and $D\\!=\\!2$, and requires $\alpha$ to exceed
a threshold value of order $1/m$ in $D\\!=\\!3$ [8]. Since we are interested
in the case of $D\\!=\\!2$, the solution for it is

$\displaystyle E-E_{\bf K}=\frac{\Lambda}{m}\cdot e^{-\frac{8\pi}{\alpha m}},$
(13)

where $\Lambda$ is a large-momentum cut-off and the non-analytic dependence on
the coupling constant $\alpha$ is clear. These results rely exclusively on the
density of states of the non-relativistic particles in $D$-dimensions and are
expected to hold for any short-range attractive potential in the continuum
[8]. One of the well-known realizations of such conditions is the celebrated
Cooper problem of the short-range attraction of electrons near the Fermi
surface [9]. The presence of the finite Fermi-momentum cutoff leads to an
effective dimensional reduction from 3D to 2D, resulting in the 2D-like
density of states and the pairing gap that closely follows Eq. (13).

For a realization of the pair-BEC, the lowest bound-state energy needs to be
below the absolute minimum of the two-particle continuum, which typically
occurs at a high-symmetry ${\bf K}$-point. In an example above, this is the
${\bf K}\\!=\\!0$ point. There, a systematic symmetry consideration of the
pairing in the higher partial-wave channels is also possible. Given the
mapping onto the one-body problem in the continuum, it can be made explicit
that with an exception of the $s$-wave channel, which is exposed above, all
other channels have nodes of their corresponding pairing potential that pass
through the minimum of the single-particle energy. This crucial feature leads
to an effective dimensional increase, as opposed to the dimensional reduction,
suggesting that the weak attraction in the continuum in 2D can only create
$s$-wave bound states, but not the ones in the higher-harmonics, which would
require the corresponding coupling strengths to exceed threshold values of
order $1/m$.

###  II.3 Pairing in the $J_{1}$–$J_{2}$ Heisenberg model

We now turn to the problem of pairing in the ferro-antiferromagnetic $S=1/2$
$J_{1}$–$J_{2}$ Heisenberg model on the square lattice with the ferromagnetic
$J_{1}\\!=\\!-1$ taken as a unit of energy and antiferromagnetic
$J_{2}\\!>\\!0$. The convenient shorthand notation ${\cal J}_{\bf q}$ in (6)
can be rewritten using standard hopping amplitudes

$\displaystyle{\cal J}_{\bf q}=-4\gamma_{\bf q}+4J_{2}\gamma^{(2)}_{\bf q},\ \
\ \mbox{with}\ \ \ \gamma_{\bf q}=\frac{1}{2}\big{(}\cos q_{x}+\cos
q_{y}\big{)},\ \ \ \mbox{and}\ \ \ \gamma^{(2)}_{\bf q}=\cos q_{x}\cos q_{y}.$
(14)

The single-magnon energy (6) for $S=1/2$ is given by

$\displaystyle\varepsilon_{{\bf q}}=H+2\big{(}1-\gamma_{\bf
q}\big{)}-2J_{2}\big{(}1-\gamma^{(2)}_{\bf q}\big{)}=\varepsilon_{{\bf
Q}}-2\gamma_{\bf q}+2J_{2}\big{(}1+\gamma^{(2)}_{\bf q}\big{)},$ (15)

where we use that for $J_{2}\\!>\\!0.5$, the minima of $\varepsilon_{{\bf q}}$
are at the non-trivial ${\bf Q}\\!=\\!(0,\pi)[(\pi,0)]$ points, with
$\varepsilon_{{\bf Q}}\\!=\\!H+2-4J_{2}$, leading to $H_{s1}\\!=\\!-2+4J_{2}$,
see (7). This is the range of $J_{2}$ that we will consider below, although
extension to the case of $J_{2}\\!<\\!0.5$, for which the minimum of
$\varepsilon_{{\bf q}}$ is at ${\bf Q}\\!=\\!0$, is easily made.

In addition to the hopping amplitude $\gamma_{\bf q}$ in (14), one can see an
immediate utility in introducing other combinations of lattice harmonics of
different symmetry

$\displaystyle\gamma^{-}_{\bf q}=\frac{1}{2}\big{(}\cos q_{x}-\cos
q_{y}\big{)},\ \ \ \gamma^{s_{xy}}_{\bf q}=[\gamma^{(2)}_{\bf q}=]\cos
q_{x}\cos q_{y},\ \ \ \mbox{and}\ \ \ \gamma^{d_{xy}}_{\bf q}=\sin q_{x}\sin
q_{y},$ (16)

as they permit an easy symmetrization of the single-magnon energies in SE (9)
and a separation of interaction potential (10) into the partial harmonics
using

$\displaystyle\gamma_{{\bf q}-{\bf p}}+\gamma_{{\bf q}-{\bf p}}=2\gamma_{\bf
p}\gamma_{\bf q}+2\gamma^{-}_{\bf p}\gamma^{-}_{\bf q},\ \ \ \ \ \
\gamma^{(2)}_{{\bf q}-{\bf p}}+\gamma^{(2)}_{{\bf q}-{\bf
p}}=2\gamma^{s_{xy}}_{\bf p}\gamma^{s_{xy}}_{\bf q}+2\gamma^{d_{xy}}_{\bf
p}\gamma^{d_{xy}}_{\bf q}.$ (17)

Introducing binding energy of the pair $2\Delta\\!>\\!0$ relative to the
absolute minimum of the two-magnon continuum via
$E\\!=\\!-2\Delta+2\varepsilon_{{\bf Q}}$, converts the energy difference in
the left-hand side of SE in (9) to

$\displaystyle E-\varepsilon_{\frac{\bf K}{2}+{\bf q}}-\varepsilon_{\frac{\bf
K}{2}-{\bf q}}=-2\Delta+4\big{(}\gamma_{\frac{\bf K}{2}}\gamma_{\bf
q}+\gamma^{-}_{\frac{\bf K}{2}}\gamma^{-}_{\bf
q}\big{)}-4J_{2}\big{(}1+\gamma^{s_{xy}}_{\frac{\bf K}{2}}\gamma^{s_{xy}}_{\bf
q}+\gamma^{d_{xy}}_{\frac{\bf K}{2}}\gamma^{d_{xy}}_{\bf q}\big{)},$ (18)

while the interaction potential (10) becomes

$\displaystyle V_{\bf K}({\bf q},{\bf p})=-8\gamma_{\bf q}\big{(}\gamma_{\bf
p}-\gamma_{\frac{\bf K}{2}}\big{)}-8\gamma^{-}_{\bf q}\big{(}\gamma^{-}_{\bf
p}-\gamma^{-}_{\frac{\bf K}{2}}\big{)}+8J_{2}\gamma^{s_{xy}}_{\bf
q}\big{(}\gamma^{s_{xy}}_{\bf p}-\gamma^{s_{xy}}_{\frac{\bf
K}{2}}\big{)}+8J_{2}\gamma^{d_{xy}}_{\bf q}\big{(}\gamma^{d_{xy}}_{\bf
p}-\gamma^{d_{xy}}_{\frac{\bf K}{2}}\big{)},$ (19)

with the terms clearly assuming the separable partial-wave decomposition
structure, although with asymmetric “shifts” of the harmonics by
$\gamma^{\gamma}_{\frac{\bf K}{2}}$ due to the hard-core components of the
original potential (10).

Figure S1: The two-magnon continuum (shaded area) and the bound state energies
(lines) along the representative ${\bf K}=[q_{x},q_{y}]$ cuts in the Brillouin
zone for (a) $J_{2}=0.6$ and (b) $J_{2}=1.0$. $M=[\pi,\pi]$ and $X=[\pi,0]$.

For the magnon-band minima at ${\bf Q}\\!=\\!(0,\pi)[(\pi,0)]$ points, there
only two distinct total momenta of their pair at which pair-BEC can occur,
${\bf K}\\!=\\!0$ and ${\bf K}\\!=\\!(\pi,\pi)$. While rigorous arguments can
be given to prove that the ${\bf K}\\!=\\!(\pi,\pi)$ sector does not support
pairing in any orbital channel in the range of $J_{2}\\!>\\!1/\pi$, we skip
them in favor of the visual demonstration of the numerical solution for the
bound state problem (9) with the energy denominator (18) and pairing potential
(19) for ${\bf K}$ along a representative cut of the Brillouin zone and for
two representative $J_{2}=0.6$ and $J_{2}=1.0$, see Fig. S1. This Figure
clearly points to the importance of the ${\bf K}\\!=\\!0$ sector in our case,
considered next.

####  II.3.1 ${\bf K}\\!=\\!0$ sector

For the ${\bf K}\\!=\\!0$ sector, $\gamma_{\frac{\bf
K}{2}}\\!=\\!\gamma^{s_{xy}}_{\frac{\bf K}{2}}\\!=\\!1$ and
$\gamma^{-}_{\frac{\bf K}{2}}\\!=\\!\gamma^{d_{xy}}_{\frac{\bf
K}{2}}\\!=\\!0$, converting the energy part of SE (18) to

$\displaystyle\Delta E_{0}({\bf q})=\big{(}E-\varepsilon_{\frac{\bf K}{2}+{\bf
q}}-\varepsilon_{\frac{\bf K}{2}-{\bf q}}\big{)}\big{|}_{{\bf
K}=0}=-2\Delta-4J_{2}\big{(}1+\gamma^{(2)}_{\bf q}\big{)}+4\gamma_{\bf
q}=-2\Delta-2\big{(}\varepsilon_{\bf q}-\varepsilon_{\bf Q}\big{)}<0,$ (20)

where the binding energy $\Delta$ is relative to the $2\varepsilon_{\bf Q}$
minimum, and the interaction potential (19) becomes,

$\displaystyle V_{0}({\bf q},{\bf p})=-8\gamma_{\bf q}\big{(}\gamma_{\bf
p}-1\big{)}-8\gamma^{-}_{\bf q}\gamma^{-}_{\bf p}+8J_{2}\gamma^{s_{xy}}_{\bf
q}\big{(}\gamma^{s_{xy}}_{\bf p}-1\big{)}+8J_{2}\gamma^{d_{xy}}_{\bf
q}\gamma^{d_{xy}}_{\bf p}=\sum_{\gamma}V_{\gamma}R_{\gamma}({\bf
q})\widetilde{R}_{\gamma}({\bf p}),$ (21)

with $\gamma\\!=\\!\\{s,d,s_{xy},d_{xy}\\}$,
$V_{\gamma}\\!=\\!\\{8,8,-8J_{2},-8J_{2}\\}$, $R_{\gamma}({\bf
q})\\!=\\!\\{\gamma_{\bf q},\gamma^{-}_{\bf q},\gamma^{s_{xy}}_{\bf
q},\gamma^{d_{xy}}_{\bf q}\\}$, and $\widetilde{R}_{\gamma}({\bf
p})\\!=\\!\\{\gamma_{\bf p}-1,\gamma^{-}_{\bf p},\gamma^{s_{xy}}_{\bf
p}-1,\gamma^{d_{xy}}_{\bf p}\\}$.

While it is tempting to conclude that this form of interaction already
guarantees decomposition of the SE in (9) into a fully orthogonal set of
algebraic equations on pairing in separate channels, the analysis provided so
far is based entirely on the lattice-harmonic expansion of the interaction
potential and magnon energies. Continuing without any symmetry analysis for a
moment, the SE (9) for the ${\bf K}\\!=\\!0$ sector can be rewritten as

$\displaystyle\psi_{0}({\bf
q})=-\frac{1}{2}\sum_{\gamma}V_{\gamma}R_{\gamma}({\bf
q})\cdot\frac{C_{\gamma}}{\Delta E_{0}({\bf q})},\ \ \ \ \mbox{with}\ \ \ \
C_{\gamma}=\frac{1}{N}\sum_{\bf p}\widetilde{R}_{\gamma}({\bf p})\psi_{0}({\bf
p}).$ (22)

Although it is clear that the SE is, indeed, reduced to a set of algebraic
equations for $\psi_{0}({\bf q})$, these equations are, or can be, coupled, as
the constants $C_{\gamma}$ in (22) can have contributions from the
$\gamma^{\prime}\\!\neq\\!\gamma$ components of $\psi_{0}({\bf q})$.

The partial waves in (21) and (22) are of the $s$, $s_{xy}$,
$d_{x^{2}-y^{2}}$, and $d_{xy}$ character, with the mirror-plane symmetries at
the ${\bf K}\\!=\\!0$ point separating $d_{x^{2}-y^{2}}$ and $d_{xy}$ partial-
wave solutions of the SE, while the two generalized $s$-wave-like solutions
are allowed to mix. Thus, Eq. (22) for the ${\bf K}\\!=\\!0$ sector decouples
into two independent algebraic equations for the $d$-waves and two coupled
equations for the $s$-waves.

The resultant equations on the pairing energy for the $d_{x^{2}-y^{2}}$ ($d$
for brevity) and $d_{xy}$ waves are given by

$\displaystyle d:\ \ \ \ \ 1=-\frac{4}{N}\sum_{\bf p}\frac{(\gamma^{-}_{\bf
p})^{2}}{\Delta E_{0}({\bf p})},\ \ \ \ \ \ \ \ \ \mbox{and}\ \ \ d_{xy}:\ \ \
\ \ 1=\frac{4J_{2}}{N}\sum_{\bf p}\frac{(\gamma^{d_{xy}}_{\bf p})^{2}}{\Delta
E_{0}({\bf p})},$ (23)

of which the $d_{xy}$-wave does not permit any solutions having negatively-
defined right-hand side. The $d$-wave equation yields a solution of our
ultimate interest that is discussed in more detail below.

The coupled $s$-wave equations can be written as

$\displaystyle C_{s}=C_{s}I_{11}+C_{s_{xy}}I_{12},\ \ \ \mbox{and}\ \ \
C_{s_{xy}}=C_{s_{xy}}I_{22}+C_{s}I_{21},\ \ \ \mbox{with}\ \ \
I_{11}=-\frac{4}{N}\sum_{\bf p}\frac{\gamma_{\bf p}(\gamma_{\bf p}-1)}{\Delta
E_{0}({\bf p})},$ (24) $\displaystyle I_{22}=\frac{4J_{2}}{N}\sum_{\bf
p}\frac{\gamma^{s_{xy}}_{\bf p}(\gamma^{s_{xy}}_{\bf p}-1)}{\Delta E_{0}({\bf
p})},\ \ \ I_{12}=\frac{4J_{2}}{N}\sum_{\bf p}\frac{\gamma^{s_{xy}}_{\bf
p}(\gamma_{\bf p}-1)}{\Delta E_{0}({\bf p})},\ \ \
I_{21}=-\frac{4}{N}\sum_{\bf p}\frac{\gamma_{\bf p}(\gamma^{s_{xy}}_{\bf
p}-1)}{\Delta E_{0}({\bf p})}.$ (25)

While, ultimately, one can simply use these expressions to numerically prove
that they are unable to yield a viable solution for the binding energy, having
an intuition on “why” would be very useful. For that, let us consider only the
diagonal $I_{\gamma\gamma}$ terms, because if they are unable to sustain a
solution on their own, the situation is unlikely to improve by the coupling
between the two channels.

For the $I_{22}$ in the $s_{xy}$-wave channel, the situation is similar to the
$d_{xy}$-wave above in that the coupling is due to a repulsive $J_{2}$ term,
but it contains a “hard-core shift” in the pairing potential, which makes the
outcome less certain analytically. However, with a minimal effort, one is able
to prove that $I_{22}$ does indeed remain negative and yields no feasible
solutions for the binding energy,

The interaction in the $s$-wave channel in $I_{11}$ is attractive
($J_{1}\\!=\\!-1\\!<\\!0$), but it has a hard-core shift in the pairing
potential, which is repulsive and is often quoted as the reason for the
suppression of the $s$-wave pairing. While this is a valid mechanism for the
lack of pairing in the pure ferromagnetic cases, which never result in the
bound states that occur below the absolute minimum of the continuum, this is
not the reason why the $s$-wave is projected out of a pairing in our case. A
counterexample is always useful. If we ignore the momentum-dependence in
$\Delta E_{0}({\bf q})$ in the denominator of the integral in $I_{11}$ in (24)
for a moment (flat magnon bands), it will yield a solution for the bound-state
energy that is actually degenerate with the one for the $d$-wave in (23). Yet,
for the true $\Delta E_{0}({\bf q})$ from (20), the $s$-wave solution is
absent. Hence, the lack of the binding is elsewhere.

The resolution is that the extended $s$-wave pairing harmonic $\gamma_{\bf
p}\\!\propto\\!(\cos q_{x}+\cos q_{y})$ in (24) has nodes of zeros passing
through the energy minima of the denominator $\Delta E_{0}({\bf p})$, which
are at the magnon energy minima ${\bf Q}\\!=\\!(0,\pi)[(\pi,0)]$, leading to a
suppressed attraction via the “dimensional increase” mechanism and giving no
bound-state solution. For the smaller values of $J_{2}\\!<\\!0.5$, the magnon
band minima migrate back to the ferromagnetic minima at ${\bf Q}\\!=\\!0$,
where the “hard-core” reasoning takes over in forbidding the ${\bf K}\\!=\\!0$
bound state in the $s$-wave channel.

Altogether, the $d$-wave ${\bf K}\\!=\\!0$ magnon bound-state in (23) is the
only candidate to have a BEC in the ferro-antiferromagnetic $J_{1}$–$J_{2}$
Heisenberg model on the square lattice. We consider it in some more details
next.

####  II.3.2 $D$-wave solution and asymptotics

We rewrite the equation on the $d$-wave pairing energy (23) in two different
ways

$\displaystyle 1=\frac{1}{2N}\sum_{\bf p}\frac{\big{(}\cos p_{x}-\cos
p_{y}\big{)}^{2}}{\Delta+2J_{2}\big{(}1+\cos p_{x}\cos
p_{y}\big{)}-\big{(}\cos p_{x}+\cos p_{y}\big{)}},\ \ \ \ \mbox{and}\ \
1=-\frac{V_{d}}{2N}\sum_{\bf p}\frac{\big{(}\gamma^{-}_{\bf
p}\big{)}^{2}}{-2\Delta-2\big{(}\varepsilon_{\bf p}-\varepsilon_{\bf
Q}\big{)}},$ (26)

with the first one being handy to use in the numerical integration that
results in the curves for $\Delta$ vs $J_{2}$ in Fig. 2(a) and
$h_{c2}\\!=\\!h_{c1}+\Delta$ in Fig. 1(a) of the main text. The second one,
with $V_{d}\\!=\\!8$ in units of $|J_{1}|$ as before, is useful to relate back
to the very beginning of our bound-state consideration, the $s$-wave pairing
in the continuum. Taking the $J_{2}\\!\gg\\!1$ limit for the denominator in
(26) and expanding near the magnon energy minima ${\bf
Q}\\!=\\!(0,\pi)[(\pi,0)]$ yields $\Delta E_{0}({\bf p})\\!=\\!-2\Delta-{\bf
p}^{2}/m_{\bf Q}$ with $1/m_{\bf Q}\\!=\\!2J_{2}$.

Figure S2: (a) Numerical results for the $d$-wave bound-state wave-function
$\psi_{d}({\bf r}_{x})$ on a semi-log plot for three representative values of
$J_{2}$ together with the asymptotic result $\psi_{d}({\bf
r})\\!=\\!Ae^{-|{\bf r}|/\xi}$, with the correlation length
$\xi\\!=\\!0.8e^{\pi J_{2}/2}$ for all three sets. (b) Same as in Fig. 1(a)
and Fig. 2(a) of the main text, with the analytical results for $\Delta$ in
(26) obtained using integration on an $L_{x}\times L_{y}$ lattice, with
effective $L_{x}=\infty$ and two values of the width, $L_{y}=8$ and
$L_{y}=12$.

Since the zero-nodes of the pairing harmonic $\gamma^{-}_{\bf p}$ are along
the $q_{x}\\!=\\!\pm q_{y}$ diagonals, perfectly missing the magnon band
minima, the coupling of the low-energy magnons is finite at ${\bf Q}$,
allowing us to approximate $(\gamma^{-}_{\bf p})^{2}\\!\approx\\!1$ near ${\bf
Q}$ and giving us all the ingredients of the “$s$-wave-like” pairing. Thus, in
the weak-coupling limit, integration in (26) yields the answer that mirrors
that of Eq. (13)

$\displaystyle 2\Delta\approx\frac{\Lambda}{m_{\bf Q}}\cdot
e^{-\frac{4\pi}{m_{\bf Q}V_{d}}}=2J_{2}\Lambda\cdot e^{-\pi J_{2}},$ (27)

with $\Lambda$ being the large-momentum cutoff and an extra factor $1/2$ in
the exponent originating from the two inequivalent minima of the magnon band
at ${\bf Q}\\!=\\!(0,\pi)$ and $(\pi,0)$. Thus, the pairing in the $d$-wave
channel in the considered problem indeed maps onto the textbook $s$-wave
solution. One of the key consequence of that is the absence of the upper limit
on $J_{2}$ for such a pairing to occur.

The other consequence is the exponential decay of the pair wave-function with
the distance, $\psi_{d}({\bf r})\\!\sim\\!e^{-|{\bf r}|/\xi}$, with the
correlation length that can be inferred from the virial theorem, ${\bf
q}_{\xi}^{2}/2m_{\bf Q}\\!\simeq\\!\Delta$, as $\xi\\!\sim\\!e^{\pi J_{2}/2}$.
In Figure S2(a), we show the real-space amplitudes of the bound-state wave-
function $\psi_{d}({\bf r})$ from (22) along the line ${\bf r}=(r_{x},0)$ for
several representative values of $J_{2}$ together with the asymptotic
exponential behavior with $\xi\\!=\\!0.8e^{\pi J_{2}/2}$. Lastly, for
$J_{2}\\!<\\!0.5$, the magnon band minimum shifts to ${\bf Q}\\!=\\!0$,
quickly terminating pairing at $J_{2}^{c}\\!=\\!0.4077593304754(1)$, as can be
inferred numerically from (26). This is the end-point of the $h_{c2}$ line in
Fig. 1(a) in zero field and $\Delta_{0}\\!=\\!0$ point in Fig. 2(a) of the
main text, see also Fig. S2(b) and its inset.

One can emulate the finite-size effects of the DMRG results for the paring
energy $\Delta$ that can be seen in Fig. 1(a) and Fig. 2(a) of the main text,
by taking the analytical integral for $\Delta$ in (26) on an $L_{x}\times
L_{y}$ lattice. The results of such an effort, with two values of $L_{y}$ and
effectively infinite $L_{x}$ are shown in Fig. S2. The results are nearly
matched by the width $L_{y}=12$, not $L_{y}=8$, which can be interpreted in
favor of the open-periodic boundary conditions employed by DMRG as giving a
better approximation of the thermodynamic limit than its nominal width would
suggest [10].

##  III multi-pair states

In this Section, we provide a comparison of the energies of the states
obtained from the DMRG calculations in the $16\\!\times\\!8$ cluster with
fixed total $S^{z}$, which gives us the energies of the multi-pair states as a
function of $J_{2}$. The energy of the $n$-magnon state is denoted as $E(n)$
and corresponds to the state with $n$ spin flips from the fully polarized
ferromagnet. In the following, we consider only even values of $n$, because
the states with odd values have consistently higher energies than the nearest
even ones, indicating that the gap to single-magnon excitations is always
present in the considered range of $J_{2}\\!>\\!0.39$ and for the low
concentration of spin flips, $n/N_{site}\ll 1$.

In Fig. S3(a) and Fig. S3(b), we show two schematic plots of the total energy
of the system $E(n)$ vs magnon number $n$ for attractive and repulsive magnon
pairs, respectively. The $n\\!=\\!0$ state is a fully polarized state and we
do not include external field into consideration yet. In either case, the
energy is lowered by the increase of the magnon density because the fully
polarized state is a highly excited state in the regime where an
antiferromagnetic state is a ground state of the system in zero field. If
magnon pairs are attractive, the (negative) energy gain by $n+2$ magnons is
always greater than the sum of such energy gains by $n$ and $2$ magnons
separately. Therefore, the $E(n)$ vs $n$ curve should be concave, as is shown
in Fig. S3(a). In this case, the system can not be stabilized at a low magnon-
pair density. Depending on the magnetic field, the ground state is either a
saturated FM state or a multi-pair state with a large concentration of
magnons, which is beyond the considered low-density regime. On the other hand,
if the magnon pairs are repulsive, the energy curve should be convex, as is
shown in Fig. S3(b), a state with a fixed magnon-pair density can be
stabilized by a field, and the magnon pair-density in the ground state will
increase continuously upon lowering of the field. In the following
consideration, it is convenient to count the energy $E(n)$ from the energy
$E(0)$ of the fully polarized state, $\bar{E}(n)\equiv E(n)-E(0)$.

Figure S3: (a) and (b) Schematic plots of the total energy $E(n)$ vs magnon
number $n$ without the magnetic field. The slope of the line connecting $E(0)$
and $E(n)$ is the multi-pair BEC instability field $h_{cn}$. The energy vs
magnon number curve is concave (convex) for attractive (repulsive)
interaction. (c) The total energy of the system with $n$ magnons $\bar{E}(n)$
relative to the energy of $\frac{n}{2}$ non-interacting pairs
$\frac{n}{2}\bar{E}(2)$ for $J_{2}=0.7$ and $J_{2}=0.5$ in a $16\times 8$
cylinder. (d) The $h$–$J_{2}$ phase diagram for the states with fixed number
of magnons $n=2,4,6,8$ in a $16\times 8$ cylinder. For $J_{2}>0.61$, magnon
pairs are repulsive and the leading instability is the single-pair-BEC
($n=2$), while for $J_{2}<0.61$ the magnon pairs are attractive, which leads
to the multi-pair instabilities with narrowing $J_{2}$ steps that resemble a
devil’s staircase.

In Fig. S3(c) we show the energy of the system with $n$ magnons, $\bar{E}(n)$,
relative to the energy of $\frac{n}{2}$ non-interacting pairs,
$\frac{n}{2}\bar{E}(2)$, for $J_{2}=0.7$ and $J_{2}=0.5$, obtained from DMRG
calculations in a $16\times 8$ cylinder. For $J_{2}=0.7$, this energy is
positive, indicating that adding another pair to the already existing ones
cost extra energy due to repulsion between pairs. On the other hand, for
$J_{2}=0.5$, this energy is negative due to attraction between pairs.

To identify which multi-pair state has the leading instability with respect to
the fully polarized state, we calculate and compare the $n$-magnon instability
field $h_{cn}$, defined as
$h_{cn}\\!=\\!\big{(}E(0)-E(n)\big{)}/n\\!=\\!-\bar{E}(n)/n$, for $n=2,4,6,8$
in a $16\times 8$ cylinder. In the schematic plots in Figs. S3(a) and S3(b),
the values of $h_{cn}$ corresponds to the discrete versions of the (negative
of the) slopes of the lines connecting $E(0)$ and $E(n)$. In the case of
attractive magnon pairs, if the multi-pair states are formed that still repel
each other, the instability field is achieved at $h_{cn}>h_{c2}$. The biggest
achievable $h_{cn}$ corresponds to the leading instability.

In Fig. S3(d), we present the $h$–$J_{2}$ phase diagram showing such
instability fields for the states with fixed number of magnons in a $16\times
8$ cylinder for $J_{2}$ from 0.5 to 0.7. To obtain this phase diagram, we
calculated the total energy of the system in magnetic field $h$ as a function
of $J_{2}$, $\widetilde{E}_{n}(h,J_{2})=E_{n}(J_{2})+nh$, where $E_{n}(J_{2})$
is the total energy of the system with $n$ magnons without magnetic field as
above. For each $(h,J_{2})$ point in this diagram, the $n$-magnon state with
the lowest $\widetilde{E}_{n}(h,J_{2})$ is chosen. The phase boundary between
the $n$-magnon state and the fully polarized state with $n=0$ corresponds to
the $h_{cn}(J_{2})$ discussed above.

For $J_{2}\gtrsim 0.61$, magnon pairs are repulsive and the leading
instability field is $h_{c2}$. The states with more pairs are stabilized at a
lower field. For $J_{2}\lesssim 0.61$, the leading instability changes from
$h_{c2}$ to $h_{c2n}$ in steps of 2, with $h_{c2n}$ increasing as $J_{2}$ is
decreased. At the same time, the $J_{2}$-ranges in which each $h_{c2n}$
dominates shrink with the increase of $n$, resembling a devil’s staircase of
an infinite sequence of transitions with the ever decreasing width. Although
the shift of the leading instability away from the $h_{c2}$ line in favor of
the multi-pair $h_{c2n}$ in our analysis is beyond doubt, the finite-size
effects prevent a definite confirmation of the devil’s staircase scenario.

In the attractive regime there are two scenarios. If multi-pair objects repel,
the system will be stabilized at $h_{cn}$. Near the boundary of the two-magnon
stability ($J_{2}\lesssim 0.61$), this enables nematic state to survive for a
certain field range below $h_{cn}$ before entering the AFM state at a lower
field. For a yet smaller $J_{2}$ ($J_{2}\lesssim 0.5$), and stronger
attraction, the multi-pair state cannot be stabilized for any finite $n$ and
the system undergoes a first-order transition to a strongly canted AFM state.
This transition is discussed in Sec. VII in more detail.

##  IV reduction of the single-magnon gap in the presence of pairs

As is discussed in the main text, the actual nematic phase occupies only a
fraction of the naïvely anticipated region between $h_{c1}$ and $h_{c2}$, with
the rest of that region taken by the AFM phase. Since $h_{c1}$ corresponds to
a closing of the single-magnon gap for the non-interacting magnons, an
expansion of the AFM phase above $h_{c1}$ suggests that the single magnon-
excitations can reduce their gap due to interaction with the existing pair-
condensate. Here, we substantiate this picture by providing quantitative data
for the reduction of the single-magnon gap based on the DMRG calculations in
the $16\\!\times\\!8$ cluster with fixed total $S^{z}$ for the representative
value of $J_{2}=0.7$.

Figure S4: (a) The total energy of the system with $n$ magnons relative to the
magnon pairing energy in a pair per magnon, $\frac{1}{2}\bar{E}_{2}$, times
$n$ vs $n$ in the $16\\!\times\\!8$ cluster with fixed total $S^{z}$ for
$J_{2}=0.7$ and $n$ from 1 to 10. Red line is the even-$n$ sector and blue
line is the odd-$n$ sector. (b) The single-magnon gap $\Delta E(n_{odd})$ vs
$n$.

First, we look at the total energy relative to the magnon pairing energy, as
is shown in Fig. S4(a). The total energy of the $n$-magnon state $\bar{E}(n)$
is the energy relative to the energy of the fully polarized state:
$\bar{E}(n)=E(n)-E(0)$. Considering it relative to the magnon pairing energy
in a single pair per magnon times $n$, $\frac{n}{2}\bar{E}(2)$, is similar to
introducing a chemical potential, so that energies of the states with a
different magnon number can be compared. One can see in Fig. S4(a), that for
the states with even values of $n$, in which all magnons form pairs, the
relative energy increases with $n$. This effect is simply due to repulsion
between magnon pairs in a finite cluster. However, for the states with the odd
$n$, in which there is one unpaired magnon, their energy decreases with $n$
despite the pair-repulsion, indicating an attraction of the single magnon to
the pairs and a reduction of the single-magnon gap.

In Fig. S4(b), we plot the single-magnon excitation gap, $\Delta\bar{E}(n)$,
defined as $\Delta\bar{E}(n)\\!=\\!(2\bar{E}(n)-\bar{E}(n-1)-\bar{E}(n+1))/2$
for $n$=odd, which also corresponds to the average vertical distance between a
blue cross and two closest red circles in Fig. S4(a). Indeed, for $n\\!=\\!1$,
the value $\Delta\bar{E}(1)$ is nothing but the magnon pairing gap $\Delta$,
see Fig. 2(a) of the main text, which separates the lowest magnon energy from
that of the pair. Generalizing it to $n\\!>\\!1$, corresponds to the energy of
an extra magnon with respect to the background energy of magnon pairs. As one
can see in Fig. S4(b), this single-magnon excitation gap reduces significantly
upon increase of $n$, implying that the magnon is attracted to the magnon
pairs already at the level of one magnon–one pair, $n\\!=\\!3$, and lowers its
energy further by interacting with multiple magnon pairs for larger $n$. The
decrease of the single-magnon gap due to interaction with the condensate of
magnon pairs results in an expansion of the AFM phase above the non-
interacting phase boundary $h_{c1}$.

##  V finite-size scaling of the nematic order parameter

To investigate the order-parameter evolution with the system size, we follow
the strategy of Ref. [10]. It suggests that the DMRG cylinders with the open-
periodic boundary conditions and the aspect ratio $L_{x}/L_{y}$ about 2.0
should cancel the leading finite-size effects for the order parameter in the
center of the cluster, closely approximating its value in the thermodynamic 2D
limit.

We use three cylinders with the aspect ratio of 2.0 and apply the nematic
pair-field of 0.1$\langle S^{-}_{i}S^{-}_{i+y}\rangle$ on both edges. The
results for the nematic order-parameter $\langle S^{-}_{i}S^{-}_{i+y}\rangle$
and for the $z$-axis magnetization $\langle S^{z}_{i}\rangle$ are shown in
Fig. S5(a) as the profiles along the width of the clusters. The
$1/L_{y}$-scaling of the nematic order parameter and magnetization in the
center of the cylinders are presented in Fig. S5(b). Both quantities show
little change, giving a strong support to the existence of the nematic phase
in 2D and providing a vindication to the strategy of Ref. [10] used in this
work, suggesting that the width-8 cylinders that we rely on in our study
approximate the 2D behavior very well.

Figure S5: Finite-size scaling of the nematic order parameter and $z$-axis
magnetization in the nematic phase for $h=0.96$ and $J_{2}=0.7$. (a) Nematic
order parameter $\langle S^{-}_{i}S^{-}_{i+y}\rangle$ (left axis) and $z$-axis
magnetization $\langle S^{z}_{i}\rangle$ (right axis) along the cylinder for
the three cluster sizes of the same aspect ratio 2.0. There is a weak
symmetry-breaking pairfield $0.1S^{-}_{i}S^{-}_{i+y}$ applied to both edges.
(b) $1/L_{y}$-scaling of $\langle S^{-}_{i}S^{-}_{i+y}\rangle$ and $\langle
S^{z}_{i}\rangle$ at the center of the cylinders in (a).

##  VI dipole and quadrupole order parameters in a Heisenberg antiferromagnet

In a multipole expansion, higher-order multipole correlations can exist merely
as a byproduct of the lower-order ones, so the nematic (quadrupolar) order
parameter should always be detectable in the presence of the conventional
dipolar order. In this section, we present a numerical demonstration of this
effect for the $J_{1}\\!>\\!0$, nearest-neighbor Heisenberg antiferromagnet in
a magnetic field. By highlighting its differences from the $J_{1}$–$J_{2}$
ferro-antiferromagnetic model that we study in the present work, we further
confirm that the existence of the pure nematic phase in the $J_{1}$–$J_{2}$
model is not a byproduct of the dipolar order.

Figure S6: (a) Spin configuration in a long-cylinder scan of the Heisenberg
$J_{1}$-AFM, with the magnetic field varying linearly across the length of the
cylinder from zero to above the saturation field. (b) Magnetization $\langle
S^{z}_{i}\rangle$ (left axis), square of the dipolar order parameter, $\langle
S^{-}_{i}\rangle^{2}$, and the quadrupolar order parameter $\langle
S^{-}_{i}S^{-}_{i+y}\rangle$, (right axis) along the cylinder in (a). (c) The
same order parameters as in (b) for the $J_{1}$–$J_{2}$ FM-AFM model for
$J_{2}=0.7$, as in Fig. 4(b) of the main text, but for a wider field range.
The shaded regions are the phases from the naive phase diagram in Fig. 1(a) of
the main text based on the values of $h_{c1}$ and $h_{c2}$. The actual AFM to
nematic phase boundary is at $h\approx 0.95$, which is way above $h_{c1}$ and
only slightly below $h_{c2}$.

In Figs. S6(a) and S6(b) we show the DMRG results for the long-cylinder scan
of the Heisenberg $J_{1}$-AFM in a magnetic field, with $J_{1}=1$ set as the
energy unit. The magnetization curve shown in Fig. S6(b) agrees very closely
with the prior spin-wave calculation [11]. One can see that the spin canting,
demonstrated by the square of the order parameter, $\langle
S^{-}_{i}\rangle^{2}$, occurs below saturation field with no sign of an
intermediate nematic phase. There is also a finite quadrupolar order
parameter, $\langle S^{-}_{i}S^{-}_{i+x(y)}\rangle$, which is detected
alongside the dipolar order, but not independently. Furthermore, the dipolar
order is always stronger than the quadrupolar one throughout our scan, with no
region where the former can dominate and indicate a possible occurrence of the
pure nematic state.

This should be contrasted with the scan of the $J_{1}$–$J_{2}$ ferro-
antiferromagnetic model for $J_{2}=0.7$, as in Fig. 4(b) of the main text, but
for a wider field range, shown in Fig. S6(c). One can observe a region near
the saturation field where the ordered moment has deviated from the full
saturation, but shows no spin canting, also indicated by the dominant
quadrupolar order, $\langle S^{-}_{i}S^{-}_{i+y}\rangle$, with the suppressed
or vanishing dipolar one, $\langle S^{-}_{i}\rangle^{2}$, in this region.
These are the distinct features of the spin nematic phase as we argue and
analyze in more detail in the main text.

##  VII first order transition from FM to AFM at $J_{2}=0.45$

Figure S7: DMRG non-scan results for $J_{2}=0.45$, see Fig. 5 in the main
text. (a) The ground state at $h=0.140$ is a strongly canted AFM state as is
indicated by the magnetization $\langle S^{z}_{i}\rangle$ (left axis), and by
the dipolar, $\langle S^{-}_{i}\rangle^{2}$, and quadrupolar, $\langle
S^{-}_{i}S^{-}_{i+y}\rangle$, order parameters (right axis). (b) The $h=0.143$
state is a fully saturated FM state. We apply small pairing and canting fields
on the left and right edges, as is described in the main text.

In this section, we provide further evidence supporting the first-order nature
of the FM to AFM transition at $J_{2}$=0.45, based on the calculations in the
clusters with fixed parameters (i.e., non-scan calculations). As is shown in
Fig. S7, for the two very close magnetic fields, $h=0.140$ and $h=0.143$, the
ground states are very different. For $h=0.140$, it is a strongly canted AFM
state, while for $h=0.143$ it is a fully saturated FM. This discontinuity in
the ground state spin configuration is a clear indication of the first-order
transition, consistent with the scan calculations, shown in Fig. 5 of the main
text.

Another general feature of the first-order transition is a possibility of the
metastable states, in which the two states that are separated by such a
transition can appear stable beyond the regions where they constitute global
energy minima. We can observe such a metastability in DMRG calculations. As
was mentioned in Sec. III, in the attractive-pair regime, the ground state is
either a saturated FM state or a strongly canted AFM state with a large magnon
density. An implication of this two drastically different outcomes is that
these states remain local minima in terms of magnon density, even if they are
not true ground states. For example, in the $h=0.143$ case (polarized FM
ground state), we can initiate DMRG calculation in the collinear Néel state,
and the cluster will be stuck in the metastable canted AFM state, showing
correlations that are similar to the one in Fig. S7(a), but having a higher
energy. On the other hand, for the repulsive-pair regime ($J_{2}>0.6$), where
all phase transitions are expected to be continuous, we do not detect such a
metastability in our DMRG calculations, in agreement with these expectations.

## References

  * Tinkham [2003] M. Tinkham,  _Group Theory and Quantum Mechanics (Dover Books on Chemistry)_ (Dover Publications, 2003).
  * Zhitomirsky and Tsunetsugu [2010] M. E. Zhitomirsky and H. Tsunetsugu, Magnon pairing in quantum spin nematic, EPL (Europhysics Letters) 92, 37001 (2010).
  * Bethe [1931] H. Bethe, Zur Theorie der Metalle, Zeitschrift für Physik 71, 205 (1931).
  * Feynman [1998] R. P. Feynman,  _Statistical Mechanics: A Set Of Lectures_ (CRC Press, Boca Raton, 1998) pp. 198–220.
  * Dyson [1956] F. J. Dyson, General theory of spin-wave interactions, Phys. Rev. 102, 1217 (1956).
  * Wortis [1963] M. Wortis, Bound States of Two Spin Waves in The Heisenberg Ferromagnet, Phys. Rev. 132, 85 (1963).
  * Mattis [2006] D. C. Mattis,  _The theory of magnetism made simple: an introduction to physical concepts and to some useful mathematical methods_ (WorId Scientific Publishing, Singapore, 2006).
  * Galitski  _et al._ [2013] V. M. Galitski, B. Karnakov, V. Galitski, and V. I. Kogan,  _Exploring Quantum Mechanics: A Collection of 700+ Solved Problems for Students, Lecturers, and Researchers_ (OUP Oxford, 2013).
  * Annett [2006] J. F. Annett,  _Superconductivity, Superfluids and Condensates_ (Oxford University Press Inc., New York, 2006).
  * White and Chernyshev [2007] S. R. White and A. L. Chernyshev, Neél Order in Square and Triangular Lattice Heisenberg Models, Phys. Rev. Lett. 99, 127004 (2007).
  * Zhitomirsky and Nikuni [1998] M. E. Zhitomirsky and T. Nikuni, Magnetization curve of a square-lattice heisenberg antiferromagnet, Phys. Rev. B 57, 5013 (1998).