††thanks: Joint first author††thanks: Joint first author

# Ferrofluid drop impacts and Rosensweig peak formation in a non-uniform
magnetic field

Amelia Cordwell  Department of Physics, The University of Auckland, Auckland
1142, New Zealand  Alex Chapple<EMAIL_ADDRESS>Department of
Physics, The University of Auckland, Auckland 1142, New Zealand  Stephen Chung
Department of Physics, The University of Auckland, Auckland 1142, New Zealand
Frederick Steven Wells  Department of Physics, The University of Auckland,
Auckland 1142, New Zealand  The MacDiarmid Institute for Advanced Materials
and Nanotechnology, New Zealand  Geoff R. Willmott<EMAIL_ADDRESS>Department of Physics, The University of Auckland, Auckland 1142, New Zealand
The MacDiarmid Institute for Advanced Materials and Nanotechnology, New
Zealand  School of Chemical Sciences, The University of Auckland, Auckland
1142, New Zealand

###### Abstract

Vertical drop impacts of ferrofluids onto glass slides in a non-uniform
magnetic field have been studied using high-speed photography. Outcomes have
been classified based on the motion of the fluid-surface contact lines, and
formation of peaks (Rosensweig instabilities) which affect the height of the
spreading drop. The largest peaks are nucleated at the edge of a spreading
drop, similarly to crown-rim instabilities in drop impacts with conventional
fluids, and remain there for an extended time. Impact Weber numbers ranged
from 18.0 to 489, and the vertical component of the $B$-field was varied
between 0 and 0.37 T at the surface by changing the vertical position of a
simple disc magnet placed below the surface. The falling drop was aligned with
the vertical cylindrical axis of the 25 mm diameter magnet, and the impacts
produced Rosensweig instabilities without splashing. At high magnetic field
strengths a stationary ring of ferrofluid forms approximately above the outer
edge of the magnet.

##  I  Introduction

Ferrofluids, which consist of ferromagnetic colloids in solution [1, 2], have
found industrial applications such as loudspeaker cooling, sealing rotating
shafts, and high-speed printing [3]. They have been used in areas as diverse
as sculpture [4], self assembly [5], medicine [3], and even digital logic [6].
Ferrofluids are interesting because various complex phenomena are produced
when they interact with applied magnetic fields [2]. Most famously,
ferrofluids form Rosensweig instabilities: a series of ordered peaks
spontaneously forms on the surface of a ferrofluid in the presence of a
magnetic field normal to the surface [1]. After peaks form they slowly arrange
into regular hexagonal or square grids, depending on the fluid’s magnetisation
history [7]. The creation of such patterns is only possible when the fluid
reaches a magnetisation above a critical value as dictated by the fluid
parameters [8]. The heights of Rosensweig peaks increase with increased
magnetic field strength [9], although in most cases the parameters of the
analytic theory describing this growth have been found empirically [7].

The impact of ferrofluids onto solid surfaces has been the subject of only a
few studies for the cases of no or uniform magnetic field [10, 11, 12, 13, 14,
15], even though the general field of liquid droplet impacts is of strong
fundamental and applied interest [16, 17, 18]. It has been found that the
spreading of an impacting ferrofluid drop is delayed by a vertical magnetic
field, but extends further than with no field applied [10]. Horizontal
magnetic fields have been found to stabilise a ferrofluid surface, delaying
the onset of instabilities [15]. Angled magnetic fields cause the fluid to
form ridges rather than hexagonally arranged peaks at high enough field
strengths [11, 12]. A study reaching relatively high magnetic field strengths
(up to $\sim$ 0.1T) showed that the normalized spreading droplet height
monotonically decreased with increased magnetic field strength for nanolitre
droplet volumes [10]. Another study [13] at similar field strengths confirmed
these findings, and demonstrated that the spreading oscillation of the droplet
in the horizontal direction can be reduced substantially with a vertical
magnetic field. The maximum spread of impacting ferrofluid droplets has been
studied in fields $\lesssim 0.3$ T [10, 14].

It is of interest to study ferrofluid droplet impacts at high enough magnetic
field strengths ($\gtrsim 0.15$ T) and with large enough droplet volumes so
that Rosensweig peaks play an important role in the spreading dynamics.
Moreover it is desirable to understand the dynamics of ferrofluid droplets in
non-uniform magnetic fields, for both fundamental and practical reasons. For
example, on-demand ferrofluid droplet formations can be created using a time-
dependent non-uniform magnetic field [19]. The field generated by a simple bar
magnet is also non-uniform. The dynamics of stationary ferrofluid droplets in
the presence of a non-uniform magnetic field have been studied [20], but the
influence of non-uniform magnetic fields on drop impacts and consequent
Rosensweig instability formation is still to be described. A ferrofluid
droplet falling prior to impact is elongated in the direction of any applied
magnetic field, taking an ellipsoidal shape [21]. In a non-uniform field this
elongation can produce asymmetries, and an extended tip forms for drops
falling vertically along the cylindrical axis of a disc magnet [22].

In this work we perform experiments and analysis to study vertical ferrofluid
drop impacts in the non-uniform magnetic field generated by a permanent disc
magnet placed below the impact surface. Our work extends to comparatively
higher magnetic field strengths and impact velocities than previous ferrofluid
drop impact studies, and we omit drops which splash or break up. We report on
the dynamic formation of Rosensweig instabilities following impact. Drop
impacts on solid surfaces for any fluid produce a spreading thin lamella, with
a raised dynamic rim at the edge of the spreading droplet [23]. So-called
‘crown’ instabilities on this rim can lead to corona splashes, where secondary
droplets are created [16, 24]. For ferrofluids, we are particularly interested
in the interplay between the formation of crown instabilities and Rosensweig
instabilities, as well as the patterning of these instabilities.

We have also found that at relatively high field strengths, ferrofluid drop
impacts can result in stationary rims of ferrofluid forming near the edge of
the magnet. Rim formations have been observed previously [23], although they
were not stationary, likely due to weaker magnetic fields.

##  II Materials and Methods

Figure 1: Schematic diagram of the experiment. A is the variable height needle
stage, B the glass slide, C the magnet on a movable stage, D the side view
camera, E the top view camera with a mirror providing an image at an oblique
angle, and F and G are LED lighting sources.

Ferrofluid droplets were dropped onto a clean glass slide placed above a
permanent neodymium disc magnet, and the resulting impacts were studied using
high-speed photography (Fig. 1). The ferrofluid (EFH1, Ferrotec) was composed
of 3 vol% Fe3O4 particles in light hydrocarbon. The material properties (based
on manufacturer’s specifications [25]) include a saturation magnetisation of
$M_{s}=44$ mT, density $\rho=1.21\times 10^{3}$ kg/m3, initial magnetic
susceptibility $\chi=2.64$, surface tension $\sigma=0.029$ N/m, and viscosity
$\eta=6$ $\mathrm{mPa}\text{\,}\mathrm{s}$.

The magnet was an axially magnetised cylindrical disc (NdFeB Grade-N38, AMF
Magnetics Australia) with diameter of $d=25$ mm and a length of $l=12.7$ mm.
It was affixed to a movable platform so that the distance between the top of
the magnet and the top of the glass slide ($h_{0}$) could be varied. The
vertical component of the magnetic flux density ($B_{z}$) was measured using a
teslameter (Group 3 DTM-133) at radial positions of $R=0,6.2$ and $12.5$ mm
for cylindrical polar co-ordinates $(R,\phi,z)$ with their origin at the
centre of the upper surface of the magnet. The resulting plot of $B_{z}$ as a
function of $h_{0}$ (equivalent to $z$) is shown in Fig. 2. A numerical model
[26] for a uniformly magnetised cylinder of the same size agreed with the on-
axis measurements for a magnetisation of $M=1.22/\mu_{0}$ T (here $\mu_{0}$ is
the permeability of free space). The accuracy of the model was reduced off-
axis, likely because the magnet is not exactly uniformly magnetised. Repeat
measurements carried out using the same magnet over a long time span ($\sim$3
years) have revealed no significant change in the on-axis field strength. It
was found that the critical value of $h_{0}$ for Rosensweig peak formation for
a ferrofluid drop on the glass slide was $\approx 40$ mm, corresponding to an
on-axis vertical flux density of around 12 mT.

Figure 2: Measured (points) and modelled (line) on-axis vertical magnetic flux
density for the 25 mm diameter disc magnet used in this study, plotted as a
function of the magnet-surface separation $h_{0}$. The model line represent a
uniformly magnetised cylinder with on-axis magnetisation of $M=1.22/\mu_{0}$
T.

Individual droplets were generated by a syringe pump and dropped from a
polyether ether ketone (PEEK) needle mounted on a movable stage. The release
point was at a variable height $L_{0}$ above the glass slide surface. The size
of the droplets could vary depending on the precise positioning of the needle
and on the magnet position, as the droplet release could be affected by the
$B$-field. Furthermore, ferrofluid droplets elongate in the direction of the
$B$-field resulting in a dynamically varying shape as the droplet falls [22].
Initial tests were performed to align the centre of the magnet with the centre
of the droplet impact. As the paramagnetic ferrofluid is attracted to areas of
higher magnetic field, the centre of a misaligned drop would shift towards the
center of the magnet approximately a minute after impact, leaving behind some
residue. The centre of the magnet could be shifted appropriately. In impact
videos, significantly misaligned droplets could be identified due to the non-
vertical drop elongation. These were excluded from analysis of height and line
dynamics (Fig. 6) but some asymmetric impacts are included and discussed in
the study of rim formation (Fig. 9).

Droplet impacts were recorded for magnet distances ranging from $h_{0}=2$ mm
to 28 mm, with drop heights ranging from $L_{0}=260$ mm to 820 mm. The glass
slide was cleaned with acetone and wiped with a paper towel in between each
impact. Two high speed cameras (a Photron Fastcam SA5 and a Photron Fastcam
Mini AX5) operating at 2000 fps (frames per second) were used to record the
impacts from side-on and top-down views respectively. The top down view was
reflected from a mirror at $\sim 45^{\circ}$ to the vertical, although a more
oblique angle was sometimes used (Fig. 1). Top-down views were not used to
make quantitative measurements, and were analysed manually where required. For
the side-view camera, impacts were back-lit using a customised LED lighting
unit, and images were spatially calibrated (typically at $\sim$25 $\mu$m per
pixel) using a fiducial reference square. 100 frames (0.05 s) were recorded
for each impact, long enough to show the initial formation of Rosensweig
peaks. In some experiments, recording continued for a longer time period to
capture the droplet settling further.

##  III Analysis

###  III.1 Image Analysis and Dimensionless Numbers

Drop impacts are most often characterised using the Weber number,

$We_{0}=\frac{\rho v^{2}D_{0}}{\sigma},$ (1)

due to the importance of competition between inertia and surface tension in
the drop outcome. Here $D_{0}$ is the diameter of the droplet, assumed to be
spherical, and $v$ is the impact velocity. However, ferrofluid drops in a
field become non-spherical, so in this study a first-principles definition of
the Weber number is used,

$We=12\frac{E_{\text{kin}}}{E_{\text{surf}}}=\frac{6mv_{c}^{2}}{\sigma
A_{s}},$ (2)

where $E_{\text{kin}}=\frac{1}{2}mv_{c}^{2}$ and $E_{\text{surf}}=\sigma
A_{s}$ are respectively the kinetic and surface energies of the drop, and
where $m$ and $A_{s}$ are respectively the drop’s mass and surface area.
Furthermore, different parts of the droplet have different velocities at
impact due to drop deformation in the non-uniform magnetic field [22].
Therefore, the center-of-mass velocity $v_{c}$ is used as the velocity scale.

Measurements of the drop’s initial geometry and velocity were made using the
frame immediately prior to impact captured by the side-view camera. Codes were
developed in Python to calculate each droplet’s volume ($V$) and surface area
($A_{s}$) by treating each row of pixels as a flat truncated cone one pixel in
height, and with the bottom and top radii calculated using the respective
widths of the imaged droplet. The measured volume was used to calculate an
equivalent diameter $D_{0}$ for each droplet, defined as the diameter of a
spherical droplet with same volume i.e. $D_{0}=2(3V$/$4\pi)^{1/3}$. Analysis
of the experimental images revealed that $v_{c}$ varied between $0.50$ and
$3.88$ m/s. The range of equivalent diameters was $0.45<D_{0}<11.6$ mm, and
the experiments covered the range $18.0<We<489$.

The value of $We$ calculated using Eq. 2 is the dimensionless number primarily
used to describe the experiments here. A comparison with other dimensionless
numbers is provided in the Supplemental Material [27], using the $D_{0}$ and
$v_{c}$ as length and velocity scales respectively. These include $We_{0}$
(Eq. 1), as well as the Reynolds ($Re=\rho v_{c}D_{0}$/$\eta$), magnetic Bond
($Bo_{m}=B_{z}^{2}D_{0}$/$\mu_{0}\sigma$) and Hartmann ($Ha=\rho
B\mathcal{M}D_{0}$/$2v_{c}\eta$, where $\mathcal{M}$ is the magnetic moment)
numbers. These dimensionless numbers are not used in the main text because the
use of $D_{0}$ (unlike the calculation of $We$) does not account for droplet
elongation.

To study the drop spreading dynamics, the contact line width ($s$, see
Supplemental Material [27]) and the maximum height of the droplet above the
surface ($H$) were also measured as a function of time after impact ($t$).
Drop impact studies often use a dimensionless spread factor, defined as the
linear extent of the spreading drop divided by the drop’s diameter prior to
impact. Here the droplets can be highly non-spherical prior to impact, so an
alternative spread factor ($\beta$) is used, in which the contact line width
is divided by the equivalent diameter,

$\beta=\frac{s}{D_{0}}.$ (3)

The height is similarly non-dimensionalized using $\delta=H$ / $D_{0}$.

Measurements taken at the point of maximum spreading can be used for
comparisons of drop impacts. However, it was found that for a ferrofluid drop
in a magnetic field the contact line width can reach an initial maximum, then
steadily grow or else decrease before growing again due to magnetic effects
(see below). Therefore we identify maximum spread with the first maximum of
the contact width after impact.

###  III.2 The Radial Field

A theoretical magnetic body force acting on the fluid in a radial direction
can be calculated. This force is caused by the radial change in the magnetic
field magnitude. For this calculation it is assumed that the fluid is at
saturation magnetisation, that the force on any volume of fluid can be
considered independently, and that the magnet is a cylinder of axially uniform
magnetisation. Under these conditions, the force on a magnetised material in
an external field,

$\displaystyle\mathbf{F}=\int_{V}\nabla(\mathbf{M}\cdot\mathbf{B})dV,$ (4)

produces a magnetic body force on a unit volume of ferrofluid at saturation
magnetisation, $\mathbf{M}=M_{s}\mathbf{B}/|B|$,

$\displaystyle\mathbf{f}=M_{s}\nabla(\mathbf{B}/|B|\cdot\mathbf{B}),$ (5)

with the radial component

$\displaystyle f_{R}=M_{s}\frac{\partial}{\partial
R}(\mathbf{B}\cdot\mathbf{B})^{1/2}.$ (6)

Figure 3 shows the radial force per unit volume of ferrofluid as a function of
$R$ at the impact surface (a distance $h_{0}$ above the magnet), computed
using Eq. 6 and the $B$-field generated from the fit to the data in Fig. 2.
When the surface is more than 7.5 mm above the magnet, there is no radial
outwards force ($f_{R}<0$ for all values of $R$). When $h_{0}<7.5$, an outward
radial force acts at radial positions from the centre of the magnet ($R=0$)
out to a point where the sign of $f_{R}$ changes. As $h_{0}$ decreases the
radial position at which $f_{R}=0$ increases, and elsewhere the magnitude of
$f_{R}$ increases. The signs of the force indicate stable equilibrium at
$f_{R}=0$, and when $h_{0}$ is small the restoring forces intensify. This
calculation suggests that ferrofluid could be pushed towards $f_{R}=0$,
potentially forming a circular rim which would become thinner and further from
the centre of the magnet with decreasing $h_{0}$.

Figure 3: Radial magnetic force per unit volume of ferrofluid at the impact
surface for an axially magnetised uniform cylinder as described in the text.
The blue line corresponds to $f_{R}=0$.

##  IV Results and Discussion

###  IV.1 Typical Results

Figure 4: Droplet impact sequences in vertical strips, labelled with the time
after impact in milliseconds. Scale bars are 1 mm. Droplet A: We = 289,
$V=7.34$ $\mu$l, $v=1.61$ m/s, $h_{0}=62.0$ mm, $B_{z}=0.007$ T. Droplet B: We
= 271, $V=6.25$ $\mu$l, $v=1.72$ m/s, $h_{0}=28.9$ mm, $B_{z}=0.025$ T.
Droplet C: We = 249, $V=6.36$ $\mu$l, $v=1.82$ m/s, $h_{0}=12.85$ mm,
$B_{z}=0.112$ T. Droplet D: We = 269, $V=4.82$ $\mu$l, $v=2.14$ m/s,
$h_{0}=5.0$ mm, $B_{z}=0.272$ T. Figure 5: For the four droplets (A-D) from
Fig. 4, plots of (a) normalized height of the droplet shape, and (b) spread
factor as defined in Eq. 3. Insets show an expanded view of the first few ms
for each figure.

Figure 4 shows four representative examples of ferrofluid droplet impacts from
this study, with their height and contact line dynamics recorded in Fig. 5.
Full videos of impacts A-D are contained in the Supplemental Material [27].
The magnetic field increases monotonically from A to D, affecting the
elongated shape at impact [22], and the centre-of-mass velocity at impact
progressively increases. The drop size decreases with $h_{0}$ because the
field also affects drop release. Droplet A is not under significant influence
of the magnetic field and simply spreads out before retracting slightly. This
deposition outcome is similar to observations of ferrofluid drop impacts onto
hydrophilic surfaces with no applied field reported previously [15]. Droplets
B-D have different dynamics. For droplet B the initial impact and spreading
produce a height and spread factor similar to droplet A (Fig. 5), despite the
elongation of the drop at impact. In general, the maximum spread factor
observed in this study at relatively low field strengths ($B_{z}\lesssim 100$
mT) follows a trend close to $\beta_{max}\propto We^{0.25}$ (see Supplemental
Material [27]), as is often observed for non-magnetic drop impacts onto solid
surfaces [28]. Unlike droplet A however, instabilities appear on the outer
edge of the rim of droplet B as it approaches maximum spread ($\approx 6$ ms),
similar in appearance to crown-rim instabilities observed for conventional
drop impacts [16, 24]. These instabilities do not generate secondary drops or
large fingers, but as the droplet retracts, they do not disappear as they
would for a conventional drop impact. Between 15 and 50 ms, the growth of
Rosensweig instabilities becomes apparent. Qualitatively, we observe that
there appears to be a correlation between the crown-like instabilities and the
Rosensweig instabilities. This suggests a mechanism in which the Rosensweig
peaks nucleate from the rim instabilities, which is reasonable given that
Rosensweig peaks are caused by a relative increase in the magnetic force
acting on small peaks on an otherwise uniform film. Once nucleated, peaks
undergo unstable growth until opposed by sufficient surface tension [1].

For droplet C, fluctuations in the rim when the drop is near maximum spread
($\approx 6$ ms) are directly observed to evolve into Rosensweig
instabilities. Near maximum spread, the peaks emerging from the expanding
lamella are more pronounced than for either droplet B or D. Over time, some of
the instabilities shift and merge into one another (see videos in the
Supplemental Material [27] and Fig. 7, below), and consequently the maximum
height appears to fluctuate. These movements reflect that the spatial
arrangement of crown-rim instabilities on a dynamically moving lamella does
not correspond to a stable (or metastable) arrangement of Rosensweig peaks.
The number of peaks and the droplet height remain constant after $\approx 40$
ms, but the peaks continue to spread apart from each other. In the contact
line dynamics (Fig. 5(b)), after maximum spread the drop initially retracts
slightly as for droplets A and B. Later droplet C spreads out, which is
discussed in Section IV.2. Droplet C produces shapes with the greatest minimum
height and lowest maximum spreading of the four droplets (Fig. 5). The height
is produced by the clear and early emergence of peaks, while it is generally
observed that the maximum spread reduces with $B_{z}$ at relatively low fields
(see Supplemental Material [27]).

For droplet D, Rosensweig instabilities form earlier than for droplet C, and
the contact line moves continuously outwards without retracting. The peak
spacing decreases monotonically between droplets B, C and D. The maximum
Rosensweig peak width in vertically non-uniform fields can be calculated as
[5, 20],

$\lambda_{c}=2\pi\sqrt{\frac{\sigma}{\rho g+M_{s}\frac{d}{dz}(B_{z})}},$ (7)

where $g$ is gravitational acceleration, and it is assumed that the magnetic
flux density is only dependent on $z$ and that the fluid is at its saturation
magnetisation. It is therefore reasonable to attribute the decreasing peak
spacing for droplets B-D to the increasing magnetic field gradient.

We note two questions that are apparent when comparing droplets A-D that are
left as suitable topics for future study. Firstly it is observed that rim
instabilities evolve into Rosensweig peaks, but it is unclear whether the
origin of the rim instabilities is always magnetic. A comparison between
droplets A and B suggests this is the case, but it is possible that the shape
and velocity of the impact for droplet B would generate rim instabilities in
the absence of a $B$-field. Here we note that the time scale for Rosensweig
instability growth decreases with increasing magnetic field, so that the
effects of magnetism are clear at maximum spread for droplets C and D.
Secondly the increasing aspect ratio of the droplets before impact may have an
important influence on the dynamics, as the upper halves of these droplets
experience conditions similar to that of a conventional droplet impacting on a
thin film. A recent study has approached the topic of elongated drop impacts,
albeit for regular ellipsoids with aspect ratios 2:1 or smaller [29]. That
study found that the scaling of maximum spreading diameter with We0.25, as
developed for spherical droplets, is preserved. Here, we do not attempt to
deconvolute the separate influences of impacting droplet shape and the
magnetic field on the spreading lamella.

###  IV.2 Impact Outcome Classification

Droplets outcomes have been classified for a large number of experiments using
measurements such as those in Fig. 5. Experiments were focused at $We$ values
and release heights below those that resulted in the droplet breaking up or
splashing, and any such cases were identified and excluded. The height
dynamics were divided into two categories based on the trend observed
following the initial spreading (up to maximum spread, or $\approx$5 ms) of
the droplet. These are heights that monotonically decrease (Decay, e.g.
droplet A), or else increase (Grow, e.g. droplets B and C). The contact line
dynamics were classified into three categories, also based on the trend
observed following the initial spreading of the droplet. Decay occurs when
$\beta$ increases to a maximum spread before decreasing (e.g. droplets A and
B), Rebound occurs when $\beta$ initially decreases after reaching maximum
spread before eventually increasing above that value (e.g. droplet C), and
Rise occurs when $\beta$ monotonically increases (e.g. droplet D).

Plots of outcomes as a function of $B_{z}$ and $We$ shown in Fig. 6 give
generalised evidence for the trends identified in Section IV.1. The trend for
droplet height (Fig. 6(a)) is notable because it is non-monotonic with respect
to field strength. The Decay classification is recorded when there is no
magnetic field present, as expected for a non-magnetic fluid on a wetting
surface. When a relatively small magnetic field is present, emergence of
Rosensweig peaks leads to an increase in the height (Grow classification), as
expected from previous observations [9]. For $B\gtrsim 0.25$ T the Decay
classification is again recorded because the Rosensweig peaks are small and
form almost immediately. The classifications can be uncertain (especially near
boundaries) because the degree of height increase or decrease can be small.
For example for droplet D, the height remains approximately constant (Fig.
5(a); classified marginally as Grow), consistent with its position on the
boundary between Decay and Grow classifications in Fig. 6(a).

The contact line dynamics in Fig. 6(b) broadly indicate a transition from
Decay to Rebound to Rise as the magnetic field strength increases. At zero
field the drops retract slightly due to the wetting properties of the fluid
following impact. With a field applied, two apparently contradictory trends
are observed: the maximum spread following impact decreases with $B_{z}$ (see
Supplemental Material [27]), while later the spreading tends to increase with
$B_{z}$, producing Rebound and Rise classifications. The former trend (but not
the latter) is supported by the calculation of radial force components (Fig.
3), which suggests that $f_{R}$ always acts inwards for $B_{z}\lesssim 0.20$ T
($h_{0}\gtrsim 7.5$ mm). The outward droplet spreading can be instead
explained by the Rosensweig peaks, which repel each other, and can therefore
force the (wetting) contact line to gradually spread. This mechanism is
consistent with the apparent dependence of the classifications on $B_{z}$ (and
hence on the vertical field which pins the peaks) in Fig. 6(b) and with the
correspondence in timing and intensity between peak formation and contact line
growth in Fig. 5(b).

To summarize Fig. 6(b), at low non-zero fields ($B_{z}\lesssim 0.04$), the
field has insufficient strength to reverse the retraction observed at zero
field. A few relatively tall Rosensweig peaks form well after the initial
spreading of the droplet, at which point $\beta$ may increase relative to the
zero-field case (Fig. 5(b)). Over a wide range of higher field strengths
($0.04\lesssim B_{z}\lesssim 0.25$), the increased pinning of the peaks
produces a spreading force that prevails following an initial small
retraction, so that the contact line width increases steadily in the long time
limit. In this range, Decay outcomes are observed for some cases at relatively
high $We$, in which the spreading and retraction driven by inertia dominates
any spreading caused by peaks. Finally, for the highest fields studied
($B_{z}\gtrsim 0.25$) a strong outward force pulls the droplet outwards
continuously post-impact. Here, outwards motion may be assisted by the radial
force (Fig. 3). These Rise outcomes can lead to formation of an outer rim, a
phenomenon discussed in Section IV.4.

Overall, Fig. 6 indicates strong correspondence between the transition from
Decay to Grow in the height dynamics, and from Rebound to Rise in the contact
line dynamics ($B_{z}\approx 0.25$). The categories are not perfectly
delineated as a function of $B_{z}$ and $We$, highlighting the difficulty of
summarising these nuanced experiments using a single dimensionless quantity as
discussed above. Differences in droplet volume and in the precise alignment of
the droplet impact point with the magnet centre are likely explanations for
the indistinct classification boundaries.

Figure 6: For ferrofluid drop impacts, classifications of (a) height and (b)
contact line dynamics as discussed in the text. Red markers indicate the
positions of experiments A - D from Figure 4, as labelled. Data points in (a)
and (b) do not exactly correspond because the imaging quality precluded a
classification in some cases.

###  IV.3 Motion of Peaks

Figure 7: Oblique top-down view of two ferrofluid drop impacts. Experimental
parameters are (b) $We$ = 220, $V=5.05$ $\mu$l, $v=2.00$ m/s, $h_{0}=23$ mm,
$B_{z}=0.041$ T, and (c) $We$ = 322, $V=4.66$ $\mu$l, $v=2.27$ m/s,
$h_{0}=16.5$ mm, $B_{z}=0.076$ T.

Figure 7 uses oblique top-down views of three droplets to explore how the
arrangement of the Rosensweig peaks is affected by the dynamics of the drop
impact. The drops in the non-zero field (Figs. 7(b) and (c)) are both
classified as “Grow” for their height dynamics and “Rebound” for their contact
line dynamics.

From impact up to $\approx$10 ms later, the lamella rim for the droplets in a
field is more prominent that for the control experiment (Fig. 7(a)), and
instabilities start to become apparent. Peaks form on the spreading rim of the
lamella, and are more widely spaced at lower field strength (Fig. 7(b)). Peaks
inside the rim form later and are clearly smaller for two reasons. Firstly,
during the spreading phase the thick rim supporting the growth of the outer
peaks contains more fluid than the thinner internal lamella. Secondly, the
instabilities in the rim form first, during spreading, potentially nucleated
by crown-rim instabilities (see Fig. 4). These peaks then have longer to grow
(and perhaps merge) to create taller peaks than peaks growing from the
lamella. The volume of fluid in each peak is less than would be supported by a
thick pool of fluid [9], so it is likely that growth of the inner peaks is
restricted by a shortage of fluid to draw upon.

Images captured approximately 1 minute after impact confirmed that, although
there was some movement of the Rosensweig peaks, configurations such as those
observed after 90 ms in Fig. 7 were reasonably stable. The peaks repel each
other and apparently reach a stable state in which fluid is not exchanged
between them. The inner peaks do not undergo sufficient merging to grow in
size significantly. The observation of larger outer and smaller inner peaks is
the inverse of the trend observed for a static droplet on a dry plate, in
which case the central peak is significantly larger [30].

###  IV.4 Rim Formation

Figure 8: Side-on (left) and top-down (right) oblique views of a ferrofluid
drop impact on to a plastic Petri dish with parameters $h_{0}$ = 1.2 mm,
$V\approx 6.2~{}\mu$l, $L_{0}$ = 290 mm. Labels denote time in ms.

In Fig. 6(b), it was established that the contact line continues to spread
outwards after the initial $\approx$5 ms following impact if the magnetic
field is strong enough (i.e. the Rise classification). To extend the study of
this spreading, some separate experiments were carried out with the magnet
placed very close to the surface ($h_{0}$ = 1.2 mm). A plastic Petri dish was
used due to the extent of the spreading drop, and so data from these
experiments were not included with the data for drops onto glass in the other
figures. Other experimental parameters are as described in Section II.

Figure 8 shows a typical example of these impacts. Rosensweig peaks form on
the edge of the droplet almost immediately upon impact, and maintain their
shape as they move radially outwards. As the peaks spread, they each leave
behind a trail of fluid, and their height decreases. Other peaks form nearer
the point of impact, and move outwards more slowly.

Eventually ($\approx 200$ ms) the outward motion is halted and the fluid at
the edge of the drop forms a circular rim, approximately following the outline
of the edge of the cylindrical magnet below the surface. This outer rim is
stationary, without any Rosensweig peaks, and is not to be confused with the
rim observed at the edge of the lamella immediately after impact. Well after
the rim has stabilised, the slow-moving peaks formed near the impact point can
move towards and coalesce with the rim, if the magnetic field is strong enough
(see also Supplemental Material [27]).

A loss of peak height in spreading ferrofluid drops has previously been
attributed to a reduction in vertical magnetic field strength [11], but here
the loss in height is likely to be due to increased magnetic field gradient
and stabilisation of the rim due to the radial field component. Rim formation
in a ferrofluid droplet has been observed previously [23], but this occurred
in a uniform vertical magnetic field created by electromagnetic coils, and the
rim was not pinned as in the present experiments.

Figure 9: Rim formation as a function of $B_{z}$ and the normalized maximum
spreading diameter. The dotted grey line corresponds to $f_{R}=0$ from Fig. 3.
Above this line the calculated force acts radially outwards, and below the
force acts inwards. Inset: a drop impact (We = 453, $V=2.79\mu$l, $v=3.23$
m/s, $h_{0}=2.26$ mm, $B=0.358$ T) pictured at $t=37$ s, showing partial rim
formation.

To further probe rim formation and the effect of impact velocity, the set of
droplet impact experiments carried out using a glass substrate for $2<h_{0}<5$
mm ($B_{z}\gtrsim 0.25$ T) was studied further. These experiments are included
in Fig. 6. Figure 9 shows the occurrence of rim formation as a function of
$B_{z}$ and the maximum spreading diameter. The droplet is considered to have
formed a rim when there is a clear discrepancy between the height of the
droplet edge and centre, as observed using both the top and side-view videos.
Some cases are labelled “unclear” due to ongoing motion of the fluid at the
end of the captured video.

Our initial hypothesis was that the transition between inwards and outwards
radial force ($f_{R}=0$, Fig. 3) would cause rim formation. A ferrofluid
element in a spreading droplet would be pushed towards that calculated radius,
and reach dynamic stability there. In agreement with this explanation, the
width of the rim (when observed) was generally thicker for smaller values of
$B_{z}$ (see Supplementary Information). However, this hypothesis does not
explain the onset of rim formation in Fig. 9, which is primarily a function of
the magnetic field strength. It appears that above some threshold value
($B_{z}\gtrsim 0.31$ T), which does not depend on the extent of horizontal
spreading, the vertical field is sufficient to pin the rim of the lamella in
place as it slows during spreading. This pinning apparently overcomes the
surface tension which acts to smooth the droplet and dissipate the rim.

It is prudent to discuss some experimental uncertainties which could affect
the conclusions here and may warrant further investigation. Firstly, the
arbitrary limit on the recording time means that the results may not represent
the equilibrium state of the rim. Secondly, many of these experiments produced
less symmetric spreading than shown in Fig. 8, as the symmetry of the
spreading was very sensitive to the position of the drop impact. When
asymmetric impact occurred, a rim could form to one side of the impact point
and spread along the shape of a circular rim (Fig. 9, inset).

Figure 10: A rim formed following drop impact of a larger amount of ferrofluid
than in Fig. 8.

In a final experiment, a drop impact was carried out using a larger volume of
ferrofluid, with the result shown in Fig. 10. A rim was clearly formed, and as
additional fluid was added to the stationary drop, the rim grew more in the
inwards direction. The radial ridges that can be observed on the inner edge of
the rim are likely to be related to those reported previously [11, 12], in
which lines rather than Rosensweig peaks appeared on the surface of a fluid
when the magnetic field was tilted by 23° to the vertical. Other workers have
observed a floating ferrofluid layer form a circular rim within a radially
inhomogeneous field [31]. These observations suggest that rim formation can be
explained by the static magnetic field, but that a dynamic event such as a
droplet impact is required to observe phenomena such as radially moving
Rosensweig peaks, and as discussed earlier can affect the final droplet shape.

##  V Conclusion

Ferrofluid drop impacts have been investigated within a non-uniform field of
the kind readily produced by a simple permanent magnet. The work has focused
on spreading droplets which form Rosensweig instabilities, and which do not
splash or break up. In general, the outcomes were more dependent on the
magnitude of the applied field than on the size, elongation, or velocity of
the impacting drops as characterised by the Weber number. The non-uniform
field is most simply characterized by the vertical component of the flux
density at the impact surface, $B_{z}$.

The motions of the fluid-surface contact line and of the spreading droplet
height were measured during the initial fast expansion to ‘maximum’ spread
($\approx$5 ms after impact), and up to around 100 ms. At $B_{z}\approx 0$,
the drops retracted slightly after reaching maximum spread, and their height
monotonically decreased. For $0<B_{z}\lesssim$ 0.25 T peaks formed so that,
after the initial fast spreading stage, the height of the spreading drop grew.
The instabilities repelled each other so that the extent of the contact line
also gradually increased. At $B_{z}\gtrsim$ 0.25 T, the greater field produced
more, smaller peaks so that there was monotonic height reduction and outwards
spreading.

The first and largest Rosensweig peaks were formed on the rim at the droplet
edge, where there is a greater volume of fluid available than in the lamella.
Although some motion and coalescence of peaks was observed over time, the
largest peaks remained at the outer edge of the drop, with smaller peaks
closer to the centre, for at least $\approx$1 s. This patterning is the
inverse of what is observed for static pools of ferrofluid. Perturbations
which evolved into the largest, outer peaks were apparent in the rim by the
time of maximum spread ($\approx$5 ms). For relatively high fields, it was
clear that these were instigated by the $B$-field. At lower fields, magnetic
perturbations are slower to form and it is possible that Rosensweig peaks are
primarily nucleated by crown-rim instabilities.

For higher fields ($B_{z}\gtrsim 0.31$ T) the peaks at the outer edge of the
spreading drop eventually coalesced to form a prominent, stationary rim near
the edge of the magnet below. The transition to rim formation is primarily a
function of the magnetic field strength, rather than being related to drop
spreading. It is suggested that the rim is pinned by the magnetic field in a
similar manner to previous observations of static ferrofluid drops.

The study of drop impacts generally has a broad parameter space, including
(for example) the size, release height, and fluid properties of the drop. For
ferrofluids this space has not been well studied, so that various interesting
extensions to the current work are possible. For example, impacts which
produce splashes can be studied, and the shape of the non-uniform $B$-field
could be varied along with its magnitude. Impacts can also be studied for a
wide range of surface properties, including different wettabilities, as well
as random or designed surface structures (roughness) at any length scale.
Ultimately there is the prospect of using drop impact to explore and control
distributions of Rosensweig peaks, thereby producing interesting patterns that
can be of practical interest.

###### Acknowledgements.

The authors wish to thank Ekaterina Zossimova and Santhosh Kumar Pandian for
minor contributions to this work.

## References

  * Rosensweig [1985] R. E. Rosensweig,  _Ferrohydrodynamics_ (Cambridge University Press, Cambridge, 1985).
  * Deb  _et al._ [2022] R. Deb, B. Sarma, and A. Dalal, Magnetowetting dynamics of sessile ferrofluid droplets: a review, Soft Matter 18, 2287 (2022).
  * Kole and Khandekar [2021] M. Kole and S. Khandekar, Engineering applications of ferrofluids: A review, Journal of Magnetism and Magnetic Materials 537, 168222 (2021).
  * Kodama [2008] S. Kodama, Dynamic ferrofluid sculpture: organic shape-changing art forms, Communications of the Association for Computing Machinery 51, 79 (2008).
  * Timonen  _et al._ [2013] J. V. Timonen, M. Latikka, L. Leibler, R. H. Ras, and O. Ikkala, Switchable static and dynamic self-assembly of magnetic droplets on superhydrophobic surfaces, Science 341, 253 (2013).
  * Katsikis  _et al._ [2015] G. Katsikis, J. S. Cybulski, and M. Prakash, Synchronous universal droplet logic and control, Nature Physics 11, 588 (2015).
  * Gollwitzer  _et al._ [2006] C. Gollwitzer, I. Rehberg, and R. Richter, Via hexagons to squares in ferrofluids: experiments on hysteretic surface transformations under variation of the normal magnetic field, Journal of Physics: Condensed Matter 18, S2643 (2006).
  * Cowley and Rosensweig [1967] M. Cowley and R. E. Rosensweig, The interfacial stability of a ferromagnetic fluid, Journal of Fluid mechanics 30, 671 (1967).
  * Friedrichs and Engel [2001] R. Friedrichs and A. Engel, Pattern and wave number selection in magnetic fluids, Phys. Rev. E 64, 021406 (2001).
  * Ahmed  _et al._ [2018a] A. Ahmed, B. A. Fleck, and P. R. Waghmare, Maximum spreading of a ferrofluid droplet under the effect of magnetic field, Physics of Fluids 30, 077102 (2018a).
  * Odenbach [2009] S. Odenbach,  _Colloidal magnetic fluids: basics, development and application of ferrofluids_ , Vol. 763 (Springer, 2009).
  * Reimann  _et al._ [2005] B. Reimann, R. Richter, H. Knieling, R. Friedrichs, and I. Rehberg, Hexagons become the secondary pattern if symmetry is broken, Phys. Rev. E 71, 055202 (2005).
  * Zhou and Jing [2019] J. Zhou and D. Jing, Effects of vertical magnetic field on impact dynamics of ferrofluid droplet onto a rigid substrate, Phys. Rev. Fluids 4, 083602 (2019).
  * Ahmed  _et al._ [2018b] A. Ahmed, A. J. Qureshi, B. A. Fleck, and P. R. Waghmare, Effects of magnetic field on the spreading dynamics of an impinging ferrofluid droplet, Journal of Colloid and Interface Science 532, 309 (2018b).
  * Sahoo  _et al._ [2021] N. Sahoo, G. Khurana, D. Samanta, and P. Dhar, Collisional ferrohydrodynamics of magnetic fluid droplets on superhydrophobic surfaces, Physics of Fluids 33, 012012 (2021).
  * Yarin [2006] A. Yarin, Drop impact dynamics: Splashing, spreading, receding, bouncing…, Annual Review of Fluid Mechanics 38, 159 (2006).
  * Josserand and Thoroddsen [2016] C. Josserand and S. T. Thoroddsen, Drop impact on a solid surface, Annu. Rev. Fluid Mech. 48, 365 (2016).
  * Marengo  _et al._ [2011] M. Marengo, C. Antonini, I. V. Roisman, and C. Tropea, Drop collisions with simple and complex surfaces, Curr. Opin. Coll. Int. Sci. 16, 292 (2011).
  * Bijarchi and Shafii [2020] M. A. Bijarchi and M. B. Shafii, Experimental investigation on the dynamics of on-demand ferrofluid drop formation under a pulse-width-modulated nonuniform magnetic field, Langmuir 36, 7724 (2020).
  * Vieu and Walter [2018] T. Vieu and C. Walter, Shape and fission instabilities of ferrofluids in non-uniform magnetic fields, Journal of Fluid Mechanics 840, 455 (2018).
  * Shi  _et al._ [2014] D. Shi, Q. Bi, Y. He, and R. Zhou, Experimental investigation on falling ferrofluid droplets in vertical magnetic fields, Experimental Thermal and Fluid Science 54, 313 (2014).
  * Wells  _et al._ [2022] F. S. Wells, M. Broom, and G. R. Willmott, Deformation of ferrofluid drops falling in a non-uniform magnetic field, unpublished (2022).
  * Sahoo  _et al._ [2022] N. Sahoo, P. Dhar, and D. Samanta, Vertical magnetic field aided droplet-impact- magnetohydrodynamics of ferrofluids, Colloids and Surfaces A: Physicochemical and Engineering Aspects 633, 127872 (2022).
  * Yarin and Weiss [1995] A. L. Yarin and D. A. Weiss, Impact of drops on solid surfaces: self-similar capillary waves, and splashing as a new type of kinematic discontinuity, Journal of Fluid Mechanics 283, 141–173 (1995).
  * [25] Ferrotec corporation, EFH1 specifications and physical properties, https://ferrofluid.ferrotec.com/products/ferrofluid-educational-fluid/efh/efh1/, accessed: 2021-12-16.
  * [26] S. M. Blinder, Magnetic field of a cylindrical bar magnet, https://demonstrations.wolfram.com/MagneticFieldOfACylindricalBarMagnet/, accessed: 2022-2-17.
  * [27] See Supplemental Material at [URL will be inserted by publisher] for drop impact videos, a comparison of dimensionless numbers, details of contact line width measurement, and further data relating to maximum spreading and rim formation .
  * Antonini  _et al._ [2012] C. Antonini, A. Amirfazli, and M. Marengo, Drop impact and wettability: From hydrophilic to superhydrophobic surfaces, Physics of Fluids 24, 102104 (2012).
  * Zhang  _et al._ [2021] X. Zhang, B. Ji, X. Liu, S. Ding, X. Wu, and J. Min, Maximum spreading and energy analysis of ellipsoidal impact droplets, Physics of Fluids 33, 052108 (2021).
  * Chen and Lo [2006] C.-Y. Chen and L.-W. Lo, Breakup of thin films of micro-magnetic drops in perpendicular fields, Journal of magnetism and magnetic materials 305, 440 (2006).
  * Bushueva  _et al._ [2011] C. Bushueva, K. Kostarev, and A. Lebedev, Evolution of a ferrofluid floating layer under the influence of an inhomogeneous magnetic field, Magnetohydrodynamics 47, 207 (2011).


# DNA melting in poor solvent

Debjyoti Majumdar<EMAIL_ADDRESS>Institute of Physics, Bhubaneswar,
Odisha 751005, India  Homi Bhabha National Institute, Training School Complex,
Anushakti Nagar, Mumbai 400094, India

###### Abstract

The melting phase diagram of a double-stranded DNA in poor solvent is studied
using the pruned and enriched Rosenbluth method on a simple cubic lattice. As
the solvent quality is changed from good to poor, there is a non-monotonic
change in the melting temperature. First-order melting transition, as in good
solvent, gives way to continuous transition and then to further broadened
transitions where the order parameter smoothly becomes zero for sufficiently
poor solvent. This change in the melting behavior is accompanied by a
continuously varying critical exponent along the melting curve, hinting at a
non-universal nature of the melting transition. Further, we show that an
unbound phase can be achieved just by changing the solvent quality.
Importantly, our results conform to the experimental findings qualitatively.

##  I Introduction

The double-stranded (ds) DNA consists of two polymer strands connected along
the contour by hydrogen bonds between the bases forming base-pairs (bps)
watson1953 . The genetic information is encoded in the sequence of these
bases. Biological processes like gene expression and replication require the
separation of these two strands and then copying the genome sequence. In
vitro, this separation of the dsDNA into two single-strands (ss) can be
carried out by changing the temperature and is known as the thermal melting
transition. In theoretical studies, the melting transition is generally
considered in the context of a linear structure of DNA, a famous example being
the Poland-Scheraga (PS) model poland1966 . The strands, being polymers by
themselves, can further undergo a significant reduction in size through the
condensation transition, either in the ds (bound) or ss (unbound) phase, when
exposed to sufficiently poor solvent. It is this condensed form, in which, DNA
is generally found in vivo, e.g., in prokaryotes comm1 .

Condensation leads to the formation of higher-order structures through loops
and domains joyeux2015 . These higher-order structures are instrumental in
modifying the melting transition due to a change in the accessibility among
the inter- and intra-strand segments of the DNA bouligand1985 ; sikorav1991 .
Further, condensation promotes long-range order $\mathcal{O}$(cellular
dimensions), essential for biological processes like morphogenesis and also
with a possible role in optimizing the renaturation rate sikorav1991 .
Moreover, since changing the temperature in vivo is not a feasible option, it
is the solvent that has to play a key role in controlling the state of the
DNA, e.g., for enzymes, changing the local micro-environment of the DNA could
be a mechanism for altering the bps stability locally cui2007 . Therefore, it
is instructive to study the DNA melting transition in light of variable
solvent quality.

For DNA, the polarity of the solvent determines the solvent quality, with all
good solvents having high polarities, e.g., water, since DNA is also a highly
polar object, owing to the negatively charged backbone. Several experimental
studies have been performed in the last four decades baldini1985 ;
hammouda2009 ; feng1999 ; potaman1980 ; rupprecht1983 ; mikhailenko2000 ;
cui2007 , studying the role of solvents in DNA melting. A common observation
is that an increase in the alcohol concentration in the DNA solution
baldini1985  or increasing the carbon content in the hydrophobic tail of the
alcohol hammouda2009 , tends to decrease and then increase the melting
temperature. While the decrease in the melting temperature is poorly
understood, the increase in the thermal stability can be attributed to the
aggregation and precipitation of the bound DNA in the presence of alcohol
potaman1980 ; rupprecht1983 ; feng1999 . It was shown that melting occurs
mainly in the region where DNA has an elongated conformation, and therefore,
there exists a significant cooperative effect between the melting and the
unfolding transition mikhailenko2000 . In Ref. cui2007 , Cui et al. found that
at room temperature, dsDNA can be denatured into ssDNA in organic poor
solvents using AFM-based single molecule force measurement experiments. This
observation was further supported by molecular dynamics simulations, where the
double-stranded form split spontaneously into two single strands when pulled
at one strand from an aqueous to a poor solvent region. On the theoretical
side, efforts were made for DNA on special lattices to incorporate the solvent
effect, e.g., in Ref. foster2009 , the zipping/unzipping transition of DNA in
variable solvent quality is studied using the Bethe lattice. Again, using a
non-linear model of the DNA, known as the Dauxois–Peyrard–Bishop model, it was
shown that the melting temperature decreases as the solvent potential
increases macedo2014 .

Although several experimental papers exist investigating the effect of solvent
quality on DNA melting, it seems that a precise reason explaining why the
melting curve behaves the way it was found in Ref. baldini1985  is lacking.
Experimental difficulties in quantifying the interactions among segments from
different phases may have hindered the development of a quantitative
understanding. This is where the role of computer simulations is important. To
the best of our knowledge, a proper account of the solvent quality, along with
the base-pairing effect in the proper dimension, using computer simulation of
lattice models, is missing in the literature. Therefore, we believe that
further work is needed for a complete understanding.

In this article, we propose a model for the condensation of a melting DNA and
show that for a flexible DNA, of finite length, with complementary base-
pairing, the temperature induced melting transition is completely modified in
a poor solvent background comm7 . Interestingly, our simple coarse-grained
lattice model, which includes solvent only implicitly, through an effective
attraction among different segments (bound and unbound), is capable of
capturing qualitatively the features found previously in the experiments
baldini1985 ; mikhailenko2000 . Further, we address several aspects absent in
the previous studies, such as bubble statistics. The finiteness of the DNA in
our model makes it relevant from an experimental standpoint, since experiments
are performed upon finite systems.

This paper is organized in the following manner: in Sec. II, we describe the
DNA model for simultaneously studying the melting and condensation phenomena.
In Sec. III, we describe the method for simulating the dsDNA on the cubic
lattice. In Sec. IV, we discuss the findings, and conclude the paper in Sec.
V.

##  II DNA model and qualitative description

In our study, we use a coarse-grained model of the DNA causo2000 , where the
two strands of the dsDNA (say strands A and B) are represented by two self-
and mutually-avoiding linear polymers, embedded on a cubic lattice in the
dilute limit. One of the ends of the DNA is fixed at the origin, while the
other end remains free to wander [Fig. 1(a)]. Each strand contains $N$
monomers, with each monomer acting as a base and capable of forming a base-
pair (bp) of energy $-E_{bp}$, with the complementary monomer (same monomer
index) along the other strand, while both the monomers occupy the same lattice
site. Apart from the bp interaction, there is an effective attraction due to
the poor solvent background, between the non-bonded inter and intra-strand
segments, with an energy $-E_{nn}$ when one lattice constant apart
grassberger1997 ; a condition referred to as the nearest-neighbor (nn). In
configurations where a contact site of the two strands (i.e., a ds segment) is
nn with another occupied lattice site, which is part of either ss or ds, we
consider only one nn interaction [Fig. 1(b)]. Such a consideration arises from
the fact that any difference in the collapse of the ds and ss segments should
arise naturally from the difference in the rigidity of them, the former being
$25$ times more rigid than the latter, in reality.

(a)

(b)

Figure 1: (Color online) (a) Schematic representation of our DNA model with nn
representing nearest-neighbor interaction due to the poor solvent (red
squiggly lines). Broken bonds enclosed within bound segments are referred to
as bubbles, and as Y-fork if at the end. (b) Possible two-step configurations
contributing to non-bonded nn and the associated energies. The total bp energy
is $-6$, $-3$ and $0$ in the unit of $E_{bp}$, for the three figures (left to
right), respectively. The total nn interaction energy is $-1$ in the unit of
$E_{nn}$ for all three cases.

The model described above contains all the necessary details required to study
the melting and collapse transition individually and also the interplay of the
two. However, for simplicity, we have neglected the rigidity of the ds
segments, along with certain other aspects of the dsDNA, such as its helical
structure and sequence heterogeneity, under the assumption that the effect of
these properties can be studied separately. Moreover, considering all of the
factors in a single simulational study would substantially limit the maximal
system size.

Our model DNA can be decomposed into three parts: bound segments, broken bonds
enclosed within bound segments known as thermal bubbles, and a Y-fork at the
free end [Fig. 1(a)]. This categorization of the DNA into three parts is
important because, as we will see, the interplay of the nn interactions among
these three segments will determine the behavior of the melting curve. A
schematic representation of our model is shown in Fig. 1(a) and the scheme for
assigning energy to the non-bonded nn in Fig. 1(b). For a typical
configuration [Fig. 1(a)], the Hamiltonian can be written as:

$\mathcal{H}=-E_{bp}\sum_{i=1}^{N}\delta_{{\bf r}_{i}^{\rm A},{\bf r}_{i}^{\rm
B}}-E_{nn}\Bigl{\\{}\sum_{i=1}^{N-2}\sum_{j=i+2}^{N}\bigl{(}(\delta_{1,\mid{\bf
r}_{i}^{\rm A}-{\bf r}_{j}^{\rm A}\mid}\\\ +\delta_{1,\mid{\bf r}_{i}^{\rm
B}-{\bf r}_{j}^{\rm B}\mid})\times(1-\delta_{{\bf r}_{i}^{\rm A},{\bf
r}_{i}^{\rm B}}\delta_{{\bf r}_{j}^{\rm A},{\bf r}_{j}^{\rm
B}})+\delta_{1,\mid{\bf r}_{i}^{\rm AB}-{\bf r}_{j}^{\rm AB}\mid}\\\
\times\delta_{{\bf r}_{i}^{\rm A},{\bf r}_{i}^{\rm B}}\delta_{{\bf r}_{j}^{\rm
A},{\bf r}_{j}^{\rm B}}\bigr{)}+\sum_{i,j=1}^{N}\delta_{1,\mid{\bf r}_{i}^{\rm
A}-{\bf r}_{j}^{\rm B}\mid}\times(1-\delta_{{\bf r}_{i}^{\rm A},{\bf
r}_{i}^{\rm B}})(1-\delta_{{\bf r}_{j}^{\rm A},{\bf r}_{j}^{\rm
B}})\Bigr{\\}},$ (1)

where ${\bf r}_{i}^{\rm A}$ is the position vector of the $i$th monomer of
strand A, AB denote ds sites which are parts of both A and B, and
$\delta_{ij}$ is the Kronecker delta. Notice that each term in Eq. 1 is
symmetric under exchange of strand tag.

Figure 2: Different possible phases and transition pathways. Is a diagonal
line possible?

The lattice polymer model, which we use to model the DNA strands, can also be
used as a model for polymer collapse in the poor solvent with appropriate
interactions, known as the interacting self-avoiding walk (ISAW)
grassberger1997 . The collapse of a polymer in three dimensions is a
continuous phase transition grassberger1997 , with mean-fieldlike critical
exponents at the $\theta$-point, e.g., the exponent $\nu=1/2$ (except for
logarithmic corrections), as for the ideal random walk, since the excluded
volume interaction is exactly canceled out by the nn interactions with the
vanishing of the second virial coefficient for infinitely long chains. In the
dilute limit, a single polymer chain in three dimensions collapses at
$T=3.717$ for $E_{nn}=1$, or equivalently at $E_{nn}/T=0.269$. To look at the
collapse transition, one can observe the fluctuation of the nn interactions.
In our case this nn interaction can be divided into three parts which
contributes differently in modifying the melting transition; among ss segments
($nn_{\text{ss}}$), among ss and ds segments ($nn_{\text{ss-ds}}$), and among
ds segments ($nn_{\text{ds}}$) only.

When we consider base-pairing and condensation together, naively, one would
expect four different types of phases : (i) condensed-bound phase, (ii)
condensed-unbound phase, (iii) coiled-unbound phase, and (iv) coiled-bound
phase [Fig. 2]. Generally, we characterize the bound (or unbound) phase by the
average number of bound bps per unit length ($n_{c}$) being
$1(\text{or}~{}0)$. The condensed phase is characterized by the scaling of the
globule size ($R$) with the contour length ($N$), i.e., $R\sim N^{1/d}$ where
$d$ is the dimension of the embedding lattice, and the usual SAW exponent in
the coiled phase $R\sim N^{\nu}$ with $\nu=0.588$ in three dimensions. Each of
these phases, as shown in Fig. 2, can be reached from the other by moving
along any axis. What is interesting if there is any pathway connecting phases
(i) to (iii) ?, which corresponds to the diagonal line in Fig. 2. This will
indicate a transition from a bound-collapsed to an unbound-coiled state. In
such a scenario, the major difference in melting will be due to the emergence
of the globule surface tension. Such a transition might be important from the
biological point of view; however, it is out of the purview of the present
model. Therefore, we mention it only as a possibility and leave it for future
studies. Note that, for a flexible DNA, the completely bound and unbound
phases are equivalent, except for a length factor of 2, when considered
concerning a change in the nn interaction energy $E_{nn}$ at a fixed
temperature $T$.

##  III Simulation algorithm

For simulation, we use a sequential Monte Carlo algorithm, called the pruned
and enriched Rosenbluth method (PERM) grassberger1997 ; causo2000 . The chains
are built in steps, with a recursive call to the perm() subroutine, adding a
monomer on top of the lastly added monomer, both the strands at once. This is
possible since we know the free sites for the next step of both the strands
and a null atmosphere comm3  would result in an immediate exit of the current
call to the subroutine. The joint possibilities arising from the individual
atmospheres is $atmos=atmos_{\rm A}\times atmos_{\rm B}$, since, at any point
of growth, the individual strands are of the same length. At $n$th step, we
calculate the following quantity

$w_{n}=\sum_{atmos}p^{k}\times q^{l},$ (2)

where $p=\exp(\beta E_{bp})$ and $q=\exp(\beta E_{nn})$ are the Boltzmann
factors for the bp and the nn interaction respectively, and $k$ and $l$ are
their respective numbers. While $k$ take values of $0$ or $1$, $l$ can take
values ranging from $0,1\cdots(2d-1)$, where $d(=3)$ is the dimension of the
underlying lattice. The final choice is made among $atmos$ number of choices
in accordance with the importance sampling scheme. The total weight of a chain
of length $N$ is $W_{N}(\beta)=\prod_{n=1}^{N}w_{n}$, from which we obtain the
estimate of the partition sum $Z_{N}(\beta)=\langle W_{N}(\beta)\rangle$, the
average being taken over the number of started tours comm4 . Pruning and
enrichment at $n$th step is performed depending upon the ratio,
$r=W_{n}/Z_{n}$. Different schemes can be employed in choosing the enrichment
threshold and the number of copies. However, we find the algorithm to be
robust in its performance. In our simulations we have used; $r>1$ as the
enrichment threshold, with $c$ different copies each with weight
$\frac{W_{n}}{c}$ and $c=min(\lfloor r\rfloor,atmos)$, present configuration
continues to grow for $r=1$, and is pruned with probability $1-r$ if $r<1$ or
continue to grow with probability $r$ but with $W_{n}=Z_{n}$ prellberg2004 .
We take averages over $10^{7}$ tours and have set $E_{bp}=k_{B}=1$ in all of
our simulations.

For monitoring lattice site occupancy, we use the method of bit maps, where
the whole cubic lattice is stored into a one-dimensional array, and, occupied
lattice sites are indexed 1(or 2) if occupied by strand A (or B), 3 if
occupied by ds, and 0 if empty. Thus, it facilitates keeping a count of the
ss-ss, ss-ds, and ds-ds nn interactions separately, with an added advantage
that the time complexity for self-avoidance check reduces significantly
compared to any other method, e.g., binary search tree. However, this added
advantage comes at the cost of a finite lattice box size, and one has to be
careful about the box boundaries. However, in our simulations, we do not
expect extended configurations, e.g., if an external force pulls the strands
along a direction, in which case linked list methods will be a better choice.

Figure 3: (Color online) Phase diagram elucidating the interplay of melting
and condensation transition. Temperature ($T$) is varied along the horizontal
axis and solvent quality ($E_{nn}$) is varied along the vertical axis, with
$E_{bp}=1$ fixed. The blue dashed line is a guide to the eye, interpolating
the melting points ($T_{m}$). To the right, we have the high-temperature
unbound phase, and the low-temperature bound phase on the left of the melting
curve. Error bars in $T_{m}$ are of the size of the data points. The numerical
determination of the melting points is discussed in Sec. IV. The diagonal
continuous red line represents the tricritical collapse for an interacting
self-avoiding walk (ISAW) which satisfies $E_{nn}/T=0.269$ at all points
grassberger1997 .

##  IV Results and discussion

In good solvent (i.e. $E_{nn}=0$), our model exhibits a first-order melting
transition at $T_{m}=0.7455$ for $E_{bp}=1$ causo2000 . The $\theta$
temperature for the collapse transition of a lattice polymer in three
dimensions is $T_{\theta}=3.717$ for $E_{nn}=1$ grassberger1997 , in unit of
Boltzmann constant ($k_{B}=1$) for both melting and collapse. Our simulations
reproduced these individual results by setting either $E_{nn}=0$ or
$E_{bp}=0$, respectively. In Fig. 3, we plot the phase diagram in the $E_{nn}$
vs. $T$ plane, with $E_{bp}=1$ fixed. Horizontally, we plot the melting
transition induced by temperature ($T$), and vertically, the collapse
transition brought about by changing the strength of the nn interaction
($E_{nn}$). We estimate the melting point temperature $(T_{m})$ and the
associated critical exponent ($\phi$) from the scaling plots of the order-
parameter $(n_{c})$, and its fluctuation ($C_{c}$) [Fig. 4]. The error bars in
$T_{m}$ are obtained from the sensitivity of the scaling plots to the $T_{m}$
value. The average binding energy per unit length $(N)$ should scale as
$n_{c}\sim N^{\phi-1}g[(T-T_{m})N^{\phi}]$, while its fluctuation should scale
as $C_{c}\sim N^{2\phi-1}h[(T-T_{m})N^{\phi}]$, where $g$ and $h$ are scaling
functions. Tuning $T_{m}$ and $\phi$ to the appropriate values would then make
the data $N$ independent, resulting in data-collapse, as in Fig. 4. The
crossover exponent $\phi$ controls the sharpness of the melting transition,
with higher values representing sharper transitions. For $\phi=1$, the
transition is categorized as first-order and continuous if $\phi<1$. In the
following paragraphs, we present our main findings, and then try to explain
them.

As we increase the nn attraction ($E_{nn}$), initially, the melting
temperature ($T_{m}$) decreases linearly with the solvent parameter $E_{nn}$
[Fig. 3]. Then, after a point $(E_{nn}\approx 0.15)$, $T_{m}$ starts
increasing, and the melting curve changes slope as it approaches the
tricritical collapse line. This behavior is qualitatively similar to what was
observed in the experiments baldini1985 ; mikhailenko2000 ; hammouda2006 .
Baldini et al. baldini1985  have further found a change in the overall nature
of the melting curve with a different solvent type. However, our model allows
change only up to an effective attraction between monomers, be it for another
solvent or a different concentration of the same solvent. The change in the
melting temperature ($T_{m}$) is further associated with a continuous decrease
in the value of the exponent $\phi$ as $E_{nn}$ increases [Fig. 4(a) inset].
This indicates a change from first-order to continuous melting transition. In
Fig. 4(b) inset, we plot the variation of the peaks of $C_{c}$ curves with
length $N$, for $E_{nn}=0.06,0.13$ and $0.17$. $\phi$ extracted from the slope
of these curves are consistent with that of the scaling plots (not shown).
This further corroborates a change in $\phi$.

The $\theta$ point for the collapse transition is obtained by fixing $T$, and
varying $E_{nn}$. We use three different $T$’s with respect to the good
solvent case; in the bound phase $(T=0.71)$, near the melting point
$(T=0.741)$, and in the unbound phase $(T=0.75)$. We found that the $\theta$
point is merely affected by the presence of the base pairing interactions.
This is evident from a comparison between the nn fluctuation in our DNA model,
with that of a single SAW of the same linear length in three dimensions, and
interacting with the same ratio of $Enn/T$. Therefore, we draw a diagonal-
straight line satisfying $E_{nn}/T=0.269$ for the collapse transition in the
phase diagram [Fig. 3]. The melting temperature corresponding to $E_{nn}=0.2$
is $T_{m}=0.742(5)$ with $\phi=0.70(8)$. Therefore, $E_{nn}/T_{m}=0.269$ which
is just over the tricritical collapse line, is the meeting point of the
melting and collapse curves.

(a)

(b)

Figure 4: (Color online) Scaling plots for (a) average number of bound bps per
unit length ($n_{c}$), and (b) its fluctuation ($C_{c}$). For (a) and (b), we
use $\phi=0.91,0.84$ and $0.70(8)$, and $T_{m}=0.743(9),0.738(9)$ and
$0.742(5)$ for different values of the nearest neighbor energy
$E_{nn}=0.03,0.15$ and $0.2$, respectively. Data shown is for chain lengths
$N=1100$ to $2000$ at an interval of $300$. (a) Inset: Crossover exponent
$\phi$ used in the scaling plots of the melting transition for different
$E_{nn}$ values. The continuous curve is a guide to the eye. (b) Inset: Log-
log plot of the $C_{c}$ peaks for different system sizes ($N$). The grey
dashed lines are fit to the points, with slope=$0.81,0.74$ and $0.63$, for
$E_{nn}=0.06,0.13$ and $0.17$, respectively.

As the melting curve crosses the collapse line (at $E_{nn}=0.2$ and
$T=0.742$), the temperature induced melting transition is gradually smeared
with the order parameter smoothly going to zero. This gradual smearing may
indicate the onset of an infinite-order transition. For finite-order phase
transitions, the finite-size scaling (FSS) is a powerful tool to extrapolate
the power-law behavior of the thermodynamic observables around the
criticality. It is performed by scaling the distance to the critical point and
the observable height as functions of the system size and the critical
exponents. However, the use of FSS is only restricted to problems with one
diverging correlation length. Therefore, as the melting curve approaches the
collapse line, with the emergence of another diverging length scale, besides
that for melting, a proper estimate of the melting temperature becomes
difficult, which limits our use of the FSS little above $E_{nn}=0.2$. Beyond
this value, we find it difficult to perform a good data collapse, apparently
because the two arms of the $C_{c}$ curves about the central peak collapses
for two different values of $T_{m}$, denoting a range of $T$ values for the
melting transition. However, such observations could be the result of poor
performance and therefore sampling in the collapsed region or may require the
usage of extended finite-size scaling (EFSS) analysis santos1981 . Notice
that, considering the system is placed at a temperature (say $T=0.741$) which
is considered to be a bound phase in the good solvent, the system can be
carried to a completely unbound state just by changing the solvent quality at
a fixed temperature, and much before the collapse transition [Fig. 3].

In a good solvent, the energetic gain over entropy leads to a bound phase,
thereby, lowering the free energy. Now, in the presence of nn attraction,
there is an added energetic gain. However, the nn interactions play a dual
role in modifying the melting transition. While $\text{nn}_{\text{ds}}$ should
help to maintain the bound phase, $\text{nn}_{\text{ss}}$ should support the
unbound phase, thereby, stabilizing the unbound phase. To supplement the fact
that the initial (i.e for small values of $E_{nn}$) decrease in $T_{m}$ is
solely due to the $\text{nn}_{\text{ss}}$, we carried out separate simulations
where we turned off the $\text{nn}_{\text{ss}}$ interactions. We found that
the melting temperature only increases, even for small $E_{nn}$ values. The
$\text{nn}_{\text{ss}}$ interactions are contributed by the broken bonds from
the thermal bubbles and the Y-fork at the end. Since, for small values of
$E_{nn}$, in the vicinity of melting point comm2 , configurations with large
bubbles that can contribute to $\text{nn}_{\text{ss}}$ should be rare, it is
the mechanism by which the Y-fork opens up that ultimately determines the
$\text{nn}_{\text{ss}}$.

(a)

(b)

Figure 5: (Color online) (a) Average number of bubbles per unit length
($n_{b}$) for chain length $N=2000$. (Inset) Semi-log plot of average bubble
size ($\ell_{b}$) along the chain, excluding the Y-fork, for different values
of $E_{nn}$ and chain length $N=2000$. For $\ell_{b}$, average is taken along
the contour length and then over many such configurations. The curves are
shifted by their corresponding $T_{m}$’s to make the melting point coincide
with zero. (b) Semi-log plot of the difference between the nn contributions
from the ss and ds segments per unit length ($N$). The negative part of the
curves for $E_{nn}=0.01$ and $0.02$ are not shown. (Inset) Fraction of broken
bonds per unit length in the Y-fork ($f_{Y}$), for different values of
$E_{nn}$.

On the other hand, with increasing $E_{nn}$, the average number of bubbles per
unit length $(n_{b})$ decreases [Fig. 5(a)], while the average bubble length
($\ell_{b}$) along the chain increases [Fig. 5(a) inset], along with a delayed
Y-fork ($f_{Y}$) opening around the melting transition [Fig. 5(b) inset]. This
picture is consistent with the melting of Gaussian strands where the majority
of the broken bonds are in the bubbles pre-melting majumdar2020 . This delay
in the Y-fork opening [Fig. 5(b) inset], along with a decreasing number and
increasing average size of bubbles, gives rise to a stabilizing effect for
sufficiently large $E_{nn}$ values ($E_{nn}\gtrsim 0.15$) upto the collapse
line [Fig. 3] majumdar2020 . Further, deep in the collapsed phase ($E_{nn}\sim
1.0$) and within the studied temperature window ($0.73<T<0.75$), there seems
to be no transition at all mikhailenko2000 , with the DNA in an unbound phase.

From the discussion in the preceding paragraph, we expect, that the slope of
the phase curve (i.e. $\frac{dE_{nn}}{dT}$) to be determined by a quantity
that changes sign as we move along the 1-dimensional critical manifold [Fig.
3]. One such quantity would be the difference between the various components
contributing to nn interactions, i.e., $\text{nn}_{\text{ss}}$ and
$\text{nn}_{\text{ds}}$ [Fig. 5(b)]. It is positive for large values of
$E_{nn}$ on the bound side of the corresponding $T_{m}$’s, and rapidly changes
sign for smaller $E_{nn}$ values, as expected for a sharper transition. From a
thermodynamic viewpoint, this is equivalent to the Clasius-Clapeyron equation,
$\frac{dE_{nn}}{dT}=\frac{L}{T\Delta nn}$, where
$L=T(S_{\text{ss}}-S_{\text{ds}})$ is the latent heat, $S$ is the entropy, and
$\Delta nn=(nn_{\text{ss}}-nn_{\text{ss-ds}}-nn_{\text{ds}})$. Below the curve
satisfying $E_{nn}/T=0.269$, polymers are in the swollen phase. Therefore,
both $nn_{\text{ds}}$ and $nn_{\text{ss}}\rightarrow 0$ for
$N\rightarrow\infty$. If there is latent heat i.e. $L\neq 0$, the slope would
be $\infty$, i.e., the phase curve is a vertical line. If $L=0$, as for
continuous transitions, then for finite $N$, $S_{\text{ss}}-S_{\text{ds}}\sim
N^{a}>0$ with $a<1$, so that the entropy difference per unit length $\delta
S/N\rightarrow 0$. So, we must have $\Delta nn\sim N^{a}$, for $N$ dependence
to cancel, and if $\Delta nn<0$ on the phase boundary then
$\frac{dE_{nn}}{dT}<0$. This inequality will change in the immediate vicinity
of the critical points for larger $E_{nn}$ values, for the slope to become
positive [Fig. 5(b)].

For an explanation of the continuous variation of the crossover exponent
$\phi$, we note that, pre-condensation, the two main length scales of the
problem is given by the two correlation lengths related to the geometry of the
polymeric strands $(\zeta_{1})$ and another related to the mean diameter of
the thermal bubbles $(\zeta_{2})$ which is expected to be related to the
average bubble length ($\ell_{b}$) causo2000 . The critical exponents $\nu$
and $\nu_{T}$, control these length scales, $\zeta_{1}$ and $\zeta_{2}$
respectively, and using which we can write causo2000 ; santos1981 ,

$\phi=\nu/\nu_{T}.$ (3)

Now, with increasing $E_{nn}$, $\nu_{T}$ increases with an increase in the
average bubble length $(\ell_{b})$ along the chain [Fig. 5(b) inset], while
$\nu=1/3$ only as the melting curve crosses the collapse line, in the
thermodynamic limit ($N\rightarrow\infty$), thereby leading to a gradual
lowering in the value of $\phi$. In other words, the introduction of effective
attraction due to poor solvent, introduces long-range interaction among
monomers which are far apart along the chain contour. This, changes the way
the correlations become infinite as $T\rightarrow T_{m}^{-}$,
$\zeta_{2}\sim(T-T_{m})^{-\nu_{T}}$, thus denoting a change in the exponent
$\nu_{T}$baxter_book .

Again, for continuous transitions, $\phi$ can be estimated from the bubble-
size distribution (bsd) at the critical point, using the expression
$P(\ell_{b})\sim\ell_{b}^{-c}$, where $c=\phi+1$ is the bubble-size exponent
(bse) carlon2002 . In Fig. 6 we have plotted the bsd at the critical points of
the melting transition for $E_{nn}=0.03,0.15,0.17,0.18,0.19$ and $0.2$, which
clearly shows that $\phi$ decreases with increasing $E_{nn}$. This,
corroborates the fact that the continuous change in $\phi$ is, in fact a
genuine change and not some artifact of the FSS originally used for estimation
of $\phi$. However, we see differences in the value of the $\phi$ obtained
from the two different methods. Interestingly, according to the estimation of
$\phi$ from bse ($c$), there should be another tricritical point along the
melting curve at $c=2$, apparently where the melting curve changes slope.

Figure 6: Bubble-size $(\ell_{b})$ distribution $(P(\ell_{b}))$ at the
critical points for $E_{nn}=0.03,0.15,0.17,0.18,0.19$ and $0.2$. Data starting
from $E_{nn}=0.15$ are shifted by factor of $2$ to $10$ in steps of $2$.
Fitting of the data is performed in the range $\ell_{b}=20$ to $200$, which
gives $c=2.3,2.07,1.99,1.94,1.89$ and $1.83$ for increasing $E_{nn}$ starting
from $0.03$ respectively.

In the study of critical phenomena, a continuous variation of the exponents,
often induced by a marginal interaction, in the renormalization group sense,
is rarely observed khan2017 ; martins2018 . For the unperturbed case, the two
universality classes for melting, using PS and alike models, are for ideal
random walk ($\phi=1/2$ and continuous melting) and self-avoiding walk
($\phi=1$ and first-order melting). The addition of solvent interaction acts
as a relevant perturbation that brings the system out from its universality
class as of the unperturbed case. However, whether the transition is non-
universal or there is an underlying weak-universality, where the critical
exponents vary, but the ratio of two critical exponents remains constant, is
to be investigated further. Such a change in the critical exponents is not new
and was observed previously upon adding the solvent effect for adsorption of
polymers on surface martins2018 .

##  V Conclusion

In conclusion, we studied the possibility of various phases resulting from the
interplay between the base-pairing and nn interaction due to poor solvent. It
is evident that the melting transition is continuous, even with excluded
volume interaction, for finite nn interaction and much before the collapse
transition. As we decrease the solvent quality, initially, there is a decrease
in the melting temperature, followed by a subsequent increase in the later
stages as the melting curve approaches the collapse line. The slope of the
melting curve is determined by a competition between the nearest-neighbor
interactions from the double-stranded and single-stranded segments. While
scaling holds in the coil region, universality does not, leading to a
continuously varying critical exponent baxter_book . We provide evidence from
three different methods to support this fact; first, from the FSS; second,
from the average size of bubbles; and third from the bubble-size distribution.
Apparently, there is a gradual crossover between the two universality classes
for DNA melting. Melting as a transition is suppressed in the deep-collapsed
region with the breakdown of FSS. However, there could be a smooth crossover
between the bound and unbound phases over a large temperature range. Apart
from that, we also see a possibility of solvent-mediated melting. Importantly,
our results are in qualitative agreement with the previous experimental
results baldini1985 ; hammouda2009 .

Including phase separation effects in the solvent could be one of the
interesting extensions of this work worth investigating. One can easily extend
this model to include different rigidity for the ss and ds segments
majumdar2021 . Adding rigidity to the ds segments can change the results
significantly since, the solvent will affect the rigid and flexible components
differently bastolla1997 . The change in the effective excluded volume
interaction should also change the bulk elastic response majumdar2020 . Note
that in a real situation, the Y-fork can exist at both ends, due to which the
slope of the melting curve might change. Our results would help set the in
vitro environmental conditions for the experiments.

Finally, it is often seen that biological systems at all scales are often
poised at criticality, which can be both statistical and dynamical.
Maintaining itself at the criticality has its advantages, e.g., high
susceptibility to external perturbations mora2011 . However, how this critical
point is searched for and maintained by these biological systems is another
puzzle since the critical region spans only around a small region (or a single
point in the thermodynamic limit) of the parameter space. The possibility that
DNA might poise itself near the critical point or if a diagonal line like in
Fig. 2 is possible is something that needs to be investigated and how it can
be advantageous for the system, given that DNA is such an important functional
unit of the cell.

## Data Availability

Data and information about the simulations are available upon reasonable
request to the author.

## Conflict of Interest

The author has no conflict of interest to declare that are relevant to the
content of this article.

## Acknowledgement

DM is indebted to Somendra M Bhattacharjee for insightful discussions.
Simulations were carried out at the SAMKHYA computing facility at IOPB.

## References

  * (1) J. D. Watson and F. H. C. Crick, Molecular Structure of Nucleic Acids: A Structure for Deoxyribose Nucleic Acid, Nature (London) 171, 737-738 (1953).
  * (2) D. Poland and H. A. Scheraga, Occurrence of a Phase Transition in Nucleic Acid Models, J. Chem. Phys. 45, 1464 (1966); D. Poland and H. A. Scheraga, Phase Transitions in One Dimension and the Helix-Coil Transition in Polyamino Acids, J. Chem. Phys. 45, 1456 (1966).
  * (3) In eukaryotes, however, the folding mechanism is rather sophisticated and requires the involvement of histomers.
  * (4) M. Joyeux, Compaction of bacterial genomic DNA: clarifying the concepts, J. Phys.: Condens. Matter 27, 383001 (2015).
  * (5) Y. Bouligand, Brisures de symétrie et morphogenése biologique, C. R. Acad. Sci. (Paris) Série Gén. La Vie des Sciences, 2, no. 2, 121-140 (1985).
  * (6) J.-L. Sikorav and G. M. Church, Complementary Recognition in Condensed DNA: Accelerated DNA Renaturation, J. Mol. Biol. 222, 1085 (1991).
  * (7) S. Cui, J. Yu, F. Kuhner, K. Schulten and H. E. Gaub, Double-Stranded DNA dissociates into single strands when dragged into a poor solvent, J. Am. Chem. Soc. 129, 14710 (2007).
  * (8) G. Baldini, H. Fu-Hua, G. Varani, L. Cordone, S. L. Fornili and G. Onori, DNA Melting Induced by Alcohols: Role of the Solvent Properties, Il Nuovo Cimento 6 D, N. 6 (1985).
  * (9) B. Hammouda, Insight into the denaturation transition of DNA, Int. J. of Bio. Macro. 45, 532–534 (2009).
  * (10) V. N. Potaman, Yu. A. Bannikov and L. S. Shlyakhtenko, Sedimentation of DNA in ethanol-water solutions within the interval of B $\rightarrow$ A transition, Nucleic Acids Res. 8, 635 (1980).
  * (11) A. Rupprecht et al., Mechanochemical study of conformational transitions and melting of Li-, Na-, K-, and CsDNA fibers in ethanol–water solutions, Biopolymers 34, 897 (1994); J. Piskur and A. Rupprecht, Aggregated DNA in ethanol solution, FEBS Lett. 375, 174 (1995).
  * (12) Y. Fang, T. S. Spisz and J. H. Hoh, Ethanol-induced structural transitions of DNA on mica, Nucleic Acids Res. 27, 1943-1949 (1999).
  * (13) S. V. Mikhailenko et al., Interplay between Folding/Unfolding and Helix/Coil Transitions in Giant DNA, Biomacromolecules 1, 597-603 (2000).
  * (14) D. P. Foster and C. Pinettes, Bethe approximation for a DNA-like self-avoiding walk model with variable solvent quality, Phys. Rev. E 79, 051108 (2009).
  * (15) D. X. Macedo, I. Guedes and E. L. Albuquerque, Thermal properties of a DNA denaturation with solvent interaction, Physica A 404, 234–241 (2014).
  * (16) Here, poor solvent refers to any finite attraction (i.e. $E_{nn}>0$) between non-bonded nearest neighbors; also see Sec. II.
  * (17) M. S. Causo, B. Coluzzi and P. Grassberger, Simple model for the DNA denaturation transition, Phys. Rev. E 62, 3958 (2000).
  * (18) P. Grassberger, Pruned-enriched Rosenbluth method: Simulations of $\theta$ polymers of chain length up to 1 000 000, Phys. Rev. E 56, 3682 (1997).
  * (19) Here, atmosphere refers to the number of possible free sites in the next step rechnitzer2008 .
  * (20) A Rechnitzer and E J Janse van Rensburg, Generalized atmospheric Rosenbluth methods (GARM), J. Phys. A: Math. Theor. 41, 442002 (2008).
  * (21) A tour refers to the total number of individual configurations, generated between two successive return to the main() function of our C program and arranged in a tree-like topology with enrichments representing the branches.
  * (22) T. Prellberg and J. Krawczyk, Flat Histogram Version of the Pruned and Enriched Rosenbluth Method, Phys. Rev. Lett. 92, 120602 (2004).
  * (23) B. Hammouda and D. Worceste, The Denaturation Transition of DNA in Mixed Solvents, Biophysical Journal, 91, 2237–2242 (2006).
  * (24) R. R. dos Santos and R. B. Stinchcombe, Finite size scaling and crossover phenomena: the XY chain in a transverse field at zero temperature, J. Phys. A 14, 2741 (1981).
  * (25) For first-order transition most of the broken bonds are in the Y-fork near melting majumdar2020 .
  * (26) D. Majumdar and S. M. Bhattacharjee, Softening of DNA near melting as disappearance of an emergent property, Phys. Rev. E 102, 032407 (2020).
  * (27) R. J. Baxter, Exactly solved models in statistical mechanics, (Academic Press, London, 1982).
  * (28) E. Carlon, E. Orlandini and A. L. Stella, Roles of Stiffness and Excluded Volume in DNA Denaturation, Phys. Rev. Lett. 88, 198101 (2002).
  * (29) N. Khan, P. Sarkar, A. Midya, P. Mandal, and P. K. Mohanty, Continuously varying critical exponents beyond weak universality, Sci. Rep. 7, 45004 (2017).
  * (30) P. H. L. Martins, J. A. Plascak and M. Bachmann, Adsorption of flexible polymer chains on a surface: Effects of different solvent conditions, J. Chem. Phys. 148, 204901 (2018).
  * (31) D. Majumdar, Elasticity of a DNA chain dotted with bubbles under force, Phys. Rev E 103, 052412 (2021).
  * (32) U. Bastolla and P. Grassberger, Phase Transitions of Single Semistiff Polymer Chains, J. Stat. Phys., 89 1061 (1997).
  * (33) T. Mora and W. Bialek, Are Biological Systems Poised at Criticality?, J Stat Phys 144, 269-302 (2011).