Parto sborone:
$ n=\left(\frac{(22854+8638\sqrt7)(34100354867927167+12888722657083728\sqrt7)^{m}-(22854-8638\sqrt7)(34100354867927167-12888722657083728\sqrt7)^{m}}{2}\right)^4 $ funziona $ \forall m\ge 0 $
Ora una dimostrazioncina...
Lemma avvelenato: Esistono infinite coppie di interi positivi $ x,y $ tali che $ x^2-7y^2=8 $ e $ 13|x $ e $ y\equiv 1\pmod 3 $
Allora... prima di tutto considero la prima condizione... è una quasi-pell... quindi non resta che trovare la soluzione minima e la soluzione minima della pell associata... sono rispettivamente (6,2) e (8,3): $ 6^2-7\cdot 4=8 $ e $ 8^2-7\cdot 9=1 $.
Dato $ x\in\mathbb{Z}[\sqrt{7}] $ definisco $ I(x) $ il coefficiente della radice di 7 e $ R(x) $ quello che resta.
Quindi tutte le soluzioni sono della forma $ R\left((6+2\sqrt{7})(8+3\sqrt{7})^n\right),I\left((6+2\sqrt{7})(8+3\sqrt{7})^n\right) $
Dimostro che se $ n\equiv 3\pmod{14} $ allora sono rispettate anche le altre 2 condizioni

Noto che $ 13|R\left((6+2\sqrt7)(8+3\sqrt{7})^3\right) $ e che $ 13|I\left((8+3\sqrt{7})^{14}\right) $ (non l'avevate notato

) da cui per induzione si ricava $ 13|R\left((6+2\sqrt{7})(8+3\sqrt{7})^{n}\right) $ per $ n\equiv 3\pmod{14} $perciò la seconda condizione è soddisfatta
Noto che per $n$ dispari $ I\left((6+2\sqrt7)(8+3\sqrt7)^n\right)\equiv 1\pmod 3 $ (quest'è più facile) quindi anche la terza e ultima condizione è soddisfatta dato che $ n\equiv 3\pmod{14} $ implica n dispari.
Sia (a,b) una soluzione di tutte le condizioni del lemma avvelenato, dimostro che $ n=b^4 $ funziona. Da questo avrei concluso perchè il lemma avvelenato mi assicura che ci sono infinite di quelle coppie (ed è facile vedere che b cambia sempre).
Prima di tutto vale $ 3|b^2+b+1 $ poichè $ b\equiv 1\pmod 3 $.
Inoltre vale $ 13|a $ quindi $ -7b^2\equiv 8\pmod{13}\rightarrow -14b^2\equiv 16\pmod{13}\Rightarrow b^2+3\equiv 0\pmod{13}\Rightarrow 13|b^2+a+3 $
E ora la magia:
$ n^2+n+1=(b^4+1)^2-b^4=(b^4+b^2+1)(b^4-b^2+1)=((b^2+1)^2-b^2)((b+3)^2-7b^2-8)= $
$ (b^2+b+1)(b^2-b+1)((b^2+3)^2-a^2)=3\cdot\left(\frac{b^2+b+1}{3}\right)(b^2-b+1)(b^2+a+3)(b^2-a+3)= $
$ 3\cdot\left(\frac{b^2+b+1}{3}\right)(b^2-b+1)\cdot 13\cdot \left(\frac{b^2+a+3}{13}\right)(b^2-a+3) $
Ed è facile notare che per b sufficientemente grande tutti i fattori finali non solo sono interi ma anche minori di $ b^2 $ e quindi se un primo divide $ n^2+n+1 $ è minore di $ b^2=\sqrt n $
E dopo la parte sborona e la dimostrazione... un poco di euristica per i più coraggiosi

Un modo stupido di partire è tentare di fare i fini cioè qualsiasi cosa che non sia fattorizzare... io c'ho provato per un pezzo prima di pensare alla fattorizzazione e non ho cavato una mazza... poi ci si illumina e si capisce che in questo problema ci deve essere un modo ganzissimo di fattorizzare... ma dio solo sa quale... la prima cosa che uno nota è che $ n^2+n+1=(n+1)^2-n $ ma allora se n fosse un quadrato avrei fattorizzato, quindi coraggiosamente piazzo n quadrato e mi ritrovo con $ (m^2+m+1)(m^2-m+1) $.
Qui arriva il panico... questi due fattori sono ancora troppo grossi, dovrei rifattorizzare di nuovo... ma insomma questi non c'è verso di farlo facilmente, l'unica cosa che viene in mente è ripiazzare m quadrato così da fattorizzarne almeno 1... se solo riuscissi a fattorizzare anche l'altro...
E qui sta l'idea ganza: finora ho fattorizzato con $ (n+1)^2-n $ ma nessuno mi impedisce di cambiare la costante 1 in un altro numero, quindi lo faccio dopo aver visto che $ m^2-m+1 $ è tosta fattorizzarlo in quel modo se m è un quadrato... allora ottengo $ (m+a)^2-(2a+1)m-(a^2-1) $ quindi se voglio fattorizzare entrambe le robe di prima devo imporre $ m,(2a+1)m-(a^2-1) $ entrambi quadrati perfetti... non resta che scegliere $a$... beh qui pell è telefonato... viene da cercare un $a$ tale che $x^2-(2a+1)y^2=a^2-1$ abbia infinite soluzioni

Ci si arma di pazienza e si nota che a=1,2 non funzionano mentre a=3 pare funzionare con grande gioia... perciò scelgo a=3 da cui esce la mitica $x^2-7y^2=8$
Questo porta alla fattorizzazione $ (y^2+y+1)(y^2-y+1)(y^2+x+3)(y^2-x+3) $
In questo modo si dimostra che infinite volte quella roba non ha fattori primi maggiori di $ (1+\epsilon)\sqrt n $ ma non mi basta... perchè il problema mi chiede proprio minori della radice... quindi resta ancora da mostrare che $y^2+y+1,y^2+x+3$ non sono primi... e qui parte la noia e la rogna

La prima idea è tentare di fattorizzare anche quelli... ma mi sono arreso quasi subito ricordando che x,y hanno dei vincoli pesanti cioè sono soluzioni di una pell e quindi non posso poi scegliere molto su di loro (ad esempio non posso chiedere che siano quadrati).
Ma allora chiunque abbia fatto un poco di problemi olimpici sa che se devo dimostrare che una roba non è prima e non si fattorizza allora tocca trovare un primo che la divide... e qua parte la caccia al tesoro... per $ y^2+y+1 $ è facile decidere... 3 divide quella roba "relativamente spesso" quindi scelgo 3 come primo e ottengo la condizione $y\equiv 1\pmod 3$
Ora tocca lavorare sul ben più tosto $y^2+x+3$ e tocca sperare nella botta di culo

Prima di tutto ho tentato di andare un poco a caso provando con i primi più facili... 2,3,5,7... non si riesce a far nulla di decente, poi un attimo prima di abbandonare ho provato a fare qualcosa di più furbo cioè imporre modulo un certo primo p (che tralascerò):
$ x^2-7y^2=8 $ e $ y^2+x+3=0 $ da cui $x^2-8=7(-x-3)\Rightarrow x^2+7x+13=0$
Perciò come minimo l'ultima equazione deve aver soluzione mod p... che fare? Beh si nota quell'intrigante 13 come termine noto... e si prova la sminchiata p=13 e 13|x e per pura magia si scopre (non è poi una gran scoperta, ma la soddisfazione non manca

) che quello unito a $x^2-7y^2=8$ implica quello che vogliamo, cioè $13|x^2+y+3$ e che non da nessun assurdo noioso causa residui quadratici controproducenti come succedeva per gli altri primi.
Ora siamo arrivati alle richieste del lemma avvelenato... non è per nulla scontato che ci siano davvero tutte le coppie cercate, pell non ci basta
Ci si arma di santa pazienza e si risolve la pell arrivando a $(6+2\sqrt{7})(8+3\sqrt{7})^n$
Da questa si vede (questo è facile davvero

) che per la condizione del $y\equiv 1\pmod 3$ è sufficiente n dispari... una gran bella soddisfazione e si ricomincia a sperare... resta quel cazzutissimo 13, che darà seri problemi.
Che fare? L'unica è iniziare a fare una smadonna di conti, magari infurbandosi

La prima cosa da controllare è che ci sia ALMENO un n dispari per cui vale la richiesta del 13... si trova dopo ben meno conti di quanto ci aspetti che n=3 funziona miracolosamente

Ora tocca mostrare che da una soluzione della richiesta se ne trova un'altra, magari più grossa... ci si pensa un poco e si trova che in effetti è equivalente a trovare un k tale che $13|I\left((8+3\sqrt{7})^k\right)$ inoltre se si vuole che vada bene anche il 3 questo k dev'essere pari... qui i conti sono abbastanza ma non troppi se si è furbi mod 13

Come fare beh tipo definisco $x_n+y_n\sqrt 7=((8+3\sqrt 7)^2)^n$ e con qualche conto si ottiene che mod 13 vale
$(x_{n+1},y_{n+1})=(x_n+2(x_n+y_n),x_n+3(x_n+y_n))$
E questo per fare i conti non è male, se ci si ricorda di lavorare mod 13... si parte da (1,0) e si vuole arrivare a (a,0)... si scopre che in "soli" 7 passaggi ci si riesce
E si arriva al miracolo di k=14... e il problema è chiuso, non resta che ricontrollare il mare di conti che si è fatto e scoprire con gioia che non ci sono errori
Infine per i sopravvissuti: è in qualche modo generalizzabile la roba sulle pell con aggiunta la condizione di divisibilità?
Di questa non so proprio nulla... ma mi pareva caruccia come domanda...
p.s. mi viene il dubbio che ci sia una soluzione terribilmente più elegante di questa... ma vabbè ha un suo che d'istruttivo questa tra tutti i contazzi

Oppure potrebbe essere che si riusciva a dire che quei 2 fattori non erano primi infinite volte senza complicarsi così tanto la vita come ho fatto io... perchè fino a lì la soluzione è carina, ma non sono riuscito a farlo

p.p.s. quale è la fonte?