# An $O(k\log{n})$ algorithm for prefix based ranked autocomplete

Dhruv Matani

(September 2, 2011)

###### Abstract

Many search engines such as Google, Bing & Yahoo! show search suggestions when
users enter search phrases on their interfaces. These suggestions are meant to
assist the user in finding what she wants quickly and also suggesting common
searches that may result in finding information that is more relevant. It also
serves the purpose of helping the user if she is not sure of what to search
for, but has a vague idea of what it is that she wants. We present an
algorithm that takes time proportional to $O(k\log{n})$, and $O(n)$ extra
space for providing the user with the top $k$ ranked suggestions out of a
corpus of $n$ possible suggestions based on the prefix of the query that she
has entered so far.

###### Contents

  1. 1 What is prefix based ranked phrase auto-complete?
  2. 2 Where is phrase auto-complete used?
  3. 3 What are the problems that phrase auto-complete solves?
  4. 4 Why does auto-complete need to be fast (responsive)?
  5. 5 Problem Statement
  6. 6 Existing approaches
    1. 6.1 Naive approach
    2. 6.2 Space intensive approach
    3. 6.3 Ternary search trees
    4. 6.4 The TASTIER approach
  7. 7 Proposed algorithm
  8. 8 Approximate matching for auto-complete
    1. 8.1 Using exact prefix-match auto-complete for approximate match auto-complete
  9. 9 Implementation & Performance Tests
  10. 10 Conclusion

##  1 What is prefix based ranked phrase auto-complete?

Given a set of $n$ strings, each having a certain weight, the problem of
finding the set of $k$ heaviest strings each of which have the prefix $q$ is
the problem of prefix based ranked phrase auto-complete.

##  2 Where is phrase auto-complete used?

Many search engines such as Google, Bing & Yahoo! show search suggestions when
users enter search phrases on their interfaces. These suggestions are meant to
assist the user in finding what she wants quickly and also suggesting common
searches that may result in finding information that is more relevant. It also
serves the purpose of helping the user if she is not sure of what to search
for, but has a vague idea of what it is that she wants.

IMDB uses search suggestions to ensure faster search results since movie
titles are mostly unique after the first few characters of typing. Users that
find it difficult to type or users on mobile devices with constrained input
methods can get to their results much faster because of auto-complete.

##  3 What are the problems that phrase auto-complete solves?

Many a times, the user can think of many keywords or phrases that can be used
to describe the concept or idea that is being searched for. Such descriptions
generally use similar meaning words. Users don’t know in advance what
vocabulary most of the literature or prior work in that field uses. Hence,
they have to painstakingly try each combination of keywords till they fine
what they are looking for.

For example, when a user is searching for auto complete, it could actually be
referred to by articles online as either autocomplete, or auto suggest, or
search suggestions, or find as you type. The user would have to try all of
them before settling on the one that returns the most relevant results. If
most of the alternatives, 3 in this case, have the same prefix, then the user
could just start typing, and the system could suggest possible completions as
new characters are entered. This greatly reduces the trial & error that the
user has to perform. Furthermore, many users aren’t even aware of entering
different search terms. With a high probability, these users will not find
what they are looking for.

##  4 Why does auto-complete need to be fast (responsive)?

The average typing speed (while composing text) is 19 words per minute[9]. The
average word length in the English language is about 9 characters[8].
Combining the two, we notice that the average time for typing a single
character is 351ms. Accounting for a network round-trip time of 200ms[7, 5],
and client processing time of 100ms, we have only about 50ms left to do the
processing at our end. The auto complete application needs to return a list of
suggestions within 50ms for it to be useful to the user. It needs to do this
for potentially every keystroke typed in by the user (there are optimizations
we can do at the user interface layer, but we won’t discuss them here).

##  5 Problem Statement

Given a list of $n$ phrases, along with a weight for each phrase and a query
prefix $q$, determine the $k$ heaviest phrases that have $q$ as their prefix.
We assume the average phrase length to be constant and shall not account for
it in the complexity calculations.

##  6 Existing approaches

We discuss three existing approaches used by various applications to provide
search suggestion for a find as you type experience on their web pages.

###  6.1 Naive approach

The naive approach involves pre-processing the input ($n$ phrases) by sorting
them in lexicographical order. Binary Search is then used to locate the
beginning and end of the candidate list of phrases. Each phrase in this
candidate list $cl$ has $q$ as its prefix.

Now, all the terms in the candidate list $cl$ are sorted in non-increasing
order by weight and the top $k$ are selected for projection.

This requires $O(log{n})$ for the binary search and $O(cl\log{cl})$ for
sorting the candidate list. This approach also requires extra space
proportional to $O(cl)$ for storing the candidate list $cl$ before sorting.
The total complexity of this approach is $O(\log{n})+O(cl\log{cl})$.

We can see that if the query prefix $q$ is short or matches many phrases, then
the candidate list $cl$ will be large and the latter factor will dominate the
complexity. We ideally want an algorithm that is not input-sensitive.

A slight variation on the above technique would be to notice that $k$ is
generally quite small (around 16), so fetching the top 16 candidates from the
candidate list of size $cl$ will cost only $O(k.cl)$ which is faster than what
we have before. This brings our runtime down to $O(\log{n})+O(k.cl)$. However,
it is still input-sensitive and that is undesirable.

###  6.2 Space intensive approach

We could maintain a lexicographically sorted list of all possible phrase
prefixes, with all the entries with the same prefix sorted by weight. This
implies sorting by two keys, namely prefix, weight. A lookup is as simple as
performing a binary search to locate the first occurrence of the query $q$ and
reading off the next $k$ entries, as long as they have the same value as the
query $q$, since they are already sorted by weight.

This technique uses extra space that is proportional to $O(string\ length)$
for each phrase in the initial set of $n$ phrases. This means that extra space
to the order of $O(n.string\ length/2)$, which roughly translates to $15n$ in
our case, is needed. We can not ignore the constant factor here since it is
very significant. For an initial corpus of $5GB$ (17 million phrases of length
30 characters each), we would land up using extra space proportional to $75GB$
in our case.

The runtime complexity for the pre-processing step is $O(15n)$ and the runtime
complexity for querying is $O(k\log{n})$.

This approach is attractive for small data sets, but starts getting very
costly in terms of memory requirements for larger data sets.

This approach is used by the redis data structure store for providing auto-
complete facilities[2].

###  6.3 Ternary search trees

We can optimize the previous approach for space by using a trie data structure
instead of duplicating each prefix of a phrase every time we encounter it.
This gives rise to the ternary tree data structure. We need to decide in
advance what the maximum value of $k$ is that we would like to support since
this method involves pre-computing a list of the $k$ heaviest completions for
a given prefix and storing them at that node. Since the actual completions can
be stored elsewhere, we incur a penalty for storing $k$ pointers to phrases at
every ternary tree node, and not a penalty for the $k$ complete phrases. We
also incur a penalty for unused trie-nodes at every step (in the
implementation that trades space for speed).

We can compute the extra space required by assuming that we would need an
average of 10 pointers (and not 16) at every ternary tree node. Since the
ternary tree structure is highly input dependent, we can’t really perform an
estimation here. We abandon this approach since it is highly input dependant
and can degenerate to requiring a lot of space.

###  6.4 The TASTIER approach

Researchers from Tsinghua University and The University of California, Irvine
have implemented a system for predicting user input and performing an auto-
complete based on partial input from the user[4]. Their technique differs from
other (and the one presented here) in that it treats the query as a set of
keywords, and support prediction of any combination of keywords that may not
be close to each other in the data. Other techniques (and the one presented
here) do prediction using sequential data to construct their model, therefore
they can only predict a word or phrase that matches sequentially with
underlying data (e.g., phrases or sentences).

However, their technique relies on:

  * •

Maintaining the previous result set for a past query from the user. This
increases memory costs on the server and leads to increased memory pressure
when serving many concurrent users

  * •

Performing a union of potentially many candidate lists (which is time
consuming)

On a corpus of 1 million entries, their system is able to answer $1$ query per
$20$ ms. The result set size here is $10$. This means that their system can
handle at best $50$ qps, which is too low for search-engine traffic.

##  7 Proposed algorithm

We present an algorithm that takes time proportional to $O(k\log{n})$, and
$O(n)$ extra space for providing the user with the top $k$ ranked suggestions
out of a corpus of $n$ possible suggestions based on the prefix of the query
that she has entered so far.

This algorithm requires us to maintain the set of phrases in a sorted array
(Figure 1) so that we can binary search over it using our query prefix q.

Figure 1: The sorted array of all completion phrases

We additionally also maintain a Segment Tree[6] (Figure 2) (which is nothing
but an Interval Tree) that stores the maximum weight of the interval at every
node. In the referenced figure, each leaf node shows the word that the node
represents, though the word is actually not stored at that node. Instead, the
number in parenthesis is stored in the node, and it denotes the weight of the
phrase that the node is associated with. Each internal node shows a range (in
square brackets) denoting the lower and upper indexes (both inclusive) in the
sorted word array that this internal node represents. Each internal node
stores the maximum weight of the nodes in the range. This weight is shown in
parenthesis in the diagram.

Figure 2: The Segment Tree keyed by phrase frequency, shown in parenthesis
after each word. The internal nodes are of the format [range](maximum key in
subtree)

While querying, we maintain an ephemeral max-heap (or max-priority queue) of
ranges that contain the maximum weighted candidate phrases. The nodes of this
max-heap are keyed by the maximum weight of phrases stored by the range of
indexes that the node represents. For example, if the node represents range
[0–20] and its key is $45$, it means that the heaviest phrase in the range
[0–20] has a weight of $45$. We always start by entering a single element
representing the entire range of interest into the max-heap and then spliting
the range containing the heaviest phrase at every step.

For the example shown in Figure 3, we show the first step in the series of
operations performed on the Segment Tree while searching for the top $4$
phrases that begin with the string b in the set shown in Figure 1. The search
starts by inserting the entire range of indices that contain the candidate
words ([0-9] in this case) into the max-heap. The top of the max-heap is then
popped and the index of the highest phrase in that range is selected. We see
that the word bacon with the weight 18 is selected and a split at the index 8
gives rise to 2 ranges, [0-7] & [9-9]. These are inserted into the max-heap
and the procedure is repeated. Figures 4, 5, & 6 show the action of the
algorithm on the sample data set.

Figure 3: Search for the maximum weighted phrase in the range [0–9]. The green
arrows show the paths taken and the red arrows show the paths not taken Figure
4: Search for the maximum weighted phrase in the range [0–7]. The green arrows
show the paths taken and the red arrows show the paths not taken Figure 5:
Search for the maximum weighted phrase in the range [0–6]. The blue arrows
show the paths and branches taken. With the help of the blue arrows, we are
able to trace the complete path of the query to the leaf nodes. The blue boses
indicate the nodes where the query for the maximum weighted node will be
restarted. Figure 6: Search for the maximum weighted phrase in the range
[0–3]. The green arrows show the paths taken and the red arrows show the paths
not taken.

Figure 7 shows the contents of the max-heap at every step of the process of
extracting the top $4$ phrases having the prefix $b$. The left column
indicates the state of the priority queue (a max heap) at every stage and the
right column indicates the output produced (phrase) at every stage along with
the index and weight (score) of that phrase. The output is generated in non-
increasing order of weight.

Figure 7:

##  8 Approximate matching for auto-complete

Researchers from Tsinghua University and The University of California, Irvine
have implemented a system for performing a fuzzy keyword match between the
user input and the indexed corpus[3]. Their technique relies on creating
multiple keyword candidate lists per keyword entered and intersecting them to
produce the output. They optimize this process by maintaining cached result
sets for previous input by the same user so that they can just trim these
lists when a new character is added.

On a data set with 1 million entries, their technique takes up to $5ms$ for
finding up to $10$ results having prefixes of length $3$ ($3$ characters in
length). This case is sure to be hit because users start off by typing at
least $2-3$ characters of their query before hoping for relevant suggestions.
In the worst case, this system can handle a load of $200$ qps. For longer
prefixes (i.e. upwards of $6$ characters in length), the query time drops
shapply to below $0.1ms$, which enables them to answer about $10,000$ qps. The
responsiveness of our technique on the other hand is independent of the number
of characters entered by the user and depends only on the corpus size. Hence,
it’s responsiveness is much more predictable, which is desirable in a
production environment.

###  8.1 Using exact prefix-match auto-complete for approximate match auto-
complete

If we pre-process the data and queries in some specific manner, we can use the
exact-match auto-complete implementation to perform approximate match auto-
complete111This need not be an approximate prefix-match auto-complete
depending on how we pre-process the data. This can be accomplished by
performing one or more of the following transformations on the phrases in the
corpus as well as the query string. The transformations must be performed in
the same order as given below:

  * •

Remove all stop-words such as a, the, have, has, of, etc …

  * •

Remove all punctuation marks, vowels, and white-spaces from the phrase

  * •

(optionally) Convert all consonants to their numeric equivalents according to
the transformations mentioned in the Soundex Algorithm222US patent 1261167, R.
C. Russell, issued 1918-04-02 & US patent 1435663, R. C. Russell, issued
1922-11-14

  * •

Collapse all repetitive runs of characters (or numbers) into a single
character (or number). e.g. rttjdddl becomes rtjdl

Now, these transformations are stored as the phrase in the auto-complete
index. The original (untransformed) phrase is stored against the transformed
phrase so that once the lookup is done, the original phrase may be returned.

When we query the auto-complete index, we perform the same transform on the
query string and find those strings that share a prefix with the query string.

In practice, we can index phrases at every stage of the transformations and
take the union of the result set after querying each index for the
corresponding query string.

##  9 Implementation & Performance Tests

The idea mentioned above has been implemented in an application called lib-
face333Available for download at https://code.google.com/p/lib-face/. lib-face
is written in the C++ programming language and uses the mongoose web-
server444mongoose is an in-process HTTP web-server available for download at
http://code.google.com/p/mongoose/ to serve requests. Once you load the corpus
into lib-face, you can query the highest weighted phrase that has a certain
prefix by sending an HTTP request to it.

Here are the results of a test run on an Amazon EC2 instance:

Operating System Ubuntu 10.04 Data Set Size (number of entries) 14,000,000
Data Set Size (bytes) 312MiB Result set size 32 Queries Per Second 6,800
Memory (Resident Size) 1,540MiB

Notes:

  1. 1.

Not all of the 1,540MiB of resident memory is used since C++’s std::vector
uses a doubling strategy to grow, which results in half the memory being
actually unused.

  2. 2.

There are CPU and network overheads when running an application on a virtual
machine v/s on real hardware. We think it’s best to benchmark on a VM rather
than real hardware since that seems to be a very common deployment paradigm
these days. Comparing with benchmarks running on real-hardware would be unfair

  3. 3.

All our benchmarks are done assuming a result set size of $32$ whereas other
groups have assumed a result set size of $10$. This will invariably reduce the
number of queries per second that we can answer.

##  10 Conclusion

The Segment Tree method we described above is used essentially for performing
a range-max query over a certain set of phrase weights. There is another
method[1] of performing range-max queries which has a runtime cost of $O(1)$
rather than $O(\log{n})$ per query. If this method is used, we can reduce the
runtime cost of our method further from $O(k\log{n})$ to $O(k\log{k})$.

If we use a suffix array to represent all the phrases in our corpus, we can
perform a match of the query with a prefix of every suffix of each phrase in
the corpus. This makes our searches even more robus without significantly
increasing our pre-processing or query time. However, we would now need to
store an entry for every starting position of a string in the suffix array in
in our RQM data structure. This would increase the size of the Segment Tree to
$O(nk)$, where $k$ is the average length of a phrase in our corpus. This will
also have an effect on the query time, taking it up to $O(\log{nk})$ from
$O(\log{n})$.

## References

  * [1] Michael A.Bender and Martin Farach-Colton, _The lca problem revisited_ , Springer-Verlag (2000), 88–94.
  * [2] Antirez, _Redis autocomplete_ , http://antirez.com/post/autocomplete-with-redis.html.
  * [3] Shengyue Ji, Shengyue Ji, Chen Li, and Jianhua Feng, _Efficient interactive fuzzy keyword search_ , ACM (2009).
  * [4] Guoliang Li, Shengyue Ji, Chen Li, and Jianhua Feng, _Efficient type-ahead search on relational data: a tastier approach_ , SIGMOD (2009).
  * [5] W. Naylor and H. Opderbeck, _Rfc619_ , 1974, http://www.faqs.org/rfcs/rfc619.html.
  * [6] _Segment tree_ , http://en.wikipedia.org/wiki/Segment_tree.
  * [7] Phillipa Sessini and Anirban Mahanti, _Observations on round-trip times of tcp connections_.
  * [8] Petr Sojka, _Notes on compound word hyphenation in tex_ , TUGboat 16 (1995), no. 3, 290–296.
  * [9] _Words per minute_ , http://en.wikipedia.org/wiki/Words_per_minute#Alphanumeric_entry.


* (126) C. Padurariu, T. Jonckheere, J. Rech, T. Martin, and D. Feinberg, Tunable pseudogaps due to nonlocal coherent transport in voltage-biased three-terminal Josephson junctions, Phys. Rev. B 95, 205437 (2017).
  * (127) J.-D. Pillet, C. Quay, P. Morfin, C. Bena, A. Levy Yeyati, and P. Joyez, Revealing the electronic structure of a carbon nanotube carrying a supercurrent, Nature Physics 6, 965 (2010).
  * (128) L. Bretheau, J.I.-J. Wang, R. Pisoni, K. Watanabe, T. Taniguchi and P. Jarillo-Herrero, Tunnelling spectroscopy of Andreev states in graphene, Nat. Phys. 13, 756 (2017).
  * (129) J. C. Cuevas, J. Hammer, J. Kopu, J. K. Viljas, and M. Eschrig, Proximity effect and multiple Andreev reflections in diffusive superconductor–normal-metal–superconductor junctions, Phys. Rev. B 73, 184505 (2006).
  * (130) M. Vanević and W. Belzig, Full counting statistics of Andreev scattering in an asymmetric chaotic cavity, Phys. Rev. B 72, 134522 (2005).
  * (131) G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion, Phys. Rev. B 25, 4515 (1982).


# Information flow-enhanced precision in collisional quantum thermometry

Taysa M. Mendonça<EMAIL_ADDRESS>Instituto de Física de São Carlos,
Universidade de São Paulo, CP 369, 13560-970, São Carlos, SP, Brazil  Centre
for Quantum Materials and Technologies, School of Mathematics and Physics,
Queen’s University, Belfast BT7 1NN, United Kingdom  Diogo O. Soares-Pinto
<EMAIL_ADDRESS>Instituto de Física de São Carlos, Universidade de São Paulo,
CP 369, 13560-970, São Carlos, SP, Brazil  Mauro Paternostro
<EMAIL_ADDRESS>Università degli Studi di Palermo, Dipartimento di
Fisica e Chimica - Emilio Segrè, via Archirafi 36, I-90123 Palermo, Italy
Centre for Quantum Materials and Technologies, School of Mathematics and
Physics, Queen’s University, Belfast BT7 1NN, United Kingdom

###### Abstract

We describe and analyze a quantum thermometer based on a multi-layered
collisional model. The proposed architecture provides significant sensitivity
even for short interaction times between the ancillae comprised in the
thermometer and the system to be probed, and a small number of information-
acquiring collisions. The assessment of the flow of information taking place
within the layered thermometer and between system and thermometer reveals that
the tuning of the mutual backflow of information has a positive influence on
the precision of thermometry, and helps unveiling the information-theoretic
mechanisms behind the working principles of the proposed architecture.

###### pacs:

Pacs

##  I Introduction

Quantum metrology uses resources such as entanglement [1, 2, 3], quantum
coherence [4, 5, 6] and squeezing [7] to enhance the precision with which
measurement protocols are able to estimate relevant physical quantities beyond
the capabilities of classical strategies. When applied to the problem of
estimating the temperature of a system, such quantum tools make quantum
thermometry [8] effective even when operating under conditions that are
classically challenging, such as when the thermometer has not yet thermalized
with the system. This makes such framework potentially useful to probe non-
equilibrium environments, which could showcase coherence and correlations.

Energetic resource value of information is an active research area in quantum
thermodynamics, they can power up, cool, or heat systems using information
instead of consuming energy. Quantum thermometry will be an integral part of
practical implementations of such quantum information engines and devices. A
particularly powerful way to model the process of estimating temperature
leverages upon the framework of collisional models [9, 10], which has been
successfully applied to describe the thermodynamics of non-equilibrium quantum
processes [11, 12, 13, 14, 15, 16].

In collisional quantum thermometry, a set of probing particles acquire
information on the temperature of the system of interest through a sequence of
short and intense, crash-like interactions. Ref. [17] explored such approach
showing how it allows for the estimation of the temperature of a quantum
system beyond the precision constraints set by thermal fluctuations [18].
Additional studies have grounded the suitability of collisional quantum
thermometry to pursue this important metrological task [15]. In this paper, we
use a collision-based approach to build a thermometric sensor that measures
the temperature of an environment indirectly, in a completely out-of-
equilibrium regime, by resorting to measurements performed on auxiliary qubits
arranged in the form of layers of chains. We go beyond state-of-the-art in
collisional quantum thermometry by alllowing explicitly for inter-layer
interactions: although this obviously increases the complexity of the problem
at hand, it also enhances the process of information-spreading that is
responsible for the performance of this metrological task. Quantitatively, we
achieve a thermometric sensitivity that is orders of magnitude larger than
that of a single-layered thermometer as proposed by [17]. We also investigated
the case where there is a complete exchange of information during the system-
ancilla interaction, a model proposed by [19]. In this case, we showed that
the same amount of information flows between the qubits, not allowing for a
gain in sensitivity.

To investigate the motivating factors for such gain in sensitivity of our
multi-layered thermometer, we studied the behavior of the information flow
using the measure of non-Markovianity proposed in Ref. [20, 21] (dubbed here
as BLP), this procedure was demonstrated in Ref.[22]. Finally, we compare the
QFI and the BLP measurement, qubit by qubit, to investigate the path taken by
information in our complete system. It was possible to determine how the
accumulation of information provides an increase in the sensitivity of our
thermometer, we also show that such accumulation can also be characterized
through Mutual Information.

##  II Description of the Model

We consider a collisional model where a qubit – embodying the system – is
initially thermalized with an environment at temperature $T$. A set of
ancillary chains is used to indirectly estimate the temperature of the
environment, as showed in the Fig. 1(a). The dynamics of the overall compound
comprises a sequence of collisions between the system and $N$ chains of $n$
ancillae each. In a schematic manner, the process unfolds according to the
following steps:
($1$) The system $S$ thermalizes with the environment through the thermal map
$\mathcal{T}_{S}=\mathrm{e}^{\mathcal{L}\tau_{SE}}$ generated by [23]

$\frac{d\rho_{S}}{dt}=\mathcal{L}(\rho_{S})=\gamma\left(\bar{n}+1\right)\mathcal{D}\left[\mathcal{S}_{-}\right]+\gamma\left(\bar{n}\right)\mathcal{D}\left[\mathcal{S}_{+}\right],$
(1)

where
$\mathcal{D}\left[L\right]=L\rho_{S}L^{{\dagger}}-\frac{1}{2}\left\\{L^{{\dagger}}L,\rho_{S}\right\\}$,
$\bar{n}=(\mathrm{e}^{\beta\Omega_{S}}-1)^{-1}$ with
$\beta=\hbar\Omega_{S}/k_{B}$ and $k_{B}$ the Boltzmann constant. Here,
$\Omega_{S}$ is the characteristic frequency of the system, $T$ is the
temperature of the environment and $\tau_{SE}$ is a time long enough for the
system qubit to be thermalized with the temperature of the environment. At the
end of this step, the qubit system will be in the Gibbs state
$\rho_{S}^{th}=\mathrm{e}^{-H_{S}/k_{B}T}/Z$, where
$H_{S}=\hbar\Omega_{S}\mathcal{S}_{z}/2$ is the Hamiltonian of the system.

The following steps involving $A_{k,j}$, i.e. the $j^{\text{th}}$ element of
the $k^{\text{th}}$ chain ($j=1,..,n$ and $k=1,..,N$), then take place: For a
set value of $j=1,..,n$, a collision between the system $S$ and $A_{{k,j}}$,
ruled by the Hamiltonian

$H_{SA_{k,j}}=\hbar
g\left(\mathcal{S}_{+}\alpha^{k,j}_{-}+\mathcal{S}_{-}\alpha^{k,j}_{+}\right)$
(2)

occurs, followed by the collision between $A_{k,j}$ and $A_{k+1,j}$ according
to the generator

$H_{A_{k,j}A_{k+1,j}}=\hbar
J\left(\alpha^{k,j}_{+}\alpha^{k+1,j}_{-}+\alpha^{k,j}_{-}\alpha^{k+1,j}_{+}\right).$
(3)

Here, $g$ ($J$) is the system-ancilla (inter-ancilla) coupling strength, which
we assume to be independent of $j$ and $k$, while the operators
$\mathcal{S}_{\pm}=(\ket{0}\pm\ket{1})_{S}/\sqrt{2}$ and
$\alpha_{\pm}^{k,j}=(\ket{0}\pm\ket{1})_{A_{k,j}}/\sqrt{2}$ refer to system
and the relevant ancilla, respectively, with $\\{\ket{0},\ket{1}\\}$ being the
computational basis of our problem (we have shifted the energy of the states
in such a way that $\mathcal{O}_{z}\ket{0}=\ket{0}$ and
$\mathcal{O}_{z}\ket{1}=-\ket{1}$).

The collisions continue so as to span all the chains (i.e. for $k=1,...,N$),
after which a thermalization process occurs again. The sequence of system-
ancilla and ancilla-ancilla collisions starts again, this time involving the
elements $j+1$ of each chain. then consider the After that, the process
restarts for the $j+1$ ancillae spanning all the $N$ chains. Fig. 1 provides a
visual depiction of the first few steps of the process, as well as its circuit
representation representing the interactions for the case of three chains
only.

The sequence of operations that describes the described interaction protocol
for a single repetition and considering a generic ancilla leads to the state

$\rho_{S,A_{1,j}\cdots
A_{N,j}}{=}\mathcal{U}_{S,A_{k,j}}\circ\left[\bigcirc^{N-1}_{k=1}\mathcal{U}_{A_{k+1,j}A_{k,j}}\right]\circ\mathcal{U}_{S,A_{k,j}}\circ\mathcal{T}_{S}(\rho_{0}),$
(4)

where
$\mathcal{U}_{S,A_{k,j}}\circ\mu=V_{S,A_{k,j}}\,\mu\,V_{S,A_{k,j}}^{\dagger}$
and
$\mathcal{U}_{A_{k,j},A_{k+1,j}}\circ\mu=V_{A_{k,j},A_{k+1,j}}\,\mu\,V_{A_{k,j},A_{k+1,j}}^{\dagger}$
for a density matrix $\mu$ and
$V_{S,A_{k,j}}=\text{exp}(-iH_{SA_{k,j}}\tau_{SA}/\hbar)$,
$V_{A_{k,j},A_{k+1,j}}=\text{exp}(-iH_{A_{k,j}A_{k+1,j}}\tau_{A}/\hbar)$. In
our analysis, we will consider values of the dimensionless system-ancilla
interaction time $g\tau_{SA}\ll{\pi/2}$, while taking the inter-ancilla one as
$J\tau_{A}=\pi/2$. Under such conditions, a full exchange of information
between the ancillae will be allowed, as $V_{A_{k,j}A_{k+1,j}}$ is the phased
SWAP gate performing the transformation
$\ket{a,a}_{A_{k,j}A_{k+1,j}}\to(-i)^{a\oplus b}\ket{b,a}_{A_{k,j}A_{k+1,j}}$
($a,b=0,1$ and $\oplus$ standing for sum modulo 2). On the other hand, the
information exchange between the system and ancilla will only be partial.

Figure 1: (a) Illustration of the scheme of interctions for multi-layered
collisional thermometry. System $S$ thermalizes with an environment in
equilibrium at temperature $T$. The inference of such parameter is made
through the sequence of collisional interactions illustrated in the main text.
The numbered steps of the scheme [from (1) to (6) in panel (a)] illustrate the
first few steps of the sequence leading to the measurements entailed by our
metrological approach to the estimate of $T$. (b) Circuit-model representation
of interactions of the steps (1)-to-(6) of our protocol. Each inter-ancilla
collision is depicted by a full phased-SWAP gate [colored in blue], while the
green-colored gates stand for is used for partial SWAPs.

##  III Analysis of the quantum Fisher Information

We make use of the formal apparatus of quantum metrology to provide a bound to
the precision of the estimation of the temperature of the environment, which
has been encoded in the state of system. In particular, using the framework of
quantum parameter estimation [24], we can upper bound to said precision
through the quantum Fisher Information (QFI)

$QFI_{k,j}(T,\rho_{A})=\underset{\Pi_{\theta,\phi}}{\text{max}}\,F_{k,j}(\Pi_{\theta,\phi},\rho_{A_{k,j}},T),$
(5)

where
$\Pi_{\theta,\phi}=\ket{\psi(\theta,\phi)}\bra{\psi(\theta,\phi)}_{A_{k,j}}$
with
$\ket{\psi(\theta,\phi)}_{A_{k,j}}=\cos\theta\ket{0}_{A_{k,j}}+e^{i\phi}\sin\theta\ket{1}_{A_{k,j}}$
($\theta\in[0,\pi]$, $\phi\in[0,2\pi]$) are the elements of a positive
operator valued measurement (POVM) on ancilla $A_{k,j}$, and
$F_{k,j}(\Pi_{\theta,\phi},\rho_{A_{k,j}},T)$ is the classical Fisher
Information of the POVM

$F_{k,j}(\Pi_{\theta,\phi},T,\rho_{A})=\iint\,d\theta\,d\phi\,p(\theta,\phi)\left(\dfrac{\partial}{\partial
T}\ln{p(\theta,\phi)}\right)^{2}.$ (6)

Here, $p(\theta,\phi)=\text{Tr}(\Pi_{\theta,\phi}\rho_{A_{k,j}})$ is the
probability to get outcome $(\theta,\phi)$ from the measurement. The QFI
implies the maximization of its classical counterpart over all possible POVMs,
and is thus optimized over the measurement strategy on the ancillae. The
precision of a probe in thermal equilibrium with a system is limited by the
Cramer-Rao bound

$\left(\Delta T\right)^{2}\geq\frac{1}{QFI(T,\rho)},$ (7)

where $(\Delta T)^{2}$ is the temperature variance. When evaluated over a
thermal state, such lower bound to the variance of a temperature estimator
(sometimes referred to as thermal Fisher information [19]) can be surpassed
through the use of quantum resources such as entanglement and coherence.

The results of the QFI measurements for $k=1,..,8$ and $j=1,..,30$ are shown
in Fig. 2(a), which shows the value of the QFI for measurements performed on
the ancillae of each chain as we vary the temperature. Specifically, we have
assumed to measure the last ancilla of each chain. As the number of layered
chains increases, the QFI grows thus lowering the upper bound entailed by the
Cramér-Rao bound and thus witnessing and increase in sensitivity. Fig. 2(b)
shows the behavior of the QFI as we vary the number of ancillae in each chain,
i.e. varying $j$, showing that, in a sufficiently layered structure, even a
few ancillae will be sufficient for our collisional thermometer to achieve
high sensitivity. Fig. 2(c) shows the maximum value of the QFI growing with
the number of layers.

Figure 2: (a) QFI versus temperature and (b) QFI versus number of ancilla on
each chain in the optimized temperature and (c) QFI Maximum versus number of
ancilla chains.

##  IV Assessing information flow and its role in the thermometric task

To understand the reason for the increased thermometric sensitivity, we will
analyze the flow of information during the various operations carried out
between system and ancillae.

We consider the system as prepared as proposed in [22]. The effect of the
system-ancilla interaction can be quantified, at every instant of time, by
considering the distinguishability between the the time-evolved form of such
initial preparations as given by the quantity

$D\left(\rho_{S_{+}A_{k,j}},\rho_{S_{-}A_{k,j}};t\right)=\|\rho_{S_{+}A_{k,j}}(t)-\rho_{S_{-}A_{k,j}}(t)\|_{1}$
(8)

with $\|\cdot\|_{1}$ the trace distance. The lack of monotonicity of
$D\left(\rho_{S_{+}A_{k,j}},\rho_{S_{-}A_{k,j}};t\right)$ signals a back-flow
of information into the system stemming from the joint system-ancillae
dynamics. Following Refs. [25, 26, 27, 28], this witnesses the non-Markovian
nature of the open dynamics undergone by $S$.

Fig. 3 reports the behavior of Eq. (8) and its derivative
$\sigma(t)=\partial_{t}D\left(\rho_{S_{+}(A_{k,j})},\rho_{S_{-}(A_{k,j})};t\right)$
when $S$ interacts with only the first ancilla of each chain. Fig. 3(a) and
(c) refer to the Markovianity in the system. Needless to say, a complementary
picture can be taken where we characterize non-Markovianity by looking at the
distinguishability between evolved states of the ancillae. The results
stemming from such a perspective are shown in Fig. 3(b) and (d). We observe a
flow of information between system and ancillae, growing with the number of
chains being considered. Moreover, there is clearly information is flowing off
the system [as $\sigma<0$ in Fig. 3(c)], and into the ancillae [whose
corresponding $\sigma$ is positive, as seen from Fig. 3(d)].

Figure 3: Measure of non-Markovianity. (a) and (b) are respectively the
results for the trace distance of the system and ancillae. (c) and (d) are the
results for the derivative of the trace distance (measure of information flow)
of the system and ancillae, respectively.

In order to gather insight in the dynamics of information flow during the
thermometric task, in the following numerical experiments we will consider
$g\tau_{SA}=J\tau_{A}=\pi/2$, i.e., a full exchange of information between
system-ancilla and ancilla-ancilla. Measurements were performed on all parties
involved int he processes ruled by $V_{S,A_{k,j}}$ and
$V_{A_{k,j},A_{k+1,j}}$. However, for a clearer visualization we plot only the
results for $j=1$ and $k=1,2$, i.e. $V_{S,A_{1,1}}$ (measuring $S$ and
$A_{1,1}$), $V_{A_{1,1},A_{2,1}}$ (measuring $A_{1,1}$ and $A_{2,1}$) and
$V_{S,A_{2,1}}$ (measuring $S$ and $A_{2,1}$). These are the operations
illustrated in steps $2$, $3$ and $4$ of Fig. 1(a). These analyses and results
does not depend on the choice of $j$ made here.

We show the behavior of the trace distance, its derivative $\sigma$ and the
QFI in Fig. 4 and Fig. 5, respectively. We use the protocol put forward in
Ref. [22], where such tools have been used to analyze the flow of information
passing through a qubit environment. We take the initial states of the system
to be $\rho_{S}^{(\pm)}(0)=\ket{\pm}\bra{\pm}_{S}$ while the initial state of
the ancillae is $\rho_{A_{k,j}}(0)=\bigotimes^{N}_{k=1}\ket{0}\bra{0}_{k}$.
The initial state of the system-ancillae compound is the uncorrelated one
$\rho_{SA}^{(\pm)}(0)=\rho_{S}^{(\pm)}(0)\otimes\rho_{A_{k,j}}(0)$.

Fig. 4(a) and Fig. 4(b) show the dynamics of the trace distance and its
derivative, respectively, for the states described above. Fig. 4(b) shows
clearly how each collision (either system-ancilla or ancilla-ancilla) results
in a pouring of information that flows according to the preferential direction

$S\to\text{chain {\hbox{k}}}\to\text{chain {\hbox{k+1}}}\to S$

for $k=1,...,N$ since the positions of the ancillae such that $k$ is odd and
$k+1$ is even. Here we can also see that information leaves a given qubit [as
witnessed by having $\sigma<0$, cf. solid line in Fig. 4(b)], and reaches the
other qubit involved in the operation [resulting in $\sigma>0$, cf. dot-
dashed, dotted, and dashed lines].

Figure 4: (a) Trace distance and (b) Information flow between the qubits
involved in interactions $2$, $3$ and $4$ in Fig. 1, the measurements are
performed for an complete exchange of information between system and ancillae
($g\tau_{SA}=\pi/2$). Solid lines show information leaving one qubit as it
arrives at another qubit, shown in dot-dashed, dashed and dotted lines. Figure
5: QFI performed on qubits involved in interactions $2$, $3$ and $4$ in Fig.
1. We consider a full exchange of information between system and ancillae
($g\tau_{SA}=\pi/2$). Measurements are performed on qubits when information is
(a) leaving and (b) entering.

We now evaluate the QFI using the same protocol described above. Figs. 5(a)
and (b) show the QFI resulting from measurements performed on the qubits from
which the information is coming out and coming in, respectively. Note that the
QFI is only defined when the information is coming in to the qubit. Fig. 5(b)
shows that the same QFI can be extracted from measurements on $S$, $A_{k,1}$
or $A_{k+1,1}$ (since $k$ is odd and $k+1$ is even), a result that is clearly
stemming from the highly symmetric interaction configuration that we have
chosen for this illustration. We see that the QFI measurement does not vary
when $g\tau_{SA}=J\tau_{A}=\pi/2$ there is no gain in sensor sensitivity, this
happens because the same amount of information is being transferred in the
full Swap operation [cf. Fig. 4(b)]. If the measured qubit is sending the
information during the full exchange interaction, the QFI measurement will not
have a defined result (Fig. 5(a)).

We now extend the analysis and return to the case of $g\tau_{SA}\ll\pi/2$,
i.e. an incomplete exchange of information between system and ancillae. In
Fig. 6 we consider the results corresponding to $j=1$ and a compound
comprising 5 chains. The QFI for the system and the ancillae are shown in
Figs. 6(a) and (b), respectively. Fig. 6(c) shows the QFI for the ancilla that
receives information during the ancilla-ancilla interaction. Since the
interactions characterizes an full Swap operation, i.e. $J\tau_{A}=\pi/2$, the
QFI measurements performed on the ancillae from which the information is
leaving the qubit follow the result shown in Fig. 5(a), thus we will not
repeat such results here.

As in Fig. 5(b), Fig. 6(a) shows that the same amount of information can be
measured in the different parts of the compound. However, Figs. 6(b) and (c)
show that there is an increase in the sensitivity of our thermometer as we
increase the number of ancillary chains. Therefore, this represents there is
an accumulation of information caused by the short interaction time between
system and ancillae $\tau_{SA}$, as is evident in the characterization of the
distinguishability in Fig. 3(d). Since we can only obtain the QFI when the
information flows into an element, we can state that there is also information
returning to the system.

Figure 6: QFI resulting from an incomplete exchange of information between
system and ancillae ($g\tau_{SA}\ll\pi/2$). The measurements needed to perform
the estimation are performed on (a) the system (during the system-ancilla
interactions); (b) the ancillae (during the system-ancilla interactions); (c)
the ancillae (during the ancilla-ancilla interactions).

The accumulation of information can also be seen when we calculate the mutual
information between $S$ and an ancillary chain. We compute the mutual
information between the elements with $j=n$ of the $k^{\text{th}}$ chain in
all interactions as

$I\left(\rho_{S}:\rho_{A_{k,n}}\right)=\mathcal{E}(\rho_{S})+\mathcal{E}(\rho_{A_{k,n}})-\mathcal{E}(\rho_{SA_{k,n}}).$
(9)

where $\mathcal{E}(\rho)=-\text{tr}(\rho\text{log}\rho)$ is the von Neumann
entropy [29]. Fig. 7 shows that as we increase the number of ancilla chains,
the greater the amount of information shared between the system and the
ancillae.

Figure 7: Mutual information between the system and the qubits of the ancilla
chains during operations $V_{S,A_{k,n}}$.

Conclusion.– We have built a model to improve and increase the sensitivity of
a quantum thermometer based on the use of a multilayered system. Even with a
short interaction time, i.e. without requesting thermalization of the system
with the thermometer, we can adjust the sensitivity of the latter by
increasing the number of its ancillary chains. We have unveiled a tight link
between the sensitivity of the estimation of temperature as provided by the
QFI and the information flow among various parts of the compound. The
accumulation of information entailed by the information flow leads to better
sensitivity (i.e. a larger QFI). This can also be seen in the increase of
mutual information with the number of chains comprised in the thermometer.

## Acknowledgments

The authors thank discussions with G. T. Landi. TMM thanks the grants No.
2021/01277-2 and No. 2022/09219-4, São Paulo Research Foundation (FAPESP).
DOSP acknowledges the support by the Brazilian funding agencies CNPq (Grant
No. 304891/2022-3), FAPESP (Grant No. 2017/03727-0), and the Brazilian
National Institute of Science and Technology of Quantum Information (INCT/IQ).
MP acknowledges the support by the European Union’s Horizon Europe EIC
Pathfinder project QuCoM (Grant Agreement No. 101046973), the Royal Society
Wolfson Fellowship (RSWF/R3/183013), the UK EPSRC (EP/T028424/1), and the
Department for the Economy Northern Ireland under the US-Ireland R&D
Partnership Programme.

## References

  * Giovannetti  _et al._ [2006] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum metrology, Phys. Rev. Lett. 96, 010401 (2006).
  * Maccone [2013] L. Maccone, Intuitive reason for the usefulness of entanglement in quantum metrology, Phys. Rev. A 88, 042109 (2013).
  * Huang  _et al._ [2016] Z. Huang, C. Macchiavello, and L. Maccone, Usefulness of entanglement-assisted quantum metrology, Phys. Rev. A 94, 012101 (2016).
  * Micadei  _et al._ [2015] K. Micadei, D. A. Rowlands, F. A. Pollock, L. C. Céleri, R. M. Serra, and K. Modi, Coherent measurements in quantum metrology, New Journal of Physics 17, 023057 (2015).
  * Pires  _et al._ [2018] D. P. Pires, I. A. Silva, E. R. deAzevedo, D. O. Soares-Pinto, and J. G. Filgueiras, Coherence orders, decoherence, and quantum metrology, Phys. Rev. A 98, 032101 (2018).
  * Castellini  _et al._ [2019] A. Castellini, R. Lo Franco, L. Lami, A. Winter, G. Adesso, and G. Compagno, Indistinguishability-enabled coherence for quantum metrology, Phys. Rev. A 100, 012308 (2019).
  * Maccone and Riccardi [2020] L. Maccone and A. Riccardi, Squeezing metrology: a unified framework, Quantum 4, 292 (2020).
  * De Pasquale and Stace [2018] A. De Pasquale and T. M. Stace, Quantum thermometry, in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 503–527.
  * Ciccarello  _et al._ [2022] F. Ciccarello, S. Lorenzo, V. Giovannetti, and G. M. Palma, Quantum collision models: Open system dynamics from repeated interactions, Physics Reports 954, 1 (2022).
  * Ciccarello [2017] F. Ciccarello, Collision models in quantum optics, Quantum Measurements and Quantum Metrology 4, 53 (2017).
  * Rodrigues  _et al._ [2019] F. L. S. Rodrigues, G. De Chiara, M. Paternostro, and G. T. Landi, Thermodynamics of weakly coherent collisional models, Phys. Rev. Lett. 123, 140601 (2019).
  * Filippov [2022] S. Filippov, Multipartite correlations in quantum collision models, Entropy 24, 10.3390/e24040508 (2022).
  * Landi [2021] G. T. Landi, Battery charging in collision models with bayesian risk strategies, Entropy 23, 10.3390/e23121627 (2021).
  * McCloskey and Paternostro [2014] R. McCloskey and M. Paternostro, Non-markovianity and system-environment correlations in a microscopic collision model, Phys. Rev. A 89, 052120 (2014).
  * O’Connor  _et al._ [2021] E. O’Connor, B. Vacchini, and S. Campbell, Stochastic collisional quantum thermometry, Entropy 23, 10.3390/e23121634 (2021).
  * Bouton  _et al._ [2020] Q. Bouton, J. Nettersheim, S. Burgardt, D. Adam, E. Lutz, and A. Widera, An endoreversible quantum heat engine driven by atomic collisions (2020), arXiv:2009.10946 [quant-ph] .
  * Seah  _et al._ [2019] S. Seah, S. Nimmrichter, D. Grimmer, J. P. Santos, V. Scarani, and G. T. Landi, Collisional quantum thermometry, Phys. Rev. Lett. 123, 180602 (2019).
  * Paris [2009a] M. G. A. Paris, Quantum estimation for quantum technology, Int. J. Quantum. Inform. 7, 125 (2009a).
  * Shu  _et al._ [2020] A. Shu, S. Seah, and V. Scarani, Surpassing the thermal Cramér-Rao bound with collisional thermometry, Phys. Rev. A 102, 042417 (2020).
  * Breuer  _et al._ [2009] H.-P. Breuer, E.-M. Laine, and J. Piilo, Measure for the degree of non-markovian behavior of quantum processes in open systems, Phys. Rev. Lett. 103, 210401 (2009).
  * Laine  _et al._ [2010] E.-M. Laine, J. Piilo, and H.-P. Breuer, Measure for the non-markovianity of quantum processes, Phys. Rev. A 81, 062115 (2010).
  * Mendonça  _et al._ [2024] T. M. Mendonça, L. C. Céleri, M. Paternostro, and D. O. Soares-Pinto, System-environment quantum information flow (2024), arXiv:2402.15483 [quant-ph] .
  * Breuer and Petruccione [2002] H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, New York, 2002).
  * Paris [2009b] M. G. A. Paris, Quantum estimation for quantum technology, Int. J. Quantum Inf 7 (2009b).
  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Colloquium: Non-markovian dynamics in open quantum systems, Rev. Mod. Phys. 88, 021002 (2016).
  * Rivas  _et al._ [2014] A. Rivas, S. F. Huelga, and M. B. Plenio, Quantum non-markovianity: characterization, quantification and detection, Reports on Progress in Physics 77, 094001 (2014).
  * Li  _et al._ [2018] L. Li, M. J. Hall, and H. M. Wiseman, Concepts of quantum non-markovianity: A hierarchy, Physics Reports 759, 1 (2018), concepts of quantum non-Markovianity: A hierarchy.
  * Chrúciński [2022] D. Chrúciński, Dynamical maps beyond markovian regime, Physics Reports 992, 1 (2022).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ , 10th ed. (Cambridge University Press, New York, 2010).


# A simple numerical method for Hele-Shaw type problems by the method of
fundamental solutions

Koya Sakakibara  Department of Applied Mathematics, Faculty of Science,
Okayama University of Science, 1-1 Ridaicho, Kita-ku, Okayama-shi, Okayama
700-0005, Japan; RIKEN iTHEMS, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan
<EMAIL_ADDRESS>,  Yusauku Shimoji  Graduate School of Science and
Technology, Meiji University, 1-1-1 Higashi-Mita, Tama-ku, Kawasaki-shi,
Kanagawa 214-8571, Japan<EMAIL_ADDRESS>and  Shigetoshi Yazaki
Department of Mathematics, School of Science and Technology, Meiji University,
1-1-1 Higashi-Mita, Tama-ku, Kawasaki-shi, Kanagawa 214-8571, Japan
<EMAIL_ADDRESS>
###### Abstract.

Hele-Shaw flows with time-dependent gaps create fingering patterns, and
magnetic fluids in Hele-Shaw cells create intriguing patterns. We propose a
simple numerical method for Hele-Shaw type problems by the method of
fundamental solutions. The method of fundamental solutions is one of the mesh-
free numerical solvers for potential problems, which provides a highly
accurate approximate solution despite its simplicity. Moreover, combining with
the asymptotic uniform distribution method, the numerical method satisfies the
volume-preserving property. We use Amano’s method to arrange the singular
points in the method of fundamental solutions. We show several numerical
results to exemplify the effectiveness of our numerical scheme.

###### Key words and phrases:

Hele-Shaw flow; magnetic fluid; time-dependent gap; the method of fundamental
solutions; Amano’s method; volume-preserving property

###### 2010 Mathematics Subject Classification:

76D27, 76E25, 65N80, 65N35, 35J05

##  1\. Introduction

The Hele-Shaw problem describes the motion of a viscous fluid in a quasi-two-
dimensional space, originating in a short paper by Henry Selby Hele-Shaw [8],
and was presented as a kind of experimental model for describing stream lines
by using viscous fluids. Let $\Omega(t)\subset\mathbb{R}^{2}$ be a bounded
region occupied by a fluid, and $\Gamma(t)=\partial\Omega(t)$ be its boundary.
Then, the Hele-Shaw problem is described as follows:

(1) $\displaystyle\begin{dcases*}\triangle p(\cdot,t)=0&in $\Omega(t)$,
$t\in(0,T)$,\\\ p(\cdot,t)=\sigma\kappa(\cdot,t)&on $\Gamma(t)$,
$t\in(0,T)$,\\\ V(\cdot,t)=-h^{2}\nabla p(\cdot,t)\cdot\bm{N}(\cdot,t)&on
$\Gamma(t)$, $t\in(0,T)$,\end{dcases*}$

where $p(\cdot,t)$ is the pressure in $\Omega(t)$, $\sigma$ is the surface
tension coefficient, $\kappa(\cdot,t)$ denotes the curvature, $V(\cdot,t)$ is
the normal velocity, $\bm{N}(\cdot,t)$ is the unit outward normal vector of
$\Gamma(t)$, and $h$ denotes the gap of the Hele-Shaw cell. The Laplace
equation for pressure $p$, the first equation in (1), comes from the Darcy’s
law

(2) $\displaystyle\bm{u}=-h^{2}\nabla p,$

and the incompressibility condition

(3) $\displaystyle\nabla\cdot\bm{u}=0,$

where $\bm{u}$ denotes the velocity field. The Dirichlet boundary condition,
the second equation in (1), expresses the surface tension. For details of
derivation, see, for example, the books [14, 6].

Figure 1. Hele-Shaw cell

There are many variants of the Hele-Shaw problem, and in particular, they have
been actively used for the experimental and mathematical study of fingering
phenomena in various scientific fields [20, 10, 9, 27, 28, 4]. Among them,
Shelley et al. formulated the Hele-Shaw problem with a time-dependent gap in
[27] (see also Fig.1):

(4) $\displaystyle\begin{dcases*}\triangle
p(\cdot,t)=\frac{\dot{h}(t)}{h(t)^{3}}&in $\Omega(t)$, $t\in(0,T)$,\\\
p(\cdot,t)=\sigma\kappa(\cdot,t)&on $\Gamma(t)$, $t\in(0,T)$,\\\
V(\cdot,t)=-h(t)^{2}\nabla p(\cdot,t)\cdot\bm{N}(\cdot,t)&on $\Gamma(t)$,
$t\in(0,T)$,\end{dcases*}$

where $h(t)$ denotes the gap of the Hele-Shaw cell depending on time $t$. In
the current situation, where the gap is time-dependent, the Darcy’s law and
the incompressibility condition are respectively given by

(5) $\displaystyle\bm{u}=-h(t)^{2}\nabla
p,\quad\nabla\cdot\bm{u}=-\frac{\dot{h}(t)}{h(t)}.$

Therefore, the equation to be satisfied by the pressure is now the Poisson
equation, not the Laplace equation. However, since a particular solution of
the Poisson equation can be obtained concretely as

(6) $\displaystyle p^{*}(\bm{x},t)=\frac{\dot{h}(t)}{4h(t)^{3}}|\bm{x}|^{2},$

by considering $p-p^{*}$ as $p$ again, (4) eventually leads to the following
system:

(7) $\displaystyle\begin{dcases*}\triangle p(\cdot,t)=0&in $\Omega(t)$,
$t\in(0,T)$\\\
p(\bm{x},t)=\sigma\kappa(\cdot,t)-\frac{\dot{h}(t)}{4h(t)^{3}}|\bm{x}|^{2}&for
$\bm{x}\in\Gamma(t)$, $t\in(0,T)$,\\\ V(\bm{x},t)=-h(t)^{2}\left(\nabla
p(\bm{x},t)+\frac{\dot{h}(t)}{2h(t)^{3}}\bm{x}\right)\cdot\bm{N}(\bm{x},t)&for
$\bm{x}\in\Gamma(t)$, $t\in(0,T)$.\end{dcases*}$

As a more advanced topic, fingering phenomena in magnetic fluids have also
been investigated [17, 19, 5, 29]. The following is a summary of the
mathematical model proposed in [29]. Darcy’s law and the incompressibility
condition are respectively given by

(8) $\displaystyle-\nabla
p-\frac{12\eta}{h(t)^{2}}\bm{u}+\frac{2M_{c}}{h(t)}\nabla\varphi_{m}=\bm{0},\quad\nabla\cdot\bm{u}=-\frac{\dot{h}(t)}{h(t)},$

where $\varphi_{m}=\varphi_{m}(\bm{x},t)$ is the magnetostatic field potential
given by

(9)
$\displaystyle\varphi_{m}(\bm{x},t)=-M_{c}\int_{\Omega(t)}\left(\frac{1}{|\bm{x}-\bm{p}|}-\frac{1}{\sqrt{|\bm{x}-\bm{p}|^{2}+h(t)^{2}}}\right)\,\mathrm{d}\bm{p},$

in which $M_{c}$ denotes the strength of the magnetization. In other words, we
obtain the following system:

(10) $\displaystyle\begin{dcases*}\triangle
p(\cdot,t)=\frac{12\eta\dot{h}(t)}{h(t)^{3}}+\frac{2M_{c}}{h(t)}\triangle\varphi_{m}(\cdot,t)&in
$\Omega(t)$, $t\in(0,T)$,\\\ p(\cdot,t)=\sigma\kappa(\cdot,t)&on $\Gamma(t)$,
$t\in(0,T)$,\\\
V(\cdot,t)=-\frac{h(t)^{2}}{12\eta}\nabla\left(p-\frac{2M_{c}}{h}\varphi_{m}\right)(\cdot,t)\cdot\bm{N}(\cdot,t)&on
$\Gamma(t)$, $t\in(0,T)$.\end{dcases*}$

By constructing a particular solution of the Poisson equation in the same way
as we derived (7), and by performing the appropriate non-dimensionalization,
we finally obtain the following model:

(16) $\displaystyle\left\\{\begin{array}[]{l}\triangle p(\cdot,t)=0\hskip
160.0pt\text{in}\ \Omega(t),\ t\in(0,T),\\\ \displaystyle
p(\cdot,t)=\kappa(\cdot,t)\left(\frac{R_{0}}{h_{0}}\right)^{1/3}\frac{\pi^{2/3}\mathrm{Bm}}{h_{*}(t)}\varphi(\cdot,t)-\left(\frac{R_{0}}{h_{0}}\right)^{2}\frac{\pi\mathrm{Ca}|\bm{x}|^{2}}{4h_{*}(t)^{2}}\\\\[8.61108pt]
\hskip 210.0pt\text{on}\ \Gamma(t),\ t\in(0,T),\\\\[4.30554pt] \displaystyle
V(\bm{x},t)=-\frac{1}{\pi\mathrm{Ca}}\left(\frac{h_{0}}{R_{0}}\right)^{2}h_{*}(t)^{2}\nabla
p(\bm{x},t)\cdot\bm{N}(\bm{x},t)-\frac{\dot{h}_{*}(t)}{2h_{*}(t)}\bm{x}\cdot\bm{N}(\bm{x},t)\\\\[8.61108pt]
\hskip 210.0pt\text{for}\ \bm{x}\in\Gamma(t),\ t\in(0,T),\end{array}\right.$

where $R_{0}$ is the radius of the initial drop, $h_{0}$ is the initial gap,
$h_{*}(t)=\exp(t)$, $\varphi(\bm{x},t)=\varphi_{m}(\bm{x},t)/M_{c}$,
$\mathrm{Bm}$ is the magnetic Bond number, and $\mathrm{Ca}$ is the capillary
number.

The magnetic fluid is one of the smart fluids. Because of the property which
the smart fluids can be handled by controlling the electric and magnetic
field, in other words, without touching directly, they are widely used in
various fields, such as electric devices, mechanical engineering, and medical
applications. In the 1960s, magnetic fluids were developed as part of space
engineering development. R. E. Rosensweig is known as a pioneer of the
investigation into magnetic fluids. His works (for instance, [3, 18]) are
cited in many papers, even though they were contributed about half of a
century ago. Rosensweig advocated the governing equation of magnetic fluid
(see [18] for more details):

(17) $\dfrac{\partial\bm{v}}{\partial
t}+\left(\bm{v}\cdot\nabla\right)\bm{v}=-\nabla p^{*}+\mu_{0}M\nabla
H+\eta\triangle\bm{v}+\rho\bm{g}.$

This paper uses the magnetic fluid Hele-Shaw model (16), which can be derived
from his equation (17).

Magnetic fluid creates intriguing patterns. For instance, many spikes are
created in three-dimensional space by magnetic fluids in the magnetic field.
The phenomena, called the spike phenomena, are well known as one of the
intriguing phenomena caused by the collaboration of magnetic fluids and the
magnetic field. On the other hand, the fluids create fingering and the so-
called labyrinth patterns in a Hele-Shaw cell. Labyrinth patterns are mainly
observed only in quasi-two-dimensional space. These two-dimensional and three-
dimensional phenomena certainly have some common features; however, the
difference of patterns shows differences between two-dimensional and three-
dimensional problems. Therefore, proposing a simple and effective numerical
method for magnetic fluid Hele-Shaw problems will significantly impact
investigations into magnetic fluids.

In order to solve these problems, various numerical methods have been proposed
so far. The most successful method is based on the boundary integral method,
originating in the celebrated paper [9] by T.Y. Hou, J.S. Lowengrub, and M.J.
Shelley. Because of the surface tension, explicit time discretization requires
tiny time increments, and implicit discretization is difficult. By
reformulating the problem differently, they overcame these difficulties and
successfully computed the time evolution using a relatively large time
increment.

The original Hele-Shaw problem is known to have the following geometric
variational structure.

  * •

Curve-shortening property: the length $|\Gamma|(t)$ of $\Gamma(t)$ decays
monotonically in time;

  * •

Area-preserving property: the area $|\Omega|(t)$ of $\Omega(t)$ is constant in
time;

  * •

Barycenter-fixed property: the barycenter of $\Omega(t)$ does not move along
time evolution.

In both the Hele-Shaw problem with time-dependent gap (7) and the Hele-Shaw
problem for magnetic fluids (16), the volume-preserving property holds; that
is, $|\Omega|(t)\cdot h(t)$ is constant in time. Therefore, it is natural to
expect that these properties also hold in some sense in numerical methods. A
method that preserves the mathematical structure of the target equation in a
discrete sense is known as a structure-preserving numerical method, or
geometric integrator, and has been the subject of many studies [7].

Recently, the first and third authors of this paper proposed a structure-
preserving numerical scheme for the Hele-Shaw problem in [24]. The essential
idea is to combine the method of fundamental solutions (MFS for short) with
the uniform distribution method (UDM for short). The MFS is a numerical method
for the Laplace equation that does not require any numerical integration and
has the remarkable property that its approximation error decays exponentially
for the number of sample points [12, 11, 2, 21]. The UDM is vital for stable
numerical computation of moving boundary problems because it keeps the
distance between adjacent nodes uniform with time evolution [15, 16, 25, 26].
Furthermore, one of the main features of the numerical scheme proposed in [24]
is that, although it is in the sense of a semi-discrete problem, the area-
preserving property is strictly achieved, and the curve-shortening and
barycenter-fixed properties hold in an asymptotic manner. In this sense, the
paper [24] is the first structure-preserving numerical method to the Hele-Shaw
problem. Other structure-preserving numerical schemes for moving boundary
problems include, for example, [13, 23].

The purpose of this paper is to construct a numerical scheme satisfying the
volume-preserving property for the Hele-Shaw problem with a time-dependent gap
(7) and the Hele-Shaw problem with magnetic fluid (16) by extending the
contents of the paper [24]. In order to solve problem (16), we need to compute
the magnetostatic potential $\varphi$, i.e., the integral over $\Omega(t)$. In
this paper, we adopt the Monte Carlo method for numerical integration. In
addition, we improve the method adopted in [24] to arrange the singular points
that determine the accuracy of the approximate solution by the MFS.

This paper is organized as follows. In section 2, we construct a numerical
scheme for problems (7) and (16). In particular, we show that the scheme
satisfies the volume-preserving property as a semi-discrete scheme. Section 3
reports numerical results for these problems and shows that the proposed
method is much simpler than the previous ones but gives similar numerical
results. Finally, in section 4, we conclude this paper by presenting future
research topics.

##  2\. Numerical scheme

We approximate the boundary curve $\Gamma(t)$ by a polygonal curve; that is,
let $\Gamma(t)=\bigcup_{i=1}^{N}\Gamma_{i}(t)$ be an $N$-sided polygonal
Jordan curve, where $\Gamma_{i}(t)$ is the $i$-th edge of $\Gamma(t)$ defined
by

(18)
$\Gamma_{i}(t)=(\bm{X}_{i-1}(t),\bm{X}_{i}(t))\mathrel{\mathop{:}}=\\{\lambda\bm{X}_{i-1}(t)+(1-\lambda)\bm{X}_{i}(t)\mid\lambda\in(0,1)\\},$

where $\bm{X}_{i}(t)$ is the $i$-th vertex of $\Gamma(t)$ ($i=1,2,\ldots,N$.).
These $N$ vertices are labeled in the anti-clockwise order throughout this
paper, and we adopt the periodic notation $\bm{X}_{0}(t)=\bm{X}_{N}(t)$ and
$\bm{X}_{N+1}(t)=\bm{X}_{1}(t)$ (see Fig.(2)).

Figure 2. Jordan polygonal curve

The motion of the polygonal curve is then described by the following evolution
law, which is a system of ordinary differential equations:

(19) $\dot{\bm{X}}_{i}(t)=V_{i}(t)\bm{N}_{i}(t)+W_{i}\bm{T}_{i}(t),\quad
i=1,2,\ldots,N,\ t\in(0,T),$

where $\bm{N}_{i}$ and $\bm{T}_{i}$ are the unit outward normal and unit
tangent vectors at $\bm{X}_{i}$, and $V_{i}$ and $W_{i}$ denote the normal and
tangential velocities at $\bm{X}_{i}$, respectively.

We here summarize our numerical scheme. The details of each step are described
in the subsequent subsections.

For a given polygonal curve $\Gamma^{(m)}$ at the $m$-th time, we compute the
right-hand side of the evolution law (19) as in the following three steps.

Step 1::

Define the unit tangent vector $\bm{T}_{i}^{(m)}$ at $\bm{X}_{i}^{(m)}$ as the
unit vector pointing in the exact middle direction of the unit tangent vectors
defined on the two adjacent edges. Then, by rotating it by $-\pi/2$, we define
the unit outward normal vector $\bm{N}_{i}^{(m)}$ at $\bm{X}_{i}^{(m)}$. We
also define the discrete curvature $\kappa_{i}^{(m)}$ on each edge as the
first variation of the perimeter.

Step 2::

The MFS solves the Dirichlet problem, (7) or (16), in the polygonal region
$\Omega^{(m)}$. For example, (7) can be formulated as in the following using
the discrete curvatures $\\{\kappa_{i}^{(m)}\\}_{i=1}^{N}$ defined in Step 1:

(20) $\displaystyle\begin{dcases*}\triangle p(\cdot,t^{(m)})=0&in
$\Omega^{(m)}$,\\\
p(\bm{x},t^{(m)})=\sigma\kappa_{i}^{(m)}-\frac{\dot{h}(t^{(m)})}{4h(t^{(m)})^{3}}|\bm{x}|^{2}&for
$\bm{x}\in\Gamma_{i}^{(m)}$, $i=1,2,\ldots,N$.\end{dcases*}$

Using the approximate solution $P^{(m)}$ by the MFS and the unit outward
normal vector $\bm{N}_{i}^{(m)}$ at $\bm{X}_{i}^{(m)}$, the normal velocity
$V_{i}^{(m)}$ can be computed. The same procedure also works for problem (16).

Step 3::

Compute the tangential velocities $\\{W_{i}^{(m)}\\}_{i=1}^{N}$ using the
uniform distribution method.

Since the right-hand side of equation (19) can be computed by the above steps,
the evolution law (19) is computed by the fourth-order Runge–Kutta method to
obtain the polygonal curve $\Gamma^{(m+1)}$ at the next time $t^{(m+1)}$.

In the following, superscripts indicating time steps are omitted so long as
there is no risk of confusion.

###  2.1. Step 1: Spatial discretization

####  2.1.1. Definitions of $\bm{T}_{i}$ and $\bm{N}_{i}$

Let $r_{i}$ be the length of $\Gamma_{i}$ and let $\bm{t}_{i}$ and
$\bm{n}_{i}$ be the unit tangent vector and unit outward normal vector defined
naturally on $\Gamma_{i}$.

Unit tangent vectors and unit outward normal vectors at vertices are not well-
defined, and there are various possible ways of defining them. In this paper,
we adopt the following definitions. Let $\varphi_{i}$ be a signed exterior
angle at $\bm{X}_{i}$; that is, the angle between two adjacent edges
$\Gamma_{i}$ and $\Gamma_{i+1}$, satisfying
$\bm{t}_{i}\cdot\bm{t}_{i+1}=\cos\varphi_{i}$. Using this $\varphi_{i}$, we
define the unit tangent vector $\bm{T}_{i}$ as follows:

(21)
$\displaystyle\bm{T}_{i}=\frac{\bm{t}_{i}+\bm{t}_{i+1}}{|\bm{t}_{i}+\bm{t}_{i+1}|}=\frac{\bm{t}_{i}+\bm{t}_{i+1}}{2\mathsf{cos}_{i}},\quad
i=1,2,\ldots,N,$

where $\mathsf{cos}_{i}\coloneqq\cos(\varphi_{i}/2)$. The unit outward normal
vector $\bm{N}_{i}$ is then defined by $\bm{N}_{i}=-\bm{T}_{i}^{\bot}$.

####  2.1.2. Definition of discrete curvature

The curvature is usually defined as satisfying the Frenet formula. However,
since each side of a polygonal curve is a line segment, the usual definition
of curvature will lead to zero curvature almost everywhere, which is not
suitable to study the motion of polygonal curves. On the other hand, the
curvature can be interpreted as the first variation of the length functional.
This paper applies this view to polygonal curves and derives an appropriate
expression for discrete curvature.

The perimeter $L$ of a polygonal curve $\Gamma$ is given by

(22) $\displaystyle L=\sum_{i=1}^{N}r_{i}.$

According to the evolution law (19), the time derivative of $L$ can be
computed as follows:

(23) $\displaystyle\dot{L}=2\sum_{i=1}^{N}V_{i}\mathsf{sin}_{i},$

where $\mathsf{sin}_{i}=\sin(\varphi_{i}/2)$. Let us assume that the following
relation holds between the normal velocities $\\{V_{i}\\}_{i=1}^{N}$ at
vertices and the normal velocities $\\{v_{i}\\}_{i=1}^{N}$ on edges:

(24) $\displaystyle V_{i}=\frac{v_{i}+v_{i+1}}{2\mathsf{cos}_{i}},\quad
i=1,2,\ldots,N,$

which is an analogy of (21). Then, the equation (23) for the time derivative
of the perimeter is

(25)
$\displaystyle\dot{L}=\sum_{i=1}^{N}\kappa_{i}v_{i}r_{i},\quad\kappa_{i}=\frac{\mathsf{tan}_{i}+\mathsf{tan}_{i-1}}{r_{i}},\quad
i=1,2,\ldots,N,$

where $\mathsf{tan}_{i}=\mathsf{sin}_{i}/\mathsf{cos}_{i}$. We define
$\kappa_{i}$ to be the discrete curvature of the $i$-th edge $\Gamma_{i}$.

###  2.2. Step 2: Computation of normal velocities by the MFS

Normal velocities are computed using the MFS, a mesh-free numerical solver for
potential problems. For example, in the Hele-Shaw problem with a time-
dependent gap (7), the boundary is approximated by a polygonal curve, and the
discrete curvature gives the curvature on edge, so the potential problem to be
solved at time $t=t^{(m)}$ is

(26) $\displaystyle\begin{dcases*}\triangle p(\cdot,t^{(m)})=0&in
$\Omega^{(m)}$,\\\
p(\bm{x},t^{(m)})=\sigma\kappa_{i}^{(m)}-\frac{\dot{h}(t^{(m)})}{4h(t^{(m)})^{3}}|\bm{x}|^{2}&for
$\bm{x}\in\Gamma_{i}^{(m)}$, $i=1,2,\ldots,N$.\end{dcases*}$

Following the paper [24], the approximate solution $P$ based on the MFS is
given by

(27) $\displaystyle P(\bm{x})=Q_{0}+\sum_{j=1}^{N}Q_{j}E_{j}(\bm{x}),\quad
E_{j}(\bm{x})=E(\bm{x}-\bm{y}_{j})-E(\bm{x}-\bm{z}_{j}),$

where $E(\bm{x})=(2\pi)^{-1}\log|\bm{x}|$ is the fundamental solution of the
Laplace operator. The points $\\{\bm{y}_{j}\\}_{j=1}^{N}$ and
$\\{\bm{z}_{j}\\}_{j=1}^{N}$ are singular points and dummy points,
respectively, chosen “suitably” from outside the region $\Omega^{(m)}$. They
play a role in controlling the accuracy of the MFS approximate solution. The
specific arrangement of these points will be discussed in subsection 3.1.
Since $P$ satisfies the Laplace equation exactly, the linear combination
coefficients $\\{Q_{j}\\}_{j=0}^{N}$ are determined by an approximate
treatment of the Dirichlet boundary conditions and the volume-preserving
property. In this study, we adopt the midpoint
$\bm{X}_{i}^{*}=(\bm{X}_{i}+\bm{X}_{i-1})/2$ of each edge as a collocation
point and approximate the Dirichlet boundary condition using the collocation
method. In other words, $\\{Q_{j}\\}_{j=0}^{N}$ are given as a solution to the
following linear system

(28) $\displaystyle
P(\bm{X}_{i}^{*})=\sigma\kappa_{i}^{(m)}-\frac{\dot{h}(t^{(m)})}{4h(t^{(m)})^{3}}|\bm{X}_{i}^{*}|^{2},\quad
i=1,2,\ldots,N,$

with a constraint

(29) $\displaystyle\sum_{j=1}^{N}Q_{j}H_{j}=0,\quad H_{j}=\sum_{i=1}^{N}\nabla
E_{j}(\bm{X}_{i}^{*})\cdot\bm{n}_{i}r_{i},\quad j=1,2,\ldots,N.$

By differentiating the approximate solution $P$ in the normal direction
$\bm{n}_{i}$, the normal velocity $v_{i}$ on the $i$-th edge can be computed
by

(30) $\displaystyle v_{i}=-h(t^{(m)})^{2}\left(\nabla
P(\bm{X}_{i}^{*})+\frac{\dot{h}(t^{(m)})}{2h(t^{(m)})}\bm{X}_{i}^{*}\right)\cdot\bm{n}_{i},\quad
i=1,2,\ldots,N.$

The great advantage of using the MFS is that the gradient $\nabla P$ can be
computed analytically. In other words, there is no need to perform numerical
differentiation. Finally, the normal velocities $\\{V_{i}\\}_{i=1}^{N}$ at
vertices can be obtained by using the relation (24).

At the end of this subsection, we clarify what constraint (29) means. To this
end, we state a little more precisely the volume-preserving property. Since
the volume $V(t)$ of the fluid region $\Omega(t)\times(0,h(t))$ is given by
$V(t)=A(t)h(t)$, where $A(t)=|\Omega|(t)$, the volume-preserving property
$\dot{V}=0$ can be rephrased as the following equation for the time evolution
of the area:

(31) $\displaystyle\dot{A}(t)=-\frac{\dot{h}(t)}{h(t)}A(t).$

According to Proposition 4 in the paper [24], the time evolution of the area
under the evolution law (19) is given by

(32)
$\displaystyle\dot{A}=\sum_{i=1}^{N}v_{i}r_{i}+\mathrm{err}_{A},\quad\mathrm{err}_{A}=\sum_{i=1}^{N}\left(W_{i}\mathsf{sin}_{i}-\frac{v_{i+1}-v_{i}}{2}\right)\frac{r_{i+1}-r_{i}}{2}.$

If the error term $\mathrm{err}_{A}$ is equal to $0$, then the normal velocity
formula (30) gives the time evolution of the area as

(33)
$\displaystyle\dot{A}(t)=-h(t)^{2}\sum_{j=1}^{N}Q_{j}H_{j}-\frac{\dot{h}(t)}{h(t)}A(t).$

Therefore, because of (31), we obtain constraint (29).

For the above argument to be consistent, the error term $\mathrm{err}_{A}$
must be equal to $0$. For this purpose, we adopt the uniform distribution
method to compute the tangential velocity. Details are given in subsection
2.3.

The same consideration can be applied to the Hele-Shaw problem for magnetic
fluids (16) to obtain an approximate solution satisfying the volume-preserving
property.

###  2.3. Step 3: Computation of tangential velocities by the UDM

As mentioned in the previous subsection, the error term must be $0$ to achieve
the volume-preserving property. Given the representation (32) of the error
term, one possibility is to apply the UDM. It is a method of moving vertices
in the tangential direction so that the length of each edge of the polygon is
equal, i.e., $r_{i}=L/N$ ($i=1,2,\ldots,N$), and is known to provide stable
numerical computations. In this study, we adopt the asymptotic UDM, which
achieves uniform distribution in an asymptotic sense. This method is robust to
the effects of numerical errors such as rounding errors, and the tangential
velocities are given as linear expressions of normal velocities. The
asymptotic UDM requires that $r_{i}(t)\rightarrow L(t)/N$ as $t\to T_{\max}$,
not $r_{i}(t)=L(t)/N$, holds, where $T_{\max}$ represents a final computation
time.

We assume that the following relations hold to derive an expression for the
tangential velocity by the asymptotic UDM:

(34) $r_{i}-\dfrac{L}{N}=\eta_{i}e^{-\mu(t)},\quad i=1,2,\ldots,N,$

where $\\{\eta_{i}\\}_{i=1}^{N}$ are parameters chosen so that their mean is
equal to $0$. The function $\mu$ controls the strength of the effect of the
uniform arrangement, and a candidate satisfies $\mu(t)\rightarrow\infty$ as
$t\to T_{\max}$. Differentiating both sides of equation (34) with respect to
$t$ yields

(35)
$\dot{r}_{i}=\dfrac{\dot{L}}{N}+\left(\dfrac{L}{N}-r_{i}\right)\omega(t),\quad
i=1,2,\ldots,N,$

where $\omega(t)=\dot{\mu}(t)$. On the other hand, if we compute the time
derivative of the length of the edge using the evolution law (19), we obtain

(36)
$\displaystyle\dot{r}_{i}=V_{i}\mathsf{sin}_{i}+V_{i-1}\mathsf{sin}_{i-1}+W_{i}\mathsf{cos}_{i}-W_{i-1}\mathsf{cos}_{i-1},\quad
i=1,2,\ldots,N.$

Combining (35) and (36), we obtain recurrence relations describing the
tangential velocity:

$W_{i}\mathsf{cos}_{i}-W_{i-1}\mathsf{cos}_{i-1}=\frac{\dot{L}}{N}+\left(\frac{L}{N}-r_{i}\right)\omega(t)-V_{i}\mathsf{sin}_{i}-V_{i-1}\mathsf{sin}_{i-1}$

for $i=2,3,\ldots,N$. Imposing the zero-average condition
$\sum_{i=1}^{N}W_{i}=0$, we obtain $N$ linearly independent equations. This
system can be solved analytically, resulting in the following expression for
the tangential velocity in the asymptotic UDM:

(37) $\displaystyle
W_{i}=\frac{\Psi_{i}+W_{1}\mathsf{cos}_{1}}{\mathsf{cos}_{i}},\quad
i=2,3,\ldots,N,\quad
W_{1}=-\frac{\sum_{i=2}^{N}\Psi_{i}/\mathsf{cos}_{i}}{\mathsf{cos}_{1}\sum_{j=1}^{N}\mathsf{cos}_{j}^{-1}},$

where

(38)
$\displaystyle\Psi_{i}=\sum_{l=2}^{i}\psi_{l},\quad\psi_{i}=-V_{i}\mathsf{sin}_{i}-V_{i-1}\mathsf{sin}_{i-1}+\frac{\dot{L}}{N}+\left(\frac{L}{N}-r_{i}\right)\omega$

for $i=2,3,\ldots,N$. Note that if the initial arrangement of vertices is
uniform, then the asymptotic UDM theoretically keeps uniformity during time
evolution.

There is another way to eliminate $\mathrm{err}_{A}$. If $W_{i}$ is defined by

(39) $\displaystyle W_{i}=\frac{v_{i+1}-v_{i}}{2\mathsf{sin}_{i}},\quad
i=1,2,\ldots,N,$

then the error term $\mathrm{err}_{A}$ also vanishes. However, in this case,
the orientation of the edges of the polygonal curve does not change in time.
Namely, the time evolution of the polygon is restricted to the crystalline
motion. This tangential velocity is useful for theoretical analysis, but it is
not suitable for numerical computation of problems with large deformation of
the domain, such as the Hele-Shaw problem. Therefore, in this study, we adopt
the asymptotic UDM given by equation (37).

##  3\. Numerical result

In this section, we show some results of our numerical computations for (7)
and (16).

###  3.1. Arrangement of points by the modified Amano’s method

In solving problems (7) or (16) by the MFS, singular points and dummy points
must be appropriately placed. In the paper [24], the singular points are
placed at a fixed distance in the unit outward normal direction from the
midpoints of the edges:

(40) $\displaystyle\bm{y}_{j}=\bm{X}_{j}^{*}+d\bm{n}_{j},\quad
j=1,2,\ldots,N.$

However, there is no mathematical background to this point placement, and
there may be better ways to arrange them.

In this study, we employ a modified Amano’s method proposed in numerical
computation of conformal mappings (see for instance [1]). Namely, the singular
points are placed according to the following relations:

(41)
$\displaystyle\bm{y}_{j}^{a}=\bm{X}_{j}^{*}+\frac{r_{a}}{2}|\bm{X}_{j+1}^{*}-\bm{X}_{j-1}^{*}|\bm{n}_{j}^{a},$

where $\bm{n}_{j}^{a}$ is a modified unit outward normal vector determined by
the positions of the three points $\bm{X}_{j-1}^{*}$, $\bm{X}_{j}^{*}$, and
$\bm{X}_{j+1}^{*}$ and is defined through

(42)
$\displaystyle\bm{n}_{j}^{a}=-\frac{(\bm{X}_{j+1}^{*}-\bm{X}_{j-1}^{*})^{\bot}}{|\bm{X}_{j+1}^{*}-\bm{X}_{j-1}^{*}|}.$

The geometric picture of the arrangement of singular points is shown in Fig.
3.

Figure 3. Arrangement of the singular point. $\bm{y}_{j}$ corresponds to the
original arrangement (40) adopted in [24], while $\bm{y}_{j}^{a}$ does to the
modified Amano’s arrangement (41).

Amano’s method can be though of as corresponding to a discrete version of
point arrangement using a conformal mapping, which is used in the theoretical
analysis of the MFS. For example, in Fig. 3, the thin line expresses the
original smooth curve. In this case, $\bm{n}_{j}^{a}$ is considered an
approximation of the unit outward normal vector at $\bm{X}_{i}^{*}$. This
perspective is correct, and it can be mathematically proved that Amano’s
method approximates the method using conformal mappings. For more details,
refer to the paper [22]. Hence, we use the modified Amano’s arrangement (41)
to place the singular points in this paper.

While the singular points have a significant impact on the accuracy of the
approximate solution, the dummy points are only introduced for the invariance
of the approximate solution and have little effect on the accuracy of the
approximate solution. Therefore, in this paper, we adopt the following simple
arrangement:

(43)
$\displaystyle\bm{z}_{j}=1000(\cos\alpha_{i},\sin\alpha_{i})^{\top},\quad\alpha_{i}=\frac{2\pi
i}{N},\quad i=1,2,\ldots,N.$

###  3.2. Hele-Shaw problem with a time-dependent gap

We show the results of numerical computations for (7). The initial curve
$\Gamma(0):[0,1]\ni u\mapsto(x_{1}(u),x_{2}(u))^{\top}\in\mathbb{R}^{2}$ is
given by

$\displaystyle x_{1}(u)=r\cos(2\pi u),\quad x_{2}=r\sin(2\pi u),$
$\displaystyle r=R_{0}+0.02(\cos(6\pi u)+\sin(14\pi u)+\cos(30\pi
u)+\sin(50\pi u)).$

We set the parameters as in Table 1.

Table 1. parameters for the Hele-Shaw problem with a time-dependent gap (7) $N=300$ | the number of vertices
---|---
$\sigma=2.0\times 10^{-4}$ | the surface tension coefficient
$\Delta t=1/(10N^{2})$ | the time-mesh size
$r_{a}=1$ | the parameter in the Amano’s method (41)
$R_{0}=1$ | the initial reference radius

The time-dependent gap $h(t)$ is given by $h(t)=\exp(t)$. In Fig.4, the time
evolution of the boundary $\Gamma(t)$ is shown at $t=0$, $0.36$, $0.73$,
$1.10$, $1.47$, $1.84$, $2.21$, and $2.58$. The shapes at each time step
almost coincide with the earlier study by Shelley et al. [27]. Note that this
is not a comparison at the same discrete time because we use a different time
discretization than in previous studies.

Fig.5 shows the time evolution of the total length, the area and the volume,
where the horizontal axis represents the time, the left vertical axis
represents the total length and the area, and the right vertical axis
represents the volume. The graph shows that the area decreases monotonically
and that the volume-preserving property is satisfied.

(a) $t=0.0$ (b) $t=0.36$ (c) $t=0.73$ (d) $t=1.10$
---
(e) $t=1.47$ (f) $t=1.84$ (g) $t=2.21$ (h) $t=2.58$
Figure 4. Results of numerical computation for Hele-Shaw flow with time-
dependent gap: the time evolution of the boundary $\Gamma(t)$ Figure 5.
Results of numerical computation for Hele-Shaw flow with time-dependent gap:
time evolution of the total length, the area and the volume

###  3.3. Hele-Shaw problem for magnetic fluids

In this subsection, we report the numerical results of the Hele-Shaw problem
for magnetic fluids (16).

####  3.3.1. Computation of the magnetostatic field potential $\varphi_{m}$ by
adaptive Monte Carlo method

To solve problem (16), we need to compute the magnetostatic field potential
$\varphi_{m}$; that is, integration on $\Omega(t)$. In this paper, we employ
the adaptive Monte Carlo method.

At every time step, the region $\mathcal{MC}(t)\supseteq\Omega(t)$ is defined
as follows:

$\displaystyle\mathcal{MC}(t)$
$\displaystyle\coloneqq\left\\{\begin{pmatrix}r(t)\cos\theta+g_{x}(t)\\\
r(t)\sin\theta+g_{y}(t)\end{pmatrix}\mathrel{}\middle|\mathrel{}r(t)\in[0,r_{\max}(t)],\
\theta\in[0,2\pi)\right\\},$ $\displaystyle r_{\mathrm{m}ax}(t)$
$\displaystyle\coloneqq\max\\{|\bm{X}_{i}(t)-\bm{g}(t)|\mid
i=1,2,\ldots,N\\},$

where $\bm{g}(t)=(g_{x}(t),g_{y}(t))^{\top}$ is the barycenter of $\Omega(t)$.
Then, $M$ sample points $\\{\bm{p}_{k}\\}_{k=1}^{M}$ are distributed in
$\mathcal{MC}$ uniformly.

Distributed sample points $\\{\bm{p}_{k}\\}_{k=1}^{M}$ are judged whether
$\bm{p}_{k}\in\Omega$ by the following algorithm. Let $a_{k}^{i}$ be the angle
between the vectors $\bm{X}_{i-1}-\bm{p}_{k}$ and $\bm{X}_{i}-\bm{p}_{k}$,
i.e.,

$a_{k}^{i}=\arccos\left(\frac{\bm{X}_{i-1}-\bm{p}_{k}}{|\bm{X}_{i-1}-\bm{p}_{k}|}\cdot\frac{\bm{X}_{i}-\bm{p}_{k}}{|\bm{X}_{i}-\bm{p}_{k}|}\right).$

If $\sum_{i=1}^{N}a_{k}^{i}=2\pi$ holds, then the sample point $\bm{p}_{k}$ is
judged to be inside $\Omega$; otherwise, it is judged outside $\Omega$.

Let $M_{\mathrm{i}n}$ be the number of points $\bm{p}_{k}$ belonging to
$\Omega$. Then the area element $dS$ of $\Omega$ is defined as $dS=A/M_{in}$.
By using this $dS$, the the value of $\varphi=\varphi_{m}/M_{c}$ at the vertex
$\bm{X}_{i}$ can be computed approximately as follows:

$\varphi(\bm{X}_{i})=-\sum_{\bm{p}_{k}\in\Omega}\left(\dfrac{1}{|\bm{X}_{i}-\bm{p}_{k}|}-\dfrac{1}{\sqrt{|\bm{X}_{i}-\bm{p}_{k}|^{2}+h^{2}}}\right)dS$

for $i=1,2,\ldots,N$.

####  3.3.2. Numerical results

We show the results of numerical computations for (16). The initial curve
$\Gamma(0):[0,1]\ni u\mapsto(x_{1}(u),x_{2}(u))^{\top}\in\mathbb{R}^{2}$ is
given by

$x_{1}(u)=r\cos(2\pi u),\quad x_{2}=r\sin(2\pi u).$

In Figs. 6 and 7, $r$ is given by

$r=R_{0}+0.02(\cos(6\pi u)+\sin(14\pi u)+\cos(30\pi u)+\sin(50\pi u)),$

and in Figs. 9 and 10, $r$ is given

$r=R_{0}+0.05(\cos(4\pi u)-\cos(10\pi u)+\cos(22\pi u)-\sin(6\pi u)+\sin(10\pi
u)).$

We set the parameters as in Table 2.

Table 2. parameters for the Hele-Shaw flow for magnetic fluids (16) $N=300$ | the number of vertices
---|---
$\Delta t=1/(10N^{2})$ | the time-mesh size
$r_{a}=1$ | the parameter in the Amano’s method (41)
$R_{0}=1$ | the initial reference radius
$h_{r}=h_{0}/R_{0}=0.25$ | the parameter introduced for simplifying the model
$M=1000$ | the number of sample points in the Monte Carlo method

In each Fig. 6, 7, 9, and 10, the time evolution of the boundary $\Gamma(t)$
is shown at $t=0$, $0.42$, $0.84$, $1.26$, $1.68$, $2.10$, $2.52$, and $2.94$.
Comparing cases $\mathrm{Bmv}=0$, $25$, and $35$, where
$\mathrm{Bmv}=(R_{0}/h_{0})^{1/3}\mathrm{Bm}$, it can be seen that the
complexity of the pattern formed by the magnetic fluid changes depending on
the strength of the magnetic fluid. Moreover, the shapes at each time almost
coincide with earlier studies by Tatulchenkov and Cebers [29]. Note that these
are not comparisons at the same discrete time because we use a different time
discretization than in previous studies.

Figs. 8 and 11 show the time evolution of the total length, the area, and the
volume, where the horizontal axis represents the time, the left vertical axis
represents the total length and the area, and the right vertical axis
represents the volume. The graph shows that the volume-preserving property
holds.

(a) $t=0.0$ (b) $t=0.42$ (c) $t=0.84$ (d) $t=1.26$
---
(e) $t=1.68$ (f) $t=2.10$ (g) $t=2.52$ (h) $t=2.94$
Figure 6. $\mathrm{Bmv}=0,\mathrm{Ca}=100$ (a) $t=0.0$ (b) $t=0.42$ (c)
$t=0.84$ (d) $t=1.26$
---
(e) $t=1.68$ (f) $t=2.10$ (g) $t=2.52$ (h) $t=2.94$
Figure 7. $\mathrm{Bmv}=25,\mathrm{Ca}=100$ Figure 8. time evolution of the
Area, the total Length and the Volume ($\mathrm{Bmv}=25,\mathrm{Ca}=100$) (a)
$t=0.0$ (b) $t=0.42$ (c) $t=0.84$ (d) $t=1.26$
---
(e) $t=1.68$ (f) $t=2.10$ (g) $t=2.52$ (h) $t=2.94$
Figure 9. $\mathrm{Bmv}=0,\mathrm{Ca}=50$ (a) $t=0.0$ (b) $t=0.42$ (c)
$t=0.84$ (d) $t=1.26$
---
(e) $t=1.68$ (f) $t=2.10$ (g) $t=2.52$ (h) $t=2.94$
Figure 10. $\mathrm{Bmv}=35,\mathrm{Ca}=50$ Figure 11. time evolution of the
Area, the total Length and the Volume ($\mathrm{Bmv}=35,\mathrm{Ca}=50$)

##  4\. Conclusion

In this paper, several Hele-Shaw problems were solved using the MFS combined
with Amano’s method and the asymptotic uniform distribution method. The shown
results coincided with preceding studies; moreover, some of the properties of
Hele-Shaw flow, for example, volume-preserving property, are precisely
satisfied. Therefore, the results ensured that our numerical method can be
adequately effective for Hele-Shaw problems. Application of the MFS in moving
boundary problems has been reported in quite a few papers; however, our
present attempt will widen the possibility of the MFS. It is expected that the
method will be applied to Hele-Shaw type problems and any other types of
potential problems.

Besides, magnetic fluid instabilities in the Hele-Shaw cell were simulated,
then the magnetic effects were indeed observed by comparing numerical results.
For large $\mathrm{Bm}$, the number of fingers of fingering patterns
increased. It is similar to the phenomena in three-dimensional space; that is,
the number of spikes increases with the high intensity of the magnetic field.
On the other hand, complex patterns such as our results do not appear in
three-dimensional space. To investigate the difference in the behavior of
magnetic fluids between two-dimensional and three-dimensional spaces
mathematically is interesting. We hope that new aspects for investigations
into magnetic fluid are found from our results.

## References

  * [1] K. Amano, D. Okano, H. Ogata, and M. Sugihara.  Numerical conformal mappings onto the linear slit domain.  Jpn. J. Ind. Appl. Math., 29(2):165–186, 2012.
  * [2] A.H. Barnett and T. Betcke.  Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains.  J. Comput. Phys., 227(14):7003–7026, 2008.
  * [3] M.D. Cowley and R.E. Rosensweig.  The interfacial stability of a ferromagnetic fluid.  J. Fluid Mech., 30:671–688, 1967.
  * [4] J. Dockery and I. Klapper.  Finger formation in biofilm layers.  SIAM J. Appl. Math., 62(3):853–869, 2001/02.
  * [5] F. Elias, C. Flament, J.-C. Bacri, and S. Neveu.  Macro-organized patterns in ferrofluid layer: Experimental studies.  Journal de Physique I, 7(5):711–728, 1997.
  * [6] B. Gustafsson and A. Vasil’ev.  Conformal and potential analysis in Hele-Shaw cells.  Advances in Mathematical Fluid Mechanics. Birkhäuser Verlag, Basel, 2006.
  * [7] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics.  Springer, Heidelberg, 2010.  Structure-preserving algorithms for ordinary differential equations, Reprint of the second (2006) edition.
  * [8] Henry S. Hele-Shaw.  The flow of water.  Nature, 58:34–36, 1898.
  * [9] T.Y. Hou, J.S. Lowengrub, and M.J. Shelley.  Removing the stiffness from interfacial flows with surface tension.  J. Comput. Phys., 114(2):312–338, 1994.
  * [10] S.D. Howison.  Fingering in Hele-Shaw cells.  J. Fluid Mech., 167:439–453, 1986.
  * [11] M. Katsurada.  Asymptotic error analysis of the charge simulation method in a Jordan region with an analytic boundary.  J. Fac. Sci. Univ. Tokyo Sect. IA Math., 37(3):635–657, 1990.
  * [12] M. Katsurada and H. Okamoto.  A mathematical study of the charge simulation method. I.  J. Fac. Sci. Univ. Tokyo Sect. IA Math., 35(3):507–518, 1988.
  * [13] T. Kemmochi.  Energy dissipative numerical schemes for gradient flows of planar curves.  BIT, 57(4):991–1017, 2017.
  * [14] H. Lamb.  Hydrodynamics.  Cambridge Mathematical Library. Cambridge University Press, Cambridge, sixth edition, 1993.  With a foreword by R. A. Caflisch [Russel E. Caflisch].
  * [15] Karol Mikula and Daniel Ševčovič.  A direct method for solving an anisotropic mean curvature flow of plane curves with an external force.  Math. Methods Appl. Sci., 27(13):1545–1565, 2004.
  * [16] Karol Mikula and Daniel Ševčovič.  Evolution of curves on a surface driven by the geodesic curvature and external force.  Appl. Anal., 85(4):345–362, 2006.
  * [17] F. Otto.  Dynamics of labyrinthine pattern formation in magnetic fluids: a mean-field theory.  Arch. Rational Mech. Anal., 141(1):63–103, 1998.
  * [18] R.E. Rosensweig.  Magnetic fluids.  Ann. Rev. Fluid Mech., 19:437–463, 1987.
  * [19] R.E. Rosensweig, M. Zahn, and R. Shumovich.  Labyrinthine instability in magnetic and dielectric fluids.  Journal of Magnetism and Magnetic Materials, 39(1-2):127–132, 1983\.
  * [20] P.G. Saffman and G. Taylor.  The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid.  Proc. Roy. Soc. London Ser. A, 245:312–329. (2 plates), 1958.
  * [21] K. Sakakibara.  Asymptotic analysis of the conventional and invariant schemes for the method of fundamental solutions applied to potential problems in doubly-connected regions.  Jpn. J. Ind. Appl. Math., 34(1):177–228, 2017.
  * [22] K. Sakakibara.  Bidirectional numerical conformal mapping based on the dipole simulation method.  Eng. Anal. Bound. Elem., 114:45–57, 2020.
  * [23] K. Sakakibara and Y. Miyatake.  A fully discrete curve-shortening polygonal evolution law for moving boundary problems.  J. Comput. Phys., 424:109857, 22, 2021.
  * [24] K. Sakakibara and S. Yazaki.  Structure-preserving numerical scheme for the one-phase Hele-Shaw problems by the method of fundamental solutions.  Comput. Math. Methods, 1(6):e1063, 25, 2019.
  * [25] Daniel Ševčovič and Shigetoshi Yazaki.  Evolution of plane curves with a curvature adjusted tangential velocity.  Jpn. J. Ind. Appl. Math., 28(3):413–442, 2011.
  * [26] Daniel Ševčovič and Shigetoshi Yazaki.  On a gradient flow of plane curves minimizing the anisoperimetric ratio.  IAENG Int. J. Appl. Math., 43(3):160–171, 2013.
  * [27] M.J. Shelley, F.-R. Tian, and K. Wlodarski.  Hele-Shaw flow and pattern formation in a time-dependent gap.  Nonlinearity, 10(6):1471–1495, 1997.
  * [28] S. Tanveer.  Surprises in viscous fingering.  J. Fluid Mech., 409:273–308, 2000.
  * [29] A. Tatulchenkov and A. Cebers.  Magnetic fluid labyrinthine instability in Hele-Shaw cell with time dependent gap.  Phys. Fluids, 20(5):054101, 10 pp., 2008.