text:
\section{Introduction}
Understanding the collapse of, fully polymeric, topologically complex objects in a bad solvent is of broad importance,
because of its relevance in the early stages of the protein folding \cite{Sadqi2003,Majumder2019}, and its connections with arresting processes or aging phenomena \cite{Majumder2017,Midya2015}.
The thermodynamic and structural properties of the macromolecular collapse in solution, occurring when the quality of the solvent is decreased below a critical value, have been exhaustively investigated over the years \cite{deGennes1979, Kamerlin2016,Grosberg1988,terashima2011,huerta2013,kang2015}, and at equilibrium, the dynamic and static behavior are well known above and below the volume phase transition \cite{Doi1986}. Comparatively, a general framework  for
the non-equilibrium aspects, as the kinetics of the collapse, is still lacking.
The development of new experimental setups for small angle X-ray scattering or single-molecule fluorescence spectroscopy allows to monitor the collapse of a single molecule \cite{Pollack2001,hofmann2012,Zhang2014}, sheding light on the aspects controlling the collapse dynamics.

Recent computational works have proposed analogies between the collapse of linear polymer chains in bad solvent \cite{Majumder2015EPL,Majumder2017,Majumder-review} and phase-ordering phenomena or coarsening systems (as foams \cite{Saint-Jalmes2006} or polymer blends \cite{Veenstra1999}), where coarsening refers to any out-of-equilibrium relaxation process involving the growth of two separated phases from an initially homogeneous mixture. In a coarsening system the characteristic length scale grows over time and exhibits key universal features like dynamic scaling and self-similarity \cite{Smerlak2018}. In the context of macromolecular collapse in bad solvent, the coarsening process initiates with the formation of small clusters of monomers along the polymeric strands. Afterwards, these clusters become stable and start to grow by withdrawing monomers from the bridges connecting them or by coalescence with other clusters. This moves forward until all the monomers pile up to a single cluster. Finally, the single cluster collapses to the ultimate equilibrium state, characterized by a compact object of melt-like density \cite{Majumder2017}.
The characteristic length scale of the coarsening process, associated to the collapse from the initial
self-avoiding conformations of the polymer strands to the late fully-collapsed state,
is given by the size of the growing aggregates of monomers.

Some theoretical approaches have been developed for describing the different scaling regimes associated to the coarsening kinetics.
In binary liquids, where the chemical potential gradient acts as the driving force \cite{Bray2002}, the domains grow with time as $\sim t^{1/3}$.
For fluids and polymers, where hydrodynamic contributions are relevant, linear growth is predicted
at late times in the viscous regime \cite{Siggia1979} .
An apparent sublinear regime for the domain growth, $\sim t^{1/2}$, has been
observed at intermediate times in simulations of a gas-liquid separating system \cite{Testard2014}, suggesting an effective interpolation between the former early and late regimes.
In a previous work by some of us \cite{Moreno2018} a scaling exponent of $\sim 0.7$ has been found at intermediate times for the growth of the domains during the deswelling of microgels in bad solvent. This accelaration (higher exponent) with respect to the system of Ref.~\cite{Testard2014} is tentatively related to the polymer connectivity of the microgels, which facilitates the merging of close clusters.

In this contribution we investigate whether the scaling behaviour found for the microgels of Ref.~\cite{Moreno2018} is a general result for other macromolecules
or if it depends significantly on the architecture of the system.
In order to quantify the growing length scale we introduce a density field representation of the macromolecules that removes artifacts arising from the local, fast
density fluctuations in the coarsening structure \cite{Testard2014}.
We establish a set of scaling laws for the time dependence of the growing domain size during the coarsening, which are independent of the solvent quality parameter in
analogy to time-temperature superposition principle.
Domain growth in microgels shows a power law, though an overshoot is found in the late stage of the collapse for the case of diamond-like networks. This unusual behavior is related to the fast late merging of the regularly distributed nucleating centers.
Power-laws are also observed for the coarsening dynamics of collapsing flexible linear chains and star polymers, though with slightly smaller exponents imputable to a lower number of nucleation centers in the absence of a network structure. Semiflexible chains present a significantly higher exponent for the domain growth.
We also analyze the fractal structure of the clusters of different sizes that are formed during the coarsening process, and find no differences between the flexible and semiflexible systems, concluding that the higher exponent for the growing length scale originates from a faster mass diffusion in the semiflexible chains induced by the stiffness.
Finally we characterize the dynamic density correlations during the collapse and relate them with the growing length scale. A common power-law is found for this relation, irrespective of the solvent quality and macromolecular architecture, which is a typical signature of the critical nature of the process.

The paper is organised as follows. In Section 2 the model and the simulation details are presented. In Section 3 a deep analysis of the kinetics of the collapse for all the investigated systems is reported. In particular, the density field construction is introduced, and data for the domain growth  and the scaling of the dynamic correlations are discussed. Section 4 summarizes our conclusions.

\section{Model and simulation details}\label{sec2}

We performed NVT molecular dynamics simulations
of star polymers with 3 and 12 arms, semiflexible and flexible linear chains, and migrogels with a disordered and a diamond-like network structure.
We used the bead-spring model
of Kremer and Grest \cite{Kremer1990} for the interactions. All the non-bonded interactions were implemented by a, purely repulsive,
Weeks-Chandler-Andersen (WCA) potential \cite{Weeks1971} that was modified in order to tune the quality of the solvent. This was achieved by introducing an attractive tail, regulated by a solvent quality parameter $\phi$ that sets the solvophobicity of the monomers \cite{soddemann2001EPJE,LoVerso2015,Gnan2017,Moreno2018}. Therefore this parameter plays the role of an effective temperature in the model. The non-bonded interactions were given by:
\begin{equation}
V_{\rm nb}(r) = \left\{
\begin{array}{lll}
V_{\rm LJ}(r)=  4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}
- \left(\frac{\sigma}{r}\right)^{6} +\frac{1}{4} \right] - \epsilon\phi &\ &\  r\leq 2^{1/6}\sigma \\
V_{\phi}(r) = \frac{1}{2} \phi\epsilon \left[\cos\left(\alpha (\frac{r}{\sigma})^2 +\beta\right)-1\right]  &\ &\  2^{1/6}\sigma < r\leq 1.5\sigma \\
 0 &\ &\ r > 1.5\sigma
\end{array}
 \right.
\label{eq:vnb}
\end{equation}
The values for the parameters  $\alpha = \pi(2.25 -2^{1/3})^{-1}$ and $\beta = 2\pi -2.25\alpha$ are chosen in order to satisfy the condition that the non-bonded potential and its first derivative are continuous both at $r=2^{1/6}\sigma$ and at the cutoff $r_{\rm c}=1.5\sigma$ \cite{soddemann2001EPJE}.
For the case $\phi=0$ the purely repulsive WCA potential is recovered and the system is in good solvent conditions.
The quality of the solvent is worsened by choosing $\phi > 0$. Bad solvent conditions are reached when $\phi$ is higher than some critical value and the system collapses. In all the studied systems the collapse transition ($\theta$-point) occurred at $\phi \sim 0.6$ (this was estimated from the maximum of the derivative of the radius o gyration vs. $\phi$).

In addition, bonded monomers interact via a finite extensible nonlinear elastic (FENE) potential, which implements the molecular connectivity. The FENE potential reads \cite{Kremer1990}:
\begin{equation}V_{\rm FENE}(r) = -\epsilon K_{\rm F}R_0^2 \ln\left[1-\left(\frac{r}{R_0\sigma}\right)^2 \right] ,
\label{eq:fene}
\end{equation}
where $K_{\rm F} = 30$ is the spring constant and $R_0 = 1.5$ is the maximum elongation. The sum of the FENE and the non-bonded potential provides a total interaction between two connected monomers showing a deep minimum at $r \sim 0.95$, which guarantees uncrossability and prevents violation of the topological constraints.
In the case of the semiflexible linear chains the bending stiffness was implemented through the
worm-like model \cite{Zierenberg2016, Kratky1949}.  Thus, the interaction for the polymer bending has the form:
\begin{equation} V_{\rm bending}(\theta) = K_{\rm s}(1-\cos\theta) ,
\label{eq:bending}
\end{equation}
where $\theta$ is the angle between two successive bonds and $K_{\rm s}$ is the strength of the bending. We used a value $K_{\rm s} = 5$.
The units of energy, length, mass and time are respectively $\epsilon$, $\sigma$, $m$ and $\tau=(\sigma^2 m/\epsilon)^{1/2}$ where $m$ is the mass of a monomer.
In the rest of the paper all the numerical values will be given in reduced units $\epsilon = \sigma = m = \tau = 1$.

We performed MD simulations at temperature $T = \epsilon/k_{\rm B} = 1.0$ (with $k_{\rm B}$ the Boltzmann constant) using a Langevin
thermostat \cite{Izaguirre2001,Smith2009}. Thus, the force experienced by the monomers is:
\begin{equation}\ddot{\textbf{r}}_{i} =  -{\bf \nabla} V(\textbf{r}_i) - \gamma \dot{\textbf{r}}_{i} + \sqrt{2 \gamma k_{\rm B} T} {\bf \zeta}(t)
,
\label{eq:langevin}
\end{equation}
where $\textbf{r}_i$ is the position vector and $V(\textbf{r}_i)$ is the total interaction potential for the monomer of index $i$. The second term of the right side of Equation~(\ref{eq:langevin}) represents viscous damping, with $\gamma$ the friction coefficient. The last term is a random uncorrelated force, $\langle \zeta^{\alpha}_i(t)\zeta^{\beta}_j(t')\rangle = \delta_{i,j}\delta_{\alpha,\beta}\delta(t-t')$ (with $\alpha,\beta$ the Cartesian components), representing the collisions with solvent particles. The Langevin thermostat acts therefore as an implicit solvent, in which every particle interacts independently with the solvent, but hydrodynamic interactions between solute particles are not considered. Their inclusion would require the use
of e.g., Lattice Boltzmann \cite{latticeboltzmann}, multi-particle collision dynamics (MPCD) \cite{Malevanets1999} or
dissipative particle dynamics \cite{Espanol1995} methods, involving a huge computational cost due to the
size of the investigated systems and the big boxes needed to avoid finite size effects. For example, in the widely used MPCD method, radii of gyration $R_{\rm g} \sim 50$ would require using boxes
of side $L_{\rm box} \sim 200$ and at least $5L_{\rm box}^3 = 4 \times 10^7$ solvent particles for both correctly implementing the hydrodynamic interactions
and avoiding significant finite-size effects \cite{Mussawisade2005}.
Still, based on previous evidence \cite{Pham2008,Nikolov2018,Camerin2018,Keidel2018} we do not expect that hydrodynamics will lead to {\it qualitatively} different results from those presented here.
We used a time step $\delta t = 0.005$, and  a friction $\gamma = 0.05$, which is high enough for good thermalization and low enough to prevent strong damping that would slow down the dynamics to time scales requiring a huge computational cost.

After generating them, equilibration of all the investigated systems (flexible and semiflexible chains, stars, disordered and diamond microgels) was performed in the limit of good solvent ($\phi=0$). To investigate the coarsening kinetics, the systems equilibrated at $\phi =0$ were quenched at infinite rate to different values of the solvent parameter $\phi$ well below the $\theta$-point, and coarsening was analyzed from the quenching instant ($t=0$) until the macromolecule reached its equilibrium fully collapsed state.

The linear chains were constructed as simple strings of beads and springs (adding the bending interaction in the case of the semiflexible chains).
We used $N=1600$ and 400 beads for the flexible and semiflexible chain, respectively.
The star polymers were constructed by linking flexible linear arms of 800 beads to a central bead.
Several protocols have been recently proposed to generate realistic models of microgels \cite{Gnan2017,Moreno2018,Minina2019,Ninarello2019,Rudyak2019,Rovigatti2019}
beyond the regular networks usually employed in the literature.
In our method we take inspiration from the synthesis of microgels in microfluidic cavities. This experimental route takes place via confinement in a droplet of  pre-existing linear polymer chains, and further inter- and intra-molecular irreversible association  \cite{Weitz2010}.
In our model each confined single chain consists of $N$ beads, and in each chain $N_{\rm r}$ of these beads are reactive (cross-linkable) groups, randomly distributed along the polymer backbone and with at least one inert group between consecutive reactive groups to prevent trivial cross-links.
In order to implement the cavity, a confining spherical, purely repulsive, LJ potential is applied to each monomer:
\begin{equation}
V_{\rm wall}(r_{\rm w}) = \left\{
\begin{array}{lll}
4\epsilon\left[\left(\frac{\sigma}{r_{\rm w}}\right)^{12}
- \left(\frac{\sigma}{r_{\rm w}}\right)^{6} +\frac{1}{4} \right] &\ &\  r_{\rm w} \leq 2^{1/6}\sigma \\
0 &\ &\ r_{\rm w} > 2^{1/6}\sigma
\end{array}
\right.
\end{equation}
where $r_{\rm w}$ represents the shortest distance from the monomer to the spherical wall.

We performed the cross-linking of $N_{\rm cha}=36$ chains in a cavity of radius $R_{\rm cav} =55$.
Each chain had $N=600$ beads, so that the finally generated microgel had 21600 beads.
The number of reactive beads in each chain was $N_{\rm r}=12$,
i.e., the fraction of cross-linked monomers in the microgel was $f = N_{\rm r}/N = 0.02$.
The number density used in the synthesis was $3N_{\rm cha}N/(4 \pi R^3_{\rm cav}) \approx  0.03$, which qualitatively
corresponds to experimental concentrations of about 30 mg/mL \cite{Kremer1990,Moreno2016JPCL,gonzalezburgos2018}.
After equilibration of the chains inside the cavity the cross-linking of the reactive groups was activated.
A permanent bond (modeled by the FENE potential) between two reactive groups was formed if:
(i) none of them was already bonded to another reactive group, and (ii) they were at a mutual capture distance $r < 1.3\sigma$.
A random choice was made in case of multiple candidates within the capture distance.
To speed up the late stage of the cross-linking process ($\leq 6 $ remaining unbonded reactive groups), a random pair was chosen from the unreacted groups
and an attractive harmonic interaction between the constituents of the pair was implemented, in order to approach them to the capture distance and form the bond.
After forming the bond the microgel was equilibrated and the procedure was repeated until full completion of the cross-linking. Then the cavity was removed to allow for swelling and equilibration of the obtained microgel. Cross-linking of 50 initial realizations, with the same former values
of $N$, $N_{\rm cha}$, $N_{\rm r}$ and  $R_{\rm cav}$ was perfomed, leading to microgels with the same number of monomers
and cross-links but topologically polydisperse \cite{Moreno2018}.
It was found that a large number of the cross-links (about 65 \%) occurred between reactive groups belonging to the same polymer chain,
forming loops. These kind of cross-links do not contribute to the
connectivity of the network and are elastically inactive. On the other hand in the diamond-like microgels no loops are present and all the cross-links
are elastically active. Therefore, for a fair comparison with the disordered ones, the diamond-like microgels were constructed with the same fraction
of cross-links (nodes) as the average number of intermolecular bonds in the disordered microgels, tuning the number of nodes so that the total mass
of both kinds of microgels was essentially the same \cite{Moreno2018}. Thus, we simulated a diamond-like microgel of $N=21615$ beads containing 78 nodes.
The diamond-microgel was generated by placing the cross-links in the nodes of a regular diamond network, and by connecting every pair of nearest-neighbour nodes through linear bead-spring rods \cite{olvera2011,winkler2014,winkler2016,Ahualli2017,Sean2018}. All the beads out of a sphere containing the selected $N$ beads were removed.

We computed the time-averaged asphericity $a$ \cite{rawdon2008} of each disordered microgel in the swollen state ($\phi=0$)
and obtained the corresponding distribution $P(a)$. For the analysis of the coarsening kinetics
we selected three disordered microgels, at the center and at the two extremes of the distribution of asphericities.
The corresponding radius of gyration, at $\phi=0$, of the selected disordered microgels is $R_{\rm g}= 49.3, 52.4$ and 64.0 for $a=0.02, 0.06$ and 0.14, respectively. In the rest of the paper these disordered microgels will be denoted in the figure legends as I-II-III from lower to higher asphericity.
The size of the other investigated systems at $\phi=0$ is  $R_{\rm g}=51.8$ (diamond-like microgel), 45.3 (3-arm star), 48.9 (12-arm star), 39.9 (flexible chain) and 26.4 (semiflexible chain).

\section{Coarsening Kinetics}\label{sec3}

The coarsening kinetics is characterized by a growing length scale. The quantitative characterization of such a length scale
can be easily affected by artifacts originating from the `structural noise' emerging during the macromolecular collapse (bridges that connect the clusters of monomers, small halls, protrusions, etc). To avoid these artifacts we introduce a smooth representation of the macromolecules through a coarse-grained density field. This method is based, originally, on the characterization of the growing length scales
in a coarsening binary Ising system \cite{Majumder2011} and was later applied in the continuous space to
a liquid-gas phase separating system \cite{Testard2011,Testard2014}. In the first case the thermal noise effects on the coarsened structure were removed by using a majority spin rule, i.e. by replacing each spin by the majority spin of its nearest neighbours. In the second case, the real particles were substituted by their local densities averaged over their nearest environment.
These averaging procedures smooth the interface corrugations, fill the smallest holes, and delete the smallest clusters in the coarsening structure. In this way the smooth density field avoids, in the calculation of the growing length scale, the effect of non-relevant minimal paths or artificial interruptions of long paths within the dense domains.
In our systems we construct our coarse-grained density field as follows:
\\
i) we divide the space into cubic cells of side $\delta$;
\\
ii) we define the local density (for each cell) by the number of monomers in a sphere:

\begin{equation}
\rho({\bf r}) = 3n({\bf r})/(4\pi r_{\rm c}^{3})
\end{equation}
where $r_{\rm c}$ is the cutoff radius of the sphere and $n({\bf r})$ is the
number of monomers at a distance $d$  $\leq$  $r_{\rm c}$ from the position $\bf r$ of the cell center;
\\
iii) we fix the value of the coarse-grained density
at ${\bf r}$ as a weighted average of the local density over the surrounding cells:
\begin{equation}
\bar{\rho}({\bf r}) = \frac{1}{8} \left[ 2\rho({\bf r}) +\sum_{\bf k}\rho({\bf r} + \delta{\bf k}) \right]
\end{equation}
where the sum is performed over the directions ${\bf k} \in \{(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)\}$.
The grid size $\delta$ and the cutoff radius $r_{\rm c}$ are chosen
in order to get a smooth density field while keeping sufficient spatial resolution in the representation of the real system.
This is achieved by using values $\delta \sim 0.5$ and $r_{\rm c} \sim 1$. In what follows we will present results for
$\delta = 0.5$ and for two specific values  $r_{\rm c} = 1.0$ and $r_{\rm c} = 1.2$.
By using a threshold value $\rho_{\rm min}$ in the density field construction the macromolecules can be seen as a coarsening biphasic system composed by 'empty' and 'filled' domains. These domains are identified
according to their low ($\bar{\rho}({\bf r}) \le \rho_{\rm min}$) or high ($\bar{\rho}({\bf r}) > \rho_{\rm min}$) local density, respectively.

\begin{figure}[pt!]
\centering
\includegraphics[width=1.0\textwidth]{snapshot-coarsening.pdf}
\caption{Snapshots of real (red beads) and density field (orange beads) coordinates of the flexible linear chain
(a), the 3-arm star (b) and the semiflexible linear chain (c). In all cases the collapse occurs at solvent parameter
$\phi=1.2$, and the parameters $r_{\rm c} = 1.0$ and $\rho_{\rm min} = 0.6$ are used to construct the density field.
The times for each panel are, from left to right:
$t = 5, 50, 126, 251, 316, 1000$ ((a) and (b)), and $t = 5, 50, 126, 398, 500, 8000$ (c).
}
\label{fig:snaps-field}
\end{figure}

Figure~\ref{fig:snaps-field}
shows snapshots of the real-coordinates (read beads) and the density field representations (orange beads) at different times during the collapse of a flexible chain (a), a 3-arm star (b) and a semiflexible chain (c), respectively.
In the real-coordinate representation all the momomers are displayed. In the density-field representation only the filled cells ($\bar{\rho}({\bf r}) > \rho_{\rm min}$) are shown.
All these snapshots  correspond
to $\phi=1.2$, a bad-solvent state well below the $\theta$-point. The final state reached at the end of the simulation (last columns on the right) is in all cases a fully collapsed macromolecule.
The collapse experienced by the flexible chain and the star polymer (panels (a) and (b)) begins with the formation of clusters of monomers along the chain, then those merge by withdrawing monomers from the bridges connecting them, in a longitudinal diffusion process.
In the stars this process also includes merging of clusters of different arms.
The collapse for the semiflexible chain (panel (c)) seems to be qualitatively different: the dense regions are better defined, as it becomes evident in the density field representation. At early times, when the conformations are close to those of $\phi=0$, there are many less clusters than in the flexible case. Indeed small fluctuations leading to local transient clustering are strongly hindered by bending stiffness. However, as time goes on and the effective monomer attraction starts to drive the collapse, the clusters quickly grow up before starting to merge into larger clusters. In contrast to the flexible systems where the chain or arm backbone can
still perform broad lateral fluctuations at early and intermediate times, in the semiflexible chains merging of the cluster proceeds along a quasi-rodlike structure during the whole process.
In the next subsections we quantify the observed similarities and differences, by analyzing the domain length distribution, the domain growth rate, as well as the dynamic correlations, during the coarsening process.

\subsection{Chord length distributions and domain growth}

The construction of the coarse-grained density field allows  to measure the distribution of the domain size from the obtained smooth biphasic structure.
We define an `interfacial cell' \cite{Moreno2018} in our system  as a filled cell with at least one adjacent empty cell.
In order to compute the characteristic lengths of the coarsened structures we use the definition of chord \cite{Testard2011,Testard2014}. This is a straight path, along one of the three $x,y,z$-directions of the grid, which is formed just by filled cells and whose two end cells are interfacial cells.
To compute the distribution of chord lengths of a given macromolecular configuration at a given time and solvent parameter $\phi$, we sampled all the existing chords by following all the possible paths along the three directions within a volume containing all the filled cells.
To improve statistics 5 random rotations of the former configuration were taken,
and the whole procedure was repeated over 5 independent realizations of the same macromolecule.

\begin{figure}[t!]
\centering
\includegraphics[width=.95\textwidth]{P_L.pdf}\hspace{0.3cm}
\caption{Normalized distribution of chord lengths at $\phi=1.2$ and different times during the collapse of the
disordered microgel of middle asphericity (a), the semiflexible chain (b), 12-arm star  (c) and 3-arm star  (d). The chord lengths are calculated for a density field
with $r_{\rm c}=1.0$ and $\rho_{\rm min}=0.6$.}
\label{fig:distlength}
\end{figure}

Figure \ref{fig:distlength} shows the normalized distributions of chord lengths, $P(L)$, at $\phi=1.2$ and different times during the collapse of the disordered microgel of middle asphericity (Figure \ref{fig:distlength}a), the semiflexible chain (\ref{fig:distlength}b), 12-arm star polymer (\ref{fig:distlength}c) and 3-arm star polymer (\ref{fig:distlength}d).
$P(L)$ shows the characteristic length distribution observed in coarsening systems  \cite{Levitz1998,Atsuko2003,Majumder2011,MajumderEPL2011,Testard2014,Majumder2017}. In particular, for earlier times $P(L)$ shows an exponential decay and extends over longer distances
as time increases, which is a consequence of the growth of the filled domains during the coarsening.
As expected, the exponential behavior saturates at long times when the fully collapsed state is reached. The flat plateau originates
from the equiprobable different straight paths that connect two points, at both sides of the outer interface, in the
fully collapsed state (where empty cells are absent). The drop from the plateau obviously reflects the finite size of the collapsed object .
Interestingly, the $P(L)$ of the semiflexible chain shows an approximate self-similar behavior over time, i.e., the decay is shifted to longer times but shows roughly the same shape, in contrast with the other systems for which the slope of the exponential is strongly time-dependent.
This is consistent with the growth and transport process of the nucleation centers anticipated in Figure \ref{fig:snaps-field}c. In fact, as shown there, the clusters are better defined and further apart, and they grow considerably before coalescing in a single cluster. Instead, in the flexible systems (panels
(a,b) of Figure \ref{fig:snaps-field}) the clusters grow and merge in a much more gradual way until the single globule is formed at late times.

We use the information on the former distributions to quantify the growing length scale characteristic of the coarsening process.
The mean domain length for a given time $t$ is obtained from the first momentum of the distribution $P$ at time $t$ , i.e.,
$L(t) = \int L' P(L'; t) dL'$.
We find the same qualitative behavior of $P(L)$ for several choices of $r_{\rm c}$ and $\rho_{\rm min}$, whereas quantitatively $L(t)$ depends on the specific
parameters used to construct the density field. Indeed, choosing, e.g.,  a lower value for $\rho_{\rm min}$ implies having a higher number of filled cells in the density field, and consequently a larger value of $L$. For this reason we quantify the domain growth rate in terms of the relative domain size, defined as:
 %
 \begin{equation}
C(t) = \dfrac{L(t)-L(0)}{L(\infty)-L(0)}
 \end{equation}
 %
In this way $C(t)$ represents a normalized mean chord length growing from zero at $t=0$ to 1 at late times in the collapsed state.
Figure \ref{fig:ct-densfield-microgel}  shows the domain growth for three different topologies of the disordered microgels (corresponding to a low, middle, and large value in the distribution of the asphericity parameter), in comparison with the diamond-like network. For each system the panel includes data of $C(t)$ for several selections of the parameters $r_{\rm c}$ and $\rho_{\rm min}$. The absolute times have been rescaled by the time $\tau_{0.5}$, defined as $C(\tau_{0.5}) =0.5$. As can be seen for all systems, after this time rescaling the different data sets of $C(t)$ nicely overlap. This demonstrates that the density field approach is consistent, as it provides a time dependence of $C(t)$ that is independent of any reasonable choice of the parameters defining the density field.
Furthermore a good overlap is observed also for different values of the solvent quality parameter ($\phi = 1.2$ and $\phi = 1.5$), indicating that the coarsening kinetics follows an effective time-temperature superposition principle.

\begin{figure}[t!]
\centering
\includegraphics[width=.95\textwidth]{plot-Ct-microgels.png}
\caption{Relative domain size $C(t)$ over time during
the collapse at $\phi = 1.2$ and $\phi = 1.5$ of three disordered microgels of small (a),
middle (b) and large asphericity (c) and of a diamond microgel (d). Different data sets correspond to different values of the parameters $r_{\rm c}$
and $\rho_{\rm min}$ used to construct the density field. The time $t$ is normalized by $\tau_{0.5}$,
defined as the time when $C = 0.5$. Straight lines in all panels are fits to a
power-law time dependence. Exponents are indicated. In panels (b) and (d) we show typical simulation snapshots at early and intermediate times.}
\label{fig:ct-densfield-microgel}
\end{figure}

As can be seen in Figure \ref{fig:ct-densfield-microgel}, the function $C(t)$ reveals that the coarsening length scale grows by following a sublinear power law $\sim t^x$. The exponents are $x \gtrsim 0.6$ for all the microgels, irrespective of the microstructural degree of disorder. The exponents obtained here for $f=0.02$ are slightly smaller than those found in our previous work \cite{Moreno2018} ($x \sim 0.7$) for a much higher degree of cross-linking ($f=0.1$).
We believe that the slightly higher value of the exponent for $f=0.1$ is just an artifact originating from the proximity of the cross-links
in such a system. Namely, since the higher density of monomers around a cross-link enhances its propensity to become a nucleating center, the much smaller distance between cross-links in the case $f=0.1$ than in $f=0.02$ may effectively accelerate the coalescence of the growing clusters,
leading to the observed faster domain growth.

Our previous work in microgels with higher $f$ suggested that the diamond network exhibited the same power-law in $C(t)$ that the disordered networks. However, for the diamond network this behaviour was developed only in a narrow time window and the results were not conclusive \cite{Moreno2018}.
The results in Figure \ref{fig:ct-densfield-microgel}d for $f=0.02$ clearly confirm the power law regime over 3 time decades.
These data also confirm the acceleration of the domain growth in the late stage of the collapse of the diamond network \cite{Moreno2018}, which is at most a marginal effect
in the disordered networks, where the power-law continues and is a good approximation until the collapsed state is reached and the ultimate plateau emerges. The strong late acceleration in the diamond network is tentatively related to the regular spatial distribution of the cross-links acting as preferential nucleating centers,
which leads to a more homogeneous collapse until all holes in the structure vanish in a short time window (compare snapshots in panels (b) and (d)),
and the filled paths in the density field experience a late sudden growth.

\begin{figure}[t!]
\centering
\includegraphics[width=.9\textwidth]{total_Ct.pdf}\hspace{0.01cm}
\vspace{0.1cm}
\caption{Relative domain size $C(t)$ over time during
the collapse at $\phi = 1.2$ of 3-arm stars (a),
12-arm stars (b), flexible chains (c) and semiflexible chains (d). Different data sets correspond to different values of the parameters $r_{\rm c}$
and $\rho_{\rm min}$ used to construct the density field.
The time $t$ is normalized by $\tau_{0.5}$, defined as in Figure~\ref{fig:ct-densfield-microgel}. Straight lines in all panels are fits to
a power-law time dependence (exponents are indicated).
}
\label{fig:ct-densfield-total}
\end{figure}

The results presented above confirm the common scaling law $C(t) \sim t^{0.6}$ for the domain growth in collapsing microgels, irrespective of their topology (disordered or regular networks). In order to search for a more general scenario, we investigate the scaling properties of the coarsening kinetics in other very different architectures, as stars with different number of arms, and flexible and semiflexible linear chains. We follow the same procedure as for the microgels, constructing the density field and analyzing the domain growth in this representation of the macromolecule. Figure~\ref{fig:ct-densfield-total} shows the corresponding results of $C(t)$ for the stars and linear chains.
We find that the systems without bending stiffness (the stars and the flexible linear chains) essentially show the same power-law
$C(t) \sim t^{0.56}$, with an exponent that is just slightly smaller than those found for the microgels. The only investigated system with bending stiffness, i.e., the semiflexible linear chains, shows a different power-law, with a clearly higher exponent $C(t) \sim t^{0.8}$.
In a first approximation it is possible to venture that the flexible nanoparticles as microgels and star polymers act in a similar fashion during most of the coarsening process: at distances shorter than the arm length and the mesh size, the arms and strands behave like linear flexible chains. Microgels show an exponent of $\sim 0.6$ for the domain growth, slightly larger than the value $\sim 0.56$ found for stars and single flexible chains. This difference might be connected with the higher concentration of nucleating centers in the microgels (due to the presence of the cross-links). Comparatively, at the local scale the semiflexible chain is influenced by the presence of the bending, which reduces its lateral fluctuations.
This reduction may promote a faster and more homogeneous aggregation of mass along the chain contour, leading to the observed steeper growth of the coarsening length scale.

\subsection{Cluster analysis}

\begin{figure}[t!]
\centering
\includegraphics[width=.6\textwidth]{plot-nclus-rclus.pdf}
\caption{Mean cluster radius vs. mean cluster population at $\phi = 1.2$, for clusters of filled cells
in the density field representation with parameters $r_{\rm c} = 1.0$ and $\rho_{\rm min} = 0.6$.
Data are shown for the disordered microgel of middle asphericity, the diamond network, the 3-arm star and the
flexible and semiflexible linear chains. The arrows indicate the approximate cluster population
at the saturation point prior the formation of the late plateau in  $N_{\rm clus}(t)$.
Lines are power-laws, exponents are indicated.}
\label{fig:n-r-clusfield}
\end{figure}

In this subsection we shed some ligth on the microscopic origin of the different scaling exponents for the
domain growth in the flexible systems and the semiflexible linear chains. We first analyze
whether stiffness plays a role on the conformational properties of the dense regions
that progressively emerge in the coarsening structure. Namely we analyze their fractal behaviour,
i.e., the power-law dependence between the size and the mass of such dense regions.
This can be done
by defining clusters of filled cells ($\rho > \rho_{\rm min}$) in the density field representation,
and determining the relation between the mean radius $R_{\rm clus}$ and mean population $N_{\rm clus}$ of the clusters.
Two filled cells
belong to a same cluster if they are adjacent, i.e., if they are connected by a vector
$\delta{\bf k} \in \{(\pm \delta,0,0),(0,\pm \delta,0),(0,0,\pm \delta)\}$.
The mean radius $R_{\rm clus}$ is just obtained as $\langle R_{\rm g}^2 \rangle^{1/2}$, with $R_{\rm g}$ the
radius of gyration of the cluster. At each time we average the former quantities over all the clusters,
obtaining the time-dependence $R_{\rm clus}(t)$ and  $N_{\rm clus}(t)$. Substitution of time provides a univoque relation $R_{\rm clus}(N_{\rm clus})$
for the fractal behavior of the clusters. Figure~\ref{fig:n-r-clusfield} shows $R_{\rm clus}(N_{\rm clus})$
for several representative flexible systems (see legend) and for the semiflexible chains. All the systems show two
different power-laws $R_{\rm clus} \sim N_{\rm clus}^\nu$ at both sides of a crossover value $N^{\ast} \sim 400$. For small clusters ($N_{\rm clus} < N^{\ast}$) an exponent $\nu \approx 0.5$ is found, which indicates approximate Gaussian statistics
of the cluster conformations. Big clusters, which are formed in the late stage of the coarsening process,
show an exponent $\nu \approx 0.25$, reflecting the expected compact structures that progressively merge until reaching the equilibrium globule.

\begin{figure}[t!]
\centering
\includegraphics[width=.9\textwidth]{plot-clusdensfield-sizeradius-microgels-new.pdf}
\caption{Normalized population (a,b) and radius (c,d) of the clusters in the density field representation of the
disordered microgel with middle asphericity (a,c) and the diamond microgel (b,d). Data correspond to a solvent quality
parameter $\phi= 1.2$, and a density field representation with parameters $r_{\rm c} = 1.0$ and $\rho_{\rm min} = 0.6$.
The arrows indicate the time scale for which $N_{\rm clus}(t) = N^{\ast} \approx 400$.
Times are rescaled by $\tau_{0.5}$ as defined in Figure~\ref{fig:ct-densfield-microgel}.}
\label{fig:nt-rt-clusfield-micro}
\end{figure}
\begin{figure}[t!]
\centering
\includegraphics[width=.9\textwidth]{plot-clusdensfield-sizeradius-linear.pdf}
\caption{Normalized population (a,b) and radius (c,d) of the clusters in the density field representation of the
flexible (a,c) and semiflexible linear chain (b,d). Data correspond to a solvent quality
parameter $\phi= 1.2$. Values ($r_{c} = 1.0$, $\rho_{\rm min} = 0.6$) and ($r_{c} = 1.2$, $\rho_{\rm min} = 0.35$) are
used to construct the density field.
The arrows indicate the time scale for which $N_{\rm clus}(t) = N^{\ast} \approx 400$.
Times are rescaled by $\tau_{0.5}$ as defined in Figure~\ref{fig:ct-densfield-microgel}.}
\label{fig:nt-rt-clusfield-linear}
\end{figure}

Figures~\ref{fig:nt-rt-clusfield-micro} and \ref{fig:nt-rt-clusfield-linear} show the time dependence of the cluster population and size for the microgels and the flexible and semiflexible linear chains.
In analogy with the normalized function for the domain growth $C(t)$, we define the normalized population and size of the
clusters as $n(t) = (N_{\rm clus}(t)- N_{\rm clus}(0))/(N_{\rm clus}(\infty)- N_{\rm clus}(0))$ and
$r(t) = (R_{\rm clus}(t)- R_{\rm clus}(0))/(R_{\rm clus}(\infty)- R_{\rm clus}(0))$, respectively, so that both $n(t)$
and $r(t)$ grow from 0 to 1.
The arrows in all panels of Figures~\ref{fig:nt-rt-clusfield-micro} and \ref{fig:nt-rt-clusfield-linear}
indicate the time scale for which $N_{\rm clus} = N^{\ast}$.
As shown in Figure~\ref{fig:n-r-clusfield}, the clusters in all systems have the same fractal behaviour.
However, their evolution in time (Figures~\ref{fig:nt-rt-clusfield-micro} and \ref{fig:nt-rt-clusfield-linear}) can
strongly depend on the system, namely there are clear differences between the flexible and semiflexible systems.
In other words, coarsening leads to mass aggregation into the same kind of clusters for all the systems, but the rate of mass aggregation
can be significantly affected by bending stiffness.
For the time scales where the mean population of the clusters is smaller than $N^{\ast}$ (clusters are still relatively small
and approximately Gaussian), we find a power-law $n(t) \sim t^{\beta}$, with exponent $\beta \approx 1.2$ for the flexible systems, and a much higher
value $\beta = 1.6$ for the semiflexible chains.
Not surprisingly, the growing length scale of the clusters follows a power-law
$r(t) \sim t^{\gamma}$ that is compatible with the same exponents found in $C(t)$ for the growing domain size
(see Figures~\ref{fig:ct-densfield-microgel} and \ref{fig:ct-densfield-total})
---though clusters and domains represent different concepts, growing length scales in the system should follow the same scaling.
The effective exponents $\beta$ and $\gamma$ for $n(t)$ and $r(t)$ in all cases are related as $\gamma/\beta \approx 0.5$.
This indeed reflects the scaling $R_{\rm clus} \sim N_{\rm clus}^{0.5}$ found in Figure~\ref{fig:n-r-clusfield}.
In the case of the flexible and semiflexible chains, and in the stars (not shown), no apparent or significant change
is observed in the power-law behavior of $n(t)$ and $r(t)$ at later times when larger clusters are formed, until the former quantities finally saturate when all clusters have merged into a single one and the systems approach the equilibrium fully collapsed state. An overshoot in $n(t)$ and $r(t)$ before saturation is however found in the microgels. This effect is specially pronounced in the diamond network,
and not surprisingly occurs roughly in the same time window as for the normalized domain length $C(t)$ (Figure~\ref{fig:ct-densfield-microgel}d).
Since $r \sim n^{\nu}$ with $\nu \leq 0.5$, the overshoot is much less pronounced in the growing size of the clusters than in the mass aggregation.

In summary, the analysis of clusters reveals that there are no significant differences in the conformations
of the dense regions that are formed during the coarsening of the different macromolecular architectures, but on the rate at which they are formed, which is essentially the same in flexible systems but can be accelerated through bending stiffness in semiflexible systems.

\subsection{Dynamic correlations}

\begin{figure}[t!]
\centering
\includegraphics[width=1\textwidth]{plot-correl-microgels-lowden.pdf}
\caption{Domain self-correlation function $P_{\rm s}(t)$ for three disordered microgels of small (a),
middle (b) and large asphericity (c), and for the diamond-like microgel (d), during the collapse at $\phi = 1.2$. The results displayed correspond to the couples of parameters ($r_{c} = 1.0$, $\rho_{\rm min} = 0.6$) and ($r_{c} = 1.2$, $\rho_{\rm min} = 0.35$)  defining the density field. In each data set the time has been rescaled by $\tau_{0.5}$ , with $\tau_{0.5}$
defined as in Figure \ref{fig:ct-densfield-microgel}. Lines are fits to power-law dependence. Exponents are indicated.}
\label{fig:correl-microgels}
\end{figure}

A well-known feature of phase separating and coarsening systems is the scaling of the dynamic correlations with the growing length scale. To test this possibility
we first define a `spin' self-correlation function $P_{\rm s}(t)$ for the density-field representation. The function is defined as \cite{Christiansen2017}:
\begin{equation}
P_{\rm s}(t) = \langle S(t)S(0) \rangle - \langle S(t)\rangle \langle S(0) \rangle
\end{equation}
where the variable $S$ is computed for each cell of the density field, and it is assigned a value 0 or 1 if the cell is empty ($\rho \leq \rho_{\rm min}$)
or filled ($\rho > \rho_{\rm min}$), respectively.

\begin{figure}[t!]
\centering
\includegraphics[width=.49\textwidth]{plot-chordlength-correl-microgels-lowden.pdf}\hspace{0.1cm}
\includegraphics[width=.49\textwidth]{plot-chordlength-correl-stars.pdf}
\caption{Domain self-correlation function $P_{\rm s}(t)$ vs. domain size $L(t)$ at the time $t$  during the
collapse, for several values of the solvent quality parameter $\phi$ and the parameters ($r_{\rm c}$, $\rho_{\rm min}$) defining the density field
(see legends). Panel (a) shows data for the diamond microgel and three disordered microgels.
Panel (b) shows data for the 3-arm and 12-arm stars.
The lines in both panels are fits to a power-law $P_{\rm s} \sim L^{-1.1}$.}
\label{fig:PvsL}
\end{figure}

In Figure \ref{fig:correl-microgels} we show the correlation function $P_{\rm s}(t)$ for the three selected disordered microgels (low, middle and high asphericity)
and for the diamond-like microgel, at solvent parameter $\Phi = 1.2$.
Data are shown for two different couples of the parameters ($r_{\rm c}$, $\rho_{\rm min}$) used to define the density field, and the times are rescaled by $\tau_{0.5}$ as defined above.
Again, the good overlap of the different data sets confirms the consistency of the density field representation for characterizing the coarsening.
The underlying complex dynamics associated to the coarsening process is reflected by the non-exponential decay of the dynamic correlations, which indeed follows a power-law,
$P_{\rm s}(t) \sim t^{-y}$. Results for all the microgels are consistent with an exponent $y \gtrsim 0.5$, irrespective of the network topology.
This is analogous to our results for the domain growth ($C(t)$, Figure~\ref{fig:ct-densfield-microgel}) in the microgels, for which the specific network topology
has no significant effect on the observed power-law.

Moreover, from the computation of $P_{\rm s}(t)$ and the average domain size $L(t)$, we can establish a direct relation between the dynamic self-correlations
and the growing characteristic length scale, just by taking at each time $t$ the corresponding values of $P_{\rm s}$ and $L$.
Figure~\ref{fig:PvsL} shows the results for this relation in the diamond-like microgel and in the three selected disordered microgels (panel (a))
and in the 3-arm and 12-arm stars (panel (b)).
The data are computed for several values of the solvent parameter $\phi$ and the parameters ($r_{\rm c}$, $\rho_{\rm min}$) used to define the density field.
Time-temperature superposition is also confirmed for this scaling relation, and
data for all flexible systems are consistent with a common power law $P_{\rm s} \sim L^{-1.1}$.
Unfortunately we could not confirm this observation for the semiflexible case, for which the power-law regime in $P_{\rm s}(L)$ developed
over a short range --- note that the size of the investigated semiflexible chains is much smaller than for the flexible systems --- and had very poor statistics. Confirmation would require to simulate much longer semiflexible chains, which would highly complicate the analysis.
Tests for much longer chains showed that they tend to collapse into rods and in an
extremely heterogeneous fashion (visiting several long-living intermediate conformations prior to the equilibrium state).

\section{Conclusions}\label{sec5}

By means of simulations, we have investigated the coarsening kinetics emerging during the collapse of several macromolecular architectures in bad solvent:
microgels with realistic (disordered) and with ideal regular (diamond) networks, star polymers of 3 and 12 arms, and linear chains.
We have also investigated the effect
of bending stiffness on the coarsening kinetics by simulating the collapse of semiflexible linear chains.
In order to remove fast fluctuations that can lead to artifacts in the characterization of the growing length scale, we have made use of a smooth
density field representation of the macromolecules. The domain growth during the coarsening follows a power-law behaviour that is independent of any
reasonable selection of the parameters used to construct the density field (grid size and threshold for defining dense regions).
The scaling behaviour is independent of the solvent quality parameter, in analogy to time-temperature superposition.
All flexible systems show approximately the same exponent for the time dependence of the coarsening length scale ($\sim t^{0.6}$). An overshoot is found in the diamond networks
in the late stage of the coarsening, which can be tentatively assigned to the regular distribution of the nodes acting as preferential
nucleating centers and their roughly simultaneous merging when the network is close to the collapsed globular state.
The length scale of coarsening shows a clearly steeper growth in the semiflexible chains ($\sim t^{0.8}$).
To elucidate the origin of this difference, we have analyzed the clusters of dense regions
formed during the coarsening in the density field representation. The clusters in all systems (flexible or semiflexible) show the same fractal behaviour, i.e.,
their size scales with the mass following the same power-laws. This suggests that
the faster growing length scale in the semiflexible chains just originates from a faster mass diffusion along the chain contour,
and not from a distinct structural feature of the aggregates formed during the coarsening process.
We have analyzed dynamic correlations, and in analogy with critical phenomena, investigated their dependence on the growing length scale $L$.
We find an apparent common power-law dependence of the correlations ($\sim L^{-1.1}$) for all the flexible systems.

This work reports, to the best of our knowledge, the first comparison of the coarsening kinetics in bad solvent
for a broad range of macromolecular architectures, and explores the role of chain stiffness.
As such, it provides a general physical scenario, and a valuable set of results for future theoretical developments
in this, still scarcely studied, fundamental problem with potential applications to e.g., protein folding.

\section*{Acknowledgements}
This research was funded by projects
PGC2018-094548-B-I00 (MCIU/AEI/FEDER, UE) and IT-1175-19 (GV, Spain).


text:
\section{Introduction}
\subsection{The conformal modulus of a doubly-connected domain}
Let $D$ be a doubly-con\-nec\-ted planar domain with nondegenerate boundary
components, then $D$ is conformally equivalent to an annulus $\{r< |z|< R\}$ (see, e.g., \cite[Chap.~V, \S~1]{goluzin1969geometric}, \cite[Sect.~3]{kuhnau2005conformal}), and its conformal modulus $m(D)$ is the number defined as follows:
$$
m(D):=\dfrac{1}{2\pi}\log\dfrac{R}{r}\,.
$$
Another approach to define the conformal modulus  is in using  the concept of extremal length of curve family (see, e.g., \cite[Chap.~1, Sect.~D]{ahlfors2006lectures}, \cite[Chap.~3,
Sect.~11]{kuhnau2005conformal}). Specifically,
$$
m(D)=\lambda(\Gamma)
$$
where $\lambda(\Gamma)$ is the extremal length of the family $\Gamma$ of curves in the domain $D$ which join the boundary components of
$D$. Moreover,
$$
m(D)={1}/{\lambda(\Gamma')}
$$
where $\Gamma'$ is the family of curves  in $D$ separating its boundary components.

Conformal moduli of doubly-connected domains are also closely related to the conformal moduli of quadrilaterals (see, e.g., \cite{ahlfors2006lectures,kuhnau2005conformal,papamichael2010numerical}). A quadrilateral
$(Q;z_1,z_2,z_3,z_4)$ is a Jordan domain $Q$ with four fixed
points (vertices) $z_k$, $1\le k\le 4$,  on its boundary; the increasing of the index $k$ corresponds to the positive bypass of the boundary $\partial Q$.

For a given quadrilateral
$\mbox{\boldmath$Q$}=(Q;z_1,z_2,z_3,z_4)$, there exists a
conformal mapping of $Q$ onto a quadrilateral $(\Pi;0, 1,1+hi,hi)$, where $\Pi=[0,1]\times[0,h]$,
$h>0$, keeping the vertices. By definition,  the conformal modulus
$m(\mbox{\boldmath$Q$})$ of $\mbox{\boldmath$Q$}$ is equal to $h$.
The conformal modulus $m(\mbox{\boldmath$Q$})$ coincides with the
extremal length of the family of curves connecting in $Q$ the boundary arcs
$z_1z_2$ and $z_3z_4$; it is also the inverse to the extremal
length of the family of curves in $Q$ joining the arcs $z_2z_3$
and $z_4z_1$. The quadrilateral
$\mbox{\boldmath$Q$}^*=(Q;z_2,z_3,z_4,z_1)$ is called conjugate to
$\mbox{\boldmath$Q$}$. It is clear that
$m(\mbox{\boldmath$Q$}^*)=m^{-1}(\mbox{\boldmath$Q$})$.
If $\mbox{\boldmath$Q$}^c=(Q^c;z_4,z_3,z_2,z_1)$, where $Q^c=\overline{\mathbb{C}}\setminus Q$, then the conformal modulus of $\mbox{\boldmath$Q$}^c$ is called the exterior modulus (see, e.g., \cite{vuorinen2013exterior}).

Over the years, there have been obtained many results on the
conformal moduli of doubly-connected domains, as well as on moduli of quadrilaterals  (see, e.g.,  \cite{dautova2018,dautova2019,dautova2021conformal, hakula2021conformal,hrv1,hrv2,kuhnau2005conformal,nass,nasser2021circular,nasyrov2015riemann,nvgiang2021,nasyrov2021moduli,vuorinen2013exterior}).

Conformal moduli of doubly-connected domains and quadrilaterals are conformally invariant. They are also quasiinvariant under quasiconformal mappings (see, e.g., \cite[Chap.~2, Sect.~A]{ahlfors2006lectures}, \cite[Chap.~3, Sect.~22]{kuhnau2005conformal}). If $f$ is an $H$-quasiconformal mapping, then
$$
\dfrac{1}{H}m(\Omega)\leq m(f(\Omega))\leq H m(\Omega)
$$
where $\Omega$ is a doubly-connected domain. The same is valid for conformal moduli of quadrilaterals.

Prof. M.~Vuorinen posed the problem of investigating the distortion of the conformal moduli of doubly-connected  domains under the action of the $H$-quasiconformal mappings,
\begin{equation}\label{cttuabg}
f_H:x+iy\mapsto Hx+iy\quad (H>1);
\end{equation}
that is the stretching along the abscissa axis with coefficient $H$. In particular, it is interesting to study its asymptotic behavior as $H\to+\infty$ (see, e.g., \cite{nasyrov2015riemann}).

In \cite{dautova2018}, \cite{dautova2019}, \cite{nasyrov2015riemann}, \cite{nvgiang2021} the results on solving the Vuorinen problem were obtained for the cases where $D$ is a rectangular frame, a symmetric or non-symmetric bounded doubly-connected domain  or unbounded doubly-connected domains, symmetric with respect to the coordinate axes. Here we study the problem by Vuorinen for the case of an unbounded doubly-connected domain of sufficiently arbitrary form, non-symmetric with respect to the coordinate axes.
\medskip

\subsection{The main result}

Let $\mathfrak{S}$ be the set of all domains $\Omega$ which satisfy the following conditions:

(i) $\Omega$ is an unbounded doubly-connected domain;

(ii) the boundary of $\Omega$ consists of the graphs of continuous functions:  functions  $g_1$, $f_1$ are given on a segment $[a, b]$ and $f_1(x)\le g_1(x)$, $x\in[a, b]$; the functions  $g_2$, $f_2$ are given on a segment $[c, d]\subset[a,b]$ and $g_2(x)\le f_2(x)<f_1(x)$, $x\in[c, d]$; moreover, $f_1(a)=g_1(a), f_1(b)=g_1(b)$; $f_2(c)=g_2(c), f_2(d)=g_2(d)$ (see. Fig.~\ref{f6}).

\begin{figure}[ht] \centering
\includegraphics[width=3.5 in]{1.eps}
\caption{A domain $\Omega$ from the set $\mathfrak{S}$.}\label{f6}
\end{figure}

We will describe the behavior of the conformal modulus of the domain $\Omega_H$ obtained from  $\Omega\in\mathfrak{S}$ by the stretching $f_H$, given by \eqref{cttuabg}, as $H\to +\infty$.

\begin{theorem1}\label{main}
If $\Omega\in\mathfrak{S}$ and $\Omega_H=f_H(\Omega)$, then
\begin{equation}\label{ctulmd1}
m(\Omega_H)\sim \dfrac{1}{\gamma H}\,, \quad H\to\infty, \quad
\mbox{\rm where }\ \gamma
=\int\limits_c^d\dfrac{dx}{f_1(x)-f_2(x)}\,.
\end{equation}
Moreover,
\begin{equation}\label{bdtmodulOmegaH}
m(\Omega_H)\leq\dfrac{1}{\gamma H}\quad  \mbox{\rm for every}\quad
H>0.
\end{equation}
\end{theorem1}

\section{Local versions of the Rad\'o convergence theorem}

In this section we will obtain a local version of the Rad\'o convergence theorem and its generalization for a sequence of doubly-connected domains with degenerate moduli (see, e.g., \cite[Theorem~3]{dautova2019}); these results will be used in the proof of Theorem~\ref{main}.

As it is well known, every conformal mapping $f$ between Jordan domains can be extended, by continuity, up to a homeomorphism of the closures of the domains (see, e.g., \cite[Chap.~2, Sect.~3]{goluzin1969geometric}); for simplicity of notation, we will also denote this homeomorphism by~$f$.

Firstly, we recall the classical Rad\'o's theorem on uniform convergence of conformal mappings. Let $E$ denote the unit disk $\{|z|<1\}$ (see, e.g., \cite[Chap.~2, Sect.~5]{goluzin1969geometric}).

\begin{theorem1}[Rad\'o]\label{thr_rado_ge}
Let $D_n$ be a sequence of Jordan domains, bounded
by curves $\Lambda_n$, which converges to a kernel with respect to a point $z_0$, and let this kernel be a Jordan domain $D$ bounded by a curve $\Lambda$. Then a sequence of continuous functions $f_n : \overline{E} \to \overline{D}_n$, mapping conformally ${E}$ onto $D_n$ and normalized by the conditions $f_n(0) = z_0$, $f_n' (0) > 0$, converges uniformly to the continuous function $f\,:\,\overline{E} \to \overline{D}$, $f(0) = z_0, f'(0) > 0$ mapping conformally $E$ onto $D$ if and only if there exist homeomorphisms between the curves $\Lambda_n$ and $\Lambda$ such that, for any $\varepsilon > 0$, there is an $N$ such that the distance between the respective points of the curves $\Lambda_n$ and $\Lambda$ is smaller than $\varepsilon$ whenever $n\geq N$.
\end{theorem1}

\begin{remark1}\label{rem1}
It is easy to prove that, under the assumptions of Theorem~\ref{thr_rado_ge},
$f^{-1}_n\rightrightarrows f^{-1}$ on every compact $K$ such that $K\subset \overline{D}_n$ for sufficiently large $n$. In particular, $f^{-1}_n\rightrightarrows f^{-1}$ on $\bigcap_{n=1}^\infty \overline{D}_n$.
\end{remark1}

First we will establish a local version of Theorem~\ref{thr_rado_ge}.

\begin{figure}[ht] \centering
\includegraphics[width=4.7 in]{2.eps}
\caption{The sequence of domains $D_n$ and its kernel.}\label{f1}
\end{figure}

\begin{theorem1}\label{thr_local_rado}
Assume that $D_n$ is a sequence of Jordan domains, bounded by curves $\Gamma_n$, which converges to a kernel with respect to a point $z_0$, and this kernel is a Jordan domain $D$ bounded by a curve $\Gamma$. Let $f_n:\overline{E}\to \overline{D}_n$ be a sequence of continuous functions, mapping conformally ${E}$ onto $D_n$ and normalized by the conditions $f_n(0) = z_0$, $
f_n'(0)>0,$ and $f:\overline{E}\to \overline{D}$ be a continuous function mapping conformally  ${E}$ onto $D$ and  satisfying the conditions $f(0) = z_0$, $f'(0)>0$. Let there exist a sequence of subarcs  $\gamma_n\subset \Gamma_n$ with endpoints $a_n$, $b_n$,  a sequence of continuous mappings $\psi_n: \gamma_n\to \Gamma$ and a point $w_0\in \Gamma$ such that:

$(i)$\ \textit{for every $n\ge 1$, $\psi_n$ is a homeomorphism of $\gamma_n\setminus\{a_n,b_n\}$ onto} $\Gamma\setminus \{w_0\}$;

$(ii)$\ for every $n\ge 1$, $\psi_n(a_n)=\psi_n(b_n)=w_0$;

$(iii)$\ for every $\varepsilon> 0$ there exists a number $N\geq 0$ such that
\begin{equation*}
\left|\psi_n(z)-z\right|<\varepsilon\quad \forall n\geq N\  \forall z\in\gamma_n;
\end{equation*}

$(iv)$\  there exists a sequence of crosscuts $\alpha_n$ of domains $D_n$ with endpoints $a_n$ and $b_n$, diameters of which tend to $0$ as $n\to\infty$.

Then $f_n$ converges locally uniformly to $f$ on $\overline{E}\setminus\{\zeta_0\}$, $\zeta_0=f^{-1}(w_0)$, and $f^{-1}_n$ converges uniformly to $f^{-1}$ on every compact subset of $\bigcap_{n=1}^\infty \overline{D}_n\setminus \{w_0\}$.
\end{theorem1}
\begin{proof}
I) It is clear that the sequence of Jordan domains $D'_n$, bounded with $\gamma_n\cup\alpha_n$, converges to $D$ as to a kernel with respect to $z_0$ (Fig.~\ref{f1}). Without loss of generality, we can assume that  $z_0\in D_n$, $n\ge 1$. Consider a sequence of conformal mappings
$$
\widetilde {f}_n: E \to D'_n
$$
with normalization $\widetilde {f}_n (0) = z_0 $, $\widetilde {f}_n '( 0)> 0 $.

\begin{lemma1}\label{wfn}
The sequence $\widetilde{f}_n$ converges uniformly to $f$ on the closed unit disk~$\overline{E}$.
\end{lemma1}

\begin{proof}
Using the conditions of the theorem, we can easily show that there exist homeomorphisms
$$
\widetilde {\psi}_n: \alpha_n \cup \gamma_n \to \Gamma,
$$
satisfying the Rad\'o condition, i.e. such that for any $\varepsilon> 0$ there exists $N\geq 0$ such that
\begin{equation*}
|\widetilde{\psi}_n(z)-z|<\varepsilon, \quad z\in\alpha_n \cup\gamma_n,
\end{equation*}
whenever $n\geq N$. By Theorem~\ref{thr_rado_ge}, we conclude that $\widetilde{f}_n$ converges uniformly to $f$ on~$\overline{E}$.
\end{proof}
\medskip

II) Here we will show that the preimages of the arcs $\alpha_n$ under the mappings $\widetilde {f}_n$ are arcs on the unit circle which are contracted to a point as $n\to\infty$.

According to the principle of the boundary correspondence under conformal mappings, every $\beta_n:=\widetilde{f}_n^{-1}(\alpha_n)$ is an arc on the unit circle, with the end points $\zeta_{n1}^a:=\widetilde{f}_n^{-1}(a_n)$ and $\zeta_{n1}^b:=\widetilde{f}_n^{-1}(b_n)$. Similarly, the points $\zeta^a_{n2}:=f^{-1}_n(a_n)$ and $\zeta^b_{n2}:=f^{-1}_n(b_n)$ are also on the boundaries of the domains $f^{-1}_n(D'_n)$.

\begin{lemma1}\label{lem1}
The diameters of $\beta_n$ tend to $0$ as $n\to\infty$. Moreover, these arcs are contracted to the  point $\zeta_0$, as $n\to\infty$.
\end{lemma1}

\begin{proof}
First we recall that by (iv) the diameters of the arcs $\alpha_n$ tend to 0, and by (ii) and (iii) their endpoints $a_n$ and $b_n$ converge to $w_0$. Therefore, $\alpha_n$ are contracted to the point $w_0$, as $n\to\infty$. Assume that $\mbox{\rm diam}(\beta_n)\not\to 0$, as $n\to\infty$.  Then there exists a subsequence $\beta_{n_k}$  of the sequence $\beta_{n}$ such that the  $\beta_{n_k}$ tends to some nondegenerate subarc $\beta$ of the unit circle. By Lemma~\ref{wfn}, we have  $\widetilde{f}_{n_k}\rightrightarrows f$ on $\overline{E}$, therefore,  $f\equiv w_0$ on $\beta$. By the boundary uniqueness theorem (see, e.g., \cite[Chap.~2, Sect.~3, Lemma~3]{goluzin1969geometric}) this implies that  $f\equiv w_0$  on the unit disk. The contradiction proves that $\mbox{\rm diam}(\beta_n)\to 0$, as $n\to\infty$.

If $\beta_n$ are not contracted to the point $\zeta_0$, then there exist a subsequence $\beta_{n_k}$ of $\beta_n$ converging to some $\zeta_1\neq \zeta_0$. But $\widetilde{f}_{n_k}$ converges uniformly to $f$ on the closed unit disk. For every sequence $\zeta_{n_k}\in\beta_{n_k}$ we have $\zeta_{n_k}\to \zeta_1$ and $f_{n_k}(\zeta_{n_k})\to w_0$, $k\to \infty$, since $f_{n_k}(\zeta_{n_k})\in \alpha_{n_k}$. But, because of uniform convergence,  $f(\zeta_1)=\lim_{k\to\infty} f_{n_k}(\zeta_{n_k})=w_0$, therefore, $\zeta_1=\zeta_0$. The contradiction proves the lemma.
\end{proof}

III) Consider the arcs $\sigma_n:=f_n^{-1}(\alpha_n)$. Denote by $\widetilde{D}_n$ the extended complex plane with the excluded arc $\beta_n$. Let $\widetilde{\sigma}_n$ be the arc symmetric to $\sigma_n$ with respect to the unit circle. Let $\widetilde{G}_n$ be the unbounded domain in the extended complex plane with the boundary $\sigma_n\cup\widetilde{\sigma}_n$.

\begin{lemma1}\label{lem4}
The sequences of domains $\widetilde{D}_n$ and $\widetilde{G}_n$ converge to the domain $\overline{\mathbb{C}}\setminus\{\zeta_0\}$, as to a kernel, as $n\to\infty$. Moreover, for every $\delta>0$ the sequences  $f_n^{- 1}\circ\widetilde{f}_n$ and $\widetilde{f}_n^{-1}\circ{f}_n$ converge locally uniformly to the identity mapping on the set $\overline{E}\setminus \{\zeta_0\}$.
\end{lemma1}

\begin{proof}
Since, by Lemma~\ref{lem1}, the arcs $\beta_n$ are contracted into the point $\zeta_0$, as $n\to+\infty$, we see that $\widetilde{D}_n$ converges to $\overline{\mathbb{C}}\setminus\{\zeta_0\}$.

Denote $g_n:=f_n^{- 1}\circ \widetilde{f}_n$.
By the symmetry principle \cite[Chap.~X, Sect.~1]{lang2003complex},  we can extend the mapping $g_n$ through the arc of the unit disk, complement to  $\beta_n$, in $\widetilde{D}_n$. The extended mapping, which we also, for simplicity of notation, denote by $g_n$, maps $\widetilde{D}_n$ onto $\widetilde{G}_n$. On the other hand, the family of the functions $g_n$ is a normal family in $\overline{\mathbb{C}}\setminus\{\zeta_0\}$. Actually, $g_n$ is bounded in $E$, therefore, in $E$ it forms a normal family. This follows that in  $\overline{\mathbb{C}}\setminus\overline{E}$ it is also normal. Thus, we only need to prove that $g_n$ is normal in a small neighborhood of every point $\omega\in\partial E$ distinct from $\zeta_0$. Let us choose a neighborhood $U$ of $\omega$ not containing the points $\zeta_0$, $0$ and $\infty$. Then in $U$ every $g_n$ does not take  the values $0$, $\infty$ and some $e^{i\phi_n}\in \partial E$. Consequently, $e^{-i\phi_n}g_n$ does not take three fixed values $0$, $1$ and $\infty$. By Montel's theorem (see, e.g., \cite[Chap.~2, Sect.~7]{goluzin1969geometric}, \cite[Chap.~7, Sect.~2]{remmert2013classical}, \cite[Chap.~2, Sect.~2]{schiff1993normal}), $e^{-i\phi_n}g_n$ is a normal family in $U$ and there exists a subsequence $e^{-i\phi_{n_k}}g_{n_k}$ converging locally uniformly in $U$ to some limit function $h$. But $\partial E$ is a compact set, therefore, passing, if necessary, to a subsequence, we can assume that the sequence $e^{i\phi_n}$ converges to some $e^{i\phi_0}\in\partial E$ in $U$. Therefore, $g_{n_k}$ converges to $e^{i\phi_0}h$. This means that the family $g_n$ is normal.

Consider any subsequence $g_{n_k}$ converging locally uniformly to some $g$. Since $g_{n_k}(0)=0$, we see that $g(0)=0$. If $g\equiv \mbox{\rm const}$, then $g\equiv 0$, but from the symmetry of $g_n$ we see that if $g_{n_k}(z)\to 0$, then $g_{n_k}(1/\overline{z})\to \infty$. The controversy proves that $g\not\equiv \mbox{\rm const}$. Then $g$ is a univalent function in $\overline{\mathbb{C}}\setminus\{\zeta_0\}$ as a nonconstant limit of univalent functions (see, e.g., \cite[Chap.~1, Sect.~2]{goluzin1969geometric}), and  $g(0)=0$, $g'(0)>0$. Moreover, $g(1/\overline{z})=1/\overline{g(z)}$. This conditions imply that $g(z)\equiv z$, $z\in \overline{\mathbb{C}}\setminus\{\zeta_0\}$. Therefore, by the generalized Carath\'eodory convergence theorem (see, e.g., \cite[Chap.~5, Sect.~5]{goluzin1969geometric}), $\widetilde{G}_n$ converges to  $\overline{\mathbb{C}}\setminus\{\zeta_0\}$ as $n\to\infty$. By the same theorem, we deduce that the sequence of inverse functions $\widetilde{f}_n^{- 1}\circ {f}_n$ converges locally uniformly in $\overline{\mathbb{C}}\setminus\{\zeta_0\}$ to the identity mapping.
\end{proof}

Now we will complete the proof of Theorem~\ref{thr_local_rado}.

From Lemmas~\ref{wfn} and \ref{lem4} we obtain that $f_n=\widetilde{f}_n\circ (\widetilde{f}_n^{- 1}\circ {f}_n)$ converges locally uniformly to $f$ in $\overline{E}\setminus \{\zeta_0\}$.
Now we will prove that $f^{-1}_n$ converges uniformly to $f^{-1}$ on every compact subset $K$ of $\bigcap_{n=1}^\infty \overline{D}_n\setminus \{w_0\}$. Assume the opposite. Then there exist $\varepsilon>0$ and a subsequence of points $z_{n_k}\in K$  such that
\begin{equation}\label{f1nk}
|f^{-1}_{n_k}(z_{n_k})-f^{-1}(z_{n_k})|\ge \varepsilon, \quad  k\ge 1.
\end{equation}
Denote $\zeta_{n_k}=f^{-1}(z_{n_k})$. Since $\zeta_{n_k}\in \overline{E}$ and $\overline{E}$ is compact, we may assume that $\zeta_{n_k}$ converges to some $\zeta^*\in \overline{E}$. Then  $z_{n_k}=f(\zeta_{n_k})\to f(\zeta^*)=:z^*$. It is cleat that $z^*\in K$, therefore, $z^*\neq w_0$. Consequently,   $\zeta^*\neq \zeta_0$. From the uniform convergence of $f^{-1}_{n_k}$ to
$f^{-1}$ on $K$ it follows that
$$f^{-1}_{n_k}(z_{n_k})=f^{-1}_{n_k}\circ f(\zeta_{n_k})\to \zeta^*.$$
Thus, we have $f^{-1}_{n_k}(z_{n_k})\to f^{-1}(z^*)=\zeta^*$, as $k\to\infty$, but this contradicts to \eqref{f1nk}

Therefore, Theorem \ref{thr_local_rado} is fully verified.
\end{proof}
\medskip

Now we will state a theorem which is a local version of Rad\'o type theorem proven in \cite[Theorem~3]{dautova2019}. In the cited theorem the case of convergence of univalent functions defined on a sequence of doubly-connected domains   converging to a simply-connected domain $G$ is considered. Here we investigate the case where the uniform convergence is violated at some point.

Let $\mathcal{D}_n$ be a sequence of non-concentric annuli, each of which is bounded by the unit circle $\{|\zeta|=1\}$ and the circle $\{|\zeta+ir_n/2|=r_n/2\}$, where $0<r_n<1$ and $\lim\limits_{n\to\infty}r_n=1$. Then, $\mathcal{D}_n$ converges to the circular lune $\mathcal{D}$ bounded by the unit circle and the circle $\{|\zeta+i/2|=1/2\}$.

Let $\mathcal{G}_n$ be a sequence of doubly-connected domains in $\mathbb{C}$, each of which is bounded by Jordan curves $\Gamma_{n1}$ and  $\Gamma_{n2}$, and $\Gamma_{n2}$ is the exterior component of $\partial \mathcal{G}_n$. Let there exist a sequence of crosscuts $\alpha_n$ of domains $\mathcal{G}_n$ with endpoints $a_n$, $b_n$, such that $\lim\limits_{n\to\infty}a_n= \lim\limits_{n\to\infty}b_n=w_0$ for some $w_0\in \mathbb{C}$ and the diameters of $\alpha_n$ tend to $0$, as $n\to\infty$. For every $n$ we denote by $\gamma_n$ a subarc of $\Gamma_{n2}$ with endpoints $a_n, b_n$ such that $\gamma_n\cup\alpha_n$ is the boundary of the Jordan domain containing $\Gamma_{n1}$.  Let the sequence $\mathcal{G}_n$ converge to a simply-connected domain $\mathcal{G}$, the boundary of which consists of two Jordan curves $\Gamma_1$ and $\Gamma_2$, and  $\Gamma_1\cap\Gamma_2=\{w_0\}$ (Fig.~\ref{f3}).

\begin{figure}[ht] \centering
\includegraphics[width=4.5 in]{3.eps}
\caption{The sequence of mappings $f_n$ and its limit.}\label{f3}
\end{figure}

\begin{theorem1}\label{thr_rado_ge1}
Let $f_n$ be a sequence of functions mapping conformally $\mathcal{D}_n$ onto $\mathcal{G}_n$, and $f_n(0)=z_{n0}\in\Gamma_{n1}$. Let the function $f$ map $\mathcal{D}$ onto $\mathcal{G}$, and $f(0)=z_0\in\Gamma_1$, where $z_0\ne w_0$. Assume that $\lim\limits_{n\to\infty}z_{n0}=z_0$, and there exist homeomorphisms
$$\psi_{n1}:\Gamma_1\to\Gamma_{n1},\,\,\, \psi_{n2}:\Gamma_2\to\Gamma_{n2}$$
such that
\begin{equation}\label{eq_rado_ge1}
\forall\varepsilon>0\ \ \exists N:\quad\forall n\geq N\quad |\psi_{n1}(t)-t|<\varepsilon,\ t\in\Gamma_1;\qquad |\psi_{n2}(s)-s|<\varepsilon,\ s\in\Gamma_2.
\end{equation}
Then, $f_n$ converges locally uniformly to $f$ on $\overline{\mathcal{D}}\setminus\{-i\}$ and $f^{-1}_n$ converges uniformly to $f^{-1}$ on every compact subset of $\bigcap_{n=1}^\infty \overline{\mathcal{G}}_n\setminus \{w_0\}$.
\end{theorem1}

To prove the theorem, we need to establish some auxiliary results.

Denote $\alpha_n':=f_n^{-1}(\alpha_n), a_n':=f_n^{-1}(a_n), b_n':=f_n^{-1}(b_n)$. According to the principle of the boundary correspondence under conformal mappings, $\alpha_n'$ is a crosscut in $\mathcal{D}_n$, with endpoints $a_n', b_n'$ on the unit circle.

\begin{lemma1}\label{lem1a}
The crosscuts $\alpha'_n$ are contracted to the point $(-i)$, as $n\to\infty$.
\end{lemma1}

The proof of Lemma~\ref{lem1a} is similar to that of Lemma~\ref{lem1}.\medskip

Let $\mathcal{D}'_n$ be the sequence of doubly-connected domains, each of which is bounded by the circle $\{|\zeta+ir_n/2|=r_n/2\}$ and $\alpha'_n\cup\alpha''_n$, where $\alpha''_n$ is the arc of the unit circle $\{|\zeta|=1\}$ with endpoints $a'_n$, $b'_n$ such that the curve $\alpha'_n\cup\alpha''_n$ contains the circle $\{|\zeta+ir_n/2|=r_n/2\}$ in its interior (Fig.~\ref{f4}).

\begin{figure}[ht] \centering
\includegraphics[width=4.5 in]{4.eps}
\caption{The sequence of mappings $f_n$ on subdomains $\mathcal{D}'_n$.}\label{f4}
\end{figure}

Let $\widetilde{\mathcal{G}}_n$ be a part of $\mathcal{G}_n$, bounded by $\gamma_n\cup\alpha_n$ and $\Gamma_{n1}$ (Fig.~\ref{f5}). From Lemma \ref{lem1a} and the conditions on $\alpha_n $, we easily deduce the following lemma.

\begin{lemma1}
The sequences of domains $\mathcal{D}_n'$ and $\widetilde{\mathcal{G}}_n$ converge to the domains $\mathcal{D}$ and $\mathcal{G}$, respectively, as $n\to\infty$.
\end{lemma1}

Denote by $\widetilde{f}_n$ a sequence of functions mapping conformally  $\widetilde{\mathcal{D}}_n$ onto $\widetilde{\mathcal{G}}_n$ such that $\widetilde{f}_n(0)=z_{n0}\in\Gamma_{n1}$; here $\widetilde{\mathcal{D}}_n$ are nonconcentric annuli bounded by the unit circle and the circles $\{|\zeta+i\widetilde{r}_n/2|=\widetilde{r}_n/2\}$, where $0<\widetilde{r}_n<1$ and $\lim\limits_{n\to\infty}\widetilde{r}_n=1$. The arc $\widetilde{\alpha}_n':=\widetilde{f}^{-1}_n(\alpha_n)$ is a subarc of the unit circle with endpoints $\widetilde{a}_n':=\widetilde{f}_n^{-1}(a_n), \widetilde{b}_n':=\widetilde{f}_n^{-1}(b_n)$. Applying the same arguments as in the proof of Lemma \ref{lem1}, we obtain that the arcs $\widetilde{\alpha}'_n$ are contracted to the point $(-i)$, as $n\to\infty$.

It is easy to see that the sequence $\widetilde{\mathcal{D}}_n$ tends to $\mathcal{D}$ as $n\to\infty$.

Let the conditions of Theorem \ref{thr_rado_ge1} hold. Then it is evident that there exist homeomorphisms $\widetilde{\psi}_{n1}:\Gamma_1\to\Gamma_{n1}$ and $\widetilde{\psi}_{n2}:\Gamma_2\to\gamma_n\cup\alpha_n$ such that
\begin{equation}\label{eq_lem_rado}
\forall\varepsilon>0\ \exists N:\quad\forall n\geq N\quad |\widetilde{\psi}_{n1}(t)-t|<\varepsilon,\ t\in\Gamma_1;\quad |\widetilde{\psi}_{n2}(s)-s|<\varepsilon,\ s\in\Gamma_2.
\end{equation}
From this fact and the generalized Rad\'o theorem \cite[Sect.~4, Theorem.~6]{dautova2018} we obtain

\begin{lemma1}\label{lem_rado}
The sequence $\widetilde{f}_n$ converges uniformly to $f$ in $\overline{\mathcal{D}}$.
\end{lemma1}

\begin{figure}[ht] \centering
\includegraphics[width=4.5 in]{5.eps}
\caption{The sequence of mappings $\widetilde{f}_n$.}\label{f5}
\end{figure}

Denote by ${\mathcal{D}}'_{n1}$ and $\widetilde{\mathcal{D}}_{n1}$ the domains which obtained from the domains $\mathcal{D}'_n$ and $\widetilde{\mathcal{D}}_n$ by gluing them with their reflections with respect to the arcs on the unit circle $|\zeta|=1$ with endpoints $a'_n, b'_n$ and $\widetilde{a}'_n, \widetilde{b}'_n$, complement to $\alpha_n$ and $\alpha'_n$ respectively. It is clear that the sequences
${\mathcal{D}}'_{n1}$ and $\widetilde{\mathcal{D}}_{n1}$  tend to the domain $\mathcal{D}_0=\{|z+1/2|>1/2, \,\text{Im} z>-1\}$ as $n\to \infty$.

Every function $f_n^{-1}\circ \widetilde{f}_n$ maps $\widetilde{\mathcal{D}}_n$ onto $\mathcal{D}'_n$.  According to the Riemann-Schwarz symmetry principle, the functions  $f_n^{-1}\circ \widetilde{f}_n$ can be extended to the domains $\widetilde{\mathcal{D}}_{n1}$; for convenience, we will also denote the extended functions by $f_n^{-1}\circ \widetilde{f}_n$; each of them  maps $\widetilde{\mathcal{D}}_{n1}$ onto
$\widetilde{\mathcal{D}}'_{n1}$.

Applying the same arguments as in the proof of Lemma \ref{lem4} we obtain

\begin{lemma1}\label{lem_rado1}
The sequence of functions $f_n^{-1}\circ\widetilde{f}_n$ converges locally uniformly to the identity mapping on $\mathcal{D}_0\setminus\{-i\}$.
\end{lemma1}

The same statement also holds for the inverse functions $\widetilde{f}_n^{-1}\circ f_n$. Then, with the use of Lemma~\ref{lem_rado}, we obtain that $f_n=\widetilde{f}_n\circ (\widetilde{f}_n^{- 1}\circ {f}_n)$ converges locally uniformly to $f$ in $\mathcal{D}_0\setminus\{-i\}$. The uniform convergence of $f_n^{-1}$ to $f^{-1}$ can be proved just in the proof of similar result from Theorem~\ref{thr_local_rado}.   This proves Theorem~\ref{thr_rado_ge1}.

\begin{remark1}\label{remfaminf}
Further we can apply Theorem~\ref{thr_rado_ge1} not only for sequences of conformal mappings but also for one parameter families of such mappings. We can also apply this theorem for the case of unbounded domains; then we understand the convergence in the spherical metric.
\end{remark1}
\section{Proofs of Theorem~\ref{main}}

First we will show that \eqref{bdtmodulOmegaH} holds.

Consider four the points $A(c, f_1(c)), B(c, f_2(c)), C(d, f_2(d)), D(d, f_1(d))$ on $\partial \Omega$. The vertical segments $AB$ and $CD$ separate $\Omega$ into two domains $\mathcal{Q}$ and $\mathcal{P}$. Let $\mathcal{Q}$ denote the bounded domain and $\mathcal{P}$ be unbounded, i.e.
$$\mathcal{Q}:=\{(x, y)\in\mathbb{C}\,:\, c\leq x\leq d, f_2(x)\leq y\leq f_1(x)\},$$
and $\mathcal{P}=\Omega\setminus \overline{\mathcal{Q}}$.

Consider the quadrilaterals
$$\mbox{\boldmath$Q$}:=\left(\mathcal{Q}; A, B, C, D\right),\quad \mbox{\boldmath$P$}:=\left(\mathcal{P}; D, C, B, A\right)$$
and their images $$
\mbox{\boldmath$Q$}_H:=f_H(\mbox{\boldmath$Q$})=\left(\mathcal{Q}_H; A_H, B_H, C_H, D_H\right),\quad  \mbox{\boldmath$P$}_H:=f_H(\mbox{\boldmath$P$})=\left(\mathcal{P}_H; D_H, C_H, B_H, A_H\right)
$$
under the mapping $f_H$. By the Gr\"otzsch inequality (see, e.g., \cite[Chap.~4, Sect.~6]{goluzin1969geometric},  \cite[Chap.~3, Sect.~7]{kuhnau2005conformal}),  we have
\begin{align}\label{m1}
m^{-1}(\Omega_H)\geq m(\mbox{\boldmath$Q$}_H) +m(\mbox{\boldmath$P$}_H).
\end{align}

From~\eqref{m1} we deduce that $m^{-1}(\Omega_H)\geq m(\mbox{\boldmath$Q$}_H)$ and from \cite[Theorem~4.2]{nvgiang2021}  it follows that
\begin{equation}\label{minv}
m(\mbox{\boldmath$Q$}_H)\ge \gamma H, \quad \gamma=\int\limits_{c
}^{d}\dfrac{dx}{f_1\left({x}\right)-f_2\left({x}\right)}.
\end{equation}
Therefore,
$m^{-1}(\Omega_H)\geq \gamma H$ and this implies~\eqref{bdtmodulOmegaH}. In particular, we have
\begin{equation}\label{momega}
m(\Omega_H)\to 0, \quad\text{ as }\quad H\to\infty.
\end{equation}

Denote by $\varphi_H$ a conformal mapping of $\Omega_H$ onto a concentric annulus
$$
\mathfrak{D}_H:= \{1 < |\zeta| < \rho_H\},
$$
such that the point $A_H$ is mapped to some point on the circle $\{|\zeta|=1\}$, say $\zeta=1$,  and $\varphi_H(B_H)$ is on the circle  $\{|\zeta|=\rho\}$.
We have
$$
m(\Omega_H)=m(\mathfrak{D}_H) = \frac{1}{2\pi}\log\rho_H,
$$
therefore, by \eqref{momega}, $\lim\limits_{H\to\infty}\rho_H=1$.
\begin{lemma1}\label{modulOmega_PQ}
We have
\begin{align}\label{ctOmega_PQ}
m^{-1}(\Omega_H)\sim m(\mbox{\boldmath $P$}_H)+m(\mbox{\boldmath $Q$}_H),\quad\mbox{\rm as }\quad H\to+\infty.
\end{align}
\end{lemma1}
\begin{proof} To show this, we will establish that for every $\varepsilon > 0$ there exists $H_0$ such that for all $H \geq H_0$ the images of the segments $A_HB_H$ and $C_HD_H$ under the map $\varphi_H$ lie in the intersections of the annulus $\mathfrak{D}_H$ and circular sectors of solution less than $\varepsilon$ centered at the origin.

First we consider the images of the segment $A_HB_H$.
The modulus is a conformal invariant, consequently, it is invariant under the shifts of the domain. Therefore, instead of $\Omega$ we can consider the domain which is obtained from $\Omega$ by translating along the real axis by the vector $(-c, 0)$. For convenience, we also denote the obtained domain by $\Omega$. In this case, we really have $c=0$,  therefore  the points $A_H=A$, $B_H=B$ and  the segment $A_HB_H=AB$ do not depend on~$H$.

Denote $r_H=\frac{2\rho_H}{1+\rho_H^2}$. It is clear that  $r_H<1$ and  $r_H\to1$ as $H\to \infty$.

The M\"obius transformation
$$
\psi_H(\zeta)=\rho_H\,\dfrac{\rho_H\zeta+i}{\zeta+i\rho_H}
$$
maps conformally the non-concentric annulus
$$
\mathfrak{G}_H:=\{\zeta\in\mathbb{C}\,:\,|\zeta+ir_H/2|>r_H/2,\, |\zeta|< 1\}.
$$
onto the annulus $\mathfrak{D}_H$.

Consider two cases.

a) Let $a<c=0$. Denote
$$
\mathfrak{G}:=\{\zeta\in\mathbb{C}\,:\,|\zeta+i/2|>1/2\},\, |\zeta|< 1\},
$$
$$
\Omega=\{(x,y):y<f_1(0)\}\setminus \{(x,f_2(0)):x\ge 0\}.
$$
By the definition of the kernel convergence, the family of domains $\Omega_H$ converges to $\Omega$ as to a kernel and also $\mathfrak{G}_H$ converges to $\mathfrak{G}$, as $H\to+\infty$.

 Now consider the conformal mappings  $\mathcal{F}_H:=\psi_H^{-1}\circ \varphi_H$. Then, by the Carath\'eodory convergence theorem~\cite[Chap.~5, Sect.~5]{goluzin1969geometric}, $\mathcal{F}_H$ converge locally uniformly in $\Omega$ to some conformal mapping $\mathcal{F}:\Omega\to\mathfrak{G}$, as $H\to+\infty$.

 Now with the help of Theorem~\ref{thr_rado_ge1}, we will show that the mappings $\mathcal{F}_H$ converge uniformly to the conformal mapping $\mathcal{F}:\Omega\to\mathfrak{G}$ on the segment $AB$. We should note that the limiting domain  $\Omega$ is not Jordan since it is a half plane cut along a ray and formally we can not employ Theorem~\ref{thr_rado_ge1} directly to the sequence $\mathcal{F}_H$ converging to $\mathcal{F}$. But we can apply to $\Omega_H$ additional conformal mappings
 $$
 g_H(z)=\delta_H\sqrt{\frac{D_H(z-A)}{D_H-z}}+c_H,
 $$
with appropriate fixed branches of the square root, where $\delta_H\to 1$ and $c_H\to0$, as $H\to\infty$.
Then the domains $g_H(\Omega_H)$ converge to a Jordan domain which is the image of $\Omega$ under the mapping $g(z)=\sqrt{z-A}$. Obviously, $g_H$ converges to $g$ uniformly on $AB$. We can choose  $\delta_H$ and $c_H$  such that $g_H(AB)\subset \overline{g(\Omega)}$. Let us fix a compact subset $K$ of  $g(\Omega)$ such that $g_H(AB)\subset \overline{g(\Omega)}$. By Theorem~\ref{thr_rado_ge1}, taking into account Remark~\ref{remfaminf}, we obtain that $\mathcal{F}_H\circ g_H^{-1}$ converge to $\mathcal{F}\circ g^{-1}$ on $K$. From this fact we conclude that $\mathcal{F}_H$ converge to $\mathcal{F}$ on $AB$.

Denote $\gamma_H=\varphi_H(AB)$.
Let us fix a neighborhood $U$ of the curve $\Gamma=\mathcal{F}(AB)$ which lies outside of some neighborhood $V$ of the point $(-i)$. Because of the uniform convergence, for sufficiently large $H$ the curves $\Gamma_H=\mathcal{F}_H(AB)$ lie in $U$. Consequently, they do not intersect $V$. We also note that for every $\delta>0$ the M{\"o}bius transforms $\psi_H$
converge uniformly to the constant $1$ on the set $|\zeta+i|\geq\delta$, as $H\to\infty$. This follows that for sufficiently large $H$ the curves $\gamma_H=\psi_H(\Gamma_H)$ lie in an arbitrarily small neighborhood of the point $1$.\smallskip

b) Let now $a=c$. This case is considered similarly to the case $a>c$; we should note that now the kernel of the sequence $\Omega_H$ equals to $\widetilde{\Omega}$ which is the complex plane with excluded rays
$\{(x,f_k(0)):x\ge 0\}$, $k=1$, $2$.

This means that for every $\varepsilon > 0$ there exists $H_0$ and $\alpha$ such that for all $H \geq H_0$ the image of the segments $A_HB_H$ under the map $\varphi_H$ lie in the sector $\alpha<\arg z<\alpha+\varepsilon$. Similarly we prove that for every $\varepsilon > 0$ there exists $H_1$ and $\beta $ such that for all $H \geq H_1$ the image of the segments $C_HD_H$ under the map $\varphi_H$ lie in the sector $\beta<\arg z<\beta +\varepsilon$. Without loss of generality we can assume that $H_0=H_1$ and $\alpha=0<\varepsilon<\beta<\beta+\varepsilon<2\pi$. Comparing the modules of $\varphi_H(\mbox{\boldmath$P$}_H)$ and $\varphi_H(\mbox{\boldmath$Q$}_H)$ with the modules of quadrilaterals lying in the annulus $\mathfrak{D}_H$ inside of the sectors bounded with the rays $\arg z = 0, \varepsilon, \beta,$ and $\beta+\varepsilon$, taking into account the conformal invariance of moduli, we conclude that for $H\ge H_0$
$$
(\beta-\varepsilon)/\log \rho_H<m(\mbox{\boldmath$Q$}_H)<(\beta+\varepsilon)/\log \rho_H,
$$
$$
(2\pi-\beta-\varepsilon)/\log \rho_H<m(\mbox{\boldmath$P$}_H)<(2\pi-\beta+\varepsilon)/\log \rho_H,
$$
Adding the above inequalities, we obtain
\begin{align}\label{modul_PH_QH_OmegaH}
(2\pi-2\varepsilon)/\log \rho_H<m(\mbox{\boldmath$Q$}_H)+m(\mbox{\boldmath$P$}_H)<(2\pi+2\varepsilon)/\log \rho_H, \quad\text{ for all }\quad H\ge H_0.
\end{align}
Because, $\varepsilon>0$ can be chosen arbitrarily small, from \eqref{modul_PH_QH_OmegaH} we get
$$
m(\mbox{\boldmath$Q$}_H)+m(\mbox{\boldmath$P$}_H)\sim 2\pi/\log\rho_H=m^{-1}(\Omega_H).
$$

The lemma is entirely proven.
\end{proof}

Now we  will estimate the behavior of the conformal modulus of $\mbox{\boldmath$P$}_H$ as $H\to+\infty$.

\begin{lemma1}\label{mdmodulP_H}
We have
\begin{align}\label{ctmodulP_H}
m(\mbox{\boldmath$P$}_H)=O(\log H),\quad \mbox{\rm as }\quad H\to+\infty.
\end{align}
\end{lemma1}
\begin{proof}
As in the proof of Lemma~\ref{modulOmega_PQ}, we shift $\Omega$ by the vector $(-c, 0)$ and keep for it and its subdomains the same notation.

Consider the quadrilateral $\mbox{\boldmath$P$}$ (Fig.~\ref{f6}).  Let us fix some real numbers $\kappa, \ell$ and positive $\varepsilon>0$ such that $f_2(c)<\kappa-\varepsilon<\kappa< f_1(c)$, $f_2(d)<\ell-\varepsilon<\ell< f_1(d)$.
If we replace the vertices of $\mbox{\boldmath$P$}$ with the points $d+i\ell$, $d+i(\ell-\varepsilon)$, $c+i(\kappa-\varepsilon)$, $c+i\kappa$,
then the modulus of the obtained quadrilateral $\mbox{\boldmath$P$}^{1}$ will be greater than the modulus of $\mbox{\boldmath$P$}$. Now we choose some rays $L_1$ and $L_2$ which come from the points $c+i\kappa$ and $d+i\ell$ and lie in $\mathcal{P}$. If we replace the side of the quadrilateral $\mbox{\boldmath$P$}^{1}$ with end points   $c+i\kappa$ and $d+i\ell$ by the union of the rays $L_1$ and $L_2$, then we obtain a quadrilateral $\mbox{\boldmath$P$}^{2}$, the modulus of which is greater than $m(\mbox{\boldmath$P$}^{1})$. Further, for sufficiently large $M$ the segment with enpoints $c+i(\kappa-M)$ and  $d+i(\kappa-M)$ is lower than the graph of $y=f_2(x)$. Replacing the boundary arc of the quadrilateral  $\mbox{\boldmath$P$}^{2}$ connecting the points $c+i\kappa$ and $d+i\ell$ and containing the graph of $y=f_2(x)$ with the three-link polygonal line $L$ with vertices at the points $c+i\kappa$, $c+i(\kappa-M)$, $d+i(\ell-M)$  and $d+i\ell$, we obtain a quadrilateral $\mbox{\boldmath$P$}^{3}$, the modulus of which is greater than  $m(\mbox{\boldmath$P$}^{2})$. We point out that the boundary of the quadrilateral consists of the rays
$L_1$, $L_2$ and the polygonal line $L$; its vertices are $d+i\ell$, $d+i(\ell-\varepsilon)$, $c+i(\kappa-\varepsilon)$, and $c+i\kappa$.

Denote by $\mbox{\boldmath$P$}^{3}_H$ the quadrilateral which is obtained from $\mbox{\boldmath$P$}^{3}$ by the stretching with coefficient $H$. It is clear that
\begin{equation}\label{mphph3}
m(\mbox{\boldmath$P$}_H)\le m(\mbox{\boldmath$P$}^3_H),
\end{equation}
therefore, to prove the lemma, we only need to prove that $m(\mbox{\boldmath$P$}^{3}_H)= O(\log H)$, $H\to+\infty$.
Let the equations of the straight lines containing $L_1$, the segment with endpoints $c+i(\kappa-M)$, $d+i(\kappa-M)$ and $L_2$ have the form $y=a_kx+b_k $, $k=1$, $2$,~$3$. Then their images under the stretching are on the lines $y=(a_k/H)x+b_k$, $1\le k\le 3$. Now we define the mapping
$$
\eta_H:x+iy\mapsto x+iv(x,y),
$$
where
\begin{align}\label{defquasi1}
v(x,y)= {\begin{cases}
y-\frac{a_1}{H}\,x-b_1,\quad \mbox{\rm if}\quad x\leq Hc,\\[0.1cm]
y- \frac{a_2}{H}\,x-b_2 -M,\quad \mbox{\rm if }\quad Hc\leq x\leq Hd,\\[0.1cm]
y- \frac{a_3}{H}\,x-b_3,\quad \mbox{\rm if}\quad x\geq Hd.
\end{cases}}
\end{align}
It is a $K$-quasiconformal homeomorphism of the plane with
$$
K=K(H)=\frac{1+k}{1-k},\quad  k=\frac{\mathfrak{a}}{\sqrt{\mathfrak{a}^2+4H^2}},\quad \mathfrak{a}:=\max_{1\le k
\le 3} |a_k|.
$$
Since $K(H)\to 1$, $H\to+\infty$, there exists $K_0\geq 1$ such that $K(H)\le K_0$

Denote by $\mbox{\boldmath$T$}_H$ the quadrilateral which is the image of $\mbox{\boldmath$P$}^{3}_H$ under the mapping $\eta_H$. Then, because of quasiinvariance of conformal modulus under quasiconformal mappings, we have
\begin{equation}\label{p3hk0}
m(\mbox{\boldmath$P$}^{3}_H)\le K(H) m(\mbox{\boldmath$T$}_H)\le  K_0 m(\mbox{\boldmath$T$}_H).
\end{equation}

Because of \eqref{mphph3} and  \eqref{p3hk0},  to prove Lemma~\ref{mdmodulP_H},
we only need to establish that
\begin{equation}\label{thoh}
m(\mbox{\boldmath$T$}_H)=O(\log H),\quad\mbox{\rm as}\quad H\to+\infty.
\end{equation}
Since conformal modulus is invariant under shifts, we can assume that $c=-d$. We note that $\mbox{\boldmath$T$}_H=(T_H; dH, dH-i\varepsilon, -dH-i\varepsilon, -dH)$ where $T_H$ is obtained from the lower half plane by removing the rectangle $[-dH, dH]\times[-M,0]$. The quadrilateral $\mbox{\boldmath$T$}_H$ is symmetric with respect to the imaginary axis. Its modulus is twice times greater than the modulus of  its right half, i.e. the quadrilateral $\mbox{\boldmath$T$}_H^1$ which is the fourth quarter of the plane with the excluded rectangle $[0, dH]\times[-M,0]$, with vertices  $dH$, $dH-i\varepsilon$, $-iM$, and $\infty$. If we replace the side of $\mbox{\boldmath$T$}_H^1$ which is the vertical ray, going down from the point $-iM$ to infinity, with the ray going from this point to the left and shift the obtained quadrilateral by the vector $-dH$ then we obtain the quadrilateral $\mbox{\boldmath$U$}_H=(U; 0, -i\varepsilon,-d_H-iM, \infty)$ the modulus of which is bigger than the modulus of $\mbox{\boldmath$T$}_H^1$. Here $U$ is the domain the boundary of which is the polygonal line consisting of the horizontal rays $\{(x, -M), x\le 0\}$, $\{(x, 0), x\ge 0\}$ and the vertical segment $[-iM, 0]$. We have
\begin{equation}\label{mthuh}
m(\mbox{\boldmath$T$}_H)\le 2m(\mbox{\boldmath$U$}_H).
\end{equation}

Now we need the following auxiliary result.

\begin{lemma1}\label{mdmodulU_H} We have
$$
m(\mbox{\boldmath$U$}_H)=O(\log H), \quad\mbox{\rm as }\quad H\to+\infty.
$$
\end{lemma1}

We note that this lemma, in fact,  is a generalization of  \cite[Sect.~4, Lemma~4.6]{nvgiang2021}.

\begin{proof} We can map conformally the lower half plane $\mathbb{L}$ onto $U$ by a function $g$ such that $g(-1) = -iM$, $g(1) = 0$, $g(\infty) = \infty$. The function $g$ has the form
\begin{equation}\label{gz}
g(\zeta)=\frac{M}{\pi}\int_1^\zeta\sqrt{\frac{\omega+1}{\omega-1}}\,d\omega
\end{equation}
where the branch of the square root is fixed such that it is positive for real $\omega>1$.
Denote by $(-R_H)$ and  $\xi$ the preimages of the points $(-dH-iM)$ and $(-i\varepsilon)$ under the mapping $g(\zeta)$, $R_H>0$, $-1<\xi<1$. From \eqref{gz} it follows that
$$
g(\zeta)\sim \frac{M}{\pi}\,\zeta,\quad\mbox{\rm as }\quad \zeta\to\infty,
$$
therefore,
$$
R_H\sim \frac{M}{\pi}\,dH.
$$
Further, the conformal modulus of $\mbox{\boldmath$U$}_H$ is equal to the conformal modulus of the quadrilateral $\mbox{\boldmath$L$}_H=(\mathbb{L}; \infty, -R_H, \xi, 0)$. In its turn, by the Riemann-Schwarz symmetry principle, $m(\mbox{\boldmath$L$}_H)$ is a half of the modulus of the doubly-connected domain $\mathcal{A}=\overline{\mathbb{C}}\setminus\left( (-\infty, -R_H]\cup [\xi, 1]\right)$, which is an image of the Teichm\"uller ring under a linear transform. Thus, we have (see, e.g., \cite[Chap.3]{ahlfors2006lectures}; \cite{avv})
$$
2\pi m(\mathcal{A})\sim \log R_H\sim\log H, \quad\mbox{\rm as }\quad H\to \infty.
$$
Since $m(\mbox{\boldmath$U$}_H)\le 2 m(\mathcal{A})$, we proved the statement of the lemma.
\end{proof}

From Lemma~\ref{mdmodulU_H} and \eqref{mthuh} it follows that \eqref{thoh} holds. This proves Lemma~\ref{mdmodulP_H}.

At last, combining \eqref{minv}, \eqref{ctOmega_PQ} and \eqref{ctmodulP_H} we obtain \eqref{ctulmd1}.
Therefore, Theorem~\ref{main} is completely proved.
\end{proof}

\section{Discussion}
There are still many interesting issues concerning investigation
of the behavior of conformal moduli of domains
under their geometric transformations. Our next task includes studying the behavior of the exterior conformal moduli of arbitrary quadrilaterals when stretched in the direction of the abscissa axis.
\section*{FUNDING}
The work of the second author is performed under the development program of the Volga Region Mathematical Center (agreement no.~075-02-2022-882).

redpajama_set_name: RedPajamaArXiv