Motility-Induced Phase Separation(MIPS)

Grupo: Bernardo Boatini e Murilo Kessler Azambuja

Introdução

A matéria ativa é um tipo de sistema fora do equilíbrio termodinâmico, onde cada "partícula" ou "agente" do sistema tem a capacidade dissipar energia na forma de forças mecânicas exercidas sobre o ambiente no qual está imerso. Esses sistemas muitas vezes podem exibir vários fenômenos novos, como movimentos coletivos quando exige-se um alinhamento das partículas (Como bio-polímeros que se auto organizam, tipo os microtúbulos que são parte do citoesqueleto celular, ou cardumes de peixe, por exemplo), ou os chamados MIPS (Motility-Induced Phase Separation) que são sistemas que apresentam uma mudança de fase física devido à interações que proíbem a ocupação simultânea de um volume do espaço por duas partículas simultâneas. Nestes sistemas, os agentes possuem a capacidede de auto-propulsão, de forma que a condição de balanço detalhado é quebrada no nível microscópico, uma vez que as partículas possuem uma direção preferencial de movimento.Estes tipos de sistema muitas vezes podem ser tratados a partir de uma abordagem hidrodinâmica, na qual são chamados de fluidos ativos, e divergem do comportamento usual de fluidos compostos de matéria usual inativa, os quais são descritos pela equação de Navier-Stokes.

Apesar desses sistemas fundamentalmente quebrarem a condição do balanço detalhado, ainda não é claro se o comportamento estatístico universal de sistemas de fluidos ativos necessariamente divergem daqueles de de sistemas em equilíbrio. A investigação de comportamentos universais nestes tipos de sistema, além de ser um interece central na física, também pode nos permitir utilizar conhecimentos já bem conhecidos (como a transição de fase em sistemas de matéria em equilíbrio termodinâmico) para descrever sistemas novos. Neste trabalho foi reproduzido o artigo onde foi investigada o comportamento crítico de um sistema com MIPS e se viu que este tipo de sistema pertence à classe de universalidade de sistemas em equilíbrio como o modelo de Ising para spin. Para tal, foi utilizada três abordagens diferentes: um modelo hidrodinâmico, uma descrição baseada em teoria de campos e a simulação de um modelo de rede hexagonal.

Modelo

O sistema geral que estamos interessados em descrever é constituído de um conjunto de partículas em um meio com friccção, i.e, em um sistema sem conservação de momentum e com interações que não permitem a ocupação simultânea de duas partículas em um mesmo sítio (FIG. 1).

FIGURA 1: Imagem retirada de . Em (a) vemos a difusão devido à rotação da direção de movimento preferencial da partícula, em (b) temos o "movimento balístico" (seta verde) e a difusão devido à translação da partícula (setas azuis) e, em (c), vemos a interação entre partículas em sítios vizinhos, impedindo a ocupação simultânea do mesmo sítio por duas partículas.

Uma forma de modelar este sistema é considerá-lo um fluido ativo polar e compressível sem interações de alinhamento entre partículas, de forma que não se espera observar movimentos coletivos. Tal sistema é descrito de forma contínua pela chamada equação de Toner-Tu.

De um ponto de vista microscópico, podemos usar uma formulação discreta baseada na teoria de campos. Neste caso, consideramos uma rede hexagonal bidimensional, com o numero de ocupação de cada sítio da rede restringido por uma constante $M$ (futuramente diremos que $M=1$, de forma que apenas uma partícula pode ocupar cada sítio).

Modelo no Continuo

As equações de Tuner-Tu são dadas por

$\partial _{t}\rho +\nabla \cdot \mathbf {g} =0,\qquad (1a)$

$\partial _{t}\mathbf {g} =-\zeta \nabla \rho -\kappa \mathbf {g} +\mu \nabla ^{2}\mathbf {g} +\mathbf {f} ,\qquad (1b)$

onde $\rho$ é o campo de densidade de massa e $\mathbf {g}$ é o campo de densidade de momentum. A equação (1a) diz respeito à conservação de massa do sistema como um todo. A equação (1b) se assemelha muito à equação de Navier-Stokes, assim podemos ver que o primeiro termo da equação caracteriza a compressibilidade do fluido, o segundo termo representa uma difusão do campo de momentum (o sinal negativo indica que o campo de momentum diminui com o passar do tempo), o terceiro termo representa a viscosidade do fluido e $\mathbf {f}$ é um ruído Gaussiano com estatísticas espaço-temporais dadas por

$\langle \mathbf {f} (t,\mathbf {r} )\rangle =0,\ \langle \mathbf {f} (t,\mathbf {r} )\mathbf {f} (t',\mathbf {r'} )\rangle =2D\delta (t-t')\delta (\mathbf {r} -\mathbf {r'} ),$

onde $D$ é a intensidade do ruído.

No mesmo artigo o autor argumenta que, devido à ausência de interações de alinhamento entre partículas, a emergência de movimento coletivo é impossível. Por este motivo, temos que o campo de densidade de momentum deve ir à zero no limite hidrodinâmico, de forma que $\mathbf {g}$ não é uma variável livre e pode ser escrita em função de $\rho$. Desta forma, obtêm que

$\partial _{t}\phi =\nabla ^{2}{\frac {\delta H}{\delta \phi }}+\nabla \mathbf {f} ,\qquad (2)$

onde $\phi (\mathbf {r} )=\rho (\mathbf {r} )-\rho _{0}$ para alguma constante $\rho _{0}$ e $H={\frac {a_{1}}{2}}\phi ^{2}+{\frac {a_{2}}{3}}\phi ^{3}+{\frac {a_{3}}{4}}\phi ^{4}+{\frac {b}{2}}{\big (}\nabla \phi {\big )}^{2}$ é o Hamiltoniano de Landau-Ginzburg. A partir da equação (2), vemos que é possível descrever o sistema a partir da manipulação de dois parâmetros principais: o campo de densidade de partículas $\rho$ (ou $\phi$, a menos de uma constante) e a intensidade $D$ do ruído $\mathbf {f}$.

Modelo no Discreto

Como já foi dito antes, no modelo discreto é utilizado uma rede bidimensional hexagonal com número de ocupação máximo $M$ para cada sítio da rede. Os agentes são divididos em seis tipos distintos de partículas ativas, cada uma se move em uma direção $\theta _{i}$ (as diagonais do hexágono) e somente vai saltar para o sítio adjacente ao longo desta direção, com uma taxa que depende da ocupação do sítio em questão. Além disso, é possível que uma partícula se transforme em outra, com uma taxa $\sigma$, correspondendo assim a um ruído rotacional de cada partícula. Utilizando a notação utilizada em , temos que a rede hexagonal é caracterizado pelo conjunto de vetores $\mathbf {R}$ e a orientação ou tipo da partícula ativa pelo índice $i$. Desta forma o número de partículas do tipo $i$ ocupando o sítio $\mathbf {R}$ da rede é denotado por $A_{\mathbf {R} }^{\theta _{i}}$. Com o formalismo desenvolvido em , temos que a teoria de campos nos fornece uma ação dada por

$S=\int dt\sum _{\mathbf {R} }{\Bigg \{}\sum _{i}{\Bigg [}{\hat {A}}_{\mathbf {R} }^{\theta _{i}}\partial _{t}A_{\mathbf {R} }^{\theta _{i}}\qquad (3a)$

$-\gamma {\Big (}M-N_{\mathbf {R} +\mathbf {e} _{\theta _{i}}}{\Big )}A_{\mathbf {R} }^{\theta _{i}}{\Big (}e^{{\hat {A}}_{\mathbf {R} +\mathbf {e} _{\theta _{i}}}^{\theta _{i}}-{\hat {A}}_{\mathbf {R} }^{\theta _{i}}}-1{\Big )}{\Bigg ]}\qquad (3b)$

$-\sum _{\langle i,j\rangle }\sigma A_{\mathbf {R} }^{\theta _{i}}{\Big (}e^{{\hat {A}}_{\mathbf {R} }^{\theta _{j}}-{\hat {A}}_{\mathbf {R} }^{\theta _{i}}}-1{\Big )}{\Bigg \}},\qquad (3c)$

onde $N_{\mathbf {R} }=\sum _{j}A_{\mathbf {R} }^{\theta _{j}}$ é o número de partículas total ocupando o sítio $\mathbf {R}$. A equação (3b), diz respeito ao termo de transição entre uma posição $\mathbf {R}$ e uma posição $\mathbf {R} +\mathbf {e} _{\theta _{i}}$ na rede hexagonal. Esta transição ocorre à uma taxa $\gamma {\Big (}M-N_{\mathbf {R} +\mathbf {e} _{\theta _{i}}}{\Big )}$, de forma que a taxa de transição para o sítio decresce com o aumento do número de partículas ocupando o mesmo. Se o número de partículas ocupando o sítio é $M$, então a taxa de transição é nula, indicando que não será possível que ocorra uma transição para este sítio. Este termo modela as interações de exclusão de volume discutidas anteriormente. Em nosso caso específico tomaremos $M=1$. Além disso, a equação (3c) corresponde à transição de tipo da partícula, i.e., à transição do movimento preferencial da partícula de uma orientação $\theta _{i}$ para uma $\theta _{j}$. Esta transição ocorre com uma taxa $\sigma$ e a notação $\langle i,j\rangle$ indica que a partícula somente trocará para uma orientação $\theta _{j}$ que seja adjascente à orientação $\theta _{i}$ inicial.

Implementação

Nessa secção serão colocados alguns detalhes da implementação da rede, das medidas e da dinâmica de Monte Carlo.

Algoritmo

A dinâmica de Monte Carlo foi implementada de forma que a cada MCS, N partículas do sistema são selecionadas aleatoriamente:

- Primeiramente a partícula decide se vai mudar seu sentido($\theta$) sorteando um valor $d$ através de uma distribuição gaussiana centrada em 0, com variância igula a intencidade do ruído $\sigma$;

- A nova direção é dada por Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \theta + d.60°} ;

- Depois disso a partícula tem probabilidade $24/30$ de andar para o sitio vizinho naquela direção(se ele estiver vazio);

- Caso a partícula não decida andar naquela direção, ela se desloca na direção de outro sitio vizinho(se ele estiver vazio) sorteado aleatoriamente(emulando o comportamento difusivo).

Redes Hexagonais

Em modelos do tipo LGCA(Lattice Gas Celular Automata) o tipo de rede, bem como as topologias usadas, podem ser determinantes para as medidas e caracterizações do estado macroscópico. A simetria hexagonal se mostra uma boa alternativa para aumentar os graus de liberdade de uma partícula numa rede, quando comparada com a simetria quadrada. Na geometria hexagonal, assim como na quadrada, todos os sítios vizinhos possuem uma mesma distancia entre si só que com 6 graus de liberdade ao invés de quatro.

FIGURA 2: Imagem adaptada de . As setas duplas em verde sinalizam os sentidos que seguem trajetórias fechadas passando por só uma borda na topologia "armchair"(a) e "zigzag" (b)

Existem duas topologias possíveis em uma rede hexagonal: "zigzag" e "armchair", e a diferença das duas está principalmente na implementação das condições de contorno periódicas. No caso da topologia em "zigzag"(FIG.2b), o espaço é representado em um paralelogramo, e existem 4 sentidos em que cruzando uma parede a partícula fecha sua trajetória sem ter que cruzar nova mente outra parede do sistema. Na topologia "armchair"(FIG.2a) por outro lado só existem 2 sentidos em que isso acontece.

Sub-Box Sampling Method

O método das subcaixas consiste em dividir o sistema em uma certa quantidade de subsistemas de interesse, permitindo observar a flutuação de densidades apenas nessas regiões.

Nesse trabalho o sistema foi dividido em 4 caixas de tamanho $L\times L$: 2 centradas no centro de massa horizontal do sistema e duas centradas na região horizontal "mais rarefeita", que consequentemente esta a uma distancia de $3L$ do centro de massa.

Quando os parâmetros permitem que haja separação de fase, duas subcaixas ficam posicionadas exatamente sobre o centro da fase liquida e a mesma coisa para a fase gasosa, permitindo uma medição mais precisa das flutuações de densidade sem considerar as regiões de fronteira.

Porém, essas flutuações se tornam especialmente relevantes no ponto crítico, onde o uso de subcaixas se torna especialmente importante na extração medidas para classificar a transição de fase do modelo em analise.

Medida 1

A primeira medida utilizada para comparar o MIPS com o Ising foi a incompressibilidade isotérmica dos subsistemas $\chi _{L}$ definida da seguinte forma

$\chi _{L}={\frac {\langle N^{2}\rangle _{L}-\langle N\rangle _{L}^{2}}{\langle N\rangle _{L}}}$

onde $\langle N\rangle _{L}$ denota a média do número de partículas contidas em um subsistema de tamanho L.

A Incompressibilidade isotérmica possui um fator de escala com o tamanho do sistema($\chi _{L}\alpha L^{\gamma /\nu }$), que é conhecido para o caso do Ising 2D ($\gamma /\nu =7/4$).

Medida 2

A segunda medida proposta é a derivada do cumulante de quarta ordem ($U_{4}$) com relação a distancia adimensional do ponto crítico ($\tau ={\frac {\sigma -\sigma _{c}}{\sigma }}$) calculada no ponto crítico($\tau =0$). O cumulante é definido da seguinte forma

$U_{4}={\frac {\langle \Delta \rho ^{2}\rangle _{L}^{2}}{\langle \Delta \rho ^{4}\rangle _{L}}}$

$\Delta \rho =\rho -\langle \rho \rangle _{L}$

Em $\langle .\rangle _{L}$ é descrito como: "[...] an average over the four sub-boxes of finite size L [...]"

O cumulante de quarta ordem possui a propriedade de ser invariante frente a escala do sistema no ponto crítico(FIG.4), e por isso a sua derivada no ponto critico se torna um fator de escala importante para comparar o estado macroscópico do MIPS e do Ising.

FIGURA 3: Imagem retirada de . Mostra o comportamento esperado do cumulante de quarta ordem para o modelo de MIPS nas proximidades da criticidade em sistemas de diferentes tamanhos.

Em $\tau =0$ a derivada do cumulante de quarta ordem, com respeito a $\tau$, possui um fator de escala $\alpha L^{1/\nu }$, que tambem pode ser encontrado no modelo de Ising($\nu =1$).

Resultados

Mostrar e explicar os Resultados

FIGURA 4:

FIGURA 5:

FIGURA 6:

FIGURA 7:

FIGURA 8: